1 /*
2 * Copyright (c) 1989, 1993
3 * The Regents of the University of California. All rights reserved.
4 *
5 * %sccs.include.redist.c%
6 */
7
8 #ifndef lint
9 static char sccsid[] = "@(#)fmod.c 8.1 (Berkeley) 06/04/93";
10 #endif /* not lint */
11
12 /* fmod.c
13 *
14 * SYNOPSIS
15 *
16 * #include <math.h>
17 * double fmod(double x, double y)
18 *
19 * DESCRIPTION
20 *
21 * The fmod function computes the floating-point remainder of x/y.
22 *
23 * RETURNS
24 *
25 * The fmod function returns the value x-i*y, for some integer i
26 * such that, if y is nonzero, the result has the same sign as x and
27 * magnitude less than the magnitude of y.
28 *
29 * On a VAX or CCI,
30 *
31 * fmod(x,0) traps/faults on floating-point divided-by-zero.
32 *
33 * On IEEE-754 conforming machines with "isnan()" primitive,
34 *
35 * fmod(x,0), fmod(INF,y) are invalid operations and NaN is returned.
36 *
37 */
38 #if !defined(vax) && !defined(tahoe)
39 extern int isnan(),finite();
40 #endif /* !defined(vax) && !defined(tahoe) */
41 extern double frexp(),ldexp(),fabs();
42
43 #ifdef TEST_FMOD
44 static double
_fmod(x,y)45 _fmod(x,y)
46 #else /* TEST_FMOD */
47 double
48 fmod(x,y)
49 #endif /* TEST_FMOD */
50 double x,y;
51 {
52 int ir,iy;
53 double r,w;
54
55 if (y == (double)0
56 #if !defined(vax) && !defined(tahoe) /* per "fmod" manual entry, SunOS 4.0 */
57 || isnan(y) || !finite(x)
58 #endif /* !defined(vax) && !defined(tahoe) */
59 )
60 return (x*y)/(x*y);
61
62 r = fabs(x);
63 y = fabs(y);
64 (void)frexp(y,&iy);
65 while (r >= y) {
66 (void)frexp(r,&ir);
67 w = ldexp(y,ir-iy);
68 r -= w <= r ? w : w*(double)0.5;
69 }
70 return x >= (double)0 ? r : -r;
71 }
72
73 #ifdef TEST_FMOD
74 extern long random();
75 extern double fmod();
76
77 #define NTEST 10000
78 #define NCASES 3
79
80 static int nfail = 0;
81
82 static void
doit(x,y)83 doit(x,y)
84 double x,y;
85 {
86 double ro = fmod(x,y),rn = _fmod(x,y);
87 if (ro != rn) {
88 (void)printf(" x = 0x%08.8x %08.8x (%24.16e)\n",x,x);
89 (void)printf(" y = 0x%08.8x %08.8x (%24.16e)\n",y,y);
90 (void)printf(" fmod = 0x%08.8x %08.8x (%24.16e)\n",ro,ro);
91 (void)printf("_fmod = 0x%08.8x %08.8x (%24.16e)\n",rn,rn);
92 (void)printf("\n");
93 }
94 }
95
main()96 main()
97 {
98 register int i,cases;
99 double x,y;
100
101 srandom(12345);
102 for (i = 0; i < NTEST; i++) {
103 x = (double)random();
104 y = (double)random();
105 for (cases = 0; cases < NCASES; cases++) {
106 switch (cases) {
107 case 0:
108 break;
109 case 1:
110 y = (double)1/y; break;
111 case 2:
112 x = (double)1/x; break;
113 default:
114 abort(); break;
115 }
116 doit(x,y);
117 doit(x,-y);
118 doit(-x,y);
119 doit(-x,-y);
120 }
121 }
122 if (nfail)
123 (void)printf("Number of failures: %d (out of a total of %d)\n",
124 nfail,NTEST*NCASES*4);
125 else
126 (void)printf("No discrepancies were found\n");
127 exit(0);
128 }
129 #endif /* TEST_FMOD */
130