libgcc/config/libbid/bid64_round_integral.c

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