]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/e_fmod.c
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, };
24 __ieee754_fmod (double x
, double y
)
26 int32_t n
, hx
, hy
, hz
, ix
, iy
, sx
, i
;
29 EXTRACT_WORDS (hx
, lx
, x
);
30 EXTRACT_WORDS (hy
, ly
, y
);
31 sx
= hx
& 0x80000000; /* sign of x */
33 hy
&= 0x7fffffff; /* |y| */
35 /* purge off exception values */
36 if ((hy
| ly
) == 0 || (hx
>= 0x7ff00000) || /* y=0,or x not finite */
37 ((hy
| ((ly
| -ly
) >> 31)) > 0x7ff00000)) /* or y is NaN */
38 return (x
* y
) / (x
* y
);
41 if ((hx
< hy
) || (lx
< ly
))
42 return x
; /* |x|<|y| return x */
44 return Zero
[(u_int32_t
) sx
>> 31]; /* |x|=|y| return x*0*/
47 /* determine ix = ilogb(x) */
48 if (__glibc_unlikely (hx
< 0x00100000)) /* subnormal x */
52 for (ix
= -1043, i
= lx
; i
> 0; i
<<= 1)
57 for (ix
= -1022, i
= (hx
<< 11); i
> 0; i
<<= 1)
62 ix
= (hx
>> 20) - 1023;
64 /* determine iy = ilogb(y) */
65 if (__glibc_unlikely (hy
< 0x00100000)) /* subnormal y */
69 for (iy
= -1043, i
= ly
; i
> 0; i
<<= 1)
74 for (iy
= -1022, i
= (hy
<< 11); i
> 0; i
<<= 1)
79 iy
= (hy
>> 20) - 1023;
81 /* set up {hx,lx}, {hy,ly} and align y to x */
82 if (__glibc_likely (ix
>= -1022))
83 hx
= 0x00100000 | (0x000fffff & hx
);
84 else /* subnormal x, shift x to normal */
89 hx
= (hx
<< n
) | (lx
>> (32 - n
));
98 if (__glibc_likely (iy
>= -1022))
99 hy
= 0x00100000 | (0x000fffff & hy
);
100 else /* subnormal y, shift y to normal */
105 hy
= (hy
<< n
) | (ly
>> (32 - n
));
119 hz
= hx
- hy
; lz
= lx
- ly
; if (lx
< ly
)
123 hx
= hx
+ hx
+ (lx
>> 31); lx
= lx
+ lx
;
127 if ((hz
| lz
) == 0) /* return sign(x)*0 */
128 return Zero
[(u_int32_t
) sx
>> 31];
129 hx
= hz
+ hz
+ (lz
>> 31); lx
= lz
+ lz
;
132 hz
= hx
- hy
; lz
= lx
- ly
; if (lx
< ly
)
139 /* convert back to floating value and restore the sign */
140 if ((hx
| lx
) == 0) /* return sign(x)*0 */
141 return Zero
[(u_int32_t
) sx
>> 31];
142 while (hx
< 0x00100000) /* normalize x */
144 hx
= hx
+ hx
+ (lx
>> 31); lx
= lx
+ lx
;
147 if (__glibc_likely (iy
>= -1022)) /* normalize output */
149 hx
= ((hx
- 0x00100000) | ((iy
+ 1023) << 20));
150 INSERT_WORDS (x
, hx
| sx
, lx
);
152 else /* subnormal output */
157 lx
= (lx
>> n
) | ((u_int32_t
) hx
<< (32 - n
));
162 lx
= (hx
<< (32 - n
)) | (lx
>> n
); hx
= sx
;
166 lx
= hx
>> (n
- 32); hx
= sx
;
168 INSERT_WORDS (x
, hx
| sx
, lx
);
169 x
*= one
; /* create necessary signal */
171 return x
; /* exact output */
173 strong_alias (__ieee754_fmod
, __fmod_finite
)