ref: 53aac62a54ebbcc1781c0000d3e2384ed038c7bb
dir: /libmath/fdlibm/e_remainder.c/
/* derived from /netlib/fdlibm */ /* @(#)e_remainder.c 1.3 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* __ieee754_remainder(x,p) * Return : * returns x REM p = x - [x/p]*p as if in infinite * precise arithmetic, where [x/p] is the (infinite bit) * integer nearest x/p (in half way case choose the even one). * Method : * Based on fmod() return x-[x/p]chopped*p exactlp. */ #include "fdlibm.h" static const double zero = 0.0; double __ieee754_remainder(double x, double p) { int hx,hp; unsigned sx,lx,lp; double p_half; hx = __HI(x); /* high word of x */ lx = __LO(x); /* low word of x */ hp = __HI(p); /* high word of p */ lp = __LO(p); /* low word of p */ sx = hx&0x80000000; hp &= 0x7fffffff; hx &= 0x7fffffff; /* purge off exception values */ if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */ if((hx>=0x7ff00000)|| /* x not finite */ ((hp>=0x7ff00000)&& /* p is NaN */ (((hp-0x7ff00000)|lp)!=0))) return (x*p)/(x*p); if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */ if (((hx-hp)|(lx-lp))==0) return zero*x; x = fabs(x); p = fabs(p); if (hp<0x00200000) { if(x+x>p) { x-=p; if(x+x>=p) x -= p; } } else { p_half = 0.5*p; if(x>p_half) { x-=p; if(x>=p_half) x -= p; } } __HI(x) ^= sx; return x; }