1 /* floating-point mod function without infinity or NaN checking */ 2 #include <math.h> 3 double fmod(double x,double y)4fmod (double x, double y) 5 { 6 int sign = 0, yexp; 7 double r, yfr; 8 9 if (y == 0) 10 return 0; 11 if (y < 0) 12 y = -y; 13 yfr = frexp (y, &yexp); 14 if (x < 0) { 15 sign = 1; 16 r = -x; 17 } else 18 r = x; 19 while (r >= y) { 20 int rexp; 21 double rfr = frexp (r, &rexp); 22 r -= ldexp (y, rexp - yexp - (rfr < yfr)); 23 } 24 if (sign) 25 r = -r; 26 return r; 27 } 28