1
2 /*******************************************************************************
3 MIT License
4 -----------
5
6 Copyright (c) 2002-2019 Advanced Micro Devices, Inc.
7
8 Permission is hereby granted, free of charge, to any person obtaining a copy
9 of this Software and associated documentaon files (the "Software"), to deal
10 in the Software without restriction, including without limitation the rights
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 copies of the Software, and to permit persons to whom the Software is
13 furnished to do so, subject to the following conditions:
14
15 The above copyright notice and this permission notice shall be included in
16 all copies or substantial portions of the Software.
17
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
24 THE SOFTWARE.
25 *******************************************************************************/
26
27 #include "libm.h"
28 #include "libm_util.h"
29
30 #define USE_NAN_WITH_FLAGS
31 #define USE_SCALEDOUBLE_3
32 #define USE_GET_FPSW_INLINE
33 #define USE_SET_FPSW_INLINE
34 #define USE_HANDLE_ERROR
35 #include "libm_inlines.h"
36 #undef USE_NAN_WITH_FLAGS
37 #undef USE_SCALEDOUBLE_3
38 #undef USE_GET_FPSW_INLINE
39 #undef USE_SET_FPSW_INLINE
40 #undef USE_HANDLE_ERROR
41
42 #if !defined(_CRTBLD_C9X)
43 #define _CRTBLD_C9X
44 #endif
45
46 #include "libm_errno.h"
47
48 /* Computes the exact product of x and y, the result being the
49 nearly doublelength number (z,zz) */
dekker_mul12(double x,double y,double * z,double * zz)50 static inline void dekker_mul12(double x, double y,
51 double *z, double *zz)
52 {
53 double hx, tx, hy, ty;
54 /* Split x into hx (head) and tx (tail). Do the same for y. */
55 unsigned long long u;
56 GET_BITS_DP64(x, u);
57 u &= 0xfffffffff8000000;
58 PUT_BITS_DP64(u, hx);
59 tx = x - hx;
60 GET_BITS_DP64(y, u);
61 u &= 0xfffffffff8000000;
62 PUT_BITS_DP64(u, hy);
63 ty = y - hy;
64 *z = x * y;
65 *zz = (((hx * hy - *z) + hx * ty) + tx * hy) + tx * ty;
66 }
67
68 #pragma function(fmod)
69 #undef _FUNCNAME
70 #if defined(COMPILING_FMOD)
fmod(double x,double y)71 double fmod(double x, double y)
72 #define _FUNCNAME "fmod"
73 #define _OPERATION OP_FMOD
74 #else
75 double remainder(double x, double y)
76 #define _FUNCNAME "remainder"
77 #define _OPERATION OP_REM
78 #endif
79 {
80 double dx, dy, scale, w, t, v, c, cc;
81 int i, ntimes, xexp, yexp;
82 unsigned long long u, ux, uy, ax, ay, todd;
83 unsigned int sw;
84
85 dx = x;
86 dy = y;
87
88
89 GET_BITS_DP64(dx, ux);
90 GET_BITS_DP64(dy, uy);
91 ax = ux & ~SIGNBIT_DP64;
92 ay = uy & ~SIGNBIT_DP64;
93 xexp = (int)((ux & EXPBITS_DP64) >> EXPSHIFTBITS_DP64);
94 yexp = (int)((uy & EXPBITS_DP64) >> EXPSHIFTBITS_DP64);
95
96 if (xexp < 1 || xexp > BIASEDEMAX_DP64 ||
97 yexp < 1 || yexp > BIASEDEMAX_DP64)
98 {
99 /* x or y is zero, denormalized, NaN or infinity */
100 if (xexp > BIASEDEMAX_DP64)
101 {
102 /* x is NaN or infinity */
103 if (ux & MANTBITS_DP64)
104 {
105 /* x is NaN */
106 return _handle_error(_FUNCNAME, _OPERATION, ux|0x0008000000000000, _DOMAIN, 0,
107 EDOM, x, y, 2);
108 }
109 else
110 {
111 /* x is infinity; result is NaN */
112 return _handle_error(_FUNCNAME, _OPERATION, INDEFBITPATT_DP64, _DOMAIN,
113 AMD_F_INVALID, EDOM, x, y, 2);
114 }
115 }
116 else if (yexp > BIASEDEMAX_DP64)
117 {
118 /* y is NaN or infinity */
119 if (uy & MANTBITS_DP64)
120 {
121 /* y is NaN */
122 return _handle_error(_FUNCNAME, _OPERATION, uy|0x0008000000000000, _DOMAIN, 0,
123 EDOM, x, y, 2);
124 }
125 else
126 {
127 #ifdef _CRTBLD_C9X
128 /* C99 return for y = +-inf is x */
129 return x;
130 #else
131 /* y is infinity; result is indefinite */
132 return _handle_error(_FUNCNAME, _OPERATION, INDEFBITPATT_DP64, _DOMAIN,
133 AMD_F_INVALID, EDOM, x, y, 2);
134 #endif
135 }
136 }
137 else if (ax == 0x0000000000000000)
138 {
139 /* x is zero */
140 if (ay == 0x0000000000000000)
141 {
142 /* y is zero */
143 return _handle_error(_FUNCNAME, _OPERATION, INDEFBITPATT_DP64, _DOMAIN,
144 AMD_F_INVALID, EDOM, x, y, 2);
145 }
146 else
147 /* C99 return for x = 0 must preserve sign */
148 return x;
149 }
150 else if (ay == 0x0000000000000000)
151 {
152 /* y is zero */
153 return _handle_error(_FUNCNAME, _OPERATION, INDEFBITPATT_DP64, _DOMAIN,
154 AMD_F_INVALID, EDOM, x, y, 2);
155 }
156
157 /* We've exhausted all other possibilities. One or both of x and
158 y must be denormalized */
159 if (xexp < 1)
160 {
161 /* x is denormalized. Figure out its exponent. */
162 u = ax;
163 while (u < IMPBIT_DP64)
164 {
165 xexp--;
166 u <<= 1;
167 }
168 }
169 if (yexp < 1)
170 {
171 /* y is denormalized. Figure out its exponent. */
172 u = ay;
173 while (u < IMPBIT_DP64)
174 {
175 yexp--;
176 u <<= 1;
177 }
178 }
179 }
180 else if (ax == ay)
181 {
182 /* abs(x) == abs(y); return zero with the sign of x */
183 PUT_BITS_DP64(ux & SIGNBIT_DP64, dx);
184 return dx;
185 }
186
187 /* Set x = abs(x), y = abs(y) */
188 PUT_BITS_DP64(ax, dx);
189 PUT_BITS_DP64(ay, dy);
190
191 if (ax < ay)
192 {
193 /* abs(x) < abs(y) */
194 #if !defined(COMPILING_FMOD)
195 if (dx > 0.5*dy)
196 dx -= dy;
197 #endif
198 return x < 0.0? -dx : dx;
199 }
200
201 /* Save the current floating-point status word. We need
202 to do this because the remainder function is always
203 exact for finite arguments, but our algorithm causes
204 the inexact flag to be raised. We therefore need to
205 restore the entry status before exiting. */
206 sw = get_fpsw_inline();
207
208 /* Set ntimes to the number of times we need to do a
209 partial remainder. If the exponent of x is an exact multiple
210 of 52 larger than the exponent of y, and the mantissa of x is
211 less than the mantissa of y, ntimes will be one too large
212 but it doesn't matter - it just means that we'll go round
213 the loop below one extra time. */
214 if (xexp <= yexp)
215 ntimes = 0;
216 else
217 ntimes = (xexp - yexp) / 52;
218
219 if (ntimes == 0)
220 {
221 w = dy;
222 scale = 1.0;
223 }
224 else
225 {
226 /* Set w = y * 2^(52*ntimes) */
227 w = scaleDouble_3(dy, ntimes * 52);
228
229 /* Set scale = 2^(-52) */
230 PUT_BITS_DP64((unsigned long long)(-52 + EXPBIAS_DP64) << EXPSHIFTBITS_DP64,
231 scale);
232 }
233
234
235 /* Each time round the loop we compute a partial remainder.
236 This is done by subtracting a large multiple of w
237 from x each time, where w is a scaled up version of y.
238 The subtraction must be performed exactly in quad
239 precision, though the result at each stage can
240 fit exactly in a double precision number. */
241 for (i = 0; i < ntimes; i++)
242 {
243 /* t is the integer multiple of w that we will subtract.
244 We use a truncated value for t.
245
246 N.B. w has been chosen so that the integer t will have
247 at most 52 significant bits. This is the amount by
248 which the exponent of the partial remainder dx gets reduced
249 every time around the loop. In theory we could use
250 53 bits in t, but the quad precision multiplication
251 routine dekker_mul12 does not allow us to do that because
252 it loses the last (106th) bit of its quad precision result. */
253
254 /* Set dx = dx - w * t, where t is equal to trunc(dx/w). */
255 t = (double)(long long)(dx / w);
256 /* At this point, t may be one too large due to
257 rounding of dx/w */
258
259 /* Compute w * t in quad precision */
260 dekker_mul12(w, t, &c, &cc);
261
262 /* Subtract w * t from dx */
263 v = dx - c;
264 dx = v + (((dx - v) - c) - cc);
265
266 /* If t was one too large, dx will be negative. Add back
267 one w */
268 /* It might be possible to speed up this loop by finding
269 a way to compute correctly truncated t directly from dx and w.
270 We would then avoid the need for this check on negative dx. */
271 if (dx < 0.0)
272 dx += w;
273
274 /* Scale w down by 2^(-52) for the next iteration */
275 w *= scale;
276 }
277
278 /* One more time */
279 /* Variable todd says whether the integer t is odd or not */
280 t = (double)(long long)(dx / w);
281 todd = ((long long)(dx / w)) & 1;
282 dekker_mul12(w, t, &c, &cc);
283 v = dx - c;
284 dx = v + (((dx - v) - c) - cc);
285 if (dx < 0.0)
286 {
287 todd = !todd;
288 dx += w;
289 }
290
291 /* At this point, dx lies in the range [0,dy) */
292 #if !defined(COMPILING_FMOD)
293 /* For the fmod function, we're done apart from setting
294 the correct sign. */
295 /* For the remainder function, we need to adjust dx
296 so that it lies in the range (-y/2, y/2] by carefully
297 subtracting w (== dy == y) if necessary. The rigmarole
298 with todd is to get the correct sign of the result
299 when x/y lies exactly half way between two integers,
300 when we need to choose the even integer. */
301 if (ay < 0x7fd0000000000000)
302 {
303 if (dx + dx > w || (todd && (dx + dx == w)))
304 dx -= w;
305 }
306 else if (dx > 0.5 * w || (todd && (dx == 0.5 * w)))
307 dx -= w;
308
309 #endif
310
311 /* **** N.B. for some reason this breaks the 32 bit version
312 of remainder when compiling with optimization. */
313 /* Restore the entry status flags */
314 set_fpsw_inline(sw);
315
316 /* Set the result sign according to input argument x */
317 return x < 0.0? -dx : dx;
318
319 }
320