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) |
150 | |