libgcc/config/libbid/bid128_round_integral.c

   1 /* Copyright (C) 2007, 2009  Free Software Foundation, Inc.
   2
   3 This file is part of GCC.
   4
   5 GCC is free software; you can redistribute it and/or modify it under
   6 the terms of the GNU General Public License as published by the Free
   7 Software Foundation; either version 3, or (at your option) any later
   8 version.
   9
  10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13 for more details.
  14
  15 Under Section 7 of GPL version 3, you are granted additional
  16 permissions described in the GCC Runtime Library Exception, version
  17 3.1, as published by the Free Software Foundation.
  18
  19 You should have received a copy of the GNU General Public License and
  20 a copy of the GCC Runtime Library Exception along with this program;
  21 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  22 <http://www.gnu.org/licenses/>.  */
  23
  24 #define BID_128RES
  25
  26 #include "bid_internal.h"
  27
  28 /*****************************************************************************
  29  *  BID128_round_integral_exact
  30  ****************************************************************************/
  31
  32 BID128_FUNCTION_ARG1 (bid128_round_integral_exact, x)
  33
  34      UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
  35      };
  36 UINT64 x_sign;
  37 UINT64 x_exp;
  38 int exp;                        // unbiased exponent
  39   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  40 UINT64 tmp64;
  41 BID_UI64DOUBLE tmp1;
  42 unsigned int x_nr_bits;
  43 int q, ind, shift;
  44 UINT128 C1;
  45 UINT256 fstar;
  46 UINT256 P256;
  47
  48   // check for NaN or Infinity
  49 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  50   // x is special
  51   if ((x.w[1] & MASK_NAN) == MASK_NAN) {        // x is NAN
  52     // if x = NaN, then res = Q (x)
  53     // check first for non-canonical NaN payload
  54     if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  55         (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
  56          (x.w[0] > 0x38c15b09ffffffffull))) {
  57       x.w[1] = x.w[1] & 0xffffc00000000000ull;
  58       x.w[0] = 0x0ull;
  59     }
  60     if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {    // x is SNAN
  61       // set invalid flag
  62       *pfpsf |= INVALID_EXCEPTION;
  63       // return quiet (x)
  64       res.w[1] = x.w[1] & 0xfc003fffffffffffull;        // clear out also G[6]-G[16]
  65       res.w[0] = x.w[0];
  66     } else {    // x is QNaN
  67       // return x
  68       res.w[1] = x.w[1] & 0xfc003fffffffffffull;        // clear out G[6]-G[16]
  69       res.w[0] = x.w[0];
  70     }
  71     BID_RETURN (res)
  72   } else {      // x is not a NaN, so it must be infinity
  73     if ((x.w[1] & MASK_SIGN) == 0x0ull) {       // x is +inf
  74       // return +inf
  75       res.w[1] = 0x7800000000000000ull;
  76       res.w[0] = 0x0000000000000000ull;
  77     } else {    // x is -inf
  78       // return -inf
  79       res.w[1] = 0xf800000000000000ull;
  80       res.w[0] = 0x0000000000000000ull;
  81     }
  82     BID_RETURN (res);
  83   }
  84 }
  85   // unpack x
  86 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  87 C1.w[1] = x.w[1] & MASK_COEFF;
  88 C1.w[0] = x.w[0];
  89
  90   // check for non-canonical values (treated as zero)
  91 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
  92   // non-canonical
  93   x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
  94   C1.w[1] = 0;  // significand high
  95   C1.w[0] = 0;  // significand low
  96 } else {        // G0_G1 != 11
  97   x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
  98   if (C1.w[1] > 0x0001ed09bead87c0ull ||
  99       (C1.w[1] == 0x0001ed09bead87c0ull
 100        && C1.w[0] > 0x378d8e63ffffffffull)) {
 101     // x is non-canonical if coefficient is larger than 10^34 -1
 102     C1.w[1] = 0;
 103     C1.w[0] = 0;
 104   } else {      // canonical
 105     ;
 106   }
 107 }
 108
 109   // test for input equal to zero
 110 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
 111   // x is 0
 112   // return 0 preserving the sign bit and the preferred exponent
 113   // of MAX(Q(x), 0)
 114   if (x_exp <= (0x1820ull << 49)) {
 115     res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
 116   } else {
 117     res.w[1] = x_sign | x_exp;
 118   }
 119   res.w[0] = 0x0000000000000000ull;
 120   BID_RETURN (res);
 121 }
 122   // x is not special and is not zero
 123
 124 switch (rnd_mode) {
 125 case ROUNDING_TO_NEAREST:
 126 case ROUNDING_TIES_AWAY:
 127   // if (exp <= -(p+1)) return 0.0
 128   if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35
 129     res.w[1] = x_sign | 0x3040000000000000ull;
 130     res.w[0] = 0x0000000000000000ull;
 131     *pfpsf |= INEXACT_EXCEPTION;
 132     BID_RETURN (res);
 133   }
 134   break;
 135 case ROUNDING_DOWN:
 136   // if (exp <= -p) return -1.0 or +0.0
 137   if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffa000000000000ull == -34
 138     if (x_sign) {
 139       // if negative, return negative 1, because we know coefficient
 140       // is non-zero (would have been caught above)
 141       res.w[1] = 0xb040000000000000ull;
 142       res.w[0] = 0x0000000000000001ull;
 143     } else {
 144       // if positive, return positive 0, because we know coefficient is
 145       // non-zero (would have been caught above)
 146       res.w[1] = 0x3040000000000000ull;
 147       res.w[0] = 0x0000000000000000ull;
 148     }
 149     *pfpsf |= INEXACT_EXCEPTION;
 150     BID_RETURN (res);
 151   }
 152   break;
 153 case ROUNDING_UP:
 154   // if (exp <= -p) return -0.0 or +1.0
 155   if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
 156     if (x_sign) {
 157       // if negative, return negative 0, because we know the coefficient
 158       // is non-zero (would have been caught above)
 159       res.w[1] = 0xb040000000000000ull;
 160       res.w[0] = 0x0000000000000000ull;
 161     } else {
 162       // if positive, return positive 1, because we know coefficient is
 163       // non-zero (would have been caught above)
 164       res.w[1] = 0x3040000000000000ull;
 165       res.w[0] = 0x0000000000000001ull;
 166     }
 167     *pfpsf |= INEXACT_EXCEPTION;
 168     BID_RETURN (res);
 169   }
 170   break;
 171 case ROUNDING_TO_ZERO:
 172   // if (exp <= -p) return -0.0 or +0.0
 173   if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
 174     res.w[1] = x_sign | 0x3040000000000000ull;
 175     res.w[0] = 0x0000000000000000ull;
 176     *pfpsf |= INEXACT_EXCEPTION;
 177     BID_RETURN (res);
 178   }
 179   break;
 180 }
 181
 182   // q = nr. of decimal digits in x
 183   //  determine first the nr. of bits in x
 184 if (C1.w[1] == 0) {
 185   if (C1.w[0] >= 0x0020000000000000ull) {       // x >= 2^53
 186     // split the 64-bit value in two 32-bit halves to avoid rounding errors
 187     if (C1.w[0] >= 0x0000000100000000ull) {     // x >= 2^32
 188       tmp1.d = (double) (C1.w[0] >> 32);        // exact conversion
 189       x_nr_bits =
 190         33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 191     } else {    // x < 2^32
 192       tmp1.d = (double) (C1.w[0]);      // exact conversion
 193       x_nr_bits =
 194         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 195     }
 196   } else {      // if x < 2^53
 197     tmp1.d = (double) C1.w[0];  // exact conversion
 198     x_nr_bits =
 199       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 200   }
 201 } else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
 202   tmp1.d = (double) C1.w[1];    // exact conversion
 203   x_nr_bits =
 204     65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 205 }
 206
 207 q = nr_digits[x_nr_bits - 1].digits;
 208 if (q == 0) {
 209   q = nr_digits[x_nr_bits - 1].digits1;
 210   if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
 211       (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
 212        C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
 213     q++;
 214 }
 215 exp = (x_exp >> 49) - 6176;
 216 if (exp >= 0) { // -exp <= 0
 217   // the argument is an integer already
 218   res.w[1] = x.w[1];
 219   res.w[0] = x.w[0];
 220   BID_RETURN (res);
 221 }
 222   // exp < 0
 223 switch (rnd_mode) {
 224 case ROUNDING_TO_NEAREST:
 225   if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
 226     // need to shift right -exp digits from the coefficient; exp will be 0
 227     ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
 228     // chop off ind digits from the lower part of C1
 229     // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
 230     tmp64 = C1.w[0];
 231     if (ind <= 19) {
 232       C1.w[0] = C1.w[0] + midpoint64[ind - 1];
 233     } else {
 234       C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
 235       C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
 236     }
 237     if (C1.w[0] < tmp64)
 238       C1.w[1]++;
 239     // calculate C* and f*
 240     // C* is actually floor(C*) in this case
 241     // C* and f* need shifting and masking, as shown by
 242     // shiftright128[] and maskhigh128[]
 243     // 1 <= x <= 34
 244     // kx = 10^(-x) = ten2mk128[ind - 1]
 245     // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 246     // the approximation of 10^(-x) was rounded up to 118 bits
 247     __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
 248     // determine the value of res and fstar
 249
 250     // determine inexactness of the rounding of C*
 251     // if (0 < f* - 1/2 < 10^(-x)) then
 252     //   the result is exact
 253     // else // if (f* - 1/2 > T*) then
 254     //   the result is inexact
 255     // Note: we are going to use ten2mk128[] instead of ten2mk128trunc[]
 256
 257     if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
 258       // redundant shift = shiftright128[ind - 1]; // shift = 0
 259       res.w[1] = P256.w[3];
 260       res.w[0] = P256.w[2];
 261       // redundant fstar.w[3] = 0;
 262       // redundant fstar.w[2] = 0;
 263       fstar.w[1] = P256.w[1];
 264       fstar.w[0] = P256.w[0];
 265       // fraction f* < 10^(-x) <=> midpoint
 266       // f* is in the right position to be compared with
 267       // 10^(-x) from ten2mk128[]
 268       // if 0 < fstar < 10^(-x), subtract 1 if odd (for rounding to even)
 269       if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
 270           ((fstar.w[1] < (ten2mk128[ind - 1].w[1]))
 271            || ((fstar.w[1] == ten2mk128[ind - 1].w[1])
 272                && (fstar.w[0] < ten2mk128[ind - 1].w[0])))) {
 273         // subract 1 to make even
 274         if (res.w[0]-- == 0) {
 275           res.w[1]--;
 276         }
 277       }
 278       if (fstar.w[1] > 0x8000000000000000ull ||
 279           (fstar.w[1] == 0x8000000000000000ull
 280            && fstar.w[0] > 0x0ull)) {
 281         // f* > 1/2 and the result may be exact
 282         tmp64 = fstar.w[1] - 0x8000000000000000ull;     // f* - 1/2
 283         if (tmp64 > ten2mk128[ind - 1].w[1] ||
 284             (tmp64 == ten2mk128[ind - 1].w[1] &&
 285              fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 286           // set the inexact flag
 287           *pfpsf |= INEXACT_EXCEPTION;
 288         }       // else the result is exact
 289       } else {  // the result is inexact; f2* <= 1/2
 290         // set the inexact flag
 291         *pfpsf |= INEXACT_EXCEPTION;
 292       }
 293     } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 294       shift = shiftright128[ind - 1];   // 3 <= shift <= 63
 295       res.w[1] = (P256.w[3] >> shift);
 296       res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
 297       // redundant fstar.w[3] = 0;
 298       fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
 299       fstar.w[1] = P256.w[1];
 300       fstar.w[0] = P256.w[0];
 301       // fraction f* < 10^(-x) <=> midpoint
 302       // f* is in the right position to be compared with
 303       // 10^(-x) from ten2mk128[]
 304       if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
 305           fstar.w[2] == 0 && (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
 306                               (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
 307                                fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
 308         // subract 1 to make even
 309         if (res.w[0]-- == 0) {
 310           res.w[1]--;
 311         }
 312       }
 313       if (fstar.w[2] > onehalf128[ind - 1] ||
 314           (fstar.w[2] == onehalf128[ind - 1]
 315            && (fstar.w[1] || fstar.w[0]))) {
 316         // f2* > 1/2 and the result may be exact
 317         // Calculate f2* - 1/2
 318         tmp64 = fstar.w[2] - onehalf128[ind - 1];
 319         if (tmp64 || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
 320             (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
 321              fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 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     } else {    // 22 <= ind - 1 <= 33
 330       shift = shiftright128[ind - 1] - 64;      // 2 <= shift <= 38
 331       res.w[1] = 0;
 332       res.w[0] = P256.w[3] >> shift;
 333       fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
 334       fstar.w[2] = P256.w[2];
 335       fstar.w[1] = P256.w[1];
 336       fstar.w[0] = P256.w[0];
 337       // fraction f* < 10^(-x) <=> midpoint
 338       // f* is in the right position to be compared with
 339       // 10^(-x) from ten2mk128[]
 340       if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
 341           fstar.w[3] == 0 && fstar.w[2] == 0 &&
 342           (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
 343            (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
 344             fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
 345         // subract 1 to make even
 346         if (res.w[0]-- == 0) {
 347           res.w[1]--;
 348         }
 349       }
 350       if (fstar.w[3] > onehalf128[ind - 1] ||
 351           (fstar.w[3] == onehalf128[ind - 1] &&
 352            (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
 353         // f2* > 1/2 and the result may be exact
 354         // Calculate f2* - 1/2
 355         tmp64 = fstar.w[3] - onehalf128[ind - 1];
 356         if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1]
 357             || (fstar.w[1] == ten2mk128[ind - 1].w[1]
 358                 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 359           // set the inexact flag
 360           *pfpsf |= INEXACT_EXCEPTION;
 361         }       // else the result is exact
 362       } else {  // the result is inexact; f2* <= 1/2
 363         // set the inexact flag
 364         *pfpsf |= INEXACT_EXCEPTION;
 365       }
 366     }
 367     res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
 368     BID_RETURN (res);
 369   } else {      // if ((q + exp) < 0) <=> q < -exp
 370     // the result is +0 or -0
 371     res.w[1] = x_sign | 0x3040000000000000ull;
 372     res.w[0] = 0x0000000000000000ull;
 373     *pfpsf |= INEXACT_EXCEPTION;
 374     BID_RETURN (res);
 375   }
 376   break;
 377 case ROUNDING_TIES_AWAY:
 378   if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
 379     // need to shift right -exp digits from the coefficient; exp will be 0
 380     ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
 381     // chop off ind digits from the lower part of C1
 382     // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
 383     tmp64 = C1.w[0];
 384     if (ind <= 19) {
 385       C1.w[0] = C1.w[0] + midpoint64[ind - 1];
 386     } else {
 387       C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
 388       C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
 389     }
 390     if (C1.w[0] < tmp64)
 391       C1.w[1]++;
 392     // calculate C* and f*
 393     // C* is actually floor(C*) in this case
 394     // C* and f* need shifting and masking, as shown by
 395     // shiftright128[] and maskhigh128[]
 396     // 1 <= x <= 34
 397     // kx = 10^(-x) = ten2mk128[ind - 1]
 398     // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 399     // the approximation of 10^(-x) was rounded up to 118 bits
 400     __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
 401     // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
 402     // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
 403     // if (0 < f* < 10^(-x)) then the result is a midpoint
 404     //   if floor(C*) is even then C* = floor(C*) - logical right
 405     //       shift; C* has p decimal digits, correct by Prop. 1)
 406     //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
 407     //       shift; C* has p decimal digits, correct by Pr. 1)
 408     // else
 409     //   C* = floor(C*) (logical right shift; C has p decimal digits,
 410     //       correct by Property 1)
 411     // n = C* * 10^(e+x)
 412
 413     // determine also the inexactness of the rounding of C*
 414     // if (0 < f* - 1/2 < 10^(-x)) then
 415     //   the result is exact
 416     // else // if (f* - 1/2 > T*) then
 417     //   the result is inexact
 418     // Note: we are going to use ten2mk128[] instead of ten2mk128trunc[]
 419     // shift right C* by Ex-128 = shiftright128[ind]
 420     if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
 421       // redundant shift = shiftright128[ind - 1]; // shift = 0
 422       res.w[1] = P256.w[3];
 423       res.w[0] = P256.w[2];
 424       // redundant fstar.w[3] = 0;
 425       // redundant fstar.w[2] = 0;
 426       fstar.w[1] = P256.w[1];
 427       fstar.w[0] = P256.w[0];
 428       if (fstar.w[1] > 0x8000000000000000ull ||
 429           (fstar.w[1] == 0x8000000000000000ull
 430            && fstar.w[0] > 0x0ull)) {
 431         // f* > 1/2 and the result may be exact
 432         tmp64 = fstar.w[1] - 0x8000000000000000ull;     // f* - 1/2
 433         if ((tmp64 > ten2mk128[ind - 1].w[1] ||
 434              (tmp64 == ten2mk128[ind - 1].w[1] &&
 435               fstar.w[0] >= ten2mk128[ind - 1].w[0]))) {
 436           // set the inexact flag
 437           *pfpsf |= INEXACT_EXCEPTION;
 438         }       // else the result is exact
 439       } else {  // the result is inexact; f2* <= 1/2
 440         // set the inexact flag
 441         *pfpsf |= INEXACT_EXCEPTION;
 442       }
 443     } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 444       shift = shiftright128[ind - 1];   // 3 <= shift <= 63
 445       res.w[1] = (P256.w[3] >> shift);
 446       res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
 447       // redundant fstar.w[3] = 0;
 448       fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
 449       fstar.w[1] = P256.w[1];
 450       fstar.w[0] = P256.w[0];
 451       if (fstar.w[2] > onehalf128[ind - 1] ||
 452           (fstar.w[2] == onehalf128[ind - 1]
 453            && (fstar.w[1] || fstar.w[0]))) {
 454         // f2* > 1/2 and the result may be exact
 455         // Calculate f2* - 1/2
 456         tmp64 = fstar.w[2] - onehalf128[ind - 1];
 457         if (tmp64 || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
 458             (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
 459              fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 460           // set the inexact flag
 461           *pfpsf |= INEXACT_EXCEPTION;
 462         }       // else the result is exact
 463       } else {  // the result is inexact; f2* <= 1/2
 464         // set the inexact flag
 465         *pfpsf |= INEXACT_EXCEPTION;
 466       }
 467     } else {    // 22 <= ind - 1 <= 33
 468       shift = shiftright128[ind - 1] - 64;      // 2 <= shift <= 38
 469       res.w[1] = 0;
 470       res.w[0] = P256.w[3] >> shift;
 471       fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
 472       fstar.w[2] = P256.w[2];
 473       fstar.w[1] = P256.w[1];
 474       fstar.w[0] = P256.w[0];
 475       if (fstar.w[3] > onehalf128[ind - 1] ||
 476           (fstar.w[3] == onehalf128[ind - 1] &&
 477            (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
 478         // f2* > 1/2 and the result may be exact
 479         // Calculate f2* - 1/2
 480         tmp64 = fstar.w[3] - onehalf128[ind - 1];
 481         if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1]
 482             || (fstar.w[1] == ten2mk128[ind - 1].w[1]
 483                 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 484           // set the inexact flag
 485           *pfpsf |= INEXACT_EXCEPTION;
 486         }       // else the result is exact
 487       } else {  // the result is inexact; f2* <= 1/2
 488         // set the inexact flag
 489         *pfpsf |= INEXACT_EXCEPTION;
 490       }
 491     }
 492     // if the result was a midpoint, it was already rounded away from zero
 493     res.w[1] |= x_sign | 0x3040000000000000ull;
 494     BID_RETURN (res);
 495   } else {      // if ((q + exp) < 0) <=> q < -exp
 496     // the result is +0 or -0
 497     res.w[1] = x_sign | 0x3040000000000000ull;
 498     res.w[0] = 0x0000000000000000ull;
 499     *pfpsf |= INEXACT_EXCEPTION;
 500     BID_RETURN (res);
 501   }
 502   break;
 503 case ROUNDING_DOWN:
 504   if ((q + exp) > 0) {  // exp < 0 and 1 <= -exp < q
 505     // need to shift right -exp digits from the coefficient; exp will be 0
 506     ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
 507     // (number of digits to be chopped off)
 508     // chop off ind digits from the lower part of C1
 509     // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
 510     // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
 511     // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
 512     // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
 513     // tmp64 = C1.w[0];
 514     // if (ind <= 19) {
 515     //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
 516     // } else {
 517     //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
 518     //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
 519     // }
 520     // if (C1.w[0] < tmp64) C1.w[1]++;
 521     // if carry-out from C1.w[0], increment C1.w[1]
 522     // calculate C* and f*
 523     // C* is actually floor(C*) in this case
 524     // C* and f* need shifting and masking, as shown by
 525     // shiftright128[] and maskhigh128[]
 526     // 1 <= x <= 34
 527     // kx = 10^(-x) = ten2mk128[ind - 1]
 528     // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 529     // the approximation of 10^(-x) was rounded up to 118 bits
 530     __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
 531     if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
 532       res.w[1] = P256.w[3];
 533       res.w[0] = P256.w[2];
 534       // redundant fstar.w[3] = 0;
 535       // redundant fstar.w[2] = 0;
 536       // redundant fstar.w[1] = P256.w[1];
 537       // redundant fstar.w[0] = P256.w[0];
 538       // fraction f* > 10^(-x) <=> inexact
 539       // f* is in the right position to be compared with
 540       // 10^(-x) from ten2mk128[]
 541       if ((P256.w[1] > ten2mk128[ind - 1].w[1])
 542           || (P256.w[1] == ten2mk128[ind - 1].w[1]
 543               && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
 544         *pfpsf |= INEXACT_EXCEPTION;
 545         // if positive, the truncated value is already the correct result
 546         if (x_sign) {   // if negative
 547           if (++res.w[0] == 0) {
 548             res.w[1]++;
 549           }
 550         }
 551       }
 552     } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 553       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
 554       res.w[1] = (P256.w[3] >> shift);
 555       res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
 556       // redundant fstar.w[3] = 0;
 557       fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
 558       fstar.w[1] = P256.w[1];
 559       fstar.w[0] = P256.w[0];
 560       // fraction f* > 10^(-x) <=> inexact
 561       // f* is in the right position to be compared with
 562       // 10^(-x) from ten2mk128[]
 563       if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
 564           (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
 565            fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 566         *pfpsf |= INEXACT_EXCEPTION;
 567         // if positive, the truncated value is already the correct result
 568         if (x_sign) {   // if negative
 569           if (++res.w[0] == 0) {
 570             res.w[1]++;
 571           }
 572         }
 573       }
 574     } else {    // 22 <= ind - 1 <= 33
 575       shift = shiftright128[ind - 1] - 64;      // 2 <= shift <= 38
 576       res.w[1] = 0;
 577       res.w[0] = P256.w[3] >> shift;
 578       fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
 579       fstar.w[2] = P256.w[2];
 580       fstar.w[1] = P256.w[1];
 581       fstar.w[0] = P256.w[0];
 582       // fraction f* > 10^(-x) <=> inexact
 583       // f* is in the right position to be compared with
 584       // 10^(-x) from ten2mk128[]
 585       if (fstar.w[3] || fstar.w[2]
 586           || fstar.w[1] > ten2mk128[ind - 1].w[1]
 587           || (fstar.w[1] == ten2mk128[ind - 1].w[1]
 588               && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 589         *pfpsf |= INEXACT_EXCEPTION;
 590         // if positive, the truncated value is already the correct result
 591         if (x_sign) {   // if negative
 592           if (++res.w[0] == 0) {
 593             res.w[1]++;
 594           }
 595         }
 596       }
 597     }
 598     res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
 599     BID_RETURN (res);
 600   } else {      // if exp < 0 and q + exp <= 0
 601     if (x_sign) {       // negative rounds down to -1.0
 602       res.w[1] = 0xb040000000000000ull;
 603       res.w[0] = 0x0000000000000001ull;
 604     } else {    // positive rpunds down to +0.0
 605       res.w[1] = 0x3040000000000000ull;
 606       res.w[0] = 0x0000000000000000ull;
 607     }
 608     *pfpsf |= INEXACT_EXCEPTION;
 609     BID_RETURN (res);
 610   }
 611   break;
 612 case ROUNDING_UP:
 613   if ((q + exp) > 0) {  // exp < 0 and 1 <= -exp < q
 614     // need to shift right -exp digits from the coefficient; exp will be 0
 615     ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
 616     // (number of digits to be chopped off)
 617     // chop off ind digits from the lower part of C1
 618     // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
 619     // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
 620     // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
 621     // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
 622     // tmp64 = C1.w[0];
 623     // if (ind <= 19) {
 624     //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
 625     // } else {
 626     //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
 627     //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
 628     // }
 629     // if (C1.w[0] < tmp64) C1.w[1]++;
 630     // if carry-out from C1.w[0], increment C1.w[1]
 631     // calculate C* and f*
 632     // C* is actually floor(C*) in this case
 633     // C* and f* need shifting and masking, as shown by
 634     // shiftright128[] and maskhigh128[]
 635     // 1 <= x <= 34
 636     // kx = 10^(-x) = ten2mk128[ind - 1]
 637     // C* = C1 * 10^(-x)
 638     // the approximation of 10^(-x) was rounded up to 118 bits
 639     __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
 640     if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
 641       res.w[1] = P256.w[3];
 642       res.w[0] = P256.w[2];
 643       // redundant fstar.w[3] = 0;
 644       // redundant fstar.w[2] = 0;
 645       // redundant fstar.w[1] = P256.w[1];
 646       // redundant fstar.w[0] = P256.w[0];
 647       // fraction f* > 10^(-x) <=> inexact
 648       // f* is in the right position to be compared with
 649       // 10^(-x) from ten2mk128[]
 650       if ((P256.w[1] > ten2mk128[ind - 1].w[1])
 651           || (P256.w[1] == ten2mk128[ind - 1].w[1]
 652               && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
 653         *pfpsf |= INEXACT_EXCEPTION;
 654         // if negative, the truncated value is already the correct result
 655         if (!x_sign) {  // if positive
 656           if (++res.w[0] == 0) {
 657             res.w[1]++;
 658           }
 659         }
 660       }
 661     } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 662       shift = shiftright128[ind - 1];   // 3 <= shift <= 63
 663       res.w[1] = (P256.w[3] >> shift);
 664       res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
 665       // redundant fstar.w[3] = 0;
 666       fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
 667       fstar.w[1] = P256.w[1];
 668       fstar.w[0] = P256.w[0];
 669       // fraction f* > 10^(-x) <=> inexact
 670       // f* is in the right position to be compared with
 671       // 10^(-x) from ten2mk128[]
 672       if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
 673           (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
 674            fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 675         *pfpsf |= INEXACT_EXCEPTION;
 676         // if negative, the truncated value is already the correct result
 677         if (!x_sign) {  // if positive
 678           if (++res.w[0] == 0) {
 679             res.w[1]++;
 680           }
 681         }
 682       }
 683     } else {    // 22 <= ind - 1 <= 33
 684       shift = shiftright128[ind - 1] - 64;      // 2 <= shift <= 38
 685       res.w[1] = 0;
 686       res.w[0] = P256.w[3] >> shift;
 687       fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
 688       fstar.w[2] = P256.w[2];
 689       fstar.w[1] = P256.w[1];
 690       fstar.w[0] = P256.w[0];
 691       // fraction f* > 10^(-x) <=> inexact
 692       // f* is in the right position to be compared with
 693       // 10^(-x) from ten2mk128[]
 694       if (fstar.w[3] || fstar.w[2]
 695           || fstar.w[1] > ten2mk128[ind - 1].w[1]
 696           || (fstar.w[1] == ten2mk128[ind - 1].w[1]
 697               && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 698         *pfpsf |= INEXACT_EXCEPTION;
 699         // if negative, the truncated value is already the correct result
 700         if (!x_sign) {  // if positive
 701           if (++res.w[0] == 0) {
 702             res.w[1]++;
 703           }
 704         }
 705       }
 706     }
 707     res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
 708     BID_RETURN (res);
 709   } else {      // if exp < 0 and q + exp <= 0
 710     if (x_sign) {       // negative rounds up to -0.0
 711       res.w[1] = 0xb040000000000000ull;
 712       res.w[0] = 0x0000000000000000ull;
 713     } else {    // positive rpunds up to +1.0
 714       res.w[1] = 0x3040000000000000ull;
 715       res.w[0] = 0x0000000000000001ull;
 716     }
 717     *pfpsf |= INEXACT_EXCEPTION;
 718     BID_RETURN (res);
 719   }
 720   break;
 721 case ROUNDING_TO_ZERO:
 722   if ((q + exp) > 0) {  // exp < 0 and 1 <= -exp < q
 723     // need to shift right -exp digits from the coefficient; exp will be 0
 724     ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
 725     // (number of digits to be chopped off)
 726     // chop off ind digits from the lower part of C1
 727     // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
 728     // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
 729     // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
 730     // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
 731     //tmp64 = C1.w[0];
 732     // if (ind <= 19) {
 733     //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
 734     // } else {
 735     //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
 736     //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
 737     // }
 738     // if (C1.w[0] < tmp64) C1.w[1]++;
 739     // if carry-out from C1.w[0], increment C1.w[1]
 740     // calculate C* and f*
 741     // C* is actually floor(C*) in this case
 742     // C* and f* need shifting and masking, as shown by
 743     // shiftright128[] and maskhigh128[]
 744     // 1 <= x <= 34
 745     // kx = 10^(-x) = ten2mk128[ind - 1]
 746     // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 747     // the approximation of 10^(-x) was rounded up to 118 bits
 748     __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
 749     if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
 750       res.w[1] = P256.w[3];
 751       res.w[0] = P256.w[2];
 752       // redundant fstar.w[3] = 0;
 753       // redundant fstar.w[2] = 0;
 754       // redundant fstar.w[1] = P256.w[1];
 755       // redundant fstar.w[0] = P256.w[0];
 756       // fraction f* > 10^(-x) <=> inexact
 757       // f* is in the right position to be compared with
 758       // 10^(-x) from ten2mk128[]
 759       if ((P256.w[1] > ten2mk128[ind - 1].w[1])
 760           || (P256.w[1] == ten2mk128[ind - 1].w[1]
 761               && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
 762         *pfpsf |= INEXACT_EXCEPTION;
 763       }
 764     } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 765       shift = shiftright128[ind - 1];   // 3 <= shift <= 63
 766       res.w[1] = (P256.w[3] >> shift);
 767       res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
 768       // redundant fstar.w[3] = 0;
 769       fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
 770       fstar.w[1] = P256.w[1];
 771       fstar.w[0] = P256.w[0];
 772       // fraction f* > 10^(-x) <=> inexact
 773       // f* is in the right position to be compared with
 774       // 10^(-x) from ten2mk128[]
 775       if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
 776           (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
 777            fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 778         *pfpsf |= INEXACT_EXCEPTION;
 779       }
 780     } else {    // 22 <= ind - 1 <= 33
 781       shift = shiftright128[ind - 1] - 64;      // 2 <= shift <= 38
 782       res.w[1] = 0;
 783       res.w[0] = P256.w[3] >> shift;
 784       fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
 785       fstar.w[2] = P256.w[2];
 786       fstar.w[1] = P256.w[1];
 787       fstar.w[0] = P256.w[0];
 788       // fraction f* > 10^(-x) <=> inexact
 789       // f* is in the right position to be compared with
 790       // 10^(-x) from ten2mk128[]
 791       if (fstar.w[3] || fstar.w[2]
 792           || fstar.w[1] > ten2mk128[ind - 1].w[1]
 793           || (fstar.w[1] == ten2mk128[ind - 1].w[1]
 794               && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
 795         *pfpsf |= INEXACT_EXCEPTION;
 796       }
 797     }
 798     res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
 799     BID_RETURN (res);
 800   } else {      // if exp < 0 and q + exp <= 0 the result is +0 or -0
 801     res.w[1] = x_sign | 0x3040000000000000ull;
 802     res.w[0] = 0x0000000000000000ull;
 803     *pfpsf |= INEXACT_EXCEPTION;
 804     BID_RETURN (res);
 805   }
 806   break;
 807 }
 808
 809 BID_RETURN (res);
 810 }
 811
 812 /*****************************************************************************
 813  *  BID128_round_integral_nearest_even
 814  ****************************************************************************/
 815
 816 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_nearest_even, x)
 817
 818      UINT128 res;
 819      UINT64 x_sign;
 820      UINT64 x_exp;
 821      int exp;                   // unbiased exponent
 822   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
 823      UINT64 tmp64;
 824      BID_UI64DOUBLE tmp1;
 825      unsigned int x_nr_bits;
 826      int q, ind, shift;
 827      UINT128 C1;
 828   // UINT128 res is C* at first - represents up to 34 decimal digits ~ 113 bits
 829      UINT256 fstar;
 830      UINT256 P256;
 831
 832   // check for NaN or Infinity
 833 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
 834     // x is special
 835 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
 836   // if x = NaN, then res = Q (x)
 837   // check first for non-canonical NaN payload
 838   if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
 839       (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
 840        (x.w[0] > 0x38c15b09ffffffffull))) {
 841     x.w[1] = x.w[1] & 0xffffc00000000000ull;
 842     x.w[0] = 0x0ull;
 843   }
 844   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
 845     // set invalid flag
 846     *pfpsf |= INVALID_EXCEPTION;
 847     // return quiet (x)
 848     res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out also G[6]-G[16]
 849     res.w[0] = x.w[0];
 850   } else {      // x is QNaN
 851     // return x
 852     res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out G[6]-G[16]
 853     res.w[0] = x.w[0];
 854   }
 855   BID_RETURN (res)
 856 } else {        // x is not a NaN, so it must be infinity
 857   if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
 858     // return +inf
 859     res.w[1] = 0x7800000000000000ull;
 860     res.w[0] = 0x0000000000000000ull;
 861   } else {      // x is -inf
 862     // return -inf
 863     res.w[1] = 0xf800000000000000ull;
 864     res.w[0] = 0x0000000000000000ull;
 865   }
 866   BID_RETURN (res);
 867 }
 868 }
 869   // unpack x
 870 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
 871 C1.w[1] = x.w[1] & MASK_COEFF;
 872 C1.w[0] = x.w[0];
 873
 874   // check for non-canonical values (treated as zero)
 875 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
 876   // non-canonical
 877   x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
 878   C1.w[1] = 0;  // significand high
 879   C1.w[0] = 0;  // significand low
 880 } else {        // G0_G1 != 11
 881   x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
 882   if (C1.w[1] > 0x0001ed09bead87c0ull ||
 883       (C1.w[1] == 0x0001ed09bead87c0ull
 884        && C1.w[0] > 0x378d8e63ffffffffull)) {
 885     // x is non-canonical if coefficient is larger than 10^34 -1
 886     C1.w[1] = 0;
 887     C1.w[0] = 0;
 888   } else {      // canonical
 889     ;
 890   }
 891 }
 892
 893   // test for input equal to zero
 894 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
 895   // x is 0
 896   // return 0 preserving the sign bit and the preferred exponent
 897   // of MAX(Q(x), 0)
 898   if (x_exp <= (0x1820ull << 49)) {
 899     res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
 900   } else {
 901     res.w[1] = x_sign | x_exp;
 902   }
 903   res.w[0] = 0x0000000000000000ull;
 904   BID_RETURN (res);
 905 }
 906   // x is not special and is not zero
 907
 908   // if (exp <= -(p+1)) return 0
 909 if (x_exp <= 0x2ffa000000000000ull) {   // 0x2ffa000000000000ull == -35
 910   res.w[1] = x_sign | 0x3040000000000000ull;
 911   res.w[0] = 0x0000000000000000ull;
 912   BID_RETURN (res);
 913 }
 914   // q = nr. of decimal digits in x
 915   //  determine first the nr. of bits in x
 916 if (C1.w[1] == 0) {
 917   if (C1.w[0] >= 0x0020000000000000ull) {       // x >= 2^53
 918     // split the 64-bit value in two 32-bit halves to avoid rounding errors
 919     if (C1.w[0] >= 0x0000000100000000ull) {     // x >= 2^32
 920       tmp1.d = (double) (C1.w[0] >> 32);        // exact conversion
 921       x_nr_bits =
 922         33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 923     } else {    // x < 2^32
 924       tmp1.d = (double) (C1.w[0]);      // exact conversion
 925       x_nr_bits =
 926         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 927     }
 928   } else {      // if x < 2^53
 929     tmp1.d = (double) C1.w[0];  // exact conversion
 930     x_nr_bits =
 931       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 932   }
 933 } else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
 934   tmp1.d = (double) C1.w[1];    // exact conversion
 935   x_nr_bits =
 936     65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 937 }
 938
 939 q = nr_digits[x_nr_bits - 1].digits;
 940 if (q == 0) {
 941   q = nr_digits[x_nr_bits - 1].digits1;
 942   if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
 943       || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
 944           C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
 945     q++;
 946 }
 947 exp = (x_exp >> 49) - 6176;
 948 if (exp >= 0) { // -exp <= 0
 949   // the argument is an integer already
 950   res.w[1] = x.w[1];
 951   res.w[0] = x.w[0];
 952   BID_RETURN (res);
 953 } else if ((q + exp) >= 0) {    // exp < 0 and 1 <= -exp <= q
 954   // need to shift right -exp digits from the coefficient; the exp will be 0
 955   ind = -exp;   // 1 <= ind <= 34; ind is a synonym for 'x'
 956   // chop off ind digits from the lower part of C1
 957   // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
 958   tmp64 = C1.w[0];
 959   if (ind <= 19) {
 960     C1.w[0] = C1.w[0] + midpoint64[ind - 1];
 961   } else {
 962     C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
 963     C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
 964   }
 965   if (C1.w[0] < tmp64)
 966     C1.w[1]++;
 967   // calculate C* and f*
 968   // C* is actually floor(C*) in this case
 969   // C* and f* need shifting and masking, as shown by
 970   // shiftright128[] and maskhigh128[]
 971   // 1 <= x <= 34
 972   // kx = 10^(-x) = ten2mk128[ind - 1]
 973   // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 974   // the approximation of 10^(-x) was rounded up to 118 bits
 975   __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
 976   // determine the value of res and fstar
 977   if (ind - 1 <= 2) {   // 0 <= ind - 1 <= 2 => shift = 0
 978     // redundant shift = shiftright128[ind - 1]; // shift = 0
 979     res.w[1] = P256.w[3];
 980     res.w[0] = P256.w[2];
 981     // redundant fstar.w[3] = 0;
 982     // redundant fstar.w[2] = 0;
 983     // redundant fstar.w[1] = P256.w[1];
 984     // redundant fstar.w[0] = P256.w[0];
 985     // fraction f* < 10^(-x) <=> midpoint
 986     // f* is in the right position to be compared with
 987     // 10^(-x) from ten2mk128[]
 988     // if 0 < fstar < 10^(-x), subtract 1 if odd (for rounding to even)
 989     if ((res.w[0] & 0x0000000000000001ull) &&   // is result odd?
 990         ((P256.w[1] < (ten2mk128[ind - 1].w[1]))
 991          || ((P256.w[1] == ten2mk128[ind - 1].w[1])
 992              && (P256.w[0] < ten2mk128[ind - 1].w[0])))) {
 993       // subract 1 to make even
 994       if (res.w[0]-- == 0) {
 995         res.w[1]--;
 996       }
 997     }
 998   } else if (ind - 1 <= 21) {   // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
 999     shift = shiftright128[ind - 1];     // 3 <= shift <= 63
1000     res.w[1] = (P256.w[3] >> shift);
1001     res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1002     // redundant fstar.w[3] = 0;
1003     fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1004     fstar.w[1] = P256.w[1];
1005     fstar.w[0] = P256.w[0];
1006     // fraction f* < 10^(-x) <=> midpoint
1007     // f* is in the right position to be compared with
1008     // 10^(-x) from ten2mk128[]
1009     if ((res.w[0] & 0x0000000000000001ull) &&   // is result odd?
1010         fstar.w[2] == 0 && (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
1011                             (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
1012                              fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
1013       // subract 1 to make even
1014       if (res.w[0]-- == 0) {
1015         res.w[1]--;
1016       }
1017     }
1018   } else {      // 22 <= ind - 1 <= 33
1019     shift = shiftright128[ind - 1] - 64;        // 2 <= shift <= 38
1020     res.w[1] = 0;
1021     res.w[0] = P256.w[3] >> shift;
1022     fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1023     fstar.w[2] = P256.w[2];
1024     fstar.w[1] = P256.w[1];
1025     fstar.w[0] = P256.w[0];
1026     // fraction f* < 10^(-x) <=> midpoint
1027     // f* is in the right position to be compared with
1028     // 10^(-x) from ten2mk128[]
1029     if ((res.w[0] & 0x0000000000000001ull) &&   // is result odd?
1030         fstar.w[3] == 0 && fstar.w[2] == 0
1031         && (fstar.w[1] < ten2mk128[ind - 1].w[1]
1032             || (fstar.w[1] == ten2mk128[ind - 1].w[1]
1033                 && fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
1034       // subract 1 to make even
1035       if (res.w[0]-- == 0) {
1036         res.w[1]--;
1037       }
1038     }
1039   }
1040   res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1041   BID_RETURN (res);
1042 } else {        // if ((q + exp) < 0) <=> q < -exp
1043   // the result is +0 or -0
1044   res.w[1] = x_sign | 0x3040000000000000ull;
1045   res.w[0] = 0x0000000000000000ull;
1046   BID_RETURN (res);
1047 }
1048 }
1049
1050 /*****************************************************************************
1051  *  BID128_round_integral_negative
1052  ****************************************************************************/
1053
1054 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_negative, x)
1055
1056      UINT128 res;
1057      UINT64 x_sign;
1058      UINT64 x_exp;
1059      int exp;                   // unbiased exponent
1060   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
1061   // (all are UINT64)
1062      BID_UI64DOUBLE tmp1;
1063      unsigned int x_nr_bits;
1064      int q, ind, shift;
1065      UINT128 C1;
1066   // UINT128 res is C* at first - represents up to 34 decimal digits ~
1067   // 113 bits
1068      UINT256 fstar;
1069      UINT256 P256;
1070
1071   // check for NaN or Infinity
1072 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1073     // x is special
1074 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1075   // if x = NaN, then res = Q (x)
1076   // check first for non-canonical NaN payload
1077   if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1078       (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1079        (x.w[0] > 0x38c15b09ffffffffull))) {
1080     x.w[1] = x.w[1] & 0xffffc00000000000ull;
1081     x.w[0] = 0x0ull;
1082   }
1083   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1084     // set invalid flag
1085     *pfpsf |= INVALID_EXCEPTION;
1086     // return quiet (x)
1087     res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out also G[6]-G[16]
1088     res.w[0] = x.w[0];
1089   } else {      // x is QNaN
1090     // return x
1091     res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out G[6]-G[16]
1092     res.w[0] = x.w[0];
1093   }
1094   BID_RETURN (res)
1095 } else {        // x is not a NaN, so it must be infinity
1096   if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1097     // return +inf
1098     res.w[1] = 0x7800000000000000ull;
1099     res.w[0] = 0x0000000000000000ull;
1100   } else {      // x is -inf
1101     // return -inf
1102     res.w[1] = 0xf800000000000000ull;
1103     res.w[0] = 0x0000000000000000ull;
1104   }
1105   BID_RETURN (res);
1106 }
1107 }
1108   // unpack x
1109 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1110 C1.w[1] = x.w[1] & MASK_COEFF;
1111 C1.w[0] = x.w[0];
1112
1113   // check for non-canonical values (treated as zero)
1114 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
1115   // non-canonical
1116   x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
1117   C1.w[1] = 0;  // significand high
1118   C1.w[0] = 0;  // significand low
1119 } else {        // G0_G1 != 11
1120   x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
1121   if (C1.w[1] > 0x0001ed09bead87c0ull ||
1122       (C1.w[1] == 0x0001ed09bead87c0ull
1123        && C1.w[0] > 0x378d8e63ffffffffull)) {
1124     // x is non-canonical if coefficient is larger than 10^34 -1
1125     C1.w[1] = 0;
1126     C1.w[0] = 0;
1127   } else {      // canonical
1128     ;
1129   }
1130 }
1131
1132   // test for input equal to zero
1133 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1134   // x is 0
1135   // return 0 preserving the sign bit and the preferred exponent
1136   // of MAX(Q(x), 0)
1137   if (x_exp <= (0x1820ull << 49)) {
1138     res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1139   } else {
1140     res.w[1] = x_sign | x_exp;
1141   }
1142   res.w[0] = 0x0000000000000000ull;
1143   BID_RETURN (res);
1144 }
1145   // x is not special and is not zero
1146
1147   // if (exp <= -p) return -1.0 or +0.0
1148 if (x_exp <= 0x2ffc000000000000ull) {   // 0x2ffc000000000000ull == -34
1149   if (x_sign) {
1150     // if negative, return negative 1, because we know the coefficient
1151     // is non-zero (would have been caught above)
1152     res.w[1] = 0xb040000000000000ull;
1153     res.w[0] = 0x0000000000000001ull;
1154   } else {
1155     // if positive, return positive 0, because we know coefficient is
1156     // non-zero (would have been caught above)
1157     res.w[1] = 0x3040000000000000ull;
1158     res.w[0] = 0x0000000000000000ull;
1159   }
1160   BID_RETURN (res);
1161 }
1162   // q = nr. of decimal digits in x
1163   // determine first the nr. of bits in x
1164 if (C1.w[1] == 0) {
1165   if (C1.w[0] >= 0x0020000000000000ull) {       // x >= 2^53
1166     // split the 64-bit value in two 32-bit halves to avoid rounding errors
1167     if (C1.w[0] >= 0x0000000100000000ull) {     // x >= 2^32
1168       tmp1.d = (double) (C1.w[0] >> 32);        // exact conversion
1169       x_nr_bits =
1170         33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1171     } else {    // x < 2^32
1172       tmp1.d = (double) (C1.w[0]);      // exact conversion
1173       x_nr_bits =
1174         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1175     }
1176   } else {      // if x < 2^53
1177     tmp1.d = (double) C1.w[0];  // exact conversion
1178     x_nr_bits =
1179       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1180   }
1181 } else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1182   tmp1.d = (double) C1.w[1];    // exact conversion
1183   x_nr_bits =
1184     65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1185 }
1186
1187 q = nr_digits[x_nr_bits - 1].digits;
1188 if (q == 0) {
1189   q = nr_digits[x_nr_bits - 1].digits1;
1190   if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1191       (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1192        C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1193     q++;
1194 }
1195 exp = (x_exp >> 49) - 6176;
1196 if (exp >= 0) { // -exp <= 0
1197   // the argument is an integer already
1198   res.w[1] = x.w[1];
1199   res.w[0] = x.w[0];
1200   BID_RETURN (res);
1201 } else if ((q + exp) > 0) {     // exp < 0 and 1 <= -exp < q
1202   // need to shift right -exp digits from the coefficient; the exp will be 0
1203   ind = -exp;   // 1 <= ind <= 34; ind is a synonym for 'x'
1204   // (number of digits to be chopped off)
1205   // chop off ind digits from the lower part of C1
1206   // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
1207   // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
1208   // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
1209   // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
1210   //tmp64 = C1.w[0];
1211   // if (ind <= 19) {
1212   //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1213   // } else {
1214   //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1215   //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1216   // }
1217   // if (C1.w[0] < tmp64) C1.w[1]++;
1218   // if carry-out from C1.w[0], increment C1.w[1]
1219   // calculate C* and f*
1220   // C* is actually floor(C*) in this case
1221   // C* and f* need shifting and masking, as shown by
1222   // shiftright128[] and maskhigh128[]
1223   // 1 <= x <= 34
1224   // kx = 10^(-x) = ten2mk128[ind - 1]
1225   // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1226   // the approximation of 10^(-x) was rounded up to 118 bits
1227   __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1228   if (ind - 1 <= 2) {   // 0 <= ind - 1 <= 2 => shift = 0
1229     res.w[1] = P256.w[3];
1230     res.w[0] = P256.w[2];
1231     // if positive, the truncated value is already the correct result
1232     if (x_sign) {       // if negative
1233       // redundant fstar.w[3] = 0;
1234       // redundant fstar.w[2] = 0;
1235       // redundant fstar.w[1] = P256.w[1];
1236       // redundant fstar.w[0] = P256.w[0];
1237       // fraction f* > 10^(-x) <=> inexact
1238       // f* is in the right position to be compared with
1239       // 10^(-x) from ten2mk128[]
1240       if ((P256.w[1] > ten2mk128[ind - 1].w[1])
1241           || (P256.w[1] == ten2mk128[ind - 1].w[1]
1242               && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
1243         if (++res.w[0] == 0) {
1244           res.w[1]++;
1245         }
1246       }
1247     }
1248   } else if (ind - 1 <= 21) {   // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1249     shift = shiftright128[ind - 1];     // 0 <= shift <= 102
1250     res.w[1] = (P256.w[3] >> shift);
1251     res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1252     // if positive, the truncated value is already the correct result
1253     if (x_sign) {       // if negative
1254       // redundant fstar.w[3] = 0;
1255       fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1256       fstar.w[1] = P256.w[1];
1257       fstar.w[0] = P256.w[0];
1258       // fraction f* > 10^(-x) <=> inexact
1259       // f* is in the right position to be compared with
1260       // 10^(-x) from ten2mk128[]
1261       if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
1262           (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
1263            fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1264         if (++res.w[0] == 0) {
1265           res.w[1]++;
1266         }
1267       }
1268     }
1269   } else {      // 22 <= ind - 1 <= 33
1270     shift = shiftright128[ind - 1] - 64;        // 2 <= shift <= 38
1271     res.w[1] = 0;
1272     res.w[0] = P256.w[3] >> shift;
1273     // if positive, the truncated value is already the correct result
1274     if (x_sign) {       // if negative
1275       fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1276       fstar.w[2] = P256.w[2];
1277       fstar.w[1] = P256.w[1];
1278       fstar.w[0] = P256.w[0];
1279       // fraction f* > 10^(-x) <=> inexact
1280       // f* is in the right position to be compared with
1281       // 10^(-x) from ten2mk128[]
1282       if (fstar.w[3] || fstar.w[2]
1283           || fstar.w[1] > ten2mk128[ind - 1].w[1]
1284           || (fstar.w[1] == ten2mk128[ind - 1].w[1]
1285               && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1286         if (++res.w[0] == 0) {
1287           res.w[1]++;
1288         }
1289       }
1290     }
1291   }
1292   res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1293   BID_RETURN (res);
1294 } else {        // if exp < 0 and q + exp <= 0
1295   if (x_sign) { // negative rounds down to -1.0
1296     res.w[1] = 0xb040000000000000ull;
1297     res.w[0] = 0x0000000000000001ull;
1298   } else {      // positive rpunds down to +0.0
1299     res.w[1] = 0x3040000000000000ull;
1300     res.w[0] = 0x0000000000000000ull;
1301   }
1302   BID_RETURN (res);
1303 }
1304 }
1305
1306 /*****************************************************************************
1307  *  BID128_round_integral_positive
1308  ****************************************************************************/
1309
1310 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_positive, x)
1311
1312      UINT128 res;
1313      UINT64 x_sign;
1314      UINT64 x_exp;
1315      int exp;                   // unbiased exponent
1316   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
1317   // (all are UINT64)
1318      BID_UI64DOUBLE tmp1;
1319      unsigned int x_nr_bits;
1320      int q, ind, shift;
1321      UINT128 C1;
1322   // UINT128 res is C* at first - represents up to 34 decimal digits ~
1323   // 113 bits
1324      UINT256 fstar;
1325      UINT256 P256;
1326
1327   // check for NaN or Infinity
1328 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1329     // x is special
1330 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1331   // if x = NaN, then res = Q (x)
1332   // check first for non-canonical NaN payload
1333   if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1334       (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1335        (x.w[0] > 0x38c15b09ffffffffull))) {
1336     x.w[1] = x.w[1] & 0xffffc00000000000ull;
1337     x.w[0] = 0x0ull;
1338   }
1339   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1340     // set invalid flag
1341     *pfpsf |= INVALID_EXCEPTION;
1342     // return quiet (x)
1343     res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out also G[6]-G[16]
1344     res.w[0] = x.w[0];
1345   } else {      // x is QNaN
1346     // return x
1347     res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out G[6]-G[16]
1348     res.w[0] = x.w[0];
1349   }
1350   BID_RETURN (res)
1351 } else {        // x is not a NaN, so it must be infinity
1352   if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1353     // return +inf
1354     res.w[1] = 0x7800000000000000ull;
1355     res.w[0] = 0x0000000000000000ull;
1356   } else {      // x is -inf
1357     // return -inf
1358     res.w[1] = 0xf800000000000000ull;
1359     res.w[0] = 0x0000000000000000ull;
1360   }
1361   BID_RETURN (res);
1362 }
1363 }
1364   // unpack x
1365 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1366 C1.w[1] = x.w[1] & MASK_COEFF;
1367 C1.w[0] = x.w[0];
1368
1369   // check for non-canonical values (treated as zero)
1370 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
1371   // non-canonical
1372   x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
1373   C1.w[1] = 0;  // significand high
1374   C1.w[0] = 0;  // significand low
1375 } else {        // G0_G1 != 11
1376   x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
1377   if (C1.w[1] > 0x0001ed09bead87c0ull ||
1378       (C1.w[1] == 0x0001ed09bead87c0ull
1379        && C1.w[0] > 0x378d8e63ffffffffull)) {
1380     // x is non-canonical if coefficient is larger than 10^34 -1
1381     C1.w[1] = 0;
1382     C1.w[0] = 0;
1383   } else {      // canonical
1384     ;
1385   }
1386 }
1387
1388   // test for input equal to zero
1389 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1390   // x is 0
1391   // return 0 preserving the sign bit and the preferred exponent
1392   // of MAX(Q(x), 0)
1393   if (x_exp <= (0x1820ull << 49)) {
1394     res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1395   } else {
1396     res.w[1] = x_sign | x_exp;
1397   }
1398   res.w[0] = 0x0000000000000000ull;
1399   BID_RETURN (res);
1400 }
1401   // x is not special and is not zero
1402
1403   // if (exp <= -p) return -0.0 or +1.0
1404 if (x_exp <= 0x2ffc000000000000ull) {   // 0x2ffc000000000000ull == -34
1405   if (x_sign) {
1406     // if negative, return negative 0, because we know the coefficient
1407     // is non-zero (would have been caught above)
1408     res.w[1] = 0xb040000000000000ull;
1409     res.w[0] = 0x0000000000000000ull;
1410   } else {
1411     // if positive, return positive 1, because we know coefficient is
1412     // non-zero (would have been caught above)
1413     res.w[1] = 0x3040000000000000ull;
1414     res.w[0] = 0x0000000000000001ull;
1415   }
1416   BID_RETURN (res);
1417 }
1418   // q = nr. of decimal digits in x
1419   // determine first the nr. of bits in x
1420 if (C1.w[1] == 0) {
1421   if (C1.w[0] >= 0x0020000000000000ull) {       // x >= 2^53
1422     // split 64-bit value in two 32-bit halves to avoid rounding errors
1423     if (C1.w[0] >= 0x0000000100000000ull) {     // x >= 2^32
1424       tmp1.d = (double) (C1.w[0] >> 32);        // exact conversion
1425       x_nr_bits =
1426         33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1427     } else {    // x < 2^32
1428       tmp1.d = (double) (C1.w[0]);      // exact conversion
1429       x_nr_bits =
1430         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1431     }
1432   } else {      // if x < 2^53
1433     tmp1.d = (double) C1.w[0];  // exact conversion
1434     x_nr_bits =
1435       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1436   }
1437 } else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1438   tmp1.d = (double) C1.w[1];    // exact conversion
1439   x_nr_bits =
1440     65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1441 }
1442
1443 q = nr_digits[x_nr_bits - 1].digits;
1444 if (q == 0) {
1445   q = nr_digits[x_nr_bits - 1].digits1;
1446   if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1447       (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1448        C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1449     q++;
1450 }
1451 exp = (x_exp >> 49) - 6176;
1452 if (exp >= 0) { // -exp <= 0
1453   // the argument is an integer already
1454   res.w[1] = x.w[1];
1455   res.w[0] = x.w[0];
1456   BID_RETURN (res);
1457 } else if ((q + exp) > 0) {     // exp < 0 and 1 <= -exp < q
1458   // need to shift right -exp digits from the coefficient; exp will be 0
1459   ind = -exp;   // 1 <= ind <= 34; ind is a synonym for 'x'
1460   // (number of digits to be chopped off)
1461   // chop off ind digits from the lower part of C1
1462   // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
1463   // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
1464   // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
1465   // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
1466   // tmp64 = C1.w[0];
1467   // if (ind <= 19) {
1468   //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1469   // } else {
1470   //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1471   //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1472   // }
1473   // if (C1.w[0] < tmp64) C1.w[1]++;
1474   // if carry-out from C1.w[0], increment C1.w[1]
1475   // calculate C* and f*
1476   // C* is actually floor(C*) in this case
1477   // C* and f* need shifting and masking, as shown by
1478   // shiftright128[] and maskhigh128[]
1479   // 1 <= x <= 34
1480   // kx = 10^(-x) = ten2mk128[ind - 1]
1481   // C* = C1 * 10^(-x)
1482   // the approximation of 10^(-x) was rounded up to 118 bits
1483   __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1484   if (ind - 1 <= 2) {   // 0 <= ind - 1 <= 2 => shift = 0
1485     res.w[1] = P256.w[3];
1486     res.w[0] = P256.w[2];
1487     // if negative, the truncated value is already the correct result
1488     if (!x_sign) {      // if positive
1489       // redundant fstar.w[3] = 0;
1490       // redundant fstar.w[2] = 0;
1491       // redundant fstar.w[1] = P256.w[1];
1492       // redundant fstar.w[0] = P256.w[0];
1493       // fraction f* > 10^(-x) <=> inexact
1494       // f* is in the right position to be compared with
1495       // 10^(-x) from ten2mk128[]
1496       if ((P256.w[1] > ten2mk128[ind - 1].w[1])
1497           || (P256.w[1] == ten2mk128[ind - 1].w[1]
1498               && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
1499         if (++res.w[0] == 0) {
1500           res.w[1]++;
1501         }
1502       }
1503     }
1504   } else if (ind - 1 <= 21) {   // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1505     shift = shiftright128[ind - 1];     // 3 <= shift <= 63
1506     res.w[1] = (P256.w[3] >> shift);
1507     res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1508     // if negative, the truncated value is already the correct result
1509     if (!x_sign) {      // if positive
1510       // redundant fstar.w[3] = 0;
1511       fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1512       fstar.w[1] = P256.w[1];
1513       fstar.w[0] = P256.w[0];
1514       // fraction f* > 10^(-x) <=> inexact
1515       // f* is in the right position to be compared with
1516       // 10^(-x) from ten2mk128[]
1517       if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
1518           (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
1519            fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1520         if (++res.w[0] == 0) {
1521           res.w[1]++;
1522         }
1523       }
1524     }
1525   } else {      // 22 <= ind - 1 <= 33
1526     shift = shiftright128[ind - 1] - 64;        // 2 <= shift <= 38
1527     res.w[1] = 0;
1528     res.w[0] = P256.w[3] >> shift;
1529     // if negative, the truncated value is already the correct result
1530     if (!x_sign) {      // if positive
1531       fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1532       fstar.w[2] = P256.w[2];
1533       fstar.w[1] = P256.w[1];
1534       fstar.w[0] = P256.w[0];
1535       // fraction f* > 10^(-x) <=> inexact
1536       // f* is in the right position to be compared with
1537       // 10^(-x) from ten2mk128[]
1538       if (fstar.w[3] || fstar.w[2]
1539           || fstar.w[1] > ten2mk128[ind - 1].w[1]
1540           || (fstar.w[1] == ten2mk128[ind - 1].w[1]
1541               && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1542         if (++res.w[0] == 0) {
1543           res.w[1]++;
1544         }
1545       }
1546     }
1547   }
1548   res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1549   BID_RETURN (res);
1550 } else {        // if exp < 0 and q + exp <= 0
1551   if (x_sign) { // negative rounds up to -0.0
1552     res.w[1] = 0xb040000000000000ull;
1553     res.w[0] = 0x0000000000000000ull;
1554   } else {      // positive rpunds up to +1.0
1555     res.w[1] = 0x3040000000000000ull;
1556     res.w[0] = 0x0000000000000001ull;
1557   }
1558   BID_RETURN (res);
1559 }
1560 }
1561
1562 /*****************************************************************************
1563  *  BID128_round_integral_zero
1564  ****************************************************************************/
1565
1566 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_zero, x)
1567
1568      UINT128 res;
1569      UINT64 x_sign;
1570      UINT64 x_exp;
1571      int exp;                   // unbiased exponent
1572   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
1573   // (all are UINT64)
1574      BID_UI64DOUBLE tmp1;
1575      unsigned int x_nr_bits;
1576      int q, ind, shift;
1577      UINT128 C1;
1578   // UINT128 res is C* at first - represents up to 34 decimal digits ~
1579   // 113 bits
1580      UINT256 P256;
1581
1582   // check for NaN or Infinity
1583 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1584     // x is special
1585 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1586   // if x = NaN, then res = Q (x)
1587   // check first for non-canonical NaN payload
1588   if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1589       (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1590        (x.w[0] > 0x38c15b09ffffffffull))) {
1591     x.w[1] = x.w[1] & 0xffffc00000000000ull;
1592     x.w[0] = 0x0ull;
1593   }
1594   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1595     // set invalid flag
1596     *pfpsf |= INVALID_EXCEPTION;
1597     // return quiet (x)
1598     res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out also G[6]-G[16]
1599     res.w[0] = x.w[0];
1600   } else {      // x is QNaN
1601     // return x
1602     res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out G[6]-G[16]
1603     res.w[0] = x.w[0];
1604   }
1605   BID_RETURN (res)
1606 } else {        // x is not a NaN, so it must be infinity
1607   if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1608     // return +inf
1609     res.w[1] = 0x7800000000000000ull;
1610     res.w[0] = 0x0000000000000000ull;
1611   } else {      // x is -inf
1612     // return -inf
1613     res.w[1] = 0xf800000000000000ull;
1614     res.w[0] = 0x0000000000000000ull;
1615   }
1616   BID_RETURN (res);
1617 }
1618 }
1619   // unpack x
1620 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1621 C1.w[1] = x.w[1] & MASK_COEFF;
1622 C1.w[0] = x.w[0];
1623
1624   // check for non-canonical values (treated as zero)
1625 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
1626   // non-canonical
1627   x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
1628   C1.w[1] = 0;  // significand high
1629   C1.w[0] = 0;  // significand low
1630 } else {        // G0_G1 != 11
1631   x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
1632   if (C1.w[1] > 0x0001ed09bead87c0ull ||
1633       (C1.w[1] == 0x0001ed09bead87c0ull
1634        && C1.w[0] > 0x378d8e63ffffffffull)) {
1635     // x is non-canonical if coefficient is larger than 10^34 -1
1636     C1.w[1] = 0;
1637     C1.w[0] = 0;
1638   } else {      // canonical
1639     ;
1640   }
1641 }
1642
1643   // test for input equal to zero
1644 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1645   // x is 0
1646   // return 0 preserving the sign bit and the preferred exponent
1647   // of MAX(Q(x), 0)
1648   if (x_exp <= (0x1820ull << 49)) {
1649     res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1650   } else {
1651     res.w[1] = x_sign | x_exp;
1652   }
1653   res.w[0] = 0x0000000000000000ull;
1654   BID_RETURN (res);
1655 }
1656   // x is not special and is not zero
1657
1658   // if (exp <= -p) return -0.0 or +0.0
1659 if (x_exp <= 0x2ffc000000000000ull) {   // 0x2ffc000000000000ull == -34
1660   res.w[1] = x_sign | 0x3040000000000000ull;
1661   res.w[0] = 0x0000000000000000ull;
1662   BID_RETURN (res);
1663 }
1664   // q = nr. of decimal digits in x
1665   // determine first the nr. of bits in x
1666 if (C1.w[1] == 0) {
1667   if (C1.w[0] >= 0x0020000000000000ull) {       // x >= 2^53
1668     // split the 64-bit value in two 32-bit halves to avoid rounding errors
1669     if (C1.w[0] >= 0x0000000100000000ull) {     // x >= 2^32
1670       tmp1.d = (double) (C1.w[0] >> 32);        // exact conversion
1671       x_nr_bits =
1672         33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1673     } else {    // x < 2^32
1674       tmp1.d = (double) (C1.w[0]);      // exact conversion
1675       x_nr_bits =
1676         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1677     }
1678   } else {      // if x < 2^53
1679     tmp1.d = (double) C1.w[0];  // exact conversion
1680     x_nr_bits =
1681       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1682   }
1683 } else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1684   tmp1.d = (double) C1.w[1];    // exact conversion
1685   x_nr_bits =
1686     65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1687 }
1688
1689 q = nr_digits[x_nr_bits - 1].digits;
1690 if (q == 0) {
1691   q = nr_digits[x_nr_bits - 1].digits1;
1692   if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1693       (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1694        C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1695     q++;
1696 }
1697 exp = (x_exp >> 49) - 6176;
1698 if (exp >= 0) { // -exp <= 0
1699   // the argument is an integer already
1700   res.w[1] = x.w[1];
1701   res.w[0] = x.w[0];
1702   BID_RETURN (res);
1703 } else if ((q + exp) > 0) {     // exp < 0 and 1 <= -exp < q
1704   // need to shift right -exp digits from the coefficient; the exp will be 0
1705   ind = -exp;   // 1 <= ind <= 34; ind is a synonym for 'x'
1706   // (number of digits to be chopped off)
1707   // chop off ind digits from the lower part of C1
1708   // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
1709   // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
1710   // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
1711   // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
1712   //tmp64 = C1.w[0];
1713   // if (ind <= 19) {
1714   //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1715   // } else {
1716   //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1717   //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1718   // }
1719   // if (C1.w[0] < tmp64) C1.w[1]++;
1720   // if carry-out from C1.w[0], increment C1.w[1]
1721   // calculate C* and f*
1722   // C* is actually floor(C*) in this case
1723   // C* and f* need shifting and masking, as shown by
1724   // shiftright128[] and maskhigh128[]
1725   // 1 <= x <= 34
1726   // kx = 10^(-x) = ten2mk128[ind - 1]
1727   // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1728   // the approximation of 10^(-x) was rounded up to 118 bits
1729   __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1730   if (ind - 1 <= 2) {   // 0 <= ind - 1 <= 2 => shift = 0
1731     res.w[1] = P256.w[3];
1732     res.w[0] = P256.w[2];
1733   } else if (ind - 1 <= 21) {   // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1734     shift = shiftright128[ind - 1];     // 3 <= shift <= 63
1735     res.w[1] = (P256.w[3] >> shift);
1736     res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1737   } else {      // 22 <= ind - 1 <= 33
1738     shift = shiftright128[ind - 1] - 64;        // 2 <= shift <= 38
1739     res.w[1] = 0;
1740     res.w[0] = P256.w[3] >> shift;
1741   }
1742   res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1743   BID_RETURN (res);
1744 } else {        // if exp < 0 and q + exp <= 0 the result is +0 or -0
1745   res.w[1] = x_sign | 0x3040000000000000ull;
1746   res.w[0] = 0x0000000000000000ull;
1747   BID_RETURN (res);
1748 }
1749 }
1750
1751 /*****************************************************************************
1752  *  BID128_round_integral_nearest_away
1753  ****************************************************************************/
1754
1755 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_nearest_away, x)
1756
1757      UINT128 res;
1758      UINT64 x_sign;
1759      UINT64 x_exp;
1760      int exp;                   // unbiased exponent
1761   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
1762   // (all are UINT64)
1763      UINT64 tmp64;
1764      BID_UI64DOUBLE tmp1;
1765      unsigned int x_nr_bits;
1766      int q, ind, shift;
1767      UINT128 C1;
1768   // UINT128 res is C* at first - represents up to 34 decimal digits ~
1769   // 113 bits
1770   // UINT256 fstar;
1771      UINT256 P256;
1772
1773   // check for NaN or Infinity
1774 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1775     // x is special
1776 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1777   // if x = NaN, then res = Q (x)
1778   // check first for non-canonical NaN payload
1779   if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1780       (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1781        (x.w[0] > 0x38c15b09ffffffffull))) {
1782     x.w[1] = x.w[1] & 0xffffc00000000000ull;
1783     x.w[0] = 0x0ull;
1784   }
1785   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1786     // set invalid flag
1787     *pfpsf |= INVALID_EXCEPTION;
1788     // return quiet (x)
1789     res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out also G[6]-G[16]
1790     res.w[0] = x.w[0];
1791   } else {      // x is QNaN
1792     // return x
1793     res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out G[6]-G[16]
1794     res.w[0] = x.w[0];
1795   }
1796   BID_RETURN (res)
1797 } else {        // x is not a NaN, so it must be infinity
1798   if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1799     // return +inf
1800     res.w[1] = 0x7800000000000000ull;
1801     res.w[0] = 0x0000000000000000ull;
1802   } else {      // x is -inf
1803     // return -inf
1804     res.w[1] = 0xf800000000000000ull;
1805     res.w[0] = 0x0000000000000000ull;
1806   }
1807   BID_RETURN (res);
1808 }
1809 }
1810   // unpack x
1811 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1812 C1.w[1] = x.w[1] & MASK_COEFF;
1813 C1.w[0] = x.w[0];
1814
1815   // check for non-canonical values (treated as zero)
1816 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
1817   // non-canonical
1818   x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
1819   C1.w[1] = 0;  // significand high
1820   C1.w[0] = 0;  // significand low
1821 } else {        // G0_G1 != 11
1822   x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
1823   if (C1.w[1] > 0x0001ed09bead87c0ull ||
1824       (C1.w[1] == 0x0001ed09bead87c0ull
1825        && C1.w[0] > 0x378d8e63ffffffffull)) {
1826     // x is non-canonical if coefficient is larger than 10^34 -1
1827     C1.w[1] = 0;
1828     C1.w[0] = 0;
1829   } else {      // canonical
1830     ;
1831   }
1832 }
1833
1834   // test for input equal to zero
1835 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1836   // x is 0
1837   // return 0 preserving the sign bit and the preferred exponent
1838   // of MAX(Q(x), 0)
1839   if (x_exp <= (0x1820ull << 49)) {
1840     res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1841   } else {
1842     res.w[1] = x_sign | x_exp;
1843   }
1844   res.w[0] = 0x0000000000000000ull;
1845   BID_RETURN (res);
1846 }
1847   // x is not special and is not zero
1848
1849   // if (exp <= -(p+1)) return 0.0
1850 if (x_exp <= 0x2ffa000000000000ull) {   // 0x2ffa000000000000ull == -35
1851   res.w[1] = x_sign | 0x3040000000000000ull;
1852   res.w[0] = 0x0000000000000000ull;
1853   BID_RETURN (res);
1854 }
1855   // q = nr. of decimal digits in x
1856   //  determine first the nr. of bits in x
1857 if (C1.w[1] == 0) {
1858   if (C1.w[0] >= 0x0020000000000000ull) {       // x >= 2^53
1859     // split the 64-bit value in two 32-bit halves to avoid rounding errors
1860     if (C1.w[0] >= 0x0000000100000000ull) {     // x >= 2^32
1861       tmp1.d = (double) (C1.w[0] >> 32);        // exact conversion
1862       x_nr_bits =
1863         33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1864     } else {    // x < 2^32
1865       tmp1.d = (double) (C1.w[0]);      // exact conversion
1866       x_nr_bits =
1867         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1868     }
1869   } else {      // if x < 2^53
1870     tmp1.d = (double) C1.w[0];  // exact conversion
1871     x_nr_bits =
1872       1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1873   }
1874 } else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1875   tmp1.d = (double) C1.w[1];    // exact conversion
1876   x_nr_bits =
1877     65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1878 }
1879
1880 q = nr_digits[x_nr_bits - 1].digits;
1881 if (q == 0) {
1882   q = nr_digits[x_nr_bits - 1].digits1;
1883   if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1884       (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1885        C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1886     q++;
1887 }
1888 exp = (x_exp >> 49) - 6176;
1889 if (exp >= 0) { // -exp <= 0
1890   // the argument is an integer already
1891   res.w[1] = x.w[1];
1892   res.w[0] = x.w[0];
1893   BID_RETURN (res);
1894 } else if ((q + exp) >= 0) {    // exp < 0 and 1 <= -exp <= q
1895   // need to shift right -exp digits from the coefficient; the exp will be 0
1896   ind = -exp;   // 1 <= ind <= 34; ind is a synonym for 'x'
1897   // chop off ind digits from the lower part of C1
1898   // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
1899   tmp64 = C1.w[0];
1900   if (ind <= 19) {
1901     C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1902   } else {
1903     C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1904     C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1905   }
1906   if (C1.w[0] < tmp64)
1907     C1.w[1]++;
1908   // calculate C* and f*
1909   // C* is actually floor(C*) in this case
1910   // C* and f* need shifting and masking, as shown by
1911   // shiftright128[] and maskhigh128[]
1912   // 1 <= x <= 34
1913   // kx = 10^(-x) = ten2mk128[ind - 1]
1914   // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1915   // the approximation of 10^(-x) was rounded up to 118 bits
1916   __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1917   // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1918   // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1919   // if (0 < f* < 10^(-x)) then the result is a midpoint
1920   //   if floor(C*) is even then C* = floor(C*) - logical right
1921   //       shift; C* has p decimal digits, correct by Prop. 1)
1922   //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
1923   //       shift; C* has p decimal digits, correct by Pr. 1)
1924   // else
1925   //   C* = floor(C*) (logical right shift; C has p decimal digits,
1926   //       correct by Property 1)
1927   // n = C* * 10^(e+x)
1928
1929   // shift right C* by Ex-128 = shiftright128[ind]
1930   if (ind - 1 <= 2) {   // 0 <= ind - 1 <= 2 => shift = 0
1931     res.w[1] = P256.w[3];
1932     res.w[0] = P256.w[2];
1933   } else if (ind - 1 <= 21) {   // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1934     shift = shiftright128[ind - 1];     // 3 <= shift <= 63
1935     res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1936     res.w[1] = (P256.w[3] >> shift);
1937   } else {      // 22 <= ind - 1 <= 33
1938     shift = shiftright128[ind - 1];     // 2 <= shift <= 38
1939     res.w[1] = 0;
1940     res.w[0] = (P256.w[3] >> (shift - 64));     // 2 <= shift - 64 <= 38
1941   }
1942   // if the result was a midpoint, it was already rounded away from zero
1943   res.w[1] |= x_sign | 0x3040000000000000ull;
1944   BID_RETURN (res);
1945 } else {        // if ((q + exp) < 0) <=> q < -exp
1946   // the result is +0 or -0
1947   res.w[1] = x_sign | 0x3040000000000000ull;
1948   res.w[0] = 0x0000000000000000ull;
1949   BID_RETURN (res);
1950 }
1951 }