libgcc/config/libbid/bid64_round_integral.c

   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_round_integral_exact
  28  ****************************************************************************/
  29
  30 #if DECIMAL_CALL_BY_REFERENCE
  31 void
  32 bid64_round_integral_exact (UINT64 * pres,
  33                             UINT64 *
  34                             px _RND_MODE_PARAM _EXC_FLAGS_PARAM
  35                             _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  36   UINT64 x = *px;
  37 #if !DECIMAL_GLOBAL_ROUNDING
  38   unsigned int rnd_mode = *prnd_mode;
  39 #endif
  40 #else
  41 UINT64
  42 bid64_round_integral_exact (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
  43                             _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  44 #endif
  45
  46   UINT64 res = 0xbaddbaddbaddbaddull;
  47   UINT64 x_sign;
  48   int exp;                      // unbiased exponent
  49   // Note: C1 represents the significand (UINT64)
  50   BID_UI64DOUBLE tmp1;
  51   int x_nr_bits;
  52   int q, ind, shift;
  53   UINT64 C1;
  54   // UINT64 res is C* at first - represents up to 16 decimal digits <= 54 bits
  55   UINT128 fstar = { {0x0ull, 0x0ull} };
  56   UINT128 P128;
  57
  58   x_sign = x & MASK_SIGN;       // 0 for positive, MASK_SIGN for negative
  59
  60   // check for NaNs and infinities
  61   if ((x & MASK_NAN) == MASK_NAN) {     // check for NaN
  62     if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
  63       x = x & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
  64     else
  65       x = x & 0xfe03ffffffffffffull;    // clear G6-G12
  66     if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
  67       // set invalid flag
  68       *pfpsf |= INVALID_EXCEPTION;
  69       // return quiet (SNaN)
  70       res = x & 0xfdffffffffffffffull;
  71     } else {    // QNaN
  72       res = x;
  73     }
  74     BID_RETURN (res);
  75   } else if ((x & MASK_INF) == MASK_INF) {      // check for Infinity
  76     res = x_sign | 0x7800000000000000ull;
  77     BID_RETURN (res);
  78   }
  79   // unpack x
  80   if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
  81     // if the steering bits are 11 (condition will be 0), then
  82     // the exponent is G[0:w+1]
  83     exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398;
  84     C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
  85     if (C1 > 9999999999999999ull) {     // non-canonical
  86       C1 = 0;
  87     }
  88   } else {      // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
  89     exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398;
  90     C1 = (x & MASK_BINARY_SIG1);
  91   }
  92
  93   // if x is 0 or non-canonical return 0 preserving the sign bit and
  94   // the preferred exponent of MAX(Q(x), 0)
  95   if (C1 == 0) {
  96     if (exp < 0)
  97       exp = 0;
  98     res = x_sign | (((UINT64) exp + 398) << 53);
  99     BID_RETURN (res);
 100   }
 101   // x is a finite non-zero number (not 0, non-canonical, or special)
 102
 103   switch (rnd_mode) {
 104   case ROUNDING_TO_NEAREST:
 105   case ROUNDING_TIES_AWAY:
 106     // return 0 if (exp <= -(p+1))
 107     if (exp <= -17) {
 108       res = x_sign | 0x31c0000000000000ull;
 109       *pfpsf |= INEXACT_EXCEPTION;
 110       BID_RETURN (res);
 111     }
 112     break;
 113   case ROUNDING_DOWN:
 114     // return 0 if (exp <= -p)
 115     if (exp <= -16) {
 116       if (x_sign) {
 117         res = 0xb1c0000000000001ull;
 118       } else {
 119         res = 0x31c0000000000000ull;
 120       }
 121       *pfpsf |= INEXACT_EXCEPTION;
 122       BID_RETURN (res);
 123     }
 124     break;
 125   case ROUNDING_UP:
 126     // return 0 if (exp <= -p)
 127     if (exp <= -16) {
 128       if (x_sign) {
 129         res = 0xb1c0000000000000ull;
 130       } else {
 131         res = 0x31c0000000000001ull;
 132       }
 133       *pfpsf |= INEXACT_EXCEPTION;
 134       BID_RETURN (res);
 135     }
 136     break;
 137   case ROUNDING_TO_ZERO:
 138     // return 0 if (exp <= -p)
 139     if (exp <= -16) {
 140       res = x_sign | 0x31c0000000000000ull;
 141       *pfpsf |= INEXACT_EXCEPTION;
 142       BID_RETURN (res);
 143     }
 144     break;
 145   }     // end switch ()
 146
 147   // q = nr. of decimal digits in x (1 <= q <= 54)
 148   //  determine first the nr. of bits in x
 149   if (C1 >= 0x0020000000000000ull) {    // x >= 2^53
 150     q = 16;
 151   } else {      // if x < 2^53
 152     tmp1.d = (double) C1;       // exact conversion
 153     x_nr_bits =
 154       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 155     q = nr_digits[x_nr_bits - 1].digits;
 156     if (q == 0) {
 157       q = nr_digits[x_nr_bits - 1].digits1;
 158       if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
 159         q++;
 160     }
 161   }
 162
 163   if (exp >= 0) {       // -exp <= 0
 164     // the argument is an integer already
 165     res = x;
 166     BID_RETURN (res);
 167   }
 168
 169   switch (rnd_mode) {
 170   case ROUNDING_TO_NEAREST:
 171     if ((q + exp) >= 0) {       // exp < 0 and 1 <= -exp <= q
 172       // need to shift right -exp digits from the coefficient; exp will be 0
 173       ind = -exp;       // 1 <= ind <= 16; ind is a synonym for 'x'
 174       // chop off ind digits from the lower part of C1
 175       // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits
 176       // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
 177       C1 = C1 + midpoint64[ind - 1];
 178       // calculate C* and f*
 179       // C* is actually floor(C*) in this case
 180       // C* and f* need shifting and masking, as shown by
 181       // shiftright128[] and maskhigh128[]
 182       // 1 <= x <= 16
 183       // kx = 10^(-x) = ten2mk64[ind - 1]
 184       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 185       // the approximation of 10^(-x) was rounded up to 64 bits
 186       __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
 187
 188       // if (0 < f* < 10^(-x)) then the result is a midpoint
 189       //   if floor(C*) is even then C* = floor(C*) - logical right
 190       //       shift; C* has p decimal digits, correct by Prop. 1)
 191       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
 192       //       shift; C* has p decimal digits, correct by Pr. 1)
 193       // else
 194       //   C* = floor(C*) (logical right shift; C has p decimal digits,
 195       //       correct by Property 1)
 196       // n = C* * 10^(e+x)
 197
 198       if (ind - 1 <= 2) {       // 0 <= ind - 1 <= 2 => shift = 0
 199         res = P128.w[1];
 200         fstar.w[1] = 0;
 201         fstar.w[0] = P128.w[0];
 202       } else if (ind - 1 <= 21) {       // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 203         shift = shiftright128[ind - 1]; // 3 <= shift <= 63
 204         res = (P128.w[1] >> shift);
 205         fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
 206         fstar.w[0] = P128.w[0];
 207       }
 208       // if (0 < f* < 10^(-x)) then the result is a midpoint
 209       // since round_to_even, subtract 1 if current result is odd
 210       if ((res & 0x0000000000000001ull) && (fstar.w[1] == 0)
 211           && (fstar.w[0] < ten2mk64[ind - 1])) {
 212         res--;
 213       }
 214       // determine inexactness of the rounding of C*
 215       // if (0 < f* - 1/2 < 10^(-x)) then
 216       //   the result is exact
 217       // else // if (f* - 1/2 > T*) then
 218       //   the result is inexact
 219       if (ind - 1 <= 2) {
 220         if (fstar.w[0] > 0x8000000000000000ull) {
 221           // f* > 1/2 and the result may be exact
 222           // fstar.w[0] - 0x8000000000000000ull is f* - 1/2
 223           if ((fstar.w[0] - 0x8000000000000000ull) > ten2mk64[ind - 1]) {
 224             // set the inexact flag
 225             *pfpsf |= INEXACT_EXCEPTION;
 226           }     // else the result is exact
 227         } else {        // the result is inexact; f2* <= 1/2
 228           // set the inexact flag
 229           *pfpsf |= INEXACT_EXCEPTION;
 230         }
 231       } else {  // if 3 <= ind - 1 <= 21
 232         if (fstar.w[1] > onehalf128[ind - 1] ||
 233             (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) {
 234           // f2* > 1/2 and the result may be exact
 235           // Calculate f2* - 1/2
 236           if (fstar.w[1] > onehalf128[ind - 1]
 237               || fstar.w[0] > ten2mk64[ind - 1]) {
 238             // set the inexact flag
 239             *pfpsf |= INEXACT_EXCEPTION;
 240           }     // else the result is exact
 241         } else {        // the result is inexact; f2* <= 1/2
 242           // set the inexact flag
 243           *pfpsf |= INEXACT_EXCEPTION;
 244         }
 245       }
 246       // set exponent to zero as it was negative before.
 247       res = x_sign | 0x31c0000000000000ull | res;
 248       BID_RETURN (res);
 249     } else {    // if exp < 0 and q + exp < 0
 250       // the result is +0 or -0
 251       res = x_sign | 0x31c0000000000000ull;
 252       *pfpsf |= INEXACT_EXCEPTION;
 253       BID_RETURN (res);
 254     }
 255     break;
 256   case ROUNDING_TIES_AWAY:
 257     if ((q + exp) >= 0) {       // exp < 0 and 1 <= -exp <= q
 258       // need to shift right -exp digits from the coefficient; exp will be 0
 259       ind = -exp;       // 1 <= ind <= 16; ind is a synonym for 'x'
 260       // chop off ind digits from the lower part of C1
 261       // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits
 262       // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
 263       C1 = C1 + midpoint64[ind - 1];
 264       // calculate C* and f*
 265       // C* is actually floor(C*) in this case
 266       // C* and f* need shifting and masking, as shown by
 267       // shiftright128[] and maskhigh128[]
 268       // 1 <= x <= 16
 269       // kx = 10^(-x) = ten2mk64[ind - 1]
 270       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 271       // the approximation of 10^(-x) was rounded up to 64 bits
 272       __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
 273
 274       // if (0 < f* < 10^(-x)) then the result is a midpoint
 275       //   C* = floor(C*) - logical right shift; C* has p decimal digits,
 276       //       correct by Prop. 1)
 277       // else
 278       //   C* = floor(C*) (logical right shift; C has p decimal digits,
 279       //       correct by Property 1)
 280       // n = C* * 10^(e+x)
 281
 282       if (ind - 1 <= 2) {       // 0 <= ind - 1 <= 2 => shift = 0
 283         res = P128.w[1];
 284         fstar.w[1] = 0;
 285         fstar.w[0] = P128.w[0];
 286       } else if (ind - 1 <= 21) {       // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 287         shift = shiftright128[ind - 1]; // 3 <= shift <= 63
 288         res = (P128.w[1] >> shift);
 289         fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
 290         fstar.w[0] = P128.w[0];
 291       }
 292       // midpoints are already rounded correctly
 293       // determine inexactness of the rounding of C*
 294       // if (0 < f* - 1/2 < 10^(-x)) then
 295       //   the result is exact
 296       // else // if (f* - 1/2 > T*) then
 297       //   the result is inexact
 298       if (ind - 1 <= 2) {
 299         if (fstar.w[0] > 0x8000000000000000ull) {
 300           // f* > 1/2 and the result may be exact
 301           // fstar.w[0] - 0x8000000000000000ull is f* - 1/2
 302           if ((fstar.w[0] - 0x8000000000000000ull) > ten2mk64[ind - 1]) {
 303             // set the inexact flag
 304             *pfpsf |= INEXACT_EXCEPTION;
 305           }     // else the result is exact
 306         } else {        // the result is inexact; f2* <= 1/2
 307           // set the inexact flag
 308           *pfpsf |= INEXACT_EXCEPTION;
 309         }
 310       } else {  // if 3 <= ind - 1 <= 21
 311         if (fstar.w[1] > onehalf128[ind - 1] ||
 312             (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) {
 313           // f2* > 1/2 and the result may be exact
 314           // Calculate f2* - 1/2
 315           if (fstar.w[1] > onehalf128[ind - 1]
 316               || fstar.w[0] > ten2mk64[ind - 1]) {
 317             // set the inexact flag
 318             *pfpsf |= INEXACT_EXCEPTION;
 319           }     // else the result is exact
 320         } else {        // the result is inexact; f2* <= 1/2
 321           // set the inexact flag
 322           *pfpsf |= INEXACT_EXCEPTION;
 323         }
 324       }
 325       // set exponent to zero as it was negative before.
 326       res = x_sign | 0x31c0000000000000ull | res;
 327       BID_RETURN (res);
 328     } else {    // if exp < 0 and q + exp < 0
 329       // the result is +0 or -0
 330       res = x_sign | 0x31c0000000000000ull;
 331       *pfpsf |= INEXACT_EXCEPTION;
 332       BID_RETURN (res);
 333     }
 334     break;
 335   case ROUNDING_DOWN:
 336     if ((q + exp) > 0) {        // exp < 0 and 1 <= -exp < q
 337       // need to shift right -exp digits from the coefficient; exp will be 0
 338       ind = -exp;       // 1 <= ind <= 16; ind is a synonym for 'x'
 339       // chop off ind digits from the lower part of C1
 340       // C1 fits in 64 bits
 341       // calculate C* and f*
 342       // C* is actually floor(C*) in this case
 343       // C* and f* need shifting and masking, as shown by
 344       // shiftright128[] and maskhigh128[]
 345       // 1 <= x <= 16
 346       // kx = 10^(-x) = ten2mk64[ind - 1]
 347       // C* = C1 * 10^(-x)
 348       // the approximation of 10^(-x) was rounded up to 64 bits
 349       __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
 350
 351       // C* = floor(C*) (logical right shift; C has p decimal digits,
 352       //       correct by Property 1)
 353       // if (0 < f* < 10^(-x)) then the result is exact
 354       // n = C* * 10^(e+x)
 355
 356       if (ind - 1 <= 2) {       // 0 <= ind - 1 <= 2 => shift = 0
 357         res = P128.w[1];
 358         fstar.w[1] = 0;
 359         fstar.w[0] = P128.w[0];
 360       } else if (ind - 1 <= 21) {       // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 361         shift = shiftright128[ind - 1]; // 3 <= shift <= 63
 362         res = (P128.w[1] >> shift);
 363         fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
 364         fstar.w[0] = P128.w[0];
 365       }
 366       // if (f* > 10^(-x)) then the result is inexact
 367       if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) {
 368         if (x_sign) {
 369           // if negative and not exact, increment magnitude
 370           res++;
 371         }
 372         *pfpsf |= INEXACT_EXCEPTION;
 373       }
 374       // set exponent to zero as it was negative before.
 375       res = x_sign | 0x31c0000000000000ull | res;
 376       BID_RETURN (res);
 377     } else {    // if exp < 0 and q + exp <= 0
 378       // the result is +0 or -1
 379       if (x_sign) {
 380         res = 0xb1c0000000000001ull;
 381       } else {
 382         res = 0x31c0000000000000ull;
 383       }
 384       *pfpsf |= INEXACT_EXCEPTION;
 385       BID_RETURN (res);
 386     }
 387     break;
 388   case ROUNDING_UP:
 389     if ((q + exp) > 0) {        // exp < 0 and 1 <= -exp < q
 390       // need to shift right -exp digits from the coefficient; exp will be 0
 391       ind = -exp;       // 1 <= ind <= 16; ind is a synonym for 'x'
 392       // chop off ind digits from the lower part of C1
 393       // C1 fits in 64 bits
 394       // calculate C* and f*
 395       // C* is actually floor(C*) in this case
 396       // C* and f* need shifting and masking, as shown by
 397       // shiftright128[] and maskhigh128[]
 398       // 1 <= x <= 16
 399       // kx = 10^(-x) = ten2mk64[ind - 1]
 400       // C* = C1 * 10^(-x)
 401       // the approximation of 10^(-x) was rounded up to 64 bits
 402       __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
 403
 404       // C* = floor(C*) (logical right shift; C has p decimal digits,
 405       //       correct by Property 1)
 406       // if (0 < f* < 10^(-x)) then the result is exact
 407       // n = C* * 10^(e+x)
 408
 409       if (ind - 1 <= 2) {       // 0 <= ind - 1 <= 2 => shift = 0
 410         res = P128.w[1];
 411         fstar.w[1] = 0;
 412         fstar.w[0] = P128.w[0];
 413       } else if (ind - 1 <= 21) {       // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 414         shift = shiftright128[ind - 1]; // 3 <= shift <= 63
 415         res = (P128.w[1] >> shift);
 416         fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
 417         fstar.w[0] = P128.w[0];
 418       }
 419       // if (f* > 10^(-x)) then the result is inexact
 420       if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) {
 421         if (!x_sign) {
 422           // if positive and not exact, increment magnitude
 423           res++;
 424         }
 425         *pfpsf |= INEXACT_EXCEPTION;
 426       }
 427       // set exponent to zero as it was negative before.
 428       res = x_sign | 0x31c0000000000000ull | res;
 429       BID_RETURN (res);
 430     } else {    // if exp < 0 and q + exp <= 0
 431       // the result is -0 or +1
 432       if (x_sign) {
 433         res = 0xb1c0000000000000ull;
 434       } else {
 435         res = 0x31c0000000000001ull;
 436       }
 437       *pfpsf |= INEXACT_EXCEPTION;
 438       BID_RETURN (res);
 439     }
 440     break;
 441   case ROUNDING_TO_ZERO:
 442     if ((q + exp) >= 0) {       // exp < 0 and 1 <= -exp <= q
 443       // need to shift right -exp digits from the coefficient; exp will be 0
 444       ind = -exp;       // 1 <= ind <= 16; ind is a synonym for 'x'
 445       // chop off ind digits from the lower part of C1
 446       // C1 fits in 127 bits
 447       // calculate C* and f*
 448       // C* is actually floor(C*) in this case
 449       // C* and f* need shifting and masking, as shown by
 450       // shiftright128[] and maskhigh128[]
 451       // 1 <= x <= 16
 452       // kx = 10^(-x) = ten2mk64[ind - 1]
 453       // C* = C1 * 10^(-x)
 454       // the approximation of 10^(-x) was rounded up to 64 bits
 455       __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
 456
 457       // C* = floor(C*) (logical right shift; C has p decimal digits,
 458       //       correct by Property 1)
 459       // if (0 < f* < 10^(-x)) then the result is exact
 460       // n = C* * 10^(e+x)
 461
 462       if (ind - 1 <= 2) {       // 0 <= ind - 1 <= 2 => shift = 0
 463         res = P128.w[1];
 464         fstar.w[1] = 0;
 465         fstar.w[0] = P128.w[0];
 466       } else if (ind - 1 <= 21) {       // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 467         shift = shiftright128[ind - 1]; // 3 <= shift <= 63
 468         res = (P128.w[1] >> shift);
 469         fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
 470         fstar.w[0] = P128.w[0];
 471       }
 472       // if (f* > 10^(-x)) then the result is inexact
 473       if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) {
 474         *pfpsf |= INEXACT_EXCEPTION;
 475       }
 476       // set exponent to zero as it was negative before.
 477       res = x_sign | 0x31c0000000000000ull | res;
 478       BID_RETURN (res);
 479     } else {    // if exp < 0 and q + exp < 0
 480       // the result is +0 or -0
 481       res = x_sign | 0x31c0000000000000ull;
 482       *pfpsf |= INEXACT_EXCEPTION;
 483       BID_RETURN (res);
 484     }
 485     break;
 486   }     // end switch ()
 487   BID_RETURN (res);
 488 }
 489
 490 /*****************************************************************************
 491  *  BID64_round_integral_nearest_even
 492  ****************************************************************************/
 493
 494 #if DECIMAL_CALL_BY_REFERENCE
 495 void
 496 bid64_round_integral_nearest_even (UINT64 * pres,
 497                                    UINT64 *
 498                                    px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 499                                    _EXC_INFO_PARAM) {
 500   UINT64 x = *px;
 501 #else
 502 UINT64
 503 bid64_round_integral_nearest_even (UINT64 x _EXC_FLAGS_PARAM
 504                                    _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
 505 #endif
 506
 507   UINT64 res = 0xbaddbaddbaddbaddull;
 508   UINT64 x_sign;
 509   int exp;                      // unbiased exponent
 510   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
 511   BID_UI64DOUBLE tmp1;
 512   int x_nr_bits;
 513   int q, ind, shift;
 514   UINT64 C1;
 515   UINT128 fstar;
 516   UINT128 P128;
 517
 518   x_sign = x & MASK_SIGN;       // 0 for positive, MASK_SIGN for negative
 519
 520   // check for NaNs and infinities
 521   if ((x & MASK_NAN) == MASK_NAN) {     // check for NaN
 522     if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
 523       x = x & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
 524     else
 525       x = x & 0xfe03ffffffffffffull;    // clear G6-G12
 526     if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
 527       // set invalid flag
 528       *pfpsf |= INVALID_EXCEPTION;
 529       // return quiet (SNaN)
 530       res = x & 0xfdffffffffffffffull;
 531     } else {    // QNaN
 532       res = x;
 533     }
 534     BID_RETURN (res);
 535   } else if ((x & MASK_INF) == MASK_INF) {      // check for Infinity
 536     res = x_sign | 0x7800000000000000ull;
 537     BID_RETURN (res);
 538   }
 539   // unpack x
 540   if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 541     // if the steering bits are 11 (condition will be 0), then
 542     // the exponent is G[0:w+1]
 543     exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398;
 544     C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 545     if (C1 > 9999999999999999ull) {     // non-canonical
 546       C1 = 0;
 547     }
 548   } else {      // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
 549     exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398;
 550     C1 = (x & MASK_BINARY_SIG1);
 551   }
 552
 553   // if x is 0 or non-canonical
 554   if (C1 == 0) {
 555     if (exp < 0)
 556       exp = 0;
 557     res = x_sign | (((UINT64) exp + 398) << 53);
 558     BID_RETURN (res);
 559   }
 560   // x is a finite non-zero number (not 0, non-canonical, or special)
 561
 562   // return 0 if (exp <= -(p+1))
 563   if (exp <= -17) {
 564     res = x_sign | 0x31c0000000000000ull;
 565     BID_RETURN (res);
 566   }
 567   // q = nr. of decimal digits in x (1 <= q <= 54)
 568   //  determine first the nr. of bits in x
 569   if (C1 >= 0x0020000000000000ull) {    // x >= 2^53
 570     q = 16;
 571   } else {      // if x < 2^53
 572     tmp1.d = (double) C1;       // exact conversion
 573     x_nr_bits =
 574       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 575     q = nr_digits[x_nr_bits - 1].digits;
 576     if (q == 0) {
 577       q = nr_digits[x_nr_bits - 1].digits1;
 578       if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
 579         q++;
 580     }
 581   }
 582
 583   if (exp >= 0) {       // -exp <= 0
 584     // the argument is an integer already
 585     res = x;
 586     BID_RETURN (res);
 587   } else if ((q + exp) >= 0) {  // exp < 0 and 1 <= -exp <= q
 588     // need to shift right -exp digits from the coefficient; the exp will be 0
 589     ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x'
 590     // chop off ind digits from the lower part of C1
 591     // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits
 592     // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
 593     C1 = C1 + midpoint64[ind - 1];
 594     // calculate C* and f*
 595     // C* is actually floor(C*) in this case
 596     // C* and f* need shifting and masking, as shown by
 597     // shiftright128[] and maskhigh128[]
 598     // 1 <= x <= 16
 599     // kx = 10^(-x) = ten2mk64[ind - 1]
 600     // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 601     // the approximation of 10^(-x) was rounded up to 64 bits
 602     __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
 603
 604     // if (0 < f* < 10^(-x)) then the result is a midpoint
 605     //   if floor(C*) is even then C* = floor(C*) - logical right
 606     //       shift; C* has p decimal digits, correct by Prop. 1)
 607     //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
 608     //       shift; C* has p decimal digits, correct by Pr. 1)
 609     // else
 610     //   C* = floor(C*) (logical right shift; C has p decimal digits,
 611     //       correct by Property 1)
 612     // n = C* * 10^(e+x)
 613
 614     if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
 615       res = P128.w[1];
 616       fstar.w[1] = 0;
 617       fstar.w[0] = P128.w[0];
 618     } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 619       shift = shiftright128[ind - 1];   // 3 <= shift <= 63
 620       res = (P128.w[1] >> shift);
 621       fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
 622       fstar.w[0] = P128.w[0];
 623     }
 624     // if (0 < f* < 10^(-x)) then the result is a midpoint
 625     // since round_to_even, subtract 1 if current result is odd
 626     if ((res & 0x0000000000000001ull) && (fstar.w[1] == 0)
 627         && (fstar.w[0] < ten2mk64[ind - 1])) {
 628       res--;
 629     }
 630     // set exponent to zero as it was negative before.
 631     res = x_sign | 0x31c0000000000000ull | res;
 632     BID_RETURN (res);
 633   } else {      // if exp < 0 and q + exp < 0
 634     // the result is +0 or -0
 635     res = x_sign | 0x31c0000000000000ull;
 636     BID_RETURN (res);
 637   }
 638 }
 639
 640 /*****************************************************************************
 641  *  BID64_round_integral_negative
 642  *****************************************************************************/
 643
 644 #if DECIMAL_CALL_BY_REFERENCE
 645 void
 646 bid64_round_integral_negative (UINT64 * pres,
 647                                UINT64 *
 648                                px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 649                                _EXC_INFO_PARAM) {
 650   UINT64 x = *px;
 651 #else
 652 UINT64
 653 bid64_round_integral_negative (UINT64 x _EXC_FLAGS_PARAM
 654                                _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
 655 #endif
 656
 657   UINT64 res = 0xbaddbaddbaddbaddull;
 658   UINT64 x_sign;
 659   int exp;                      // unbiased exponent
 660   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
 661   BID_UI64DOUBLE tmp1;
 662   int x_nr_bits;
 663   int q, ind, shift;
 664   UINT64 C1;
 665   // UINT64 res is C* at first - represents up to 34 decimal digits ~ 113 bits
 666   UINT128 fstar;
 667   UINT128 P128;
 668
 669   x_sign = x & MASK_SIGN;       // 0 for positive, MASK_SIGN for negative
 670
 671   // check for NaNs and infinities
 672   if ((x & MASK_NAN) == MASK_NAN) {     // check for NaN
 673     if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
 674       x = x & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
 675     else
 676       x = x & 0xfe03ffffffffffffull;    // clear G6-G12
 677     if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
 678       // set invalid flag
 679       *pfpsf |= INVALID_EXCEPTION;
 680       // return quiet (SNaN)
 681       res = x & 0xfdffffffffffffffull;
 682     } else {    // QNaN
 683       res = x;
 684     }
 685     BID_RETURN (res);
 686   } else if ((x & MASK_INF) == MASK_INF) {      // check for Infinity
 687     res = x_sign | 0x7800000000000000ull;
 688     BID_RETURN (res);
 689   }
 690   // unpack x
 691   if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 692     // if the steering bits are 11 (condition will be 0), then
 693     // the exponent is G[0:w+1]
 694     exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398;
 695     C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 696     if (C1 > 9999999999999999ull) {     // non-canonical
 697       C1 = 0;
 698     }
 699   } else {      // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
 700     exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398;
 701     C1 = (x & MASK_BINARY_SIG1);
 702   }
 703
 704   // if x is 0 or non-canonical
 705   if (C1 == 0) {
 706     if (exp < 0)
 707       exp = 0;
 708     res = x_sign | (((UINT64) exp + 398) << 53);
 709     BID_RETURN (res);
 710   }
 711   // x is a finite non-zero number (not 0, non-canonical, or special)
 712
 713   // return 0 if (exp <= -p)
 714   if (exp <= -16) {
 715     if (x_sign) {
 716       res = 0xb1c0000000000001ull;
 717     } else {
 718       res = 0x31c0000000000000ull;
 719     }
 720     BID_RETURN (res);
 721   }
 722   // q = nr. of decimal digits in x (1 <= q <= 54)
 723   //  determine first the nr. of bits in x
 724   if (C1 >= 0x0020000000000000ull) {    // x >= 2^53
 725     q = 16;
 726   } else {      // if x < 2^53
 727     tmp1.d = (double) C1;       // exact conversion
 728     x_nr_bits =
 729       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 730     q = nr_digits[x_nr_bits - 1].digits;
 731     if (q == 0) {
 732       q = nr_digits[x_nr_bits - 1].digits1;
 733       if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
 734         q++;
 735     }
 736   }
 737
 738   if (exp >= 0) {       // -exp <= 0
 739     // the argument is an integer already
 740     res = x;
 741     BID_RETURN (res);
 742   } else if ((q + exp) > 0) {   // exp < 0 and 1 <= -exp < q
 743     // need to shift right -exp digits from the coefficient; the exp will be 0
 744     ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x'
 745     // chop off ind digits from the lower part of C1
 746     // C1 fits in 64 bits
 747     // calculate C* and f*
 748     // C* is actually floor(C*) in this case
 749     // C* and f* need shifting and masking, as shown by
 750     // shiftright128[] and maskhigh128[]
 751     // 1 <= x <= 16
 752     // kx = 10^(-x) = ten2mk64[ind - 1]
 753     // C* = C1 * 10^(-x)
 754     // the approximation of 10^(-x) was rounded up to 64 bits
 755     __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
 756
 757     // C* = floor(C*) (logical right shift; C has p decimal digits,
 758     //       correct by Property 1)
 759     // if (0 < f* < 10^(-x)) then the result is exact
 760     // n = C* * 10^(e+x)
 761
 762     if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
 763       res = P128.w[1];
 764       fstar.w[1] = 0;
 765       fstar.w[0] = P128.w[0];
 766     } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 767       shift = shiftright128[ind - 1];   // 3 <= shift <= 63
 768       res = (P128.w[1] >> shift);
 769       fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
 770       fstar.w[0] = P128.w[0];
 771     }
 772     // if (f* > 10^(-x)) then the result is inexact
 773     if (x_sign
 774         && ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1]))) {
 775       // if negative and not exact, increment magnitude
 776       res++;
 777     }
 778     // set exponent to zero as it was negative before.
 779     res = x_sign | 0x31c0000000000000ull | res;
 780     BID_RETURN (res);
 781   } else {      // if exp < 0 and q + exp <= 0
 782     // the result is +0 or -1
 783     if (x_sign) {
 784       res = 0xb1c0000000000001ull;
 785     } else {
 786       res = 0x31c0000000000000ull;
 787     }
 788     BID_RETURN (res);
 789   }
 790 }
 791
 792 /*****************************************************************************
 793  *  BID64_round_integral_positive
 794  ****************************************************************************/
 795
 796 #if DECIMAL_CALL_BY_REFERENCE
 797 void
 798 bid64_round_integral_positive (UINT64 * pres,
 799                                UINT64 *
 800                                px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 801                                _EXC_INFO_PARAM) {
 802   UINT64 x = *px;
 803 #else
 804 UINT64
 805 bid64_round_integral_positive (UINT64 x _EXC_FLAGS_PARAM
 806                                _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
 807 #endif
 808
 809   UINT64 res = 0xbaddbaddbaddbaddull;
 810   UINT64 x_sign;
 811   int exp;                      // unbiased exponent
 812   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
 813   BID_UI64DOUBLE tmp1;
 814   int x_nr_bits;
 815   int q, ind, shift;
 816   UINT64 C1;
 817   // UINT64 res is C* at first - represents up to 34 decimal digits ~ 113 bits
 818   UINT128 fstar;
 819   UINT128 P128;
 820
 821   x_sign = x & MASK_SIGN;       // 0 for positive, MASK_SIGN for negative
 822
 823   // check for NaNs and infinities
 824   if ((x & MASK_NAN) == MASK_NAN) {     // check for NaN
 825     if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
 826       x = x & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
 827     else
 828       x = x & 0xfe03ffffffffffffull;    // clear G6-G12
 829     if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
 830       // set invalid flag
 831       *pfpsf |= INVALID_EXCEPTION;
 832       // return quiet (SNaN)
 833       res = x & 0xfdffffffffffffffull;
 834     } else {    // QNaN
 835       res = x;
 836     }
 837     BID_RETURN (res);
 838   } else if ((x & MASK_INF) == MASK_INF) {      // check for Infinity
 839     res = x_sign | 0x7800000000000000ull;
 840     BID_RETURN (res);
 841   }
 842   // unpack x
 843   if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 844     // if the steering bits are 11 (condition will be 0), then
 845     // the exponent is G[0:w+1]
 846     exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398;
 847     C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 848     if (C1 > 9999999999999999ull) {     // non-canonical
 849       C1 = 0;
 850     }
 851   } else {      // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
 852     exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398;
 853     C1 = (x & MASK_BINARY_SIG1);
 854   }
 855
 856   // if x is 0 or non-canonical
 857   if (C1 == 0) {
 858     if (exp < 0)
 859       exp = 0;
 860     res = x_sign | (((UINT64) exp + 398) << 53);
 861     BID_RETURN (res);
 862   }
 863   // x is a finite non-zero number (not 0, non-canonical, or special)
 864
 865   // return 0 if (exp <= -p)
 866   if (exp <= -16) {
 867     if (x_sign) {
 868       res = 0xb1c0000000000000ull;
 869     } else {
 870       res = 0x31c0000000000001ull;
 871     }
 872     BID_RETURN (res);
 873   }
 874   // q = nr. of decimal digits in x (1 <= q <= 54)
 875   //  determine first the nr. of bits in x
 876   if (C1 >= 0x0020000000000000ull) {    // x >= 2^53
 877     q = 16;
 878   } else {      // if x < 2^53
 879     tmp1.d = (double) C1;       // exact conversion
 880     x_nr_bits =
 881       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 882     q = nr_digits[x_nr_bits - 1].digits;
 883     if (q == 0) {
 884       q = nr_digits[x_nr_bits - 1].digits1;
 885       if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
 886         q++;
 887     }
 888   }
 889
 890   if (exp >= 0) {       // -exp <= 0
 891     // the argument is an integer already
 892     res = x;
 893     BID_RETURN (res);
 894   } else if ((q + exp) > 0) {   // exp < 0 and 1 <= -exp < q
 895     // need to shift right -exp digits from the coefficient; the exp will be 0
 896     ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x'
 897     // chop off ind digits from the lower part of C1
 898     // C1 fits in 64 bits
 899     // calculate C* and f*
 900     // C* is actually floor(C*) in this case
 901     // C* and f* need shifting and masking, as shown by
 902     // shiftright128[] and maskhigh128[]
 903     // 1 <= x <= 16
 904     // kx = 10^(-x) = ten2mk64[ind - 1]
 905     // C* = C1 * 10^(-x)
 906     // the approximation of 10^(-x) was rounded up to 64 bits
 907     __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
 908
 909     // C* = floor(C*) (logical right shift; C has p decimal digits,
 910     //       correct by Property 1)
 911     // if (0 < f* < 10^(-x)) then the result is exact
 912     // n = C* * 10^(e+x)
 913
 914     if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
 915       res = P128.w[1];
 916       fstar.w[1] = 0;
 917       fstar.w[0] = P128.w[0];
 918     } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 919       shift = shiftright128[ind - 1];   // 3 <= shift <= 63
 920       res = (P128.w[1] >> shift);
 921       fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
 922       fstar.w[0] = P128.w[0];
 923     }
 924     // if (f* > 10^(-x)) then the result is inexact
 925     if (!x_sign
 926         && ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1]))) {
 927       // if positive and not exact, increment magnitude
 928       res++;
 929     }
 930     // set exponent to zero as it was negative before.
 931     res = x_sign | 0x31c0000000000000ull | res;
 932     BID_RETURN (res);
 933   } else {      // if exp < 0 and q + exp <= 0
 934     // the result is -0 or +1
 935     if (x_sign) {
 936       res = 0xb1c0000000000000ull;
 937     } else {
 938       res = 0x31c0000000000001ull;
 939     }
 940     BID_RETURN (res);
 941   }
 942 }
 943
 944 /*****************************************************************************
 945  *  BID64_round_integral_zero
 946  ****************************************************************************/
 947
 948 #if DECIMAL_CALL_BY_REFERENCE
 949 void
 950 bid64_round_integral_zero (UINT64 * pres,
 951                            UINT64 *
 952                            px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 953                            _EXC_INFO_PARAM) {
 954   UINT64 x = *px;
 955 #else
 956 UINT64
 957 bid64_round_integral_zero (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 958                            _EXC_INFO_PARAM) {
 959 #endif
 960
 961   UINT64 res = 0xbaddbaddbaddbaddull;
 962   UINT64 x_sign;
 963   int exp;                      // unbiased exponent
 964   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
 965   BID_UI64DOUBLE tmp1;
 966   int x_nr_bits;
 967   int q, ind, shift;
 968   UINT64 C1;
 969   // UINT64 res is C* at first - represents up to 34 decimal digits ~ 113 bits
 970   UINT128 P128;
 971
 972   x_sign = x & MASK_SIGN;       // 0 for positive, MASK_SIGN for negative
 973
 974   // check for NaNs and infinities
 975   if ((x & MASK_NAN) == MASK_NAN) {     // check for NaN
 976     if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
 977       x = x & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
 978     else
 979       x = x & 0xfe03ffffffffffffull;    // clear G6-G12
 980     if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
 981       // set invalid flag
 982       *pfpsf |= INVALID_EXCEPTION;
 983       // return quiet (SNaN)
 984       res = x & 0xfdffffffffffffffull;
 985     } else {    // QNaN
 986       res = x;
 987     }
 988     BID_RETURN (res);
 989   } else if ((x & MASK_INF) == MASK_INF) {      // check for Infinity
 990     res = x_sign | 0x7800000000000000ull;
 991     BID_RETURN (res);
 992   }
 993   // unpack x
 994   if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 995     // if the steering bits are 11 (condition will be 0), then
 996     // the exponent is G[0:w+1]
 997     exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398;
 998     C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 999     if (C1 > 9999999999999999ull) {     // non-canonical
1000       C1 = 0;
1001     }
1002   } else {      // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
1003     exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398;
1004     C1 = (x & MASK_BINARY_SIG1);
1005   }
1006
1007   // if x is 0 or non-canonical
1008   if (C1 == 0) {
1009     if (exp < 0)
1010       exp = 0;
1011     res = x_sign | (((UINT64) exp + 398) << 53);
1012     BID_RETURN (res);
1013   }
1014   // x is a finite non-zero number (not 0, non-canonical, or special)
1015
1016   // return 0 if (exp <= -p)
1017   if (exp <= -16) {
1018     res = x_sign | 0x31c0000000000000ull;
1019     BID_RETURN (res);
1020   }
1021   // q = nr. of decimal digits in x (1 <= q <= 54)
1022   //  determine first the nr. of bits in x
1023   if (C1 >= 0x0020000000000000ull) {    // x >= 2^53
1024     q = 16;
1025   } else {      // if x < 2^53
1026     tmp1.d = (double) C1;       // exact conversion
1027     x_nr_bits =
1028       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1029     q = nr_digits[x_nr_bits - 1].digits;
1030     if (q == 0) {
1031       q = nr_digits[x_nr_bits - 1].digits1;
1032       if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1033         q++;
1034     }
1035   }
1036
1037   if (exp >= 0) {       // -exp <= 0
1038     // the argument is an integer already
1039     res = x;
1040     BID_RETURN (res);
1041   } else if ((q + exp) >= 0) {  // exp < 0 and 1 <= -exp <= q
1042     // need to shift right -exp digits from the coefficient; the exp will be 0
1043     ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x'
1044     // chop off ind digits from the lower part of C1
1045     // C1 fits in 127 bits
1046     // calculate C* and f*
1047     // C* is actually floor(C*) in this case
1048     // C* and f* need shifting and masking, as shown by
1049     // shiftright128[] and maskhigh128[]
1050     // 1 <= x <= 16
1051     // kx = 10^(-x) = ten2mk64[ind - 1]
1052     // C* = C1 * 10^(-x)
1053     // the approximation of 10^(-x) was rounded up to 64 bits
1054     __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
1055
1056     // C* = floor(C*) (logical right shift; C has p decimal digits,
1057     //       correct by Property 1)
1058     // if (0 < f* < 10^(-x)) then the result is exact
1059     // n = C* * 10^(e+x)
1060
1061     if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
1062       res = P128.w[1];
1063       // redundant fstar.w[1] = 0;
1064       // redundant fstar.w[0] = P128.w[0];
1065     } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1066       shift = shiftright128[ind - 1];   // 3 <= shift <= 63
1067       res = (P128.w[1] >> shift);
1068       // redundant fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
1069       // redundant fstar.w[0] = P128.w[0];
1070     }
1071     // if (f* > 10^(-x)) then the result is inexact
1072     // if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind-1])){
1073     //   // redundant
1074     // }
1075     // set exponent to zero as it was negative before.
1076     res = x_sign | 0x31c0000000000000ull | res;
1077     BID_RETURN (res);
1078   } else {      // if exp < 0 and q + exp < 0
1079     // the result is +0 or -0
1080     res = x_sign | 0x31c0000000000000ull;
1081     BID_RETURN (res);
1082   }
1083 }
1084
1085 /*****************************************************************************
1086  *  BID64_round_integral_nearest_away
1087  ****************************************************************************/
1088
1089 #if DECIMAL_CALL_BY_REFERENCE
1090 void
1091 bid64_round_integral_nearest_away (UINT64 * pres,
1092                                    UINT64 *
1093                                    px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1094                                    _EXC_INFO_PARAM) {
1095   UINT64 x = *px;
1096 #else
1097 UINT64
1098 bid64_round_integral_nearest_away (UINT64 x _EXC_FLAGS_PARAM
1099                                    _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
1100 #endif
1101
1102   UINT64 res = 0xbaddbaddbaddbaddull;
1103   UINT64 x_sign;
1104   int exp;                      // unbiased exponent
1105   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1106   BID_UI64DOUBLE tmp1;
1107   int x_nr_bits;
1108   int q, ind, shift;
1109   UINT64 C1;
1110   UINT128 P128;
1111
1112   x_sign = x & MASK_SIGN;       // 0 for positive, MASK_SIGN for negative
1113
1114   // check for NaNs and infinities
1115   if ((x & MASK_NAN) == MASK_NAN) {     // check for NaN
1116     if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
1117       x = x & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
1118     else
1119       x = x & 0xfe03ffffffffffffull;    // clear G6-G12
1120     if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
1121       // set invalid flag
1122       *pfpsf |= INVALID_EXCEPTION;
1123       // return quiet (SNaN)
1124       res = x & 0xfdffffffffffffffull;
1125     } else {    // QNaN
1126       res = x;
1127     }
1128     BID_RETURN (res);
1129   } else if ((x & MASK_INF) == MASK_INF) {      // check for Infinity
1130     res = x_sign | 0x7800000000000000ull;
1131     BID_RETURN (res);
1132   }
1133   // unpack x
1134   if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
1135     // if the steering bits are 11 (condition will be 0), then
1136     // the exponent is G[0:w+1]
1137     exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398;
1138     C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
1139     if (C1 > 9999999999999999ull) {     // non-canonical
1140       C1 = 0;
1141     }
1142   } else {      // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
1143     exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398;
1144     C1 = (x & MASK_BINARY_SIG1);
1145   }
1146
1147   // if x is 0 or non-canonical
1148   if (C1 == 0) {
1149     if (exp < 0)
1150       exp = 0;
1151     res = x_sign | (((UINT64) exp + 398) << 53);
1152     BID_RETURN (res);
1153   }
1154   // x is a finite non-zero number (not 0, non-canonical, or special)
1155
1156   // return 0 if (exp <= -(p+1))
1157   if (exp <= -17) {
1158     res = x_sign | 0x31c0000000000000ull;
1159     BID_RETURN (res);
1160   }
1161   // q = nr. of decimal digits in x (1 <= q <= 54)
1162   //  determine first the nr. of bits in x
1163   if (C1 >= 0x0020000000000000ull) {    // x >= 2^53
1164     q = 16;
1165   } else {      // if x < 2^53
1166     tmp1.d = (double) C1;       // exact conversion
1167     x_nr_bits =
1168       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1169     q = nr_digits[x_nr_bits - 1].digits;
1170     if (q == 0) {
1171       q = nr_digits[x_nr_bits - 1].digits1;
1172       if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1173         q++;
1174     }
1175   }
1176
1177   if (exp >= 0) {       // -exp <= 0
1178     // the argument is an integer already
1179     res = x;
1180     BID_RETURN (res);
1181   } else if ((q + exp) >= 0) {  // exp < 0 and 1 <= -exp <= q
1182     // need to shift right -exp digits from the coefficient; the exp will be 0
1183     ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x'
1184     // chop off ind digits from the lower part of C1
1185     // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits
1186     // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
1187     C1 = C1 + midpoint64[ind - 1];
1188     // calculate C* and f*
1189     // C* is actually floor(C*) in this case
1190     // C* and f* need shifting and masking, as shown by
1191     // shiftright128[] and maskhigh128[]
1192     // 1 <= x <= 16
1193     // kx = 10^(-x) = ten2mk64[ind - 1]
1194     // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1195     // the approximation of 10^(-x) was rounded up to 64 bits
1196     __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
1197
1198     // if (0 < f* < 10^(-x)) then the result is a midpoint
1199     //   C* = floor(C*) - logical right shift; C* has p decimal digits,
1200     //       correct by Prop. 1)
1201     // else
1202     //   C* = floor(C*) (logical right shift; C has p decimal digits,
1203     //       correct by Property 1)
1204     // n = C* * 10^(e+x)
1205
1206     if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
1207       res = P128.w[1];
1208     } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1209       shift = shiftright128[ind - 1];   // 3 <= shift <= 63
1210       res = (P128.w[1] >> shift);
1211     }
1212     // midpoints are already rounded correctly
1213     // set exponent to zero as it was negative before.
1214     res = x_sign | 0x31c0000000000000ull | res;
1215     BID_RETURN (res);
1216   } else {      // if exp < 0 and q + exp < 0
1217     // the result is +0 or -0
1218     res = x_sign | 0x31c0000000000000ull;
1219     BID_RETURN (res);
1220   }
1221 }