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