1 /* @(#)e_fmod.c 5.1 93/09/24 */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
10 * ====================================================
12 * $NetBSD: e_fmod.c,v 1.11 2002/05/26 22:01:49 wiz Exp $
13 * $DragonFly: src/lib/libm/src/e_fmod.c,v 1.1 2005/07/26 21:15:20 joerg Exp $
18 * Return x mod y in exact arithmetic
19 * Method: shift and subtract
23 #include "math_private.h"
25 static const double one
= 1.0, Zero
[] = {0.0, -0.0,};
28 fmod(double x
, double y
)
30 int32_t n
,hx
,hy
,hz
,ix
,iy
,sx
,i
;
33 EXTRACT_WORDS(hx
,lx
,x
);
34 EXTRACT_WORDS(hy
,ly
,y
);
35 sx
= hx
&0x80000000; /* sign of x */
37 hy
&= 0x7fffffff; /* |y| */
39 /* purge off exception values */
40 if((hy
|ly
)==0||(hx
>=0x7ff00000)|| /* y=0,or x not finite */
41 ((hy
|((ly
|-ly
)>>31))>0x7ff00000)) /* or y is NaN */
44 if((hx
<hy
)||(lx
<ly
)) return x
; /* |x|<|y| return x */
46 return Zero
[(u_int32_t
)sx
>>31]; /* |x|=|y| return x*0*/
49 /* determine ix = ilogb(x) */
50 if(hx
<0x00100000) { /* subnormal x */
52 for (ix
= -1043, i
=lx
; i
>0; i
<<=1) ix
-=1;
54 for (ix
= -1022,i
=(hx
<<11); i
>0; i
<<=1) ix
-=1;
56 } else ix
= (hx
>>20)-1023;
58 /* determine iy = ilogb(y) */
59 if(hy
<0x00100000) { /* subnormal y */
61 for (iy
= -1043, i
=ly
; i
>0; i
<<=1) iy
-=1;
63 for (iy
= -1022,i
=(hy
<<11); i
>0; i
<<=1) iy
-=1;
65 } else iy
= (hy
>>20)-1023;
67 /* set up {hx,lx}, {hy,ly} and align y to x */
69 hx
= 0x00100000|(0x000fffff&hx
);
70 else { /* subnormal x, shift x to normal */
73 hx
= (hx
<<n
)|(lx
>>(32-n
));
81 hy
= 0x00100000|(0x000fffff&hy
);
82 else { /* subnormal y, shift y to normal */
85 hy
= (hy
<<n
)|(ly
>>(32-n
));
96 hz
=hx
-hy
;lz
=lx
-ly
; if(lx
<ly
) hz
-= 1;
97 if(hz
<0){hx
= hx
+hx
+(lx
>>31); lx
= lx
+lx
;}
99 if((hz
|lz
)==0) /* return sign(x)*0 */
100 return Zero
[(u_int32_t
)sx
>>31];
101 hx
= hz
+hz
+(lz
>>31); lx
= lz
+lz
;
104 hz
=hx
-hy
;lz
=lx
-ly
; if(lx
<ly
) hz
-= 1;
105 if(hz
>=0) {hx
=hz
;lx
=lz
;}
107 /* convert back to floating value and restore the sign */
108 if((hx
|lx
)==0) /* return sign(x)*0 */
109 return Zero
[(u_int32_t
)sx
>>31];
110 while(hx
<0x00100000) { /* normalize x */
111 hx
= hx
+hx
+(lx
>>31); lx
= lx
+lx
;
114 if(iy
>= -1022) { /* normalize output */
115 hx
= ((hx
-0x00100000)|((iy
+1023)<<20));
116 INSERT_WORDS(x
,hx
|sx
,lx
);
117 } else { /* subnormal output */
120 lx
= (lx
>>n
)|((u_int32_t
)hx
<<(32-n
));
123 lx
= (hx
<<(32-n
))|(lx
>>n
); hx
= sx
;
125 lx
= hx
>>(n
-32); hx
= sx
;
127 INSERT_WORDS(x
,hx
|sx
,lx
);
128 x
*= one
; /* create necessary signal */
130 return x
; /* exact output */