libgcc/config/libbid/bid64_next.c

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