00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #if defined(LIBM_SCCS) && !defined(lint)
00014 static char rcsid[] = "$NetBSD: e_remainder.c,v 1.8 1995/05/10 20:46:05 jtc Exp $";
00015 #endif
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "math.h"
00027 #include "mathP.h"
00028
00029 #ifdef __STDC__
00030 static const double zero = 0.0;
00031 #else
00032 static double zero = 0.0;
00033 #endif
00034
00035
00036 #ifdef __STDC__
00037 double __ieee754_remainder(double x, double p)
00038 #else
00039 double __ieee754_remainder(x,p)
00040 double x,p;
00041 #endif
00042 {
00043 int32_t hx,hp;
00044 u_int32_t sx,lx,lp;
00045 double p_half;
00046
00047 EXTRACT_WORDS(hx,lx,x);
00048 EXTRACT_WORDS(hp,lp,p);
00049 sx = hx&0x80000000;
00050 hp &= 0x7fffffff;
00051 hx &= 0x7fffffff;
00052
00053
00054 if((hp|lp)==0) return (x*p)/(x*p);
00055 if((hx>=0x7ff00000)||
00056 ((hp>=0x7ff00000)&&
00057 (((hp-0x7ff00000)|lp)!=0)))
00058 return (x*p)/(x*p);
00059
00060
00061 if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);
00062 if (((hx-hp)|(lx-lp))==0) return zero*x;
00063 x = fabs(x);
00064 p = fabs(p);
00065 if (hp<0x00200000) {
00066 if(x+x>p) {
00067 x-=p;
00068 if(x+x>=p) x -= p;
00069 }
00070 } else {
00071 p_half = 0.5*p;
00072 if(x>p_half) {
00073 x-=p;
00074 if(x>=p_half) x -= p;
00075 }
00076 }
00077 GET_HIGH_WORD(hx,x);
00078 SET_HIGH_WORD(x,hx^sx);
00079 return x;
00080 }