1 /* $NetBSD: n_fmod.c,v 1.11 2018/11/09 10:19:47 martin 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. Neither the name of the University nor the names of its contributors 15 * may be used to endorse or promote products derived from this software 16 * without specific prior written permission. 17 * 18 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 19 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 20 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 21 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 22 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 23 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 24 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 25 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 26 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 27 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 28 * SUCH DAMAGE. 29 */ 30 31 #ifndef lint 32 #if 0 33 static char sccsid[] = "@(#)fmod.c 8.1 (Berkeley) 6/4/93"; 34 #endif 35 #endif /* not lint */ 36 37 #include "mathimpl.h" 38 39 /* fmod.c 40 * 41 * SYNOPSIS 42 * 43 * #include <math.h> 44 * double fmod(double x, double y) 45 * 46 * DESCRIPTION 47 * 48 * The fmod function computes the floating-point remainder of x/y. 49 * 50 * RETURNS 51 * 52 * The fmod function returns the value x-i*y, for some integer i 53 * such that, if y is nonzero, the result has the same sign as x and 54 * magnitude less than the magnitude of y. 55 * 56 * On a VAX or CCI, 57 * 58 * fmod(x,0) traps/faults on floating-point divided-by-zero. 59 * 60 * On IEEE-754 conforming machines with "isnan()" primitive, 61 * 62 * fmod(x,0), fmod(INF,y) are invalid operations and NaN is returned. 63 * 64 */ 65 #if !defined(__vax__) && !defined(tahoe) 66 extern int isnan(),finite(); 67 #endif /* !defined(__vax__) && !defined(tahoe) */ 68 69 #if DBL_MANT_DIG == LDBL_MANT_DIG && DBL_MIN_EXP == LDBL_MIN_EXP \ 70 && DBL_MIN_EXP == LDBL_MIN_EXP 71 #ifdef __weak_alias 72 __weak_alias(fmodl, fmod); 73 __weak_alias(modfl, fmod); 74 #endif 75 #endif 76 77 #ifdef TEST_FMOD 78 static double 79 _fmod(double x, double y) 80 #else /* TEST_FMOD */ 81 double 82 fmod(double x, double y) 83 #endif /* TEST_FMOD */ 84 { 85 int ir,iy; 86 double r,w; 87 88 if (y == (double)0 89 #if !defined(__vax__) && !defined(tahoe) /* per "fmod" manual entry, SunOS 4.0 */ 90 || isnan(y) || !finite(x) 91 #endif /* !defined(__vax__) && !defined(tahoe) */ 92 ) 93 return (x*y)/(x*y); 94 95 r = fabs(x); 96 y = fabs(y); 97 (void)frexp(y,&iy); 98 while (r >= y) { 99 (void)frexp(r,&ir); 100 w = ldexp(y,ir-iy); 101 r -= w <= r ? w : w*(double)0.5; 102 } 103 return x >= (double)0 ? r : -r; 104 } 105 106 float 107 fmodf(float x, float y) 108 { 109 return fmod(x, y); 110 } 111 112 #ifdef TEST_FMOD 113 #include <math.h> 114 #include <stdio.h> 115 #include <stdlib.h> 116 117 #define NTEST 10000 118 #define NCASES 3 119 120 static int nfail = 0; 121 static void 122 prf(const char *s, double d) 123 { 124 union { 125 double d; 126 unsigned long long u; 127 } x; 128 x.d = d; 129 printf("%s = %#016.16llx (%24.16e)\n", s, x.u, x.d); 130 } 131 132 static void 133 doit(double x, double y) 134 { 135 double ro = fmod(x,y),rn = _fmod(x,y); 136 if (ro != rn) { 137 prf(" x ", x); 138 prf(" y ", y); 139 prf(" fmod", ro); 140 prf("_fmod", rn); 141 (void)printf("\n"); 142 } 143 } 144 145 int 146 main(int argc, char **argv) 147 { 148 int i,cases; 149 double x,y; 150 151 srandom(12345); 152 for (i = 0; i < NTEST; i++) { 153 x = (double)random(); 154 y = (double)random(); 155 for (cases = 0; cases < NCASES; cases++) { 156 switch (cases) { 157 case 0: 158 break; 159 case 1: 160 y = (double)1/y; break; 161 case 2: 162 x = (double)1/x; break; 163 default: 164 abort(); break; 165 } 166 doit(x,y); 167 doit(x,-y); 168 doit(-x,y); 169 doit(-x,-y); 170 } 171 } 172 if (nfail) 173 (void)printf("Number of failures: %d (out of a total of %d)\n", 174 nfail,NTEST*NCASES*4); 175 else 176 (void)printf("No discrepancies were found\n"); 177 exit(0); 178 } 179 #endif /* TEST_FMOD */ 180