1 /* Copyright (C) 2007-2023 Free Software Foundation, Inc.
2
3 This file is part of GCC.
4
5 GCC is free software; you can redistribute it and/or modify it under
6 the terms of the GNU General Public License as published by the Free
7 Software Foundation; either version 3, or (at your option) any later
8 version.
9
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 for more details.
14
15 Under Section 7 of GPL version 3, you are granted additional
16 permissions described in the GCC Runtime Library Exception, version
17 3.1, as published by the Free Software Foundation.
18
19 You should have received a copy of the GNU General Public License and
20 a copy of the GCC Runtime Library Exception along with this program;
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22 <http://www.gnu.org/licenses/>. */
23
24 #include "bid_internal.h"
25
26 /*****************************************************************************
27 * BID64 nextup
28 ****************************************************************************/
29
30 #if DECIMAL_CALL_BY_REFERENCE
31 void
32 bid64_nextup (UINT64 * pres,
33 UINT64 *
34 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
35 UINT64 x = *px;
36 #else
37 UINT64
38 bid64_nextup (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
39 _EXC_INFO_PARAM) {
40 #endif
41
42 UINT64 res;
43 UINT64 x_sign;
44 UINT64 x_exp;
45 BID_UI64DOUBLE tmp1;
46 int x_nr_bits;
47 int q1, ind;
48 UINT64 C1; // C1 represents x_signif (UINT64)
49
50 // check for NaNs and infinities
51 if ((x & MASK_NAN) == MASK_NAN) { // check for NaN
52 if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
53 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
54 else
55 x = x & 0xfe03ffffffffffffull; // clear G6-G12
56 if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
57 // set invalid flag
58 *pfpsf |= INVALID_EXCEPTION;
59 // return quiet (SNaN)
60 res = x & 0xfdffffffffffffffull;
61 } else { // QNaN
62 res = x;
63 }
64 BID_RETURN (res);
65 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
66 if (!(x & 0x8000000000000000ull)) { // x is +inf
67 res = 0x7800000000000000ull;
68 } else { // x is -inf
69 res = 0xf7fb86f26fc0ffffull; // -MAXFP = -999...99 * 10^emax
70 }
71 BID_RETURN (res);
72 }
73 // unpack the argument
74 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
75 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
76 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
77 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
78 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
79 if (C1 > 9999999999999999ull) { // non-canonical
80 x_exp = 0;
81 C1 = 0;
82 }
83 } else {
84 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
85 C1 = x & MASK_BINARY_SIG1;
86 }
87
88 // check for zeros (possibly from non-canonical values)
89 if (C1 == 0x0ull) {
90 // x is 0
91 res = 0x0000000000000001ull; // MINFP = 1 * 10^emin
92 } else { // x is not special and is not zero
93 if (x == 0x77fb86f26fc0ffffull) {
94 // x = +MAXFP = 999...99 * 10^emax
95 res = 0x7800000000000000ull; // +inf
96 } else if (x == 0x8000000000000001ull) {
97 // x = -MINFP = 1...99 * 10^emin
98 res = 0x8000000000000000ull; // -0
99 } else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
100 // can add/subtract 1 ulp to the significand
101
102 // Note: we could check here if x >= 10^16 to speed up the case q1 =16
103 // q1 = nr. of decimal digits in x (1 <= q1 <= 54)
104 // determine first the nr. of bits in x
105 if (C1 >= MASK_BINARY_OR2) { // x >= 2^53
106 // split the 64-bit value in two 32-bit halves to avoid rounding errors
107 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
108 tmp1.d = (double) (C1 >> 32); // exact conversion
109 x_nr_bits =
110 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
111 } else { // x < 2^32
112 tmp1.d = (double) C1; // exact conversion
113 x_nr_bits =
114 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
115 }
116 } else { // if x < 2^53
117 tmp1.d = (double) C1; // exact conversion
118 x_nr_bits =
119 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
120 }
121 q1 = nr_digits[x_nr_bits - 1].digits;
122 if (q1 == 0) {
123 q1 = nr_digits[x_nr_bits - 1].digits1;
124 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
125 q1++;
126 }
127 // if q1 < P16 then pad the significand with zeros
128 if (q1 < P16) {
129 if (x_exp > (UINT64) (P16 - q1)) {
130 ind = P16 - q1; // 1 <= ind <= P16 - 1
131 // pad with P16 - q1 zeros, until exponent = emin
132 // C1 = C1 * 10^ind
133 C1 = C1 * ten2k64[ind];
134 x_exp = x_exp - ind;
135 } else { // pad with zeros until the exponent reaches emin
136 ind = x_exp;
137 C1 = C1 * ten2k64[ind];
138 x_exp = EXP_MIN;
139 }
140 }
141 if (!x_sign) { // x > 0
142 // add 1 ulp (add 1 to the significand)
143 C1++;
144 if (C1 == 0x002386f26fc10000ull) { // if C1 = 10^16
145 C1 = 0x00038d7ea4c68000ull; // C1 = 10^15
146 x_exp++;
147 }
148 // Ok, because MAXFP = 999...99 * 10^emax was caught already
149 } else { // x < 0
150 // subtract 1 ulp (subtract 1 from the significand)
151 C1--;
152 if (C1 == 0x00038d7ea4c67fffull && x_exp != 0) { // if C1 = 10^15 - 1
153 C1 = 0x002386f26fc0ffffull; // C1 = 10^16 - 1
154 x_exp--;
155 }
156 }
157 // assemble the result
158 // if significand has 54 bits
159 if (C1 & MASK_BINARY_OR2) {
160 res =
161 x_sign | (x_exp << 51) | MASK_STEERING_BITS | (C1 &
162 MASK_BINARY_SIG2);
163 } else { // significand fits in 53 bits
164 res = x_sign | (x_exp << 53) | C1;
165 }
166 } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
167 } // end x is not special and is not zero
168 BID_RETURN (res);
169 }
170
171 /*****************************************************************************
172 * BID64 nextdown
173 ****************************************************************************/
174
175 #if DECIMAL_CALL_BY_REFERENCE
176 void
177 bid64_nextdown (UINT64 * pres,
178 UINT64 *
179 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
180 UINT64 x = *px;
181 #else
182 UINT64
183 bid64_nextdown (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
184 _EXC_INFO_PARAM) {
185 #endif
186
187 UINT64 res;
188 UINT64 x_sign;
189 UINT64 x_exp;
190 BID_UI64DOUBLE tmp1;
191 int x_nr_bits;
192 int q1, ind;
193 UINT64 C1; // C1 represents x_signif (UINT64)
194
195 // check for NaNs and infinities
196 if ((x & MASK_NAN) == MASK_NAN) { // check for NaN
197 if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
198 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
199 else
200 x = x & 0xfe03ffffffffffffull; // clear G6-G12
201 if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
202 // set invalid flag
203 *pfpsf |= INVALID_EXCEPTION;
204 // return quiet (SNaN)
205 res = x & 0xfdffffffffffffffull;
206 } else { // QNaN
207 res = x;
208 }
209 BID_RETURN (res);
210 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
211 if (x & 0x8000000000000000ull) { // x is -inf
212 res = 0xf800000000000000ull;
213 } else { // x is +inf
214 res = 0x77fb86f26fc0ffffull; // +MAXFP = +999...99 * 10^emax
215 }
216 BID_RETURN (res);
217 }
218 // unpack the argument
219 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
220 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
221 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
222 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
223 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
224 if (C1 > 9999999999999999ull) { // non-canonical
225 x_exp = 0;
226 C1 = 0;
227 }
228 } else {
229 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
230 C1 = x & MASK_BINARY_SIG1;
231 }
232
233 // check for zeros (possibly from non-canonical values)
234 if (C1 == 0x0ull) {
235 // x is 0
236 res = 0x8000000000000001ull; // -MINFP = -1 * 10^emin
237 } else { // x is not special and is not zero
238 if (x == 0xf7fb86f26fc0ffffull) {
239 // x = -MAXFP = -999...99 * 10^emax
240 res = 0xf800000000000000ull; // -inf
241 } else if (x == 0x0000000000000001ull) {
242 // x = +MINFP = 1...99 * 10^emin
243 res = 0x0000000000000000ull; // -0
244 } else { // -MAXFP + 1ulp <= x <= -MINFP OR MINFP + 1 ulp <= x <= MAXFP
245 // can add/subtract 1 ulp to the significand
246
247 // Note: we could check here if x >= 10^16 to speed up the case q1 =16
248 // q1 = nr. of decimal digits in x (1 <= q1 <= 16)
249 // determine first the nr. of bits in x
250 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
251 // split the 64-bit value in two 32-bit halves to avoid
252 // rounding errors
253 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
254 tmp1.d = (double) (C1 >> 32); // exact conversion
255 x_nr_bits =
256 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
257 } else { // x < 2^32
258 tmp1.d = (double) C1; // exact conversion
259 x_nr_bits =
260 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
261 }
262 } else { // if x < 2^53
263 tmp1.d = (double) C1; // exact conversion
264 x_nr_bits =
265 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
266 }
267 q1 = nr_digits[x_nr_bits - 1].digits;
268 if (q1 == 0) {
269 q1 = nr_digits[x_nr_bits - 1].digits1;
270 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
271 q1++;
272 }
273 // if q1 < P16 then pad the significand with zeros
274 if (q1 < P16) {
275 if (x_exp > (UINT64) (P16 - q1)) {
276 ind = P16 - q1; // 1 <= ind <= P16 - 1
277 // pad with P16 - q1 zeros, until exponent = emin
278 // C1 = C1 * 10^ind
279 C1 = C1 * ten2k64[ind];
280 x_exp = x_exp - ind;
281 } else { // pad with zeros until the exponent reaches emin
282 ind = x_exp;
283 C1 = C1 * ten2k64[ind];
284 x_exp = EXP_MIN;
285 }
286 }
287 if (x_sign) { // x < 0
288 // add 1 ulp (add 1 to the significand)
289 C1++;
290 if (C1 == 0x002386f26fc10000ull) { // if C1 = 10^16
291 C1 = 0x00038d7ea4c68000ull; // C1 = 10^15
292 x_exp++;
293 // Ok, because -MAXFP = -999...99 * 10^emax was caught already
294 }
295 } else { // x > 0
296 // subtract 1 ulp (subtract 1 from the significand)
297 C1--;
298 if (C1 == 0x00038d7ea4c67fffull && x_exp != 0) { // if C1 = 10^15 - 1
299 C1 = 0x002386f26fc0ffffull; // C1 = 10^16 - 1
300 x_exp--;
301 }
302 }
303 // assemble the result
304 // if significand has 54 bits
305 if (C1 & MASK_BINARY_OR2) {
306 res =
307 x_sign | (x_exp << 51) | MASK_STEERING_BITS | (C1 &
308 MASK_BINARY_SIG2);
309 } else { // significand fits in 53 bits
310 res = x_sign | (x_exp << 53) | C1;
311 }
312 } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
313 } // end x is not special and is not zero
314 BID_RETURN (res);
315 }
316
317 /*****************************************************************************
318 * BID64 nextafter
319 ****************************************************************************/
320
321 #if DECIMAL_CALL_BY_REFERENCE
322 void
323 bid64_nextafter (UINT64 * pres, UINT64 * px,
324 UINT64 *
325 py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
326 UINT64 x = *px;
327 UINT64 y = *py;
328 #else
329 UINT64
330 bid64_nextafter (UINT64 x,
331 UINT64 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
332 _EXC_INFO_PARAM) {
333 #endif
334
335 UINT64 res;
336 UINT64 tmp1, tmp2;
337 FPSC tmp_fpsf = 0; // dummy fpsf for calls to comparison functions
338 int res1, res2;
339
340 // check for NaNs or infinities
341 if (((x & MASK_SPECIAL) == MASK_SPECIAL) ||
342 ((y & MASK_SPECIAL) == MASK_SPECIAL)) {
343 // x is NaN or infinity or y is NaN or infinity
344
345 if ((x & MASK_NAN) == MASK_NAN) { // x is NAN
346 if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
347 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
348 else
349 x = x & 0xfe03ffffffffffffull; // clear G6-G12
350 if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNAN
351 // set invalid flag
352 *pfpsf |= INVALID_EXCEPTION;
353 // return quiet (x)
354 res = x & 0xfdffffffffffffffull;
355 } else { // x is QNaN
356 if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN
357 // set invalid flag
358 *pfpsf |= INVALID_EXCEPTION;
359 }
360 // return x
361 res = x;
362 }
363 BID_RETURN (res);
364 } else if ((y & MASK_NAN) == MASK_NAN) { // y is NAN
365 if ((y & 0x0003ffffffffffffull) > 999999999999999ull)
366 y = y & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
367 else
368 y = y & 0xfe03ffffffffffffull; // clear G6-G12
369 if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN
370 // set invalid flag
371 *pfpsf |= INVALID_EXCEPTION;
372 // return quiet (y)
373 res = y & 0xfdffffffffffffffull;
374 } else { // y is QNaN
375 // return y
376 res = y;
377 }
378 BID_RETURN (res);
379 } else { // at least one is infinity
380 if ((x & MASK_ANY_INF) == MASK_INF) { // x = inf
381 x = x & (MASK_SIGN | MASK_INF);
382 }
383 if ((y & MASK_ANY_INF) == MASK_INF) { // y = inf
384 y = y & (MASK_SIGN | MASK_INF);
385 }
386 }
387 }
388 // neither x nor y is NaN
389
390 // if not infinity, check for non-canonical values x (treated as zero)
391 if ((x & MASK_ANY_INF) != MASK_INF) { // x != inf
392 // unpack x
393 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
394 // if the steering bits are 11 (condition will be 0), then
395 // the exponent is G[0:w+1]
396 if (((x & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
397 9999999999999999ull) {
398 // non-canonical
399 x = (x & MASK_SIGN) | ((x & MASK_BINARY_EXPONENT2) << 2);
400 }
401 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) x is unch.
402 ; // canonical
403 }
404 }
405 // no need to check for non-canonical y
406
407 // neither x nor y is NaN
408 tmp_fpsf = *pfpsf; // save fpsf
409 #if DECIMAL_CALL_BY_REFERENCE
410 bid64_quiet_equal (&res1, px,
411 py _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
412 bid64_quiet_greater (&res2, px,
413 py _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
414 #else
415 res1 =
416 bid64_quiet_equal (x,
417 y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
418 res2 =
419 bid64_quiet_greater (x,
420 y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
421 #endif
422 *pfpsf = tmp_fpsf; // restore fpsf
423 if (res1) { // x = y
424 // return x with the sign of y
425 res = (y & 0x8000000000000000ull) | (x & 0x7fffffffffffffffull);
426 } else if (res2) { // x > y
427 #if DECIMAL_CALL_BY_REFERENCE
428 bid64_nextdown (&res,
429 px _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
430 #else
431 res =
432 bid64_nextdown (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
433 #endif
434 } else { // x < y
435 #if DECIMAL_CALL_BY_REFERENCE
436 bid64_nextup (&res, px _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
437 #else
438 res = bid64_nextup (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
439 #endif
440 }
441 // if the operand x is finite but the result is infinite, signal
442 // overflow and inexact
443 if (((x & MASK_INF) != MASK_INF) && ((res & MASK_INF) == MASK_INF)) {
444 // set the inexact flag
445 *pfpsf |= INEXACT_EXCEPTION;
446 // set the overflow flag
447 *pfpsf |= OVERFLOW_EXCEPTION;
448 }
449 // if the result is in (-10^emin, 10^emin), and is different from the
450 // operand x, signal underflow and inexact
451 tmp1 = 0x00038d7ea4c68000ull; // +100...0[16] * 10^emin
452 tmp2 = res & 0x7fffffffffffffffull;
453 tmp_fpsf = *pfpsf; // save fpsf
454 #if DECIMAL_CALL_BY_REFERENCE
455 bid64_quiet_greater (&res1, &tmp1,
456 &tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
457 _EXC_INFO_ARG);
458 bid64_quiet_not_equal (&res2, &x,
459 &res _EXC_FLAGS_ARG _EXC_MASKS_ARG
460 _EXC_INFO_ARG);
461 #else
462 res1 =
463 bid64_quiet_greater (tmp1,
464 tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
465 _EXC_INFO_ARG);
466 res2 =
467 bid64_quiet_not_equal (x,
468 res _EXC_FLAGS_ARG _EXC_MASKS_ARG
469 _EXC_INFO_ARG);
470 #endif
471 *pfpsf = tmp_fpsf; // restore fpsf
472 if (res1 && res2) {
473 // if (bid64_quiet_greater (tmp1, tmp2, &tmp_fpsf) &&
474 // bid64_quiet_not_equal (x, res, &tmp_fpsf)) {
475 // set the inexact flag
476 *pfpsf |= INEXACT_EXCEPTION;
477 // set the underflow flag
478 *pfpsf |= UNDERFLOW_EXCEPTION;
479 }
480 BID_RETURN (res);
481 }