1 /* e_fmodl.c -- long double version of e_fmod.c.
4 * ====================================================
5 * Copyright (C) 1993, 2011 by Sun Microsystems, Inc. All rights reserved.
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
11 * ====================================================
15 * __ieee754_fmodl(x,y)
16 * Return x mod y in exact arithmetic
17 * Method: shift and subtract
21 #include <math_private.h>
22 #include <libm-alias-finite.h>
24 static const _Float128 one
= 1.0, Zero
[] = {0.0, -0.0,};
27 __ieee754_fmodl (_Float128 x
, _Float128 y
)
29 int64_t n
,hx
,hy
,hz
,ix
,iy
,sx
,i
;
32 GET_LDOUBLE_WORDS64(hx
,lx
,x
);
33 GET_LDOUBLE_WORDS64(hy
,ly
,y
);
34 sx
= hx
&0x8000000000000000ULL
; /* sign of x */
36 hy
&= 0x7fffffffffffffffLL
; /* |y| */
38 /* purge off exception values */
39 if((hy
|ly
)==0||(hx
>=0x7fff000000000000LL
)|| /* y=0,or x not finite */
40 ((hy
|((ly
|-ly
)>>63))>0x7fff000000000000LL
)) /* or y is NaN */
43 if((hx
<hy
)||(lx
<ly
)) return x
; /* |x|<|y| return x */
45 return Zero
[(uint64_t)sx
>>63]; /* |x|=|y| return x*0*/
48 /* determine ix = ilogb(x) */
49 if(hx
<0x0001000000000000LL
) { /* subnormal x */
51 for (ix
= -16431, i
=lx
; i
>0; i
<<=1) ix
-=1;
53 for (ix
= -16382, i
=hx
<<15; i
>0; i
<<=1) ix
-=1;
55 } else ix
= (hx
>>48)-0x3fff;
57 /* determine iy = ilogb(y) */
58 if(hy
<0x0001000000000000LL
) { /* subnormal y */
60 for (iy
= -16431, i
=ly
; i
>0; i
<<=1) iy
-=1;
62 for (iy
= -16382, i
=hy
<<15; i
>0; i
<<=1) iy
-=1;
64 } else iy
= (hy
>>48)-0x3fff;
66 /* set up {hx,lx}, {hy,ly} and align y to x */
68 hx
= 0x0001000000000000LL
|(0x0000ffffffffffffLL
&hx
);
69 else { /* subnormal x, shift x to normal */
72 hx
= (hx
<<n
)|(lx
>>(64-n
));
80 hy
= 0x0001000000000000LL
|(0x0000ffffffffffffLL
&hy
);
81 else { /* subnormal y, shift y to normal */
84 hy
= (hy
<<n
)|(ly
>>(64-n
));
95 hz
=hx
-hy
;lz
=lx
-ly
; if(lx
<ly
) hz
-= 1;
96 if(hz
<0){hx
= hx
+hx
+(lx
>>63); lx
= lx
+lx
;}
98 if((hz
|lz
)==0) /* return sign(x)*0 */
99 return Zero
[(uint64_t)sx
>>63];
100 hx
= hz
+hz
+(lz
>>63); lx
= lz
+lz
;
103 hz
=hx
-hy
;lz
=lx
-ly
; if(lx
<ly
) hz
-= 1;
104 if(hz
>=0) {hx
=hz
;lx
=lz
;}
106 /* convert back to floating value and restore the sign */
107 if((hx
|lx
)==0) /* return sign(x)*0 */
108 return Zero
[(uint64_t)sx
>>63];
109 while(hx
<0x0001000000000000LL
) { /* normalize x */
110 hx
= hx
+hx
+(lx
>>63); lx
= lx
+lx
;
113 if(iy
>= -16382) { /* normalize output */
114 hx
= ((hx
-0x0001000000000000LL
)|((iy
+16383)<<48));
115 SET_LDOUBLE_WORDS64(x
,hx
|sx
,lx
);
116 } else { /* subnormal output */
119 lx
= (lx
>>n
)|((uint64_t)hx
<<(64-n
));
122 lx
= (hx
<<(64-n
))|(lx
>>n
); hx
= sx
;
124 lx
= hx
>>(n
-64); hx
= sx
;
126 SET_LDOUBLE_WORDS64(x
,hx
|sx
,lx
);
127 x
*= one
; /* create necessary signal */
129 return x
; /* exact output */
131 libm_alias_finite (__ieee754_fmodl
, __fmodl
)