1 /* Copyright (C) 2007-2022 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
bid64_nextup(UINT64 * pres,UINT64 * px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)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 }
482