1 /* $NetBSD: n_fmod.c,v 1.3 1998/10/20 02:26:11 matt Exp $ */ 2 /* 3 * Copyright (c) 1989, 1993 4 * The Regents of the University of California. All rights reserved. 5 * 6 * Redistribution and use in source and binary forms, with or without 7 * modification, are permitted provided that the following conditions 8 * are met: 9 * 1. Redistributions of source code must retain the above copyright 10 * notice, this list of conditions and the following disclaimer. 11 * 2. Redistributions in binary form must reproduce the above copyright 12 * notice, this list of conditions and the following disclaimer in the 13 * documentation and/or other materials provided with the distribution. 14 * 3. All advertising materials mentioning features or use of this software 15 * must display the following acknowledgement: 16 * This product includes software developed by the University of 17 * California, Berkeley and its contributors. 18 * 4. Neither the name of the University nor the names of its contributors 19 * may be used to endorse or promote products derived from this software 20 * without specific prior written permission. 21 * 22 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 23 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 24 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 25 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 26 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 27 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 28 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 29 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 30 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 31 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 32 * SUCH DAMAGE. 33 */ 34 35 #ifndef lint 36 #if 0 37 static char sccsid[] = "@(#)fmod.c 8.1 (Berkeley) 6/4/93"; 38 #endif 39 #endif /* not lint */ 40 41 #include "mathimpl.h" 42 43 /* fmod.c 44 * 45 * SYNOPSIS 46 * 47 * #include <math.h> 48 * double fmod(double x, double y) 49 * 50 * DESCRIPTION 51 * 52 * The fmod function computes the floating-point remainder of x/y. 53 * 54 * RETURNS 55 * 56 * The fmod function returns the value x-i*y, for some integer i 57 * such that, if y is nonzero, the result has the same sign as x and 58 * magnitude less than the magnitude of y. 59 * 60 * On a VAX or CCI, 61 * 62 * fmod(x,0) traps/faults on floating-point divided-by-zero. 63 * 64 * On IEEE-754 conforming machines with "isnan()" primitive, 65 * 66 * fmod(x,0), fmod(INF,y) are invalid operations and NaN is returned. 67 * 68 */ 69 #if !defined(__vax__) && !defined(tahoe) 70 extern int isnan(),finite(); 71 #endif /* !defined(__vax__) && !defined(tahoe) */ 72 73 #ifdef TEST_FMOD 74 static double 75 _fmod(x,y) 76 #else /* TEST_FMOD */ 77 double 78 fmod(x,y) 79 #endif /* TEST_FMOD */ 80 double x,y; 81 { 82 int ir,iy; 83 double r,w; 84 85 if (y == (double)0 86 #if !defined(__vax__) && !defined(tahoe) /* per "fmod" manual entry, SunOS 4.0 */ 87 || isnan(y) || !finite(x) 88 #endif /* !defined(__vax__) && !defined(tahoe) */ 89 ) 90 return (x*y)/(x*y); 91 92 r = fabs(x); 93 y = fabs(y); 94 (void)frexp(y,&iy); 95 while (r >= y) { 96 (void)frexp(r,&ir); 97 w = ldexp(y,ir-iy); 98 r -= w <= r ? w : w*(double)0.5; 99 } 100 return x >= (double)0 ? r : -r; 101 } 102 103 #ifdef TEST_FMOD 104 extern long random(); 105 extern double fmod(); 106 107 #define NTEST 10000 108 #define NCASES 3 109 110 static int nfail = 0; 111 112 static void 113 doit(x,y) 114 double x,y; 115 { 116 double ro = fmod(x,y),rn = _fmod(x,y); 117 if (ro != rn) { 118 (void)printf(" x = 0x%08.8x %08.8x (%24.16e)\n",x,x); 119 (void)printf(" y = 0x%08.8x %08.8x (%24.16e)\n",y,y); 120 (void)printf(" fmod = 0x%08.8x %08.8x (%24.16e)\n",ro,ro); 121 (void)printf("_fmod = 0x%08.8x %08.8x (%24.16e)\n",rn,rn); 122 (void)printf("\n"); 123 } 124 } 125 126 main() 127 { 128 register int i,cases; 129 double x,y; 130 131 srandom(12345); 132 for (i = 0; i < NTEST; i++) { 133 x = (double)random(); 134 y = (double)random(); 135 for (cases = 0; cases < NCASES; cases++) { 136 switch (cases) { 137 case 0: 138 break; 139 case 1: 140 y = (double)1/y; break; 141 case 2: 142 x = (double)1/x; break; 143 default: 144 abort(); break; 145 } 146 doit(x,y); 147 doit(x,-y); 148 doit(-x,y); 149 doit(-x,-y); 150 } 151 } 152 if (nfail) 153 (void)printf("Number of failures: %d (out of a total of %d)\n", 154 nfail,NTEST*NCASES*4); 155 else 156 (void)printf("No discrepancies were found\n"); 157 exit(0); 158 } 159 #endif /* TEST_FMOD */ 160