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