]>
Commit | Line | Data |
---|---|---|
f7eac6eb RM |
1 | /* |
2 | * ==================================================== | |
3 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | |
4 | * | |
5 | * Developed at SunPro, a Sun Microsystems, Inc. business. | |
6 | * Permission to use, copy, modify, and distribute this | |
0ac5ae23 | 7 | * software is freely granted, provided that this notice |
f7eac6eb RM |
8 | * is preserved. |
9 | * ==================================================== | |
10 | */ | |
11 | ||
0ac5ae23 | 12 | /* |
f7eac6eb RM |
13 | * __ieee754_fmod(x,y) |
14 | * Return x mod y in exact arithmetic | |
15 | * Method: shift and subtract | |
16 | */ | |
17 | ||
1ed0291c RH |
18 | #include <math.h> |
19 | #include <math_private.h> | |
f7eac6eb | 20 | |
c5d5d574 | 21 | static const double one = 1.0, Zero[] = { 0.0, -0.0, }; |
f7eac6eb | 22 | |
0ac5ae23 UD |
23 | double |
24 | __ieee754_fmod (double x, double y) | |
f7eac6eb | 25 | { |
c5d5d574 OB |
26 | int32_t n, hx, hy, hz, ix, iy, sx, i; |
27 | u_int32_t lx, ly, lz; | |
f7eac6eb | 28 | |
c5d5d574 OB |
29 | EXTRACT_WORDS (hx, lx, x); |
30 | EXTRACT_WORDS (hy, ly, y); | |
31 | sx = hx & 0x80000000; /* sign of x */ | |
32 | hx ^= sx; /* |x| */ | |
33 | hy &= 0x7fffffff; /* |y| */ | |
f7eac6eb | 34 | |
c5d5d574 OB |
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); | |
39 | if (hx <= hy) | |
40 | { | |
41 | if ((hx < hy) || (lx < ly)) | |
42 | return x; /* |x|<|y| return x */ | |
43 | if (lx == ly) | |
44 | return Zero[(u_int32_t) sx >> 31]; /* |x|=|y| return x*0*/ | |
45 | } | |
f7eac6eb | 46 | |
c5d5d574 | 47 | /* determine ix = ilogb(x) */ |
a1ffb40e | 48 | if (__glibc_unlikely (hx < 0x00100000)) /* subnormal x */ |
c5d5d574 OB |
49 | { |
50 | if (hx == 0) | |
51 | { | |
52 | for (ix = -1043, i = lx; i > 0; i <<= 1) | |
53 | ix -= 1; | |
54 | } | |
55 | else | |
56 | { | |
57 | for (ix = -1022, i = (hx << 11); i > 0; i <<= 1) | |
58 | ix -= 1; | |
59 | } | |
60 | } | |
61 | else | |
62 | ix = (hx >> 20) - 1023; | |
f7eac6eb | 63 | |
c5d5d574 | 64 | /* determine iy = ilogb(y) */ |
a1ffb40e | 65 | if (__glibc_unlikely (hy < 0x00100000)) /* subnormal y */ |
c5d5d574 OB |
66 | { |
67 | if (hy == 0) | |
68 | { | |
69 | for (iy = -1043, i = ly; i > 0; i <<= 1) | |
70 | iy -= 1; | |
71 | } | |
72 | else | |
73 | { | |
74 | for (iy = -1022, i = (hy << 11); i > 0; i <<= 1) | |
75 | iy -= 1; | |
76 | } | |
77 | } | |
78 | else | |
79 | iy = (hy >> 20) - 1023; | |
f7eac6eb | 80 | |
c5d5d574 | 81 | /* set up {hx,lx}, {hy,ly} and align y to x */ |
a1ffb40e | 82 | if (__glibc_likely (ix >= -1022)) |
c5d5d574 OB |
83 | hx = 0x00100000 | (0x000fffff & hx); |
84 | else /* subnormal x, shift x to normal */ | |
85 | { | |
86 | n = -1022 - ix; | |
87 | if (n <= 31) | |
88 | { | |
89 | hx = (hx << n) | (lx >> (32 - n)); | |
90 | lx <<= n; | |
91 | } | |
92 | else | |
93 | { | |
94 | hx = lx << (n - 32); | |
95 | lx = 0; | |
96 | } | |
97 | } | |
a1ffb40e | 98 | if (__glibc_likely (iy >= -1022)) |
c5d5d574 OB |
99 | hy = 0x00100000 | (0x000fffff & hy); |
100 | else /* subnormal y, shift y to normal */ | |
101 | { | |
102 | n = -1022 - iy; | |
103 | if (n <= 31) | |
104 | { | |
105 | hy = (hy << n) | (ly >> (32 - n)); | |
106 | ly <<= n; | |
f7eac6eb | 107 | } |
c5d5d574 OB |
108 | else |
109 | { | |
110 | hy = ly << (n - 32); | |
111 | ly = 0; | |
f7eac6eb | 112 | } |
c5d5d574 | 113 | } |
f7eac6eb | 114 | |
c5d5d574 OB |
115 | /* fix point fmod */ |
116 | n = ix - iy; | |
117 | while (n--) | |
118 | { | |
119 | hz = hx - hy; lz = lx - ly; if (lx < ly) | |
120 | hz -= 1; | |
121 | if (hz < 0) | |
122 | { | |
123 | hx = hx + hx + (lx >> 31); lx = lx + lx; | |
f7eac6eb | 124 | } |
c5d5d574 OB |
125 | else |
126 | { | |
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; | |
130 | } | |
131 | } | |
132 | hz = hx - hy; lz = lx - ly; if (lx < ly) | |
133 | hz -= 1; | |
134 | if (hz >= 0) | |
135 | { | |
136 | hx = hz; lx = lz; | |
137 | } | |
f7eac6eb | 138 | |
c5d5d574 OB |
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 */ | |
143 | { | |
144 | hx = hx + hx + (lx >> 31); lx = lx + lx; | |
145 | iy -= 1; | |
146 | } | |
a1ffb40e | 147 | if (__glibc_likely (iy >= -1022)) /* normalize output */ |
c5d5d574 OB |
148 | { |
149 | hx = ((hx - 0x00100000) | ((iy + 1023) << 20)); | |
150 | INSERT_WORDS (x, hx | sx, lx); | |
151 | } | |
152 | else /* subnormal output */ | |
153 | { | |
154 | n = -1022 - iy; | |
155 | if (n <= 20) | |
156 | { | |
157 | lx = (lx >> n) | ((u_int32_t) hx << (32 - n)); | |
158 | hx >>= n; | |
159 | } | |
160 | else if (n <= 31) | |
161 | { | |
162 | lx = (hx << (32 - n)) | (lx >> n); hx = sx; | |
f7eac6eb | 163 | } |
c5d5d574 OB |
164 | else |
165 | { | |
166 | lx = hx >> (n - 32); hx = sx; | |
f7eac6eb | 167 | } |
c5d5d574 OB |
168 | INSERT_WORDS (x, hx | sx, lx); |
169 | x *= one; /* create necessary signal */ | |
170 | } | |
171 | return x; /* exact output */ | |
f7eac6eb | 172 | } |
0ac5ae23 | 173 | strong_alias (__ieee754_fmod, __fmod_finite) |