xref: /csrg-svn/lib/libm/common_source/fmod.c (revision 38203)
1 /*
2  * Copyright (c) 1989 The Regents of the University of California.
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms are permitted
6  * provided that the above copyright notice and this paragraph are
7  * duplicated in all such forms and that any documentation,
8  * advertising materials, and other materials related to such
9  * distribution and use acknowledge that the software was developed
10  * by the University of California, Berkeley.  The name of the
11  * University may not be used to endorse or promote products derived
12  * from this software without specific prior written permission.
13  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
14  * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
15  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
16  */
17 
18 #ifndef lint
19 static char sccsid[] = "@(#)fmod.c	5.1 (Berkeley) 05/30/89";
20 #endif /* not lint */
21 
22 /* fmod.c
23  *
24  * SYNOPSIS
25  *
26  *    #include <math.h>
27  *    double fmod(double x, double y)
28  *
29  * DESCRIPTION
30  *
31  *    The fmod function computes the floating-point remainder of x/y.
32  *
33  * RETURNS
34  *
35  *    The fmod function returns the value x-i*y, for some integer i
36  * such that, if y is nonzero, the result has the same sign as x and
37  * magnitude less than the magnitude of y.
38  *
39  * On a VAX or CCI,
40  *
41  *    fmod(x,0) traps/faults on floating-point divided-by-zero.
42  *
43  * On IEEE-754 conforming machines with "isnan()" primitive,
44  *
45  *    fmod(x,0), fmod(INF,y) are invalid operations and NaN is returned.
46  *
47  */
48 #if !defined(vax) && !defined(tahoe)
49 extern int isnan(),finite();
50 #endif	/* !defined(vax) && !defined(tahoe) */
51 extern double frexp(),ldexp(),fabs();
52 
53 #ifdef TEST_FMOD
54 static double
55 _fmod(x,y)
56 #else	/* TEST_FMOD */
57 double
58 fmod(x,y)
59 #endif	/* TEST_FMOD */
60 double x,y;
61 {
62 	int ir,iy;
63 	double r,w;
64 
65 	if (y == (double)0
66 #if !defined(vax) && !defined(tahoe)	/* per "fmod" manual entry, SunOS 4.0 */
67 		|| isnan(y) || !finite(x)
68 #endif	/* !defined(vax) && !defined(tahoe) */
69 	    )
70 	    return (x*y)/(x*y);
71 
72 	r = fabs(x);
73 	y = fabs(y);
74 	(void)frexp(y,&iy);
75 	while (r >= y) {
76 		(void)frexp(r,&ir);
77 		w = ldexp(y,ir-iy);
78 		r -= w <= r ? w : w*(double)0.5;
79 	}
80 	return x >= (double)0 ? r : -r;
81 }
82 
83 #ifdef TEST_FMOD
84 extern long random();
85 extern double fmod();
86 
87 #define	NTEST	10000
88 #define	NCASES	3
89 
90 static int nfail = 0;
91 
92 static void
93 doit(x,y)
94 double x,y;
95 {
96 	double ro = fmod(x,y),rn = _fmod(x,y);
97 	if (ro != rn) {
98 		(void)printf(" x    = 0x%08.8x %08.8x (%24.16e)\n",x,x);
99 		(void)printf(" y    = 0x%08.8x %08.8x (%24.16e)\n",y,y);
100 		(void)printf(" fmod = 0x%08.8x %08.8x (%24.16e)\n",ro,ro);
101 		(void)printf("_fmod = 0x%08.8x %08.8x (%24.16e)\n",rn,rn);
102 		(void)printf("\n");
103 	}
104 }
105 
106 main()
107 {
108 	register int i,cases;
109 	double x,y;
110 
111 	srandom(12345);
112 	for (i = 0; i < NTEST; i++) {
113 		x = (double)random();
114 		y = (double)random();
115 		for (cases = 0; cases < NCASES; cases++) {
116 			switch (cases) {
117 			case 0:
118 				break;
119 			case 1:
120 				y = (double)1/y; break;
121 			case 2:
122 				x = (double)1/x; break;
123 			default:
124 				abort(); break;
125 			}
126 			doit(x,y);
127 			doit(x,-y);
128 			doit(-x,y);
129 			doit(-x,-y);
130 		}
131 	}
132 	if (nfail)
133 		(void)printf("Number of failures: %d (out of a total of %d)\n",
134 			nfail,NTEST*NCASES*4);
135 	else
136 		(void)printf("No discrepancies were found\n");
137 	exit(0);
138 }
139 #endif	/* TEST_FMOD */
140