1 /* e_fmodl.c -- long double version of e_fmod.c.
2 */
3 /*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
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
10 * is preserved.
11 * ====================================================
12 */
13
14 /*
15 * __ieee754_fmodl(x,y)
16 * Return x mod y in exact arithmetic
17 * Method: shift and subtract
18 */
19
20 #include <math.h>
21 #include <math_private.h>
22 #include <ieee754.h>
23 #include <libm-alias-finite.h>
24
25 static const long double one = 1.0, Zero[] = {0.0, -0.0,};
26
27 long double
28 __ieee754_fmodl (long double x, long double y)
29 {
30 int64_t hx, hy, hz, sx, sy;
31 uint64_t lx, ly, lz;
32 int n, ix, iy;
33 double xhi, xlo, yhi, ylo;
34
35 ldbl_unpack (x, &xhi, &xlo);
36 EXTRACT_WORDS64 (hx, xhi);
37 EXTRACT_WORDS64 (lx, xlo);
38 ldbl_unpack (y, &yhi, &ylo);
39 EXTRACT_WORDS64 (hy, yhi);
40 EXTRACT_WORDS64 (ly, ylo);
41 sx = hx&0x8000000000000000ULL; /* sign of x */
42 hx ^= sx; /* |x| */
43 sy = hy&0x8000000000000000ULL; /* sign of y */
44 hy ^= sy; /* |y| */
45
46 /* purge off exception values */
47 if(__builtin_expect(hy==0 ||
48 (hx>=0x7ff0000000000000LL)|| /* y=0,or x not finite */
49 (hy>0x7ff0000000000000LL),0)) /* or y is NaN */
50 return (x*y)/(x*y);
51 if (__glibc_unlikely (hx <= hy))
52 {
53 /* If |x| < |y| return x. */
54 if (hx < hy)
55 return x;
56 /* At this point the absolute value of the high doubles of
57 x and y must be equal. */
58 if ((lx & 0x7fffffffffffffffLL) == 0
59 && (ly & 0x7fffffffffffffffLL) == 0)
60 /* Both low parts are zero. The result should be an
61 appropriately signed zero, but the subsequent logic
62 could treat them as unequal, depending on the signs
63 of the low parts. */
64 return Zero[(uint64_t) sx >> 63];
65 /* If the low double of y is the same sign as the high
66 double of y (ie. the low double increases |y|)... */
67 if (((ly ^ sy) & 0x8000000000000000LL) == 0
68 /* ... then a different sign low double to high double
69 for x or same sign but lower magnitude... */
70 && (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy))
71 /* ... means |x| < |y|. */
72 return x;
73 /* If the low double of x differs in sign to the high
74 double of x (ie. the low double decreases |x|)... */
75 if (((lx ^ sx) & 0x8000000000000000LL) != 0
76 /* ... then a different sign low double to high double
77 for y with lower magnitude (we've already caught
78 the same sign for y case above)... */
79 && (int64_t) (lx ^ sx) > (int64_t) (ly ^ sy))
80 /* ... means |x| < |y|. */
81 return x;
82 /* If |x| == |y| return x*0. */
83 if ((lx ^ sx) == (ly ^ sy))
84 return Zero[(uint64_t) sx >> 63];
85 }
86
87 /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
88 bit mantissa so the following operations will give the correct
89 result. */
90 ldbl_extract_mantissa(&hx, &lx, &ix, x);
91 ldbl_extract_mantissa(&hy, &ly, &iy, y);
92
93 if (__glibc_unlikely (ix == -IEEE754_DOUBLE_BIAS))
94 {
95 /* subnormal x, shift x to normal. */
96 while ((hx & (1LL << 48)) == 0)
97 {
98 hx = (hx << 1) | (lx >> 63);
99 lx = lx << 1;
100 ix -= 1;
101 }
102 }
103
104 if (__glibc_unlikely (iy == -IEEE754_DOUBLE_BIAS))
105 {
106 /* subnormal y, shift y to normal. */
107 while ((hy & (1LL << 48)) == 0)
108 {
109 hy = (hy << 1) | (ly >> 63);
110 ly = ly << 1;
111 iy -= 1;
112 }
113 }
114
115 /* fix point fmod */
116 n = ix - iy;
117 while(n--) {
118 hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
119 if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
120 else {
121 if((hz|lz)==0) /* return sign(x)*0 */
122 return Zero[(uint64_t)sx>>63];
123 hx = hz+hz+(lz>>63); lx = lz+lz;
124 }
125 }
126 hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
127 if(hz>=0) {hx=hz;lx=lz;}
128
129 /* convert back to floating value and restore the sign */
130 if((hx|lx)==0) /* return sign(x)*0 */
131 return Zero[(uint64_t)sx>>63];
132 while(hx<0x0001000000000000LL) { /* normalize x */
133 hx = hx+hx+(lx>>63); lx = lx+lx;
134 iy -= 1;
135 }
136 if(__builtin_expect(iy>= -1022,0)) { /* normalize output */
137 x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
138 } else { /* subnormal output */
139 n = -1022 - iy;
140 /* We know 1 <= N <= 52, and that there are no nonzero
141 bits in places below 2^-1074. */
142 lx = (lx >> n) | ((uint64_t) hx << (64 - n));
143 hx >>= n;
144 x = ldbl_insert_mantissa((sx>>63), -1023, hx, lx);
145 x *= one; /* create necessary signal */
146 }
147 return x; /* exact output */
148 }
149 libm_alias_finite (__ieee754_fmodl, __fmodl)