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