libgcc/config/libbid/bid64_round_integral.c

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