1 /* e_fmodl.c -- long double version of e_fmod.c.
2 * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
5 * ====================================================
6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 * Developed at SunPro, a Sun Microsystems, Inc. business.
9 * Permission to use, copy, modify, and distribute this
10 * software is freely granted, provided that this notice
12 * ====================================================
16 * __ieee754_fmodl(x,y)
17 * Return x mod y in exact arithmetic
18 * Method: shift and subtract
22 #include "math_private.h"
25 static const long double one
= 1.0, Zero
[] = {0.0, -0.0,};
27 static long double one
= 1.0, Zero
[] = {0.0, -0.0,};
31 long double __ieee754_fmodl(long double x
, long double y
)
33 long double __ieee754_fmodl(x
,y
)
37 int64_t n
,hx
,hy
,hz
,ix
,iy
,sx
,i
;
40 GET_LDOUBLE_WORDS64(hx
,lx
,x
);
41 GET_LDOUBLE_WORDS64(hy
,ly
,y
);
42 sx
= hx
&0x8000000000000000ULL
; /* sign of x */
44 hy
&= 0x7fffffffffffffffLL
; /* |y| */
46 /* purge off exception values */
47 if((hy
|ly
)==0||(hx
>=0x7fff000000000000LL
)|| /* y=0,or x not finite */
48 ((hy
|((ly
|-ly
)>>63))>0x7fff000000000000LL
)) /* or y is NaN */
51 if((hx
<hy
)||(lx
<ly
)) return x
; /* |x|<|y| return x */
53 return Zero
[(u_int64_t
)sx
>>63]; /* |x|=|y| return x*0*/
56 /* determine ix = ilogb(x) */
57 if(hx
<0x0001000000000000LL
) { /* subnormal x */
59 for (ix
= -16431, i
=lx
; i
>0; i
<<=1) ix
-=1;
61 for (ix
= -16382, i
=hx
<<15; i
>0; i
<<=1) ix
-=1;
63 } else ix
= (hx
>>48)-0x3fff;
65 /* determine iy = ilogb(y) */
66 if(hy
<0x0001000000000000LL
) { /* subnormal y */
68 for (iy
= -16431, i
=ly
; i
>0; i
<<=1) iy
-=1;
70 for (iy
= -16382, i
=hy
<<15; i
>0; i
<<=1) iy
-=1;
72 } else iy
= (hy
>>48)-0x3fff;
74 /* set up {hx,lx}, {hy,ly} and align y to x */
76 hx
= 0x0001000000000000LL
|(0x0000ffffffffffffLL
&hx
);
77 else { /* subnormal x, shift x to normal */
80 hx
= (hx
<<n
)|(lx
>>(64-n
));
88 hy
= 0x0001000000000000LL
|(0x0000ffffffffffffLL
&hy
);
89 else { /* subnormal y, shift y to normal */
92 hy
= (hy
<<n
)|(ly
>>(64-n
));
103 hz
=hx
-hy
;lz
=lx
-ly
; if(lx
<ly
) hz
-= 1;
104 if(hz
<0){hx
= hx
+hx
+(lx
>>63); lx
= lx
+lx
;}
106 if((hz
|lz
)==0) /* return sign(x)*0 */
107 return Zero
[(u_int64_t
)sx
>>63];
108 hx
= hz
+hz
+(lz
>>63); lx
= lz
+lz
;
111 hz
=hx
-hy
;lz
=lx
-ly
; if(lx
<ly
) hz
-= 1;
112 if(hz
>=0) {hx
=hz
;lx
=lz
;}
114 /* convert back to floating value and restore the sign */
115 if((hx
|lx
)==0) /* return sign(x)*0 */
116 return Zero
[(u_int64_t
)sx
>>63];
117 while(hx
<0x0001000000000000LL
) { /* normalize x */
118 hx
= hx
+hx
+(lx
>>63); lx
= lx
+lx
;
121 if(iy
>= -16382) { /* normalize output */
122 hx
= ((hx
-0x0001000000000000LL
)|((iy
+16383)<<48));
123 SET_LDOUBLE_WORDS64(x
,hx
|sx
,lx
);
124 } else { /* subnormal output */
127 lx
= (lx
>>n
)|((u_int64_t
)hx
<<(64-n
));
130 lx
= (hx
<<(64-n
))|(lx
>>n
); hx
= sx
;
132 lx
= hx
>>(n
-64); hx
= sx
;
134 SET_LDOUBLE_WORDS64(x
,hx
|sx
,lx
);
135 x
*= one
; /* create necessary signal */
137 return x
; /* exact output */