2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
9 * ====================================================
14 * Return x mod y in exact arithmetic
15 * Method: shift and subtract
19 #include "math_private.h"
21 static const double one
= 1.0, Zero
[] = {0.0, -0.0,};
23 double __ieee754_fmod(double x
, double y
)
25 int32_t n
,hx
,hy
,hz
,ix
,iy
,sx
,i
;
28 EXTRACT_WORDS(hx
,lx
,x
);
29 EXTRACT_WORDS(hy
,ly
,y
);
30 sx
= hx
&0x80000000; /* sign of x */
32 hy
&= 0x7fffffff; /* |y| */
34 /* purge off exception values */
35 if((hy
|ly
)==0||(hx
>=0x7ff00000)|| /* y=0,or x not finite */
36 ((hy
|((ly
|-ly
)>>31))>0x7ff00000)) /* or y is NaN */
39 if((hx
<hy
)||(lx
<ly
)) return x
; /* |x|<|y| return x */
41 return Zero
[(u_int32_t
)sx
>>31]; /* |x|=|y| return x*0*/
44 /* determine ix = ilogb(x) */
45 if(hx
<0x00100000) { /* subnormal x */
47 for (ix
= -1043, i
=lx
; i
>0; i
<<=1) ix
-=1;
49 for (ix
= -1022,i
=(hx
<<11); i
>0; i
<<=1) ix
-=1;
51 } else ix
= (hx
>>20)-1023;
53 /* determine iy = ilogb(y) */
54 if(hy
<0x00100000) { /* subnormal y */
56 for (iy
= -1043, i
=ly
; i
>0; i
<<=1) iy
-=1;
58 for (iy
= -1022,i
=(hy
<<11); i
>0; i
<<=1) iy
-=1;
60 } else iy
= (hy
>>20)-1023;
62 /* set up {hx,lx}, {hy,ly} and align y to x */
64 hx
= 0x00100000|(0x000fffff&hx
);
65 else { /* subnormal x, shift x to normal */
68 hx
= (hx
<<n
)|(lx
>>(32-n
));
76 hy
= 0x00100000|(0x000fffff&hy
);
77 else { /* subnormal y, shift y to normal */
80 hy
= (hy
<<n
)|(ly
>>(32-n
));
91 hz
=hx
-hy
;lz
=lx
-ly
; if(lx
<ly
) hz
-= 1;
92 if(hz
<0){hx
= hx
+hx
+(lx
>>31); lx
= lx
+lx
;}
94 if((hz
|lz
)==0) /* return sign(x)*0 */
95 return Zero
[(u_int32_t
)sx
>>31];
96 hx
= hz
+hz
+(lz
>>31); lx
= lz
+lz
;
99 hz
=hx
-hy
;lz
=lx
-ly
; if(lx
<ly
) hz
-= 1;
100 if(hz
>=0) {hx
=hz
;lx
=lz
;}
102 /* convert back to floating value and restore the sign */
103 if((hx
|lx
)==0) /* return sign(x)*0 */
104 return Zero
[(u_int32_t
)sx
>>31];
105 while(hx
<0x00100000) { /* normalize x */
106 hx
= hx
+hx
+(lx
>>31); lx
= lx
+lx
;
109 if(iy
>= -1022) { /* normalize output */
110 hx
= ((hx
-0x00100000)|((iy
+1023)<<20));
111 INSERT_WORDS(x
,hx
|sx
,lx
);
112 } else { /* subnormal output */
115 lx
= (lx
>>n
)|((u_int32_t
)hx
<<(32-n
));
118 lx
= (hx
<<(32-n
))|(lx
>>n
); hx
= sx
;
120 lx
= hx
>>(n
-32); hx
= sx
;
122 INSERT_WORDS(x
,hx
|sx
,lx
);
123 x
*= one
; /* create necessary signal */
125 return x
; /* exact output */
128 strong_alias(__ieee754_fmod
, fmod
)
129 libm_hidden_def(fmod
)