/* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, 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 "math.h" #include "math_private.h" static const double zero = 0.0; double attribute_hidden __ieee754_remainder(double x, double p) { int32_t hx,hp; u_int32_t sx,lx,lp; double p_half; EXTRACT_WORDS(hx,lx,x); EXTRACT_WORDS(hp,lp,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; } } GET_HIGH_WORD(hx,x); SET_HIGH_WORD(x,hx^sx); return x; } /* * wrapper remainder(x,p) */ #ifndef _IEEE_LIBM double remainder(double x, double y) { double z = __ieee754_remainder(x, y); if (_LIB_VERSION == _IEEE_ || isnan(y)) return z; if (y == 0.0) return __kernel_standard(x, y, 28); /* remainder(x,0) */ return z; } strong_alias(remainder, drem) #else strong_alias(__ieee754_remainder, remainder) strong_alias(__ieee754_remainder, drem) #endif libm_hidden_def(remainder)