| /* | |
| * ==================================================== | |
| * 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. | |
| * ==================================================== | |
| */ | |
| #include "bsd_private.h" | |
| double bsd__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 ((long double)x*p)/((long double)x*p); | |
| if (hp<=0x7fdfffff) x = bsd__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); | |
| if ((hx&0x7fffffff)==0) hx = 0; | |
| SET_HIGH_WORD(x,hx^sx); | |
| return x; | |
| } |