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