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