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 }