libgcc/config/libbid/bid128_to_uint64.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 #include "bid_internal.h"
  25
  26 /*****************************************************************************
  27  *  BID128_to_uint64_rnint
  28  ****************************************************************************/
  29
  30 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
  31                                           bid128_to_uint64_rnint, x)
  32
  33      UINT64 res;
  34      UINT64 x_sign;
  35      UINT64 x_exp;
  36      int exp;                   // unbiased exponent
  37   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  38      UINT64 tmp64;
  39      BID_UI64DOUBLE tmp1;
  40      unsigned int x_nr_bits;
  41      int q, ind, shift;
  42      UINT128 C1, C;
  43      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  44      UINT256 fstar;
  45      UINT256 P256;
  46
  47   // unpack x
  48 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  49 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
  50 C1.w[1] = x.w[1] & MASK_COEFF;
  51 C1.w[0] = x.w[0];
  52
  53   // check for NaN or Infinity
  54 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  55     // x is special
  56 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
  57   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
  58     // set invalid flag
  59     *pfpsf |= INVALID_EXCEPTION;
  60     // return Integer Indefinite
  61     res = 0x8000000000000000ull;
  62   } else {      // x is QNaN
  63     // set invalid flag
  64     *pfpsf |= INVALID_EXCEPTION;
  65     // return Integer Indefinite
  66     res = 0x8000000000000000ull;
  67   }
  68   BID_RETURN (res);
  69 } else {        // x is not a NaN, so it must be infinity
  70   if (!x_sign) {        // x is +inf
  71     // set invalid flag
  72     *pfpsf |= INVALID_EXCEPTION;
  73     // return Integer Indefinite
  74     res = 0x8000000000000000ull;
  75   } else {      // x is -inf
  76     // set invalid flag
  77     *pfpsf |= INVALID_EXCEPTION;
  78     // return Integer Indefinite
  79     res = 0x8000000000000000ull;
  80   }
  81   BID_RETURN (res);
  82 }
  83 }
  84   // check for non-canonical values (after the check for special values)
  85 if ((C1.w[1] > 0x0001ed09bead87c0ull)
  86     || (C1.w[1] == 0x0001ed09bead87c0ull
  87         && (C1.w[0] > 0x378d8e63ffffffffull))
  88     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  89   res = 0x0000000000000000ull;
  90   BID_RETURN (res);
  91 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  92   // x is 0
  93   res = 0x0000000000000000ull;
  94   BID_RETURN (res);
  95 } else {        // x is not special and is not zero
  96
  97   // q = nr. of decimal digits in x
  98   //  determine first the nr. of bits in x
  99   if (C1.w[1] == 0) {
 100     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
 101       // split the 64-bit value in two 32-bit halves to avoid rounding errors
 102       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
 103         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
 104         x_nr_bits =
 105           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 106       } else {  // x < 2^32
 107         tmp1.d = (double) (C1.w[0]);    // exact conversion
 108         x_nr_bits =
 109           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 110       }
 111     } else {    // if x < 2^53
 112       tmp1.d = (double) C1.w[0];        // exact conversion
 113       x_nr_bits =
 114         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 115     }
 116   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
 117     tmp1.d = (double) C1.w[1];  // exact conversion
 118     x_nr_bits =
 119       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 120   }
 121   q = nr_digits[x_nr_bits - 1].digits;
 122   if (q == 0) {
 123     q = nr_digits[x_nr_bits - 1].digits1;
 124     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
 125         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
 126             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
 127       q++;
 128   }
 129   exp = (x_exp >> 49) - 6176;
 130
 131   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
 132     // set invalid flag
 133     *pfpsf |= INVALID_EXCEPTION;
 134     // return Integer Indefinite
 135     res = 0x8000000000000000ull;
 136     BID_RETURN (res);
 137   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
 138     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
 139     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
 140     // the cases that do not fit are identified here; the ones that fit
 141     // fall through and will be handled with other cases further,
 142     // under '1 <= q + exp <= 20'
 143     if (x_sign) {       // if n < 0 and q + exp = 20
 144       // if n < -1/2 then n cannot be converted to uint64 with RN
 145       // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 1/2
 146       // <=> 0.c(0)c(1)...c(q-1) * 10^21 > 0x05, 1<=q<=34
 147       // <=> C * 10^(21-q) > 0x05, 1<=q<=34
 148       if (q == 21) {
 149         // C > 5
 150         if (C1.w[1] != 0 || C1.w[0] > 0x05ull) {
 151           // set invalid flag
 152           *pfpsf |= INVALID_EXCEPTION;
 153           // return Integer Indefinite
 154           res = 0x8000000000000000ull;
 155           BID_RETURN (res);
 156         }
 157         // else cases that can be rounded to 64-bit unsigned int fall through
 158         // to '1 <= q + exp <= 20'
 159       } else {
 160         // if 1 <= q <= 20
 161         //   C * 10^(21-q) > 5 is true because C >= 1 and 10^(21-q) >= 10
 162         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
 163         //   C > 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
 164         // set invalid flag
 165         *pfpsf |= INVALID_EXCEPTION;
 166         // return Integer Indefinite
 167         res = 0x8000000000000000ull;
 168         BID_RETURN (res);
 169       }
 170     } else {    // if n > 0 and q + exp = 20
 171       // if n >= 2^64 - 1/2 then n is too large
 172       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
 173       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
 174       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
 175       // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
 176       if (q == 1) {
 177         // C * 10^20 >= 0x9fffffffffffffffb
 178         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
 179         if (C.w[1] > 0x09 || (C.w[1] == 0x09
 180                               && C.w[0] >= 0xfffffffffffffffbull)) {
 181           // set invalid flag
 182           *pfpsf |= INVALID_EXCEPTION;
 183           // return Integer Indefinite
 184           res = 0x8000000000000000ull;
 185           BID_RETURN (res);
 186         }
 187         // else cases that can be rounded to a 64-bit int fall through
 188         // to '1 <= q + exp <= 20'
 189       } else if (q <= 19) {
 190         // C * 10^(21-q) >= 0x9fffffffffffffffb
 191         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
 192         if (C.w[1] > 0x09 || (C.w[1] == 0x09
 193                               && C.w[0] >= 0xfffffffffffffffbull)) {
 194           // set invalid flag
 195           *pfpsf |= INVALID_EXCEPTION;
 196           // return Integer Indefinite
 197           res = 0x8000000000000000ull;
 198           BID_RETURN (res);
 199         }
 200         // else cases that can be rounded to a 64-bit int fall through
 201         // to '1 <= q + exp <= 20'
 202       } else if (q == 20) {
 203         // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
 204         C.w[0] = C1.w[0] + C1.w[0];
 205         C.w[1] = C1.w[1] + C1.w[1];
 206         if (C.w[0] < C1.w[0])
 207           C.w[1]++;
 208         if (C.w[1] > 0x01 || (C.w[1] == 0x01
 209                               && C.w[0] >= 0xffffffffffffffffull)) {
 210           // set invalid flag
 211           *pfpsf |= INVALID_EXCEPTION;
 212           // return Integer Indefinite
 213           res = 0x8000000000000000ull;
 214           BID_RETURN (res);
 215         }
 216         // else cases that can be rounded to a 64-bit int fall through
 217         // to '1 <= q + exp <= 20'
 218       } else if (q == 21) {
 219         // C >= 0x9fffffffffffffffb
 220         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
 221                                && C1.w[0] >= 0xfffffffffffffffbull)) {
 222           // set invalid flag
 223           *pfpsf |= INVALID_EXCEPTION;
 224           // return Integer Indefinite
 225           res = 0x8000000000000000ull;
 226           BID_RETURN (res);
 227         }
 228         // else cases that can be rounded to a 64-bit int fall through
 229         // to '1 <= q + exp <= 20'
 230       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
 231         // C  >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
 232         C.w[1] = 0x09;
 233         C.w[0] = 0xfffffffffffffffbull;
 234         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
 235         if (C1.w[1] > C.w[1]
 236             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
 237           // set invalid flag
 238           *pfpsf |= INVALID_EXCEPTION;
 239           // return Integer Indefinite
 240           res = 0x8000000000000000ull;
 241           BID_RETURN (res);
 242         }
 243         // else cases that can be rounded to a 64-bit int fall through
 244         // to '1 <= q + exp <= 20'
 245       }
 246     }
 247   }
 248   // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
 249   // Note: some of the cases tested for above fall through to this point
 250   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
 251     // return 0
 252     res = 0x0000000000000000ull;
 253     BID_RETURN (res);
 254   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
 255     // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
 256     //   res = 0
 257     // else if x > 0
 258     //   res = +1
 259     // else // if x < 0
 260     //   invalid exc
 261     ind = q - 1;
 262     if (ind <= 18) {    // 0 <= ind <= 18
 263       if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
 264         res = 0x0000000000000000ull;    // return 0
 265       } else if (!x_sign) {     // n > 0
 266         res = 0x00000001;       // return +1
 267       } else {
 268         res = 0x8000000000000000ull;
 269         *pfpsf |= INVALID_EXCEPTION;
 270       }
 271     } else {    // 19 <= ind <= 33
 272       if ((C1.w[1] < midpoint128[ind - 19].w[1])
 273           || ((C1.w[1] == midpoint128[ind - 19].w[1])
 274               && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
 275         res = 0x0000000000000000ull;    // return 0
 276       } else if (!x_sign) {     // n > 0
 277         res = 0x00000001;       // return +1
 278       } else {
 279         res = 0x8000000000000000ull;
 280         *pfpsf |= INVALID_EXCEPTION;
 281       }
 282     }
 283   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
 284     // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
 285     // to nearest to a 64-bit unsigned signed integer
 286     if (x_sign) {       // x <= -1
 287       // set invalid flag
 288       *pfpsf |= INVALID_EXCEPTION;
 289       // return Integer Indefinite
 290       res = 0x8000000000000000ull;
 291       BID_RETURN (res);
 292     }
 293     // 1 <= x < 2^64-1/2 so x can be rounded
 294     // to nearest to a 64-bit unsigned integer
 295     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
 296       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
 297       // chop off ind digits from the lower part of C1
 298       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
 299       tmp64 = C1.w[0];
 300       if (ind <= 19) {
 301         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
 302       } else {
 303         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
 304         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
 305       }
 306       if (C1.w[0] < tmp64)
 307         C1.w[1]++;
 308       // calculate C* and f*
 309       // C* is actually floor(C*) in this case
 310       // C* and f* need shifting and masking, as shown by
 311       // shiftright128[] and maskhigh128[]
 312       // 1 <= x <= 33
 313       // kx = 10^(-x) = ten2mk128[ind - 1]
 314       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 315       // the approximation of 10^(-x) was rounded up to 118 bits
 316       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
 317       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
 318         Cstar.w[1] = P256.w[3];
 319         Cstar.w[0] = P256.w[2];
 320         fstar.w[3] = 0;
 321         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
 322         fstar.w[1] = P256.w[1];
 323         fstar.w[0] = P256.w[0];
 324       } else {  // 22 <= ind - 1 <= 33
 325         Cstar.w[1] = 0;
 326         Cstar.w[0] = P256.w[3];
 327         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
 328         fstar.w[2] = P256.w[2];
 329         fstar.w[1] = P256.w[1];
 330         fstar.w[0] = P256.w[0];
 331       }
 332       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
 333       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
 334       // if (0 < f* < 10^(-x)) then the result is a midpoint
 335       //   if floor(C*) is even then C* = floor(C*) - logical right
 336       //       shift; C* has p decimal digits, correct by Prop. 1)
 337       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
 338       //       shift; C* has p decimal digits, correct by Pr. 1)
 339       // else
 340       //   C* = floor(C*) (logical right shift; C has p decimal digits,
 341       //       correct by Property 1)
 342       // n = C* * 10^(e+x)
 343
 344       // shift right C* by Ex-128 = shiftright128[ind]
 345       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
 346       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
 347         Cstar.w[0] =
 348           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
 349         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
 350       } else {  // 22 <= ind - 1 <= 33
 351         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
 352       }
 353       // if the result was a midpoint it was rounded away from zero, so
 354       // it will need a correction
 355       // check for midpoints
 356       if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
 357           && (fstar.w[1] || fstar.w[0])
 358           && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
 359               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
 360                   && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
 361         // the result is a midpoint; round to nearest
 362         if (Cstar.w[0] & 0x01) {        // Cstar.w[0] is odd; MP in [EVEN, ODD]
 363           // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
 364           Cstar.w[0]--; // Cstar.w[0] is now even
 365         }       // else MP in [ODD, EVEN]
 366       }
 367       res = Cstar.w[0]; // the result is positive
 368     } else if (exp == 0) {
 369       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
 370       // res = C (exact)
 371       res = C1.w[0];
 372     } else {
 373       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
 374       // res = C * 10^exp (exact) - must fit in 64 bits
 375       res = C1.w[0] * ten2k64[exp];
 376     }
 377   }
 378 }
 379
 380 BID_RETURN (res);
 381 }
 382
 383 /*****************************************************************************
 384  *  BID128_to_uint64_xrnint
 385  ****************************************************************************/
 386
 387 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
 388                                           bid128_to_uint64_xrnint, x)
 389
 390      UINT64 res;
 391      UINT64 x_sign;
 392      UINT64 x_exp;
 393      int exp;                   // unbiased exponent
 394   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
 395      UINT64 tmp64, tmp64A;
 396      BID_UI64DOUBLE tmp1;
 397      unsigned int x_nr_bits;
 398      int q, ind, shift;
 399      UINT128 C1, C;
 400      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
 401      UINT256 fstar;
 402      UINT256 P256;
 403
 404   // unpack x
 405 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
 406 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
 407 C1.w[1] = x.w[1] & MASK_COEFF;
 408 C1.w[0] = x.w[0];
 409
 410   // check for NaN or Infinity
 411 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
 412     // x is special
 413 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
 414   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
 415     // set invalid flag
 416     *pfpsf |= INVALID_EXCEPTION;
 417     // return Integer Indefinite
 418     res = 0x8000000000000000ull;
 419   } else {      // x is QNaN
 420     // set invalid flag
 421     *pfpsf |= INVALID_EXCEPTION;
 422     // return Integer Indefinite
 423     res = 0x8000000000000000ull;
 424   }
 425   BID_RETURN (res);
 426 } else {        // x is not a NaN, so it must be infinity
 427   if (!x_sign) {        // x is +inf
 428     // set invalid flag
 429     *pfpsf |= INVALID_EXCEPTION;
 430     // return Integer Indefinite
 431     res = 0x8000000000000000ull;
 432   } else {      // x is -inf
 433     // set invalid flag
 434     *pfpsf |= INVALID_EXCEPTION;
 435     // return Integer Indefinite
 436     res = 0x8000000000000000ull;
 437   }
 438   BID_RETURN (res);
 439 }
 440 }
 441   // check for non-canonical values (after the check for special values)
 442 if ((C1.w[1] > 0x0001ed09bead87c0ull)
 443     || (C1.w[1] == 0x0001ed09bead87c0ull
 444         && (C1.w[0] > 0x378d8e63ffffffffull))
 445     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
 446   res = 0x0000000000000000ull;
 447   BID_RETURN (res);
 448 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
 449   // x is 0
 450   res = 0x0000000000000000ull;
 451   BID_RETURN (res);
 452 } else {        // x is not special and is not zero
 453
 454   // q = nr. of decimal digits in x
 455   //  determine first the nr. of bits in x
 456   if (C1.w[1] == 0) {
 457     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
 458       // split the 64-bit value in two 32-bit halves to avoid rounding errors
 459       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
 460         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
 461         x_nr_bits =
 462           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 463       } else {  // x < 2^32
 464         tmp1.d = (double) (C1.w[0]);    // exact conversion
 465         x_nr_bits =
 466           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 467       }
 468     } else {    // if x < 2^53
 469       tmp1.d = (double) C1.w[0];        // exact conversion
 470       x_nr_bits =
 471         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 472     }
 473   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
 474     tmp1.d = (double) C1.w[1];  // exact conversion
 475     x_nr_bits =
 476       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 477   }
 478   q = nr_digits[x_nr_bits - 1].digits;
 479   if (q == 0) {
 480     q = nr_digits[x_nr_bits - 1].digits1;
 481     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
 482         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
 483             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
 484       q++;
 485   }
 486   exp = (x_exp >> 49) - 6176;
 487
 488   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
 489     // set invalid flag
 490     *pfpsf |= INVALID_EXCEPTION;
 491     // return Integer Indefinite
 492     res = 0x8000000000000000ull;
 493     BID_RETURN (res);
 494   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
 495     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
 496     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
 497     // the cases that do not fit are identified here; the ones that fit
 498     // fall through and will be handled with other cases further,
 499     // under '1 <= q + exp <= 20'
 500     if (x_sign) {       // if n < 0 and q + exp = 20
 501       // if n < -1/2 then n cannot be converted to uint64 with RN
 502       // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 1/2
 503       // <=> 0.c(0)c(1)...c(q-1) * 10^21 > 0x05, 1<=q<=34
 504       // <=> C * 10^(21-q) > 0x05, 1<=q<=34
 505       if (q == 21) {
 506         // C > 5
 507         if (C1.w[1] != 0 || C1.w[0] > 0x05ull) {
 508           // set invalid flag
 509           *pfpsf |= INVALID_EXCEPTION;
 510           // return Integer Indefinite
 511           res = 0x8000000000000000ull;
 512           BID_RETURN (res);
 513         }
 514         // else cases that can be rounded to 64-bit unsigned int fall through
 515         // to '1 <= q + exp <= 20'
 516       } else {
 517         // if 1 <= q <= 20
 518         //   C * 10^(21-q) > 5 is true because C >= 1 and 10^(21-q) >= 10
 519         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
 520         //   C > 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
 521         // set invalid flag
 522         *pfpsf |= INVALID_EXCEPTION;
 523         // return Integer Indefinite
 524         res = 0x8000000000000000ull;
 525         BID_RETURN (res);
 526       }
 527     } else {    // if n > 0 and q + exp = 20
 528       // if n >= 2^64 - 1/2 then n is too large
 529       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
 530       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
 531       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
 532       // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
 533       if (q == 1) {
 534         // C * 10^20 >= 0x9fffffffffffffffb
 535         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
 536         if (C.w[1] > 0x09 || (C.w[1] == 0x09
 537                               && C.w[0] >= 0xfffffffffffffffbull)) {
 538           // set invalid flag
 539           *pfpsf |= INVALID_EXCEPTION;
 540           // return Integer Indefinite
 541           res = 0x8000000000000000ull;
 542           BID_RETURN (res);
 543         }
 544         // else cases that can be rounded to a 64-bit int fall through
 545         // to '1 <= q + exp <= 20'
 546       } else if (q <= 19) {
 547         // C * 10^(21-q) >= 0x9fffffffffffffffb
 548         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
 549         if (C.w[1] > 0x09 || (C.w[1] == 0x09
 550                               && C.w[0] >= 0xfffffffffffffffbull)) {
 551           // set invalid flag
 552           *pfpsf |= INVALID_EXCEPTION;
 553           // return Integer Indefinite
 554           res = 0x8000000000000000ull;
 555           BID_RETURN (res);
 556         }
 557         // else cases that can be rounded to a 64-bit int fall through
 558         // to '1 <= q + exp <= 20'
 559       } else if (q == 20) {
 560         // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
 561         C.w[0] = C1.w[0] + C1.w[0];
 562         C.w[1] = C1.w[1] + C1.w[1];
 563         if (C.w[0] < C1.w[0])
 564           C.w[1]++;
 565         if (C.w[1] > 0x01 || (C.w[1] == 0x01
 566                               && C.w[0] >= 0xffffffffffffffffull)) {
 567           // set invalid flag
 568           *pfpsf |= INVALID_EXCEPTION;
 569           // return Integer Indefinite
 570           res = 0x8000000000000000ull;
 571           BID_RETURN (res);
 572         }
 573         // else cases that can be rounded to a 64-bit int fall through
 574         // to '1 <= q + exp <= 20'
 575       } else if (q == 21) {
 576         // C >= 0x9fffffffffffffffb
 577         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
 578                                && C1.w[0] >= 0xfffffffffffffffbull)) {
 579           // set invalid flag
 580           *pfpsf |= INVALID_EXCEPTION;
 581           // return Integer Indefinite
 582           res = 0x8000000000000000ull;
 583           BID_RETURN (res);
 584         }
 585         // else cases that can be rounded to a 64-bit int fall through
 586         // to '1 <= q + exp <= 20'
 587       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
 588         // C  >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
 589         C.w[1] = 0x09;
 590         C.w[0] = 0xfffffffffffffffbull;
 591         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
 592         if (C1.w[1] > C.w[1]
 593             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
 594           // set invalid flag
 595           *pfpsf |= INVALID_EXCEPTION;
 596           // return Integer Indefinite
 597           res = 0x8000000000000000ull;
 598           BID_RETURN (res);
 599         }
 600         // else cases that can be rounded to a 64-bit int fall through
 601         // to '1 <= q + exp <= 20'
 602       }
 603     }
 604   }
 605   // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
 606   // Note: some of the cases tested for above fall through to this point
 607   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
 608     // set inexact flag
 609     *pfpsf |= INEXACT_EXCEPTION;
 610     // return 0
 611     res = 0x0000000000000000ull;
 612     BID_RETURN (res);
 613   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
 614     // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
 615     //   res = 0
 616     // else if x > 0
 617     //   res = +1
 618     // else // if x < 0
 619     //   invalid exc
 620     ind = q - 1;
 621     if (ind <= 18) {    // 0 <= ind <= 18
 622       if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
 623         res = 0x0000000000000000ull;    // return 0
 624       } else if (!x_sign) {     // n > 0
 625         res = 0x00000001;       // return +1
 626       } else {
 627         res = 0x8000000000000000ull;
 628         *pfpsf |= INVALID_EXCEPTION;
 629         BID_RETURN (res);
 630       }
 631     } else {    // 19 <= ind <= 33
 632       if ((C1.w[1] < midpoint128[ind - 19].w[1])
 633           || ((C1.w[1] == midpoint128[ind - 19].w[1])
 634               && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
 635         res = 0x0000000000000000ull;    // return 0
 636       } else if (!x_sign) {     // n > 0
 637         res = 0x00000001;       // return +1
 638       } else {
 639         res = 0x8000000000000000ull;
 640         *pfpsf |= INVALID_EXCEPTION;
 641         BID_RETURN (res);
 642       }
 643     }
 644     // set inexact flag
 645     *pfpsf |= INEXACT_EXCEPTION;
 646   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
 647     // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
 648     // to nearest to a 64-bit unsigned signed integer
 649     if (x_sign) {       // x <= -1
 650       // set invalid flag
 651       *pfpsf |= INVALID_EXCEPTION;
 652       // return Integer Indefinite
 653       res = 0x8000000000000000ull;
 654       BID_RETURN (res);
 655     }
 656     // 1 <= x < 2^64-1/2 so x can be rounded
 657     // to nearest to a 64-bit unsigned integer
 658     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
 659       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
 660       // chop off ind digits from the lower part of C1
 661       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
 662       tmp64 = C1.w[0];
 663       if (ind <= 19) {
 664         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
 665       } else {
 666         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
 667         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
 668       }
 669       if (C1.w[0] < tmp64)
 670         C1.w[1]++;
 671       // calculate C* and f*
 672       // C* is actually floor(C*) in this case
 673       // C* and f* need shifting and masking, as shown by
 674       // shiftright128[] and maskhigh128[]
 675       // 1 <= x <= 33
 676       // kx = 10^(-x) = ten2mk128[ind - 1]
 677       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 678       // the approximation of 10^(-x) was rounded up to 118 bits
 679       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
 680       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
 681         Cstar.w[1] = P256.w[3];
 682         Cstar.w[0] = P256.w[2];
 683         fstar.w[3] = 0;
 684         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
 685         fstar.w[1] = P256.w[1];
 686         fstar.w[0] = P256.w[0];
 687       } else {  // 22 <= ind - 1 <= 33
 688         Cstar.w[1] = 0;
 689         Cstar.w[0] = P256.w[3];
 690         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
 691         fstar.w[2] = P256.w[2];
 692         fstar.w[1] = P256.w[1];
 693         fstar.w[0] = P256.w[0];
 694       }
 695       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
 696       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
 697       // if (0 < f* < 10^(-x)) then the result is a midpoint
 698       //   if floor(C*) is even then C* = floor(C*) - logical right
 699       //       shift; C* has p decimal digits, correct by Prop. 1)
 700       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
 701       //       shift; C* has p decimal digits, correct by Pr. 1)
 702       // else
 703       //   C* = floor(C*) (logical right shift; C has p decimal digits,
 704       //       correct by Property 1)
 705       // n = C* * 10^(e+x)
 706
 707       // shift right C* by Ex-128 = shiftright128[ind]
 708       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
 709       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
 710         Cstar.w[0] =
 711           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
 712         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
 713       } else {  // 22 <= ind - 1 <= 33
 714         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
 715       }
 716       // determine inexactness of the rounding of C*
 717       // if (0 < f* - 1/2 < 10^(-x)) then
 718       //   the result is exact
 719       // else // if (f* - 1/2 > T*) then
 720       //   the result is inexact
 721       if (ind - 1 <= 2) {
 722         if (fstar.w[1] > 0x8000000000000000ull ||
 723             (fstar.w[1] == 0x8000000000000000ull
 724              && fstar.w[0] > 0x0ull)) {
 725           // f* > 1/2 and the result may be exact
 726           tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
 727           if (tmp64 > ten2mk128trunc[ind - 1].w[1]
 728               || (tmp64 == ten2mk128trunc[ind - 1].w[1]
 729                   && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
 730             // set the inexact flag
 731             *pfpsf |= INEXACT_EXCEPTION;
 732           }     // else the result is exact
 733         } else {        // the result is inexact; f2* <= 1/2
 734           // set the inexact flag
 735           *pfpsf |= INEXACT_EXCEPTION;
 736         }
 737       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
 738         if (fstar.w[3] > 0x0 ||
 739             (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
 740             (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
 741              (fstar.w[1] || fstar.w[0]))) {
 742           // f2* > 1/2 and the result may be exact
 743           // Calculate f2* - 1/2
 744           tmp64 = fstar.w[2] - onehalf128[ind - 1];
 745           tmp64A = fstar.w[3];
 746           if (tmp64 > fstar.w[2])
 747             tmp64A--;
 748           if (tmp64A || tmp64
 749               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
 750               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
 751                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
 752             // set the inexact flag
 753             *pfpsf |= INEXACT_EXCEPTION;
 754           }     // else the result is exact
 755         } else {        // the result is inexact; f2* <= 1/2
 756           // set the inexact flag
 757           *pfpsf |= INEXACT_EXCEPTION;
 758         }
 759       } else {  // if 22 <= ind <= 33
 760         if (fstar.w[3] > onehalf128[ind - 1] ||
 761             (fstar.w[3] == onehalf128[ind - 1] &&
 762              (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
 763           // f2* > 1/2 and the result may be exact
 764           // Calculate f2* - 1/2
 765           tmp64 = fstar.w[3] - onehalf128[ind - 1];
 766           if (tmp64 || fstar.w[2]
 767               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
 768               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
 769                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
 770             // set the inexact flag
 771             *pfpsf |= INEXACT_EXCEPTION;
 772           }     // else the result is exact
 773         } else {        // the result is inexact; f2* <= 1/2
 774           // set the inexact flag
 775           *pfpsf |= INEXACT_EXCEPTION;
 776         }
 777       }
 778
 779       // if the result was a midpoint it was rounded away from zero, so
 780       // it will need a correction
 781       // check for midpoints
 782       if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
 783           && (fstar.w[1] || fstar.w[0])
 784           && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
 785               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
 786                   && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
 787         // the result is a midpoint; round to nearest
 788         if (Cstar.w[0] & 0x01) {        // Cstar.w[0] is odd; MP in [EVEN, ODD]
 789           // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
 790           Cstar.w[0]--; // Cstar.w[0] is now even
 791         }       // else MP in [ODD, EVEN]
 792       }
 793       res = Cstar.w[0]; // the result is positive
 794     } else if (exp == 0) {
 795       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
 796       // res = C (exact)
 797       res = C1.w[0];
 798     } else {
 799       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
 800       // res = C * 10^exp (exact) - must fit in 64 bits
 801       res = C1.w[0] * ten2k64[exp];
 802     }
 803   }
 804 }
 805
 806 BID_RETURN (res);
 807 }
 808
 809 /*****************************************************************************
 810  *  BID128_to_uint64_floor
 811  ****************************************************************************/
 812
 813 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
 814                                           bid128_to_uint64_floor, x)
 815
 816      UINT64 res;
 817      UINT64 x_sign;
 818      UINT64 x_exp;
 819      int exp;                   // unbiased exponent
 820   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
 821      BID_UI64DOUBLE tmp1;
 822      unsigned int x_nr_bits;
 823      int q, ind, shift;
 824      UINT128 C1, C;
 825      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
 826      UINT256 P256;
 827
 828   // unpack x
 829 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
 830 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
 831 C1.w[1] = x.w[1] & MASK_COEFF;
 832 C1.w[0] = x.w[0];
 833
 834   // check for NaN or Infinity
 835 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
 836     // x is special
 837 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
 838   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
 839     // set invalid flag
 840     *pfpsf |= INVALID_EXCEPTION;
 841     // return Integer Indefinite
 842     res = 0x8000000000000000ull;
 843   } else {      // x is QNaN
 844     // set invalid flag
 845     *pfpsf |= INVALID_EXCEPTION;
 846     // return Integer Indefinite
 847     res = 0x8000000000000000ull;
 848   }
 849   BID_RETURN (res);
 850 } else {        // x is not a NaN, so it must be infinity
 851   if (!x_sign) {        // x is +inf
 852     // set invalid flag
 853     *pfpsf |= INVALID_EXCEPTION;
 854     // return Integer Indefinite
 855     res = 0x8000000000000000ull;
 856   } else {      // x is -inf
 857     // set invalid flag
 858     *pfpsf |= INVALID_EXCEPTION;
 859     // return Integer Indefinite
 860     res = 0x8000000000000000ull;
 861   }
 862   BID_RETURN (res);
 863 }
 864 }
 865   // check for non-canonical values (after the check for special values)
 866 if ((C1.w[1] > 0x0001ed09bead87c0ull)
 867     || (C1.w[1] == 0x0001ed09bead87c0ull
 868         && (C1.w[0] > 0x378d8e63ffffffffull))
 869     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
 870   res = 0x0000000000000000ull;
 871   BID_RETURN (res);
 872 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
 873   // x is 0
 874   res = 0x0000000000000000ull;
 875   BID_RETURN (res);
 876 } else {        // x is not special and is not zero
 877
 878   // if n < 0 then n cannot be converted to uint64 with RM
 879   if (x_sign) { // if n < 0 and q + exp = 20
 880     // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 0
 881     // set invalid flag
 882     *pfpsf |= INVALID_EXCEPTION;
 883     // return Integer Indefinite
 884     res = 0x8000000000000000ull;
 885     BID_RETURN (res);
 886   }
 887   // q = nr. of decimal digits in x
 888   //  determine first the nr. of bits in x
 889   if (C1.w[1] == 0) {
 890     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
 891       // split the 64-bit value in two 32-bit halves to avoid rounding errors
 892       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
 893         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
 894         x_nr_bits =
 895           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 896       } else {  // x < 2^32
 897         tmp1.d = (double) (C1.w[0]);    // exact conversion
 898         x_nr_bits =
 899           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 900       }
 901     } else {    // if x < 2^53
 902       tmp1.d = (double) C1.w[0];        // exact conversion
 903       x_nr_bits =
 904         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 905     }
 906   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
 907     tmp1.d = (double) C1.w[1];  // exact conversion
 908     x_nr_bits =
 909       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 910   }
 911   q = nr_digits[x_nr_bits - 1].digits;
 912   if (q == 0) {
 913     q = nr_digits[x_nr_bits - 1].digits1;
 914     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
 915         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
 916             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
 917       q++;
 918   }
 919   exp = (x_exp >> 49) - 6176;
 920   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
 921     // set invalid flag
 922     *pfpsf |= INVALID_EXCEPTION;
 923     // return Integer Indefinite
 924     res = 0x8000000000000000ull;
 925     BID_RETURN (res);
 926   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
 927     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
 928     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
 929     // the cases that do not fit are identified here; the ones that fit
 930     // fall through and will be handled with other cases further,
 931     // under '1 <= q + exp <= 20'
 932     // if n > 0 and q + exp = 20
 933     // if n >= 2^64 then n is too large
 934     // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
 935     // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
 936     // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
 937     // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
 938     if (q == 1) {
 939       // C * 10^20 >= 0xa0000000000000000
 940       __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);    // 10^20 * C
 941       if (C.w[1] >= 0x0a) {
 942         // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
 943         // set invalid flag
 944         *pfpsf |= INVALID_EXCEPTION;
 945         // return Integer Indefinite
 946         res = 0x8000000000000000ull;
 947         BID_RETURN (res);
 948       }
 949       // else cases that can be rounded to a 64-bit int fall through
 950       // to '1 <= q + exp <= 20'
 951     } else if (q <= 19) {
 952       // C * 10^(21-q) >= 0xa0000000000000000
 953       __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
 954       if (C.w[1] >= 0x0a) {
 955         // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
 956         // set invalid flag
 957         *pfpsf |= INVALID_EXCEPTION;
 958         // return Integer Indefinite
 959         res = 0x8000000000000000ull;
 960         BID_RETURN (res);
 961       }
 962       // else cases that can be rounded to a 64-bit int fall through
 963       // to '1 <= q + exp <= 20'
 964     } else if (q == 20) {
 965       // C >= 0x10000000000000000
 966       if (C1.w[1] >= 0x01) {
 967         // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
 968         // set invalid flag
 969         *pfpsf |= INVALID_EXCEPTION;
 970         // return Integer Indefinite
 971         res = 0x8000000000000000ull;
 972         BID_RETURN (res);
 973       }
 974       // else cases that can be rounded to a 64-bit int fall through
 975       // to '1 <= q + exp <= 20'
 976     } else if (q == 21) {
 977       // C >= 0xa0000000000000000
 978       if (C1.w[1] >= 0x0a) {
 979         // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
 980         // set invalid flag
 981         *pfpsf |= INVALID_EXCEPTION;
 982         // return Integer Indefinite
 983         res = 0x8000000000000000ull;
 984         BID_RETURN (res);
 985       }
 986       // else cases that can be rounded to a 64-bit int fall through
 987       // to '1 <= q + exp <= 20'
 988     } else {    // if 22 <= q <= 34 => 1 <= q - 21 <= 13
 989       // C  >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
 990       C.w[1] = 0x0a;
 991       C.w[0] = 0x0000000000000000ull;
 992       __mul_128x64_to_128 (C, ten2k64[q - 21], C);
 993       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
 994         // set invalid flag
 995         *pfpsf |= INVALID_EXCEPTION;
 996         // return Integer Indefinite
 997         res = 0x8000000000000000ull;
 998         BID_RETURN (res);
 999       }
1000       // else cases that can be rounded to a 64-bit int fall through
1001       // to '1 <= q + exp <= 20'
1002     }
1003   }
1004   // n is not too large to be converted to int64 if 0 <= n < 2^64
1005   // Note: some of the cases tested for above fall through to this point
1006   if ((q + exp) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
1007     // return 0
1008     res = 0x0000000000000000ull;
1009     BID_RETURN (res);
1010   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1011     // 1 <= x < 2^64 so x can be rounded
1012     // down to a 64-bit unsigned signed integer
1013     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1014       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1015       // chop off ind digits from the lower part of C1
1016       // C1 fits in 127 bits
1017       // calculate C* and f*
1018       // C* is actually floor(C*) in this case
1019       // C* and f* need shifting and masking, as shown by
1020       // shiftright128[] and maskhigh128[]
1021       // 1 <= x <= 33
1022       // kx = 10^(-x) = ten2mk128[ind - 1]
1023       // C* = C1 * 10^(-x)
1024       // the approximation of 10^(-x) was rounded up to 118 bits
1025       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1026       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1027         Cstar.w[1] = P256.w[3];
1028         Cstar.w[0] = P256.w[2];
1029       } else {  // 22 <= ind - 1 <= 33
1030         Cstar.w[1] = 0;
1031         Cstar.w[0] = P256.w[3];
1032       }
1033       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1034       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1035       // C* = floor(C*) (logical right shift; C has p decimal digits,
1036       //     correct by Property 1)
1037       // n = C* * 10^(e+x)
1038
1039       // shift right C* by Ex-128 = shiftright128[ind]
1040       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
1041       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1042         Cstar.w[0] =
1043           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1044         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1045       } else {  // 22 <= ind - 1 <= 33
1046         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
1047       }
1048       res = Cstar.w[0]; // the result is positive
1049     } else if (exp == 0) {
1050       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
1051       // res = C (exact)
1052       res = C1.w[0];
1053     } else {
1054       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
1055       // res = C * 10^exp (exact) - must fit in 64 bits
1056       res = C1.w[0] * ten2k64[exp];
1057     }
1058   }
1059 }
1060
1061 BID_RETURN (res);
1062 }
1063
1064 /*****************************************************************************
1065  *  BID128_to_uint64_xfloor
1066  ****************************************************************************/
1067
1068 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
1069                                           bid128_to_uint64_xfloor, x)
1070
1071      UINT64 res;
1072      UINT64 x_sign;
1073      UINT64 x_exp;
1074      int exp;                   // unbiased exponent
1075   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1076      BID_UI64DOUBLE tmp1;
1077      unsigned int x_nr_bits;
1078      int q, ind, shift;
1079      UINT128 C1, C;
1080      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
1081      UINT256 fstar;
1082      UINT256 P256;
1083
1084   // unpack x
1085 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1086 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
1087 C1.w[1] = x.w[1] & MASK_COEFF;
1088 C1.w[0] = x.w[0];
1089
1090   // check for NaN or Infinity
1091 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1092     // x is special
1093 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1094   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1095     // set invalid flag
1096     *pfpsf |= INVALID_EXCEPTION;
1097     // return Integer Indefinite
1098     res = 0x8000000000000000ull;
1099   } else {      // x is QNaN
1100     // set invalid flag
1101     *pfpsf |= INVALID_EXCEPTION;
1102     // return Integer Indefinite
1103     res = 0x8000000000000000ull;
1104   }
1105   BID_RETURN (res);
1106 } else {        // x is not a NaN, so it must be infinity
1107   if (!x_sign) {        // x is +inf
1108     // set invalid flag
1109     *pfpsf |= INVALID_EXCEPTION;
1110     // return Integer Indefinite
1111     res = 0x8000000000000000ull;
1112   } else {      // x is -inf
1113     // set invalid flag
1114     *pfpsf |= INVALID_EXCEPTION;
1115     // return Integer Indefinite
1116     res = 0x8000000000000000ull;
1117   }
1118   BID_RETURN (res);
1119 }
1120 }
1121   // check for non-canonical values (after the check for special values)
1122 if ((C1.w[1] > 0x0001ed09bead87c0ull)
1123     || (C1.w[1] == 0x0001ed09bead87c0ull
1124         && (C1.w[0] > 0x378d8e63ffffffffull))
1125     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1126   res = 0x0000000000000000ull;
1127   BID_RETURN (res);
1128 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1129   // x is 0
1130   res = 0x0000000000000000ull;
1131   BID_RETURN (res);
1132 } else {        // x is not special and is not zero
1133
1134   // if n < 0 then n cannot be converted to uint64 with RM
1135   if (x_sign) { // if n < 0 and q + exp = 20
1136     // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 0
1137     // set invalid flag
1138     *pfpsf |= INVALID_EXCEPTION;
1139     // return Integer Indefinite
1140     res = 0x8000000000000000ull;
1141     BID_RETURN (res);
1142   }
1143   // q = nr. of decimal digits in x
1144   //  determine first the nr. of bits in x
1145   if (C1.w[1] == 0) {
1146     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
1147       // split the 64-bit value in two 32-bit halves to avoid rounding errors
1148       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
1149         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
1150         x_nr_bits =
1151           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1152       } else {  // x < 2^32
1153         tmp1.d = (double) (C1.w[0]);    // exact conversion
1154         x_nr_bits =
1155           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1156       }
1157     } else {    // if x < 2^53
1158       tmp1.d = (double) C1.w[0];        // exact conversion
1159       x_nr_bits =
1160         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1161     }
1162   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1163     tmp1.d = (double) C1.w[1];  // exact conversion
1164     x_nr_bits =
1165       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1166   }
1167   q = nr_digits[x_nr_bits - 1].digits;
1168   if (q == 0) {
1169     q = nr_digits[x_nr_bits - 1].digits1;
1170     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1171         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1172             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1173       q++;
1174   }
1175   exp = (x_exp >> 49) - 6176;
1176   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1177     // set invalid flag
1178     *pfpsf |= INVALID_EXCEPTION;
1179     // return Integer Indefinite
1180     res = 0x8000000000000000ull;
1181     BID_RETURN (res);
1182   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1183     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1184     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1185     // the cases that do not fit are identified here; the ones that fit
1186     // fall through and will be handled with other cases further,
1187     // under '1 <= q + exp <= 20'
1188     // if n > 0 and q + exp = 20
1189     // if n >= 2^64 then n is too large
1190     // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
1191     // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
1192     // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
1193     // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
1194     if (q == 1) {
1195       // C * 10^20 >= 0xa0000000000000000
1196       __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);    // 10^20 * C
1197       if (C.w[1] >= 0x0a) {
1198         // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
1199         // set invalid flag
1200         *pfpsf |= INVALID_EXCEPTION;
1201         // return Integer Indefinite
1202         res = 0x8000000000000000ull;
1203         BID_RETURN (res);
1204       }
1205       // else cases that can be rounded to a 64-bit int fall through
1206       // to '1 <= q + exp <= 20'
1207     } else if (q <= 19) {
1208       // C * 10^(21-q) >= 0xa0000000000000000
1209       __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
1210       if (C.w[1] >= 0x0a) {
1211         // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
1212         // set invalid flag
1213         *pfpsf |= INVALID_EXCEPTION;
1214         // return Integer Indefinite
1215         res = 0x8000000000000000ull;
1216         BID_RETURN (res);
1217       }
1218       // else cases that can be rounded to a 64-bit int fall through
1219       // to '1 <= q + exp <= 20'
1220     } else if (q == 20) {
1221       // C >= 0x10000000000000000
1222       if (C1.w[1] >= 0x01) {
1223         // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
1224         // set invalid flag
1225         *pfpsf |= INVALID_EXCEPTION;
1226         // return Integer Indefinite
1227         res = 0x8000000000000000ull;
1228         BID_RETURN (res);
1229       }
1230       // else cases that can be rounded to a 64-bit int fall through
1231       // to '1 <= q + exp <= 20'
1232     } else if (q == 21) {
1233       // C >= 0xa0000000000000000
1234       if (C1.w[1] >= 0x0a) {
1235         // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
1236         // set invalid flag
1237         *pfpsf |= INVALID_EXCEPTION;
1238         // return Integer Indefinite
1239         res = 0x8000000000000000ull;
1240         BID_RETURN (res);
1241       }
1242       // else cases that can be rounded to a 64-bit int fall through
1243       // to '1 <= q + exp <= 20'
1244     } else {    // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1245       // C  >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
1246       C.w[1] = 0x0a;
1247       C.w[0] = 0x0000000000000000ull;
1248       __mul_128x64_to_128 (C, ten2k64[q - 21], C);
1249       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
1250         // set invalid flag
1251         *pfpsf |= INVALID_EXCEPTION;
1252         // return Integer Indefinite
1253         res = 0x8000000000000000ull;
1254         BID_RETURN (res);
1255       }
1256       // else cases that can be rounded to a 64-bit int fall through
1257       // to '1 <= q + exp <= 20'
1258     }
1259   }
1260   // n is not too large to be converted to int64 if 0 <= n < 2^64
1261   // Note: some of the cases tested for above fall through to this point
1262   if ((q + exp) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
1263     // set inexact flag
1264     *pfpsf |= INEXACT_EXCEPTION;
1265     // return 0
1266     res = 0x0000000000000000ull;
1267     BID_RETURN (res);
1268   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1269     // 1 <= x < 2^64 so x can be rounded
1270     // down to a 64-bit unsigned signed integer
1271     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1272       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1273       // chop off ind digits from the lower part of C1
1274       // C1 fits in 127 bits
1275       // calculate C* and f*
1276       // C* is actually floor(C*) in this case
1277       // C* and f* need shifting and masking, as shown by
1278       // shiftright128[] and maskhigh128[]
1279       // 1 <= x <= 33
1280       // kx = 10^(-x) = ten2mk128[ind - 1]
1281       // C* = C1 * 10^(-x)
1282       // the approximation of 10^(-x) was rounded up to 118 bits
1283       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1284       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1285         Cstar.w[1] = P256.w[3];
1286         Cstar.w[0] = P256.w[2];
1287         fstar.w[3] = 0;
1288         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1289         fstar.w[1] = P256.w[1];
1290         fstar.w[0] = P256.w[0];
1291       } else {  // 22 <= ind - 1 <= 33
1292         Cstar.w[1] = 0;
1293         Cstar.w[0] = P256.w[3];
1294         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1295         fstar.w[2] = P256.w[2];
1296         fstar.w[1] = P256.w[1];
1297         fstar.w[0] = P256.w[0];
1298       }
1299       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1300       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1301       // C* = floor(C*) (logical right shift; C has p decimal digits,
1302       //     correct by Property 1)
1303       // n = C* * 10^(e+x)
1304
1305       // shift right C* by Ex-128 = shiftright128[ind]
1306       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
1307       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1308         Cstar.w[0] =
1309           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1310         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1311       } else {  // 22 <= ind - 1 <= 33
1312         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
1313       }
1314       // determine inexactness of the rounding of C*
1315       // if (0 < f* < 10^(-x)) then
1316       //   the result is exact
1317       // else // if (f* > T*) then
1318       //   the result is inexact
1319       if (ind - 1 <= 2) {
1320         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1] ||
1321             (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] &&
1322              fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1323           // set the inexact flag
1324           *pfpsf |= INEXACT_EXCEPTION;
1325         }       // else the result is exact
1326       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
1327         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1328             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1329                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1330           // set the inexact flag
1331           *pfpsf |= INEXACT_EXCEPTION;
1332         }       // else the result is exact
1333       } else {  // if 22 <= ind <= 33
1334         if (fstar.w[3] || fstar.w[2]
1335             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1336             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1337                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1338           // set the inexact flag
1339           *pfpsf |= INEXACT_EXCEPTION;
1340         }       // else the result is exact
1341       }
1342
1343       res = Cstar.w[0]; // the result is positive
1344     } else if (exp == 0) {
1345       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
1346       // res = C (exact)
1347       res = C1.w[0];
1348     } else {
1349       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
1350       // res = C * 10^exp (exact) - must fit in 64 bits
1351       res = C1.w[0] * ten2k64[exp];
1352     }
1353   }
1354 }
1355
1356 BID_RETURN (res);
1357 }
1358
1359 /*****************************************************************************
1360  *  BID128_to_uint64_ceil
1361  ****************************************************************************/
1362
1363 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64, bid128_to_uint64_ceil,
1364                                           x)
1365
1366      UINT64 res;
1367      UINT64 x_sign;
1368      UINT64 x_exp;
1369      int exp;                   // unbiased exponent
1370   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1371      BID_UI64DOUBLE tmp1;
1372      unsigned int x_nr_bits;
1373      int q, ind, shift;
1374      UINT128 C1, C;
1375      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
1376      UINT256 fstar;
1377      UINT256 P256;
1378
1379   // unpack x
1380 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1381 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
1382 C1.w[1] = x.w[1] & MASK_COEFF;
1383 C1.w[0] = x.w[0];
1384
1385   // check for NaN or Infinity
1386 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1387     // x is special
1388 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1389   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1390     // set invalid flag
1391     *pfpsf |= INVALID_EXCEPTION;
1392     // return Integer Indefinite
1393     res = 0x8000000000000000ull;
1394   } else {      // x is QNaN
1395     // set invalid flag
1396     *pfpsf |= INVALID_EXCEPTION;
1397     // return Integer Indefinite
1398     res = 0x8000000000000000ull;
1399   }
1400   BID_RETURN (res);
1401 } else {        // x is not a NaN, so it must be infinity
1402   if (!x_sign) {        // x is +inf
1403     // set invalid flag
1404     *pfpsf |= INVALID_EXCEPTION;
1405     // return Integer Indefinite
1406     res = 0x8000000000000000ull;
1407   } else {      // x is -inf
1408     // set invalid flag
1409     *pfpsf |= INVALID_EXCEPTION;
1410     // return Integer Indefinite
1411     res = 0x8000000000000000ull;
1412   }
1413   BID_RETURN (res);
1414 }
1415 }
1416   // check for non-canonical values (after the check for special values)
1417 if ((C1.w[1] > 0x0001ed09bead87c0ull)
1418     || (C1.w[1] == 0x0001ed09bead87c0ull
1419         && (C1.w[0] > 0x378d8e63ffffffffull))
1420     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1421   res = 0x0000000000000000ull;
1422   BID_RETURN (res);
1423 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1424   // x is 0
1425   res = 0x0000000000000000ull;
1426   BID_RETURN (res);
1427 } else {        // x is not special and is not zero
1428
1429   // q = nr. of decimal digits in x
1430   //  determine first the nr. of bits in x
1431   if (C1.w[1] == 0) {
1432     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
1433       // split the 64-bit value in two 32-bit halves to avoid rounding errors
1434       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
1435         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
1436         x_nr_bits =
1437           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1438       } else {  // x < 2^32
1439         tmp1.d = (double) (C1.w[0]);    // exact conversion
1440         x_nr_bits =
1441           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1442       }
1443     } else {    // if x < 2^53
1444       tmp1.d = (double) C1.w[0];        // exact conversion
1445       x_nr_bits =
1446         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1447     }
1448   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1449     tmp1.d = (double) C1.w[1];  // exact conversion
1450     x_nr_bits =
1451       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1452   }
1453   q = nr_digits[x_nr_bits - 1].digits;
1454   if (q == 0) {
1455     q = nr_digits[x_nr_bits - 1].digits1;
1456     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1457         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1458             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1459       q++;
1460   }
1461   exp = (x_exp >> 49) - 6176;
1462   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1463     // set invalid flag
1464     *pfpsf |= INVALID_EXCEPTION;
1465     // return Integer Indefinite
1466     res = 0x8000000000000000ull;
1467     BID_RETURN (res);
1468   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1469     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1470     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1471     // the cases that do not fit are identified here; the ones that fit
1472     // fall through and will be handled with other cases further,
1473     // under '1 <= q + exp <= 20'
1474     if (x_sign) {       // if n < 0 and q + exp = 20
1475       // if n <= -1 then n cannot be converted to uint64 with RZ
1476       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
1477       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
1478       // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
1479       if (q == 21) {
1480         // C >= a
1481         if (C1.w[1] != 0 || C1.w[0] >= 0x0aull) {
1482           // set invalid flag
1483           *pfpsf |= INVALID_EXCEPTION;
1484           // return Integer Indefinite
1485           res = 0x8000000000000000ull;
1486           BID_RETURN (res);
1487         }
1488         // else cases that can be rounded to 64-bit unsigned int fall through
1489         // to '1 <= q + exp <= 20'
1490       } else {
1491         // if 1 <= q <= 20
1492         //   C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
1493         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1494         //  C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
1495         // set invalid flag
1496         *pfpsf |= INVALID_EXCEPTION;
1497         // return Integer Indefinite
1498         res = 0x8000000000000000ull;
1499         BID_RETURN (res);
1500       }
1501     } else {    // if n > 0 and q + exp = 20
1502       // if n > 2^64 - 1 then n is too large
1503       // <=> c(0)c(1)...c(19).c(20)...c(q-1) > 2^64 - 1
1504       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 > 2^64 - 1
1505       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 > 10 * (2^64 - 1)
1506       // <=> C * 10^(21-q) > 0x9fffffffffffffff6, 1<=q<=34
1507       if (q == 1) {
1508         // C * 10^20 > 0x9fffffffffffffff6
1509         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
1510         if (C.w[1] > 0x09 || (C.w[1] == 0x09
1511                               && C.w[0] > 0xfffffffffffffff6ull)) {
1512           // set invalid flag
1513           *pfpsf |= INVALID_EXCEPTION;
1514           // return Integer Indefinite
1515           res = 0x8000000000000000ull;
1516           BID_RETURN (res);
1517         }
1518         // else cases that can be rounded to a 64-bit int fall through
1519         // to '1 <= q + exp <= 20'
1520       } else if (q <= 19) {
1521         // C * 10^(21-q) > 0x9fffffffffffffff6
1522         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
1523         if (C.w[1] > 0x09 || (C.w[1] == 0x09
1524                               && C.w[0] > 0xfffffffffffffff6ull)) {
1525           // set invalid flag
1526           *pfpsf |= INVALID_EXCEPTION;
1527           // return Integer Indefinite
1528           res = 0x8000000000000000ull;
1529           BID_RETURN (res);
1530         }
1531         // else cases that can be rounded to a 64-bit int fall through
1532         // to '1 <= q + exp <= 20'
1533       } else if (q == 20) {
1534         // C > 0xffffffffffffffff
1535         if (C1.w[1]) {
1536           // set invalid flag
1537           *pfpsf |= INVALID_EXCEPTION;
1538           // return Integer Indefinite
1539           res = 0x8000000000000000ull;
1540           BID_RETURN (res);
1541         }
1542         // else cases that can be rounded to a 64-bit int fall through
1543         // to '1 <= q + exp <= 20'
1544       } else if (q == 21) {
1545         // C > 0x9fffffffffffffff6
1546         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
1547                                && C1.w[0] > 0xfffffffffffffff6ull)) {
1548           // set invalid flag
1549           *pfpsf |= INVALID_EXCEPTION;
1550           // return Integer Indefinite
1551           res = 0x8000000000000000ull;
1552           BID_RETURN (res);
1553         }
1554         // else cases that can be rounded to a 64-bit int fall through
1555         // to '1 <= q + exp <= 20'
1556       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1557         // C  > 10^(q-21) * 0x9fffffffffffffff6 max 44 bits x 68 bits
1558         C.w[1] = 0x09;
1559         C.w[0] = 0xfffffffffffffff6ull;
1560         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
1561         if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
1562           // set invalid flag
1563           *pfpsf |= INVALID_EXCEPTION;
1564           // return Integer Indefinite
1565           res = 0x8000000000000000ull;
1566           BID_RETURN (res);
1567         }
1568         // else cases that can be rounded to a 64-bit int fall through
1569         // to '1 <= q + exp <= 20'
1570       }
1571     }
1572   }
1573   // n is not too large to be converted to int64 if -1 < n <= 2^64 - 1
1574   // Note: some of the cases tested for above fall through to this point
1575   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1576     // return 0 or 1
1577     if (x_sign)
1578       res = 0x0000000000000000ull;
1579     else
1580       res = 0x0000000000000001ull;
1581     BID_RETURN (res);
1582   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1583     // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
1584     // to zero to a 64-bit unsigned signed integer
1585     if (x_sign) {       // x <= -1
1586       // set invalid flag
1587       *pfpsf |= INVALID_EXCEPTION;
1588       // return Integer Indefinite
1589       res = 0x8000000000000000ull;
1590       BID_RETURN (res);
1591     }
1592     // 1 <= x <= 2^64 - 1 so x can be rounded
1593     // to zero to a 64-bit unsigned integer
1594     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1595       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1596       // chop off ind digits from the lower part of C1
1597       // C1 fits in 127 bits
1598       // calculate C* and f*
1599       // C* is actually floor(C*) in this case
1600       // C* and f* need shifting and masking, as shown by
1601       // shiftright128[] and maskhigh128[]
1602       // 1 <= x <= 33
1603       // kx = 10^(-x) = ten2mk128[ind - 1]
1604       // C* = C1 * 10^(-x)
1605       // the approximation of 10^(-x) was rounded up to 118 bits
1606       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1607       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1608         Cstar.w[1] = P256.w[3];
1609         Cstar.w[0] = P256.w[2];
1610         fstar.w[3] = 0;
1611         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1612         fstar.w[1] = P256.w[1];
1613         fstar.w[0] = P256.w[0];
1614       } else {  // 22 <= ind - 1 <= 33
1615         Cstar.w[1] = 0;
1616         Cstar.w[0] = P256.w[3];
1617         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1618         fstar.w[2] = P256.w[2];
1619         fstar.w[1] = P256.w[1];
1620         fstar.w[0] = P256.w[0];
1621       }
1622       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1623       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1624       // C* = floor(C*) (logical right shift; C has p decimal digits,
1625       //     correct by Property 1)
1626       // n = C* * 10^(e+x)
1627
1628       // shift right C* by Ex-128 = shiftright128[ind]
1629       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
1630       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1631         Cstar.w[0] =
1632           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1633         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1634       } else {  // 22 <= ind - 1 <= 33
1635         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
1636       }
1637       // if the result is positive and inexact, need to add 1 to it
1638
1639       // determine inexactness of the rounding of C*
1640       // if (0 < f* < 10^(-x)) then
1641       //   the result is exact
1642       // else // if (f* > T*) then
1643       //   the result is inexact
1644       if (ind - 1 <= 2) {
1645         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1646             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1647                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1648           if (!x_sign) {        // positive and inexact
1649             Cstar.w[0]++;
1650             if (Cstar.w[0] == 0x0)
1651               Cstar.w[1]++;
1652           }
1653         }       // else the result is exact
1654       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
1655         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1656             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1657                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1658           if (!x_sign) {        // positive and inexact
1659             Cstar.w[0]++;
1660             if (Cstar.w[0] == 0x0)
1661               Cstar.w[1]++;
1662           }
1663         }       // else the result is exact
1664       } else {  // if 22 <= ind <= 33
1665         if (fstar.w[3] || fstar.w[2]
1666             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1667             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1668                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1669           if (!x_sign) {        // positive and inexact
1670             Cstar.w[0]++;
1671             if (Cstar.w[0] == 0x0)
1672               Cstar.w[1]++;
1673           }
1674         }       // else the result is exact
1675       }
1676
1677       res = Cstar.w[0]; // the result is positive
1678     } else if (exp == 0) {
1679       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
1680       // res = C (exact)
1681       res = C1.w[0];
1682     } else {
1683       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
1684       // res = C * 10^exp (exact) - must fit in 64 bits
1685       res = C1.w[0] * ten2k64[exp];
1686     }
1687   }
1688 }
1689
1690 BID_RETURN (res);
1691 }
1692
1693 /*****************************************************************************
1694  *  BID128_to_uint64_xceil
1695  ****************************************************************************/
1696
1697 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
1698                                           bid128_to_uint64_xceil, x)
1699
1700      UINT64 res;
1701      UINT64 x_sign;
1702      UINT64 x_exp;
1703      int exp;                   // unbiased exponent
1704   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1705      BID_UI64DOUBLE tmp1;
1706      unsigned int x_nr_bits;
1707      int q, ind, shift;
1708      UINT128 C1, C;
1709      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
1710      UINT256 fstar;
1711      UINT256 P256;
1712
1713   // unpack x
1714 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1715 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
1716 C1.w[1] = x.w[1] & MASK_COEFF;
1717 C1.w[0] = x.w[0];
1718
1719   // check for NaN or Infinity
1720 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1721     // x is special
1722 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1723   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1724     // set invalid flag
1725     *pfpsf |= INVALID_EXCEPTION;
1726     // return Integer Indefinite
1727     res = 0x8000000000000000ull;
1728   } else {      // x is QNaN
1729     // set invalid flag
1730     *pfpsf |= INVALID_EXCEPTION;
1731     // return Integer Indefinite
1732     res = 0x8000000000000000ull;
1733   }
1734   BID_RETURN (res);
1735 } else {        // x is not a NaN, so it must be infinity
1736   if (!x_sign) {        // x is +inf
1737     // set invalid flag
1738     *pfpsf |= INVALID_EXCEPTION;
1739     // return Integer Indefinite
1740     res = 0x8000000000000000ull;
1741   } else {      // x is -inf
1742     // set invalid flag
1743     *pfpsf |= INVALID_EXCEPTION;
1744     // return Integer Indefinite
1745     res = 0x8000000000000000ull;
1746   }
1747   BID_RETURN (res);
1748 }
1749 }
1750   // check for non-canonical values (after the check for special values)
1751 if ((C1.w[1] > 0x0001ed09bead87c0ull)
1752     || (C1.w[1] == 0x0001ed09bead87c0ull
1753         && (C1.w[0] > 0x378d8e63ffffffffull))
1754     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1755   res = 0x0000000000000000ull;
1756   BID_RETURN (res);
1757 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1758   // x is 0
1759   res = 0x0000000000000000ull;
1760   BID_RETURN (res);
1761 } else {        // x is not special and is not zero
1762
1763   // q = nr. of decimal digits in x
1764   //  determine first the nr. of bits in x
1765   if (C1.w[1] == 0) {
1766     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
1767       // split the 64-bit value in two 32-bit halves to avoid rounding errors
1768       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
1769         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
1770         x_nr_bits =
1771           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1772       } else {  // x < 2^32
1773         tmp1.d = (double) (C1.w[0]);    // exact conversion
1774         x_nr_bits =
1775           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1776       }
1777     } else {    // if x < 2^53
1778       tmp1.d = (double) C1.w[0];        // exact conversion
1779       x_nr_bits =
1780         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1781     }
1782   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1783     tmp1.d = (double) C1.w[1];  // exact conversion
1784     x_nr_bits =
1785       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1786   }
1787   q = nr_digits[x_nr_bits - 1].digits;
1788   if (q == 0) {
1789     q = nr_digits[x_nr_bits - 1].digits1;
1790     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1791         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1792             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1793       q++;
1794   }
1795   exp = (x_exp >> 49) - 6176;
1796   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1797     // set invalid flag
1798     *pfpsf |= INVALID_EXCEPTION;
1799     // return Integer Indefinite
1800     res = 0x8000000000000000ull;
1801     BID_RETURN (res);
1802   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1803     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1804     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1805     // the cases that do not fit are identified here; the ones that fit
1806     // fall through and will be handled with other cases further,
1807     // under '1 <= q + exp <= 20'
1808     if (x_sign) {       // if n < 0 and q + exp = 20
1809       // if n <= -1 then n cannot be converted to uint64 with RZ
1810       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
1811       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
1812       // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
1813       if (q == 21) {
1814         // C >= a
1815         if (C1.w[1] != 0 || C1.w[0] >= 0x0aull) {
1816           // set invalid flag
1817           *pfpsf |= INVALID_EXCEPTION;
1818           // return Integer Indefinite
1819           res = 0x8000000000000000ull;
1820           BID_RETURN (res);
1821         }
1822         // else cases that can be rounded to 64-bit unsigned int fall through
1823         // to '1 <= q + exp <= 20'
1824       } else {
1825         // if 1 <= q <= 20
1826         //   C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
1827         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1828         //  C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
1829         // set invalid flag
1830         *pfpsf |= INVALID_EXCEPTION;
1831         // return Integer Indefinite
1832         res = 0x8000000000000000ull;
1833         BID_RETURN (res);
1834       }
1835     } else {    // if n > 0 and q + exp = 20
1836       // if n > 2^64 - 1 then n is too large
1837       // <=> c(0)c(1)...c(19).c(20)...c(q-1) > 2^64 - 1
1838       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 > 2^64 - 1
1839       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 > 10 * (2^64 - 1)
1840       // <=> C * 10^(21-q) > 0x9fffffffffffffff6, 1<=q<=34
1841       if (q == 1) {
1842         // C * 10^20 > 0x9fffffffffffffff6
1843         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
1844         if (C.w[1] > 0x09 || (C.w[1] == 0x09
1845                               && C.w[0] > 0xfffffffffffffff6ull)) {
1846           // set invalid flag
1847           *pfpsf |= INVALID_EXCEPTION;
1848           // return Integer Indefinite
1849           res = 0x8000000000000000ull;
1850           BID_RETURN (res);
1851         }
1852         // else cases that can be rounded to a 64-bit int fall through
1853         // to '1 <= q + exp <= 20'
1854       } else if (q <= 19) {
1855         // C * 10^(21-q) > 0x9fffffffffffffff6
1856         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
1857         if (C.w[1] > 0x09 || (C.w[1] == 0x09
1858                               && C.w[0] > 0xfffffffffffffff6ull)) {
1859           // set invalid flag
1860           *pfpsf |= INVALID_EXCEPTION;
1861           // return Integer Indefinite
1862           res = 0x8000000000000000ull;
1863           BID_RETURN (res);
1864         }
1865         // else cases that can be rounded to a 64-bit int fall through
1866         // to '1 <= q + exp <= 20'
1867       } else if (q == 20) {
1868         // C > 0xffffffffffffffff
1869         if (C1.w[1]) {
1870           // set invalid flag
1871           *pfpsf |= INVALID_EXCEPTION;
1872           // return Integer Indefinite
1873           res = 0x8000000000000000ull;
1874           BID_RETURN (res);
1875         }
1876         // else cases that can be rounded to a 64-bit int fall through
1877         // to '1 <= q + exp <= 20'
1878       } else if (q == 21) {
1879         // C > 0x9fffffffffffffff6
1880         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
1881                                && C1.w[0] > 0xfffffffffffffff6ull)) {
1882           // set invalid flag
1883           *pfpsf |= INVALID_EXCEPTION;
1884           // return Integer Indefinite
1885           res = 0x8000000000000000ull;
1886           BID_RETURN (res);
1887         }
1888         // else cases that can be rounded to a 64-bit int fall through
1889         // to '1 <= q + exp <= 20'
1890       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1891         // C  > 10^(q-21) * 0x9fffffffffffffff6 max 44 bits x 68 bits
1892         C.w[1] = 0x09;
1893         C.w[0] = 0xfffffffffffffff6ull;
1894         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
1895         if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
1896           // set invalid flag
1897           *pfpsf |= INVALID_EXCEPTION;
1898           // return Integer Indefinite
1899           res = 0x8000000000000000ull;
1900           BID_RETURN (res);
1901         }
1902         // else cases that can be rounded to a 64-bit int fall through
1903         // to '1 <= q + exp <= 20'
1904       }
1905     }
1906   }
1907   // n is not too large to be converted to int64 if -1 < n <= 2^64 - 1
1908   // Note: some of the cases tested for above fall through to this point
1909   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1910     // set inexact flag
1911     *pfpsf |= INEXACT_EXCEPTION;
1912     // return 0 or 1
1913     if (x_sign)
1914       res = 0x0000000000000000ull;
1915     else
1916       res = 0x0000000000000001ull;
1917     BID_RETURN (res);
1918   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1919     // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
1920     // to zero to a 64-bit unsigned signed integer
1921     if (x_sign) {       // x <= -1
1922       // set invalid flag
1923       *pfpsf |= INVALID_EXCEPTION;
1924       // return Integer Indefinite
1925       res = 0x8000000000000000ull;
1926       BID_RETURN (res);
1927     }
1928     // 1 <= x <= 2^64 - 1 so x can be rounded
1929     // to zero to a 64-bit unsigned integer
1930     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1931       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1932       // chop off ind digits from the lower part of C1
1933       // C1 fits in 127 bits
1934       // calculate C* and f*
1935       // C* is actually floor(C*) in this case
1936       // C* and f* need shifting and masking, as shown by
1937       // shiftright128[] and maskhigh128[]
1938       // 1 <= x <= 33
1939       // kx = 10^(-x) = ten2mk128[ind - 1]
1940       // C* = C1 * 10^(-x)
1941       // the approximation of 10^(-x) was rounded up to 118 bits
1942       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1943       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1944         Cstar.w[1] = P256.w[3];
1945         Cstar.w[0] = P256.w[2];
1946         fstar.w[3] = 0;
1947         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1948         fstar.w[1] = P256.w[1];
1949         fstar.w[0] = P256.w[0];
1950       } else {  // 22 <= ind - 1 <= 33
1951         Cstar.w[1] = 0;
1952         Cstar.w[0] = P256.w[3];
1953         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1954         fstar.w[2] = P256.w[2];
1955         fstar.w[1] = P256.w[1];
1956         fstar.w[0] = P256.w[0];
1957       }
1958       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1959       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1960       // C* = floor(C*) (logical right shift; C has p decimal digits,
1961       //     correct by Property 1)
1962       // n = C* * 10^(e+x)
1963
1964       // shift right C* by Ex-128 = shiftright128[ind]
1965       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
1966       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1967         Cstar.w[0] =
1968           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1969         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1970       } else {  // 22 <= ind - 1 <= 33
1971         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
1972       }
1973       // if the result is positive and inexact, need to add 1 to it
1974
1975       // determine inexactness of the rounding of C*
1976       // if (0 < f* < 10^(-x)) then
1977       //   the result is exact
1978       // else // if (f* > T*) then
1979       //   the result is inexact
1980       if (ind - 1 <= 2) {
1981         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1982             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1983                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1984           if (!x_sign) {        // positive and inexact
1985             Cstar.w[0]++;
1986             if (Cstar.w[0] == 0x0)
1987               Cstar.w[1]++;
1988           }
1989           // set the inexact flag
1990           *pfpsf |= INEXACT_EXCEPTION;
1991         }       // else the result is exact
1992       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
1993         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1994             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1995                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1996           if (!x_sign) {        // positive and inexact
1997             Cstar.w[0]++;
1998             if (Cstar.w[0] == 0x0)
1999               Cstar.w[1]++;
2000           }
2001           // set the inexact flag
2002           *pfpsf |= INEXACT_EXCEPTION;
2003         }       // else the result is exact
2004       } else {  // if 22 <= ind <= 33
2005         if (fstar.w[3] || fstar.w[2]
2006             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2007             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2008                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2009           if (!x_sign) {        // positive and inexact
2010             Cstar.w[0]++;
2011             if (Cstar.w[0] == 0x0)
2012               Cstar.w[1]++;
2013           }
2014           // set the inexact flag
2015           *pfpsf |= INEXACT_EXCEPTION;
2016         }       // else the result is exact
2017       }
2018
2019       res = Cstar.w[0]; // the result is positive
2020     } else if (exp == 0) {
2021       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2022       // res = C (exact)
2023       res = C1.w[0];
2024     } else {
2025       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2026       // res = C * 10^exp (exact) - must fit in 64 bits
2027       res = C1.w[0] * ten2k64[exp];
2028     }
2029   }
2030 }
2031
2032 BID_RETURN (res);
2033 }
2034
2035 /*****************************************************************************
2036  *  BID128_to_uint64_int
2037  ****************************************************************************/
2038
2039 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64, bid128_to_uint64_int,
2040                                           x)
2041
2042      UINT64 res;
2043      UINT64 x_sign;
2044      UINT64 x_exp;
2045      int exp;                   // unbiased exponent
2046   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2047      BID_UI64DOUBLE tmp1;
2048      unsigned int x_nr_bits;
2049      int q, ind, shift;
2050      UINT128 C1, C;
2051      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
2052      UINT256 P256;
2053
2054   // unpack x
2055 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
2056 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
2057 C1.w[1] = x.w[1] & MASK_COEFF;
2058 C1.w[0] = x.w[0];
2059
2060   // check for NaN or Infinity
2061 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2062     // x is special
2063 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
2064   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
2065     // set invalid flag
2066     *pfpsf |= INVALID_EXCEPTION;
2067     // return Integer Indefinite
2068     res = 0x8000000000000000ull;
2069   } else {      // x is QNaN
2070     // set invalid flag
2071     *pfpsf |= INVALID_EXCEPTION;
2072     // return Integer Indefinite
2073     res = 0x8000000000000000ull;
2074   }
2075   BID_RETURN (res);
2076 } else {        // x is not a NaN, so it must be infinity
2077   if (!x_sign) {        // x is +inf
2078     // set invalid flag
2079     *pfpsf |= INVALID_EXCEPTION;
2080     // return Integer Indefinite
2081     res = 0x8000000000000000ull;
2082   } else {      // x is -inf
2083     // set invalid flag
2084     *pfpsf |= INVALID_EXCEPTION;
2085     // return Integer Indefinite
2086     res = 0x8000000000000000ull;
2087   }
2088   BID_RETURN (res);
2089 }
2090 }
2091   // check for non-canonical values (after the check for special values)
2092 if ((C1.w[1] > 0x0001ed09bead87c0ull)
2093     || (C1.w[1] == 0x0001ed09bead87c0ull
2094         && (C1.w[0] > 0x378d8e63ffffffffull))
2095     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2096   res = 0x0000000000000000ull;
2097   BID_RETURN (res);
2098 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2099   // x is 0
2100   res = 0x0000000000000000ull;
2101   BID_RETURN (res);
2102 } else {        // x is not special and is not zero
2103
2104   // q = nr. of decimal digits in x
2105   //  determine first the nr. of bits in x
2106   if (C1.w[1] == 0) {
2107     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
2108       // split the 64-bit value in two 32-bit halves to avoid rounding errors
2109       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
2110         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
2111         x_nr_bits =
2112           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2113       } else {  // x < 2^32
2114         tmp1.d = (double) (C1.w[0]);    // exact conversion
2115         x_nr_bits =
2116           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2117       }
2118     } else {    // if x < 2^53
2119       tmp1.d = (double) C1.w[0];        // exact conversion
2120       x_nr_bits =
2121         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2122     }
2123   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2124     tmp1.d = (double) C1.w[1];  // exact conversion
2125     x_nr_bits =
2126       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2127   }
2128   q = nr_digits[x_nr_bits - 1].digits;
2129   if (q == 0) {
2130     q = nr_digits[x_nr_bits - 1].digits1;
2131     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2132         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2133             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2134       q++;
2135   }
2136   exp = (x_exp >> 49) - 6176;
2137
2138   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
2139     // set invalid flag
2140     *pfpsf |= INVALID_EXCEPTION;
2141     // return Integer Indefinite
2142     res = 0x8000000000000000ull;
2143     BID_RETURN (res);
2144   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
2145     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2146     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
2147     // the cases that do not fit are identified here; the ones that fit
2148     // fall through and will be handled with other cases further,
2149     // under '1 <= q + exp <= 20'
2150     if (x_sign) {       // if n < 0 and q + exp = 20
2151       // if n <= -1 then n cannot be converted to uint64 with RZ
2152       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
2153       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
2154       // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
2155       if (q == 21) {
2156         // C >= a
2157         if (C1.w[1] != 0 || C1.w[0] >= 0x0aull) {
2158           // set invalid flag
2159           *pfpsf |= INVALID_EXCEPTION;
2160           // return Integer Indefinite
2161           res = 0x8000000000000000ull;
2162           BID_RETURN (res);
2163         }
2164         // else cases that can be rounded to 64-bit unsigned int fall through
2165         // to '1 <= q + exp <= 20'
2166       } else {
2167         // if 1 <= q <= 20
2168         //   C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
2169         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2170         //  C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
2171         // set invalid flag
2172         *pfpsf |= INVALID_EXCEPTION;
2173         // return Integer Indefinite
2174         res = 0x8000000000000000ull;
2175         BID_RETURN (res);
2176       }
2177     } else {    // if n > 0 and q + exp = 20
2178       // if n >= 2^64 then n is too large
2179       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
2180       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
2181       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
2182       // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
2183       if (q == 1) {
2184         // C * 10^20 >= 0xa0000000000000000
2185         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
2186         if (C.w[1] >= 0x0a) {
2187           // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2188           // set invalid flag
2189           *pfpsf |= INVALID_EXCEPTION;
2190           // return Integer Indefinite
2191           res = 0x8000000000000000ull;
2192           BID_RETURN (res);
2193         }
2194         // else cases that can be rounded to a 64-bit int fall through
2195         // to '1 <= q + exp <= 20'
2196       } else if (q <= 19) {
2197         // C * 10^(21-q) >= 0xa0000000000000000
2198         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
2199         if (C.w[1] >= 0x0a) {
2200           // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2201           // set invalid flag
2202           *pfpsf |= INVALID_EXCEPTION;
2203           // return Integer Indefinite
2204           res = 0x8000000000000000ull;
2205           BID_RETURN (res);
2206         }
2207         // else cases that can be rounded to a 64-bit int fall through
2208         // to '1 <= q + exp <= 20'
2209       } else if (q == 20) {
2210         // C >= 0x10000000000000000
2211         if (C1.w[1] >= 0x01) {
2212           // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
2213           // set invalid flag
2214           *pfpsf |= INVALID_EXCEPTION;
2215           // return Integer Indefinite
2216           res = 0x8000000000000000ull;
2217           BID_RETURN (res);
2218         }
2219         // else cases that can be rounded to a 64-bit int fall through
2220         // to '1 <= q + exp <= 20'
2221       } else if (q == 21) {
2222         // C >= 0xa0000000000000000
2223         if (C1.w[1] >= 0x0a) {
2224           // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
2225           // set invalid flag
2226           *pfpsf |= INVALID_EXCEPTION;
2227           // return Integer Indefinite
2228           res = 0x8000000000000000ull;
2229           BID_RETURN (res);
2230         }
2231         // else cases that can be rounded to a 64-bit int fall through
2232         // to '1 <= q + exp <= 20'
2233       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2234         // C  >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
2235         C.w[1] = 0x0a;
2236         C.w[0] = 0x0000000000000000ull;
2237         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
2238         if (C1.w[1] > C.w[1]
2239             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2240           // set invalid flag
2241           *pfpsf |= INVALID_EXCEPTION;
2242           // return Integer Indefinite
2243           res = 0x8000000000000000ull;
2244           BID_RETURN (res);
2245         }
2246         // else cases that can be rounded to a 64-bit int fall through
2247         // to '1 <= q + exp <= 20'
2248       }
2249     }
2250   }
2251   // n is not too large to be converted to int64 if -1 < n < 2^64
2252   // Note: some of the cases tested for above fall through to this point
2253   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
2254     // return 0
2255     res = 0x0000000000000000ull;
2256     BID_RETURN (res);
2257   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
2258     // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
2259     // to zero to a 64-bit unsigned signed integer
2260     if (x_sign) {       // x <= -1
2261       // set invalid flag
2262       *pfpsf |= INVALID_EXCEPTION;
2263       // return Integer Indefinite
2264       res = 0x8000000000000000ull;
2265       BID_RETURN (res);
2266     }
2267     // 1 <= x < 2^64 so x can be rounded
2268     // to zero to a 64-bit unsigned integer
2269     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
2270       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
2271       // chop off ind digits from the lower part of C1
2272       // C1 fits in 127 bits
2273       // calculate C* and f*
2274       // C* is actually floor(C*) in this case
2275       // C* and f* need shifting and masking, as shown by
2276       // shiftright128[] and maskhigh128[]
2277       // 1 <= x <= 33
2278       // kx = 10^(-x) = ten2mk128[ind - 1]
2279       // C* = C1 * 10^(-x)
2280       // the approximation of 10^(-x) was rounded up to 118 bits
2281       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2282       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2283         Cstar.w[1] = P256.w[3];
2284         Cstar.w[0] = P256.w[2];
2285       } else {  // 22 <= ind - 1 <= 33
2286         Cstar.w[1] = 0;
2287         Cstar.w[0] = P256.w[3];
2288       }
2289       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2290       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2291       // C* = floor(C*) (logical right shift; C has p decimal digits,
2292       //     correct by Property 1)
2293       // n = C* * 10^(e+x)
2294
2295       // shift right C* by Ex-128 = shiftright128[ind]
2296       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
2297       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2298         Cstar.w[0] =
2299           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2300         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2301       } else {  // 22 <= ind - 1 <= 33
2302         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
2303       }
2304       res = Cstar.w[0]; // the result is positive
2305     } else if (exp == 0) {
2306       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2307       // res = C (exact)
2308       res = C1.w[0];
2309     } else {
2310       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2311       // res = C * 10^exp (exact) - must fit in 64 bits
2312       res = C1.w[0] * ten2k64[exp];
2313     }
2314   }
2315 }
2316
2317 BID_RETURN (res);
2318 }
2319
2320 /*****************************************************************************
2321  *  BID128_to_uint64_xint
2322  ****************************************************************************/
2323
2324 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64, bid128_to_uint64_xint,
2325                                           x)
2326
2327      UINT64 res;
2328      UINT64 x_sign;
2329      UINT64 x_exp;
2330      int exp;                   // unbiased exponent
2331   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2332      BID_UI64DOUBLE tmp1;
2333      unsigned int x_nr_bits;
2334      int q, ind, shift;
2335      UINT128 C1, C;
2336      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
2337      UINT256 fstar;
2338      UINT256 P256;
2339
2340   // unpack x
2341 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
2342 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
2343 C1.w[1] = x.w[1] & MASK_COEFF;
2344 C1.w[0] = x.w[0];
2345
2346   // check for NaN or Infinity
2347 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2348     // x is special
2349 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
2350   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
2351     // set invalid flag
2352     *pfpsf |= INVALID_EXCEPTION;
2353     // return Integer Indefinite
2354     res = 0x8000000000000000ull;
2355   } else {      // x is QNaN
2356     // set invalid flag
2357     *pfpsf |= INVALID_EXCEPTION;
2358     // return Integer Indefinite
2359     res = 0x8000000000000000ull;
2360   }
2361   BID_RETURN (res);
2362 } else {        // x is not a NaN, so it must be infinity
2363   if (!x_sign) {        // x is +inf
2364     // set invalid flag
2365     *pfpsf |= INVALID_EXCEPTION;
2366     // return Integer Indefinite
2367     res = 0x8000000000000000ull;
2368   } else {      // x is -inf
2369     // set invalid flag
2370     *pfpsf |= INVALID_EXCEPTION;
2371     // return Integer Indefinite
2372     res = 0x8000000000000000ull;
2373   }
2374   BID_RETURN (res);
2375 }
2376 }
2377   // check for non-canonical values (after the check for special values)
2378 if ((C1.w[1] > 0x0001ed09bead87c0ull)
2379     || (C1.w[1] == 0x0001ed09bead87c0ull
2380         && (C1.w[0] > 0x378d8e63ffffffffull))
2381     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2382   res = 0x0000000000000000ull;
2383   BID_RETURN (res);
2384 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2385   // x is 0
2386   res = 0x0000000000000000ull;
2387   BID_RETURN (res);
2388 } else {        // x is not special and is not zero
2389
2390   // q = nr. of decimal digits in x
2391   //  determine first the nr. of bits in x
2392   if (C1.w[1] == 0) {
2393     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
2394       // split the 64-bit value in two 32-bit halves to avoid rounding errors
2395       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
2396         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
2397         x_nr_bits =
2398           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2399       } else {  // x < 2^32
2400         tmp1.d = (double) (C1.w[0]);    // exact conversion
2401         x_nr_bits =
2402           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2403       }
2404     } else {    // if x < 2^53
2405       tmp1.d = (double) C1.w[0];        // exact conversion
2406       x_nr_bits =
2407         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2408     }
2409   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2410     tmp1.d = (double) C1.w[1];  // exact conversion
2411     x_nr_bits =
2412       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2413   }
2414   q = nr_digits[x_nr_bits - 1].digits;
2415   if (q == 0) {
2416     q = nr_digits[x_nr_bits - 1].digits1;
2417     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2418         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2419             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2420       q++;
2421   }
2422   exp = (x_exp >> 49) - 6176;
2423   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
2424     // set invalid flag
2425     *pfpsf |= INVALID_EXCEPTION;
2426     // return Integer Indefinite
2427     res = 0x8000000000000000ull;
2428     BID_RETURN (res);
2429   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
2430     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2431     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
2432     // the cases that do not fit are identified here; the ones that fit
2433     // fall through and will be handled with other cases further,
2434     // under '1 <= q + exp <= 20'
2435     if (x_sign) {       // if n < 0 and q + exp = 20
2436       // if n <= -1 then n cannot be converted to uint64 with RZ
2437       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
2438       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
2439       // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
2440       if (q == 21) {
2441         // C >= a
2442         if (C1.w[1] != 0 || C1.w[0] >= 0x0aull) {
2443           // set invalid flag
2444           *pfpsf |= INVALID_EXCEPTION;
2445           // return Integer Indefinite
2446           res = 0x8000000000000000ull;
2447           BID_RETURN (res);
2448         }
2449         // else cases that can be rounded to 64-bit unsigned int fall through
2450         // to '1 <= q + exp <= 20'
2451       } else {
2452         // if 1 <= q <= 20
2453         //   C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
2454         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2455         //  C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
2456         // set invalid flag
2457         *pfpsf |= INVALID_EXCEPTION;
2458         // return Integer Indefinite
2459         res = 0x8000000000000000ull;
2460         BID_RETURN (res);
2461       }
2462     } else {    // if n > 0 and q + exp = 20
2463       // if n >= 2^64 then n is too large
2464       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
2465       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
2466       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
2467       // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
2468       if (q == 1) {
2469         // C * 10^20 >= 0xa0000000000000000
2470         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
2471         if (C.w[1] >= 0x0a) {
2472           // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2473           // set invalid flag
2474           *pfpsf |= INVALID_EXCEPTION;
2475           // return Integer Indefinite
2476           res = 0x8000000000000000ull;
2477           BID_RETURN (res);
2478         }
2479         // else cases that can be rounded to a 64-bit int fall through
2480         // to '1 <= q + exp <= 20'
2481       } else if (q <= 19) {
2482         // C * 10^(21-q) >= 0xa0000000000000000
2483         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
2484         if (C.w[1] >= 0x0a) {
2485           // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2486           // set invalid flag
2487           *pfpsf |= INVALID_EXCEPTION;
2488           // return Integer Indefinite
2489           res = 0x8000000000000000ull;
2490           BID_RETURN (res);
2491         }
2492         // else cases that can be rounded to a 64-bit int fall through
2493         // to '1 <= q + exp <= 20'
2494       } else if (q == 20) {
2495         // C >= 0x10000000000000000
2496         if (C1.w[1] >= 0x01) {
2497           // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
2498           // set invalid flag
2499           *pfpsf |= INVALID_EXCEPTION;
2500           // return Integer Indefinite
2501           res = 0x8000000000000000ull;
2502           BID_RETURN (res);
2503         }
2504         // else cases that can be rounded to a 64-bit int fall through
2505         // to '1 <= q + exp <= 20'
2506       } else if (q == 21) {
2507         // C >= 0xa0000000000000000
2508         if (C1.w[1] >= 0x0a) {
2509           // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
2510           // set invalid flag
2511           *pfpsf |= INVALID_EXCEPTION;
2512           // return Integer Indefinite
2513           res = 0x8000000000000000ull;
2514           BID_RETURN (res);
2515         }
2516         // else cases that can be rounded to a 64-bit int fall through
2517         // to '1 <= q + exp <= 20'
2518       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2519         // C  >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
2520         C.w[1] = 0x0a;
2521         C.w[0] = 0x0000000000000000ull;
2522         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
2523         if (C1.w[1] > C.w[1]
2524             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2525           // set invalid flag
2526           *pfpsf |= INVALID_EXCEPTION;
2527           // return Integer Indefinite
2528           res = 0x8000000000000000ull;
2529           BID_RETURN (res);
2530         }
2531         // else cases that can be rounded to a 64-bit int fall through
2532         // to '1 <= q + exp <= 20'
2533       }
2534     }
2535   }
2536   // n is not too large to be converted to int64 if -1 < n < 2^64
2537   // Note: some of the cases tested for above fall through to this point
2538   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
2539     // set inexact flag
2540     *pfpsf |= INEXACT_EXCEPTION;
2541     // return 0
2542     res = 0x0000000000000000ull;
2543     BID_RETURN (res);
2544   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
2545     // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
2546     // to zero to a 64-bit unsigned signed integer
2547     if (x_sign) {       // x <= -1
2548       // set invalid flag
2549       *pfpsf |= INVALID_EXCEPTION;
2550       // return Integer Indefinite
2551       res = 0x8000000000000000ull;
2552       BID_RETURN (res);
2553     }
2554     // 1 <= x < 2^64 so x can be rounded
2555     // to zero to a 64-bit unsigned integer
2556     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
2557       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
2558       // chop off ind digits from the lower part of C1
2559       // C1 fits in 127 bits
2560       // calculate C* and f*
2561       // C* is actually floor(C*) in this case
2562       // C* and f* need shifting and masking, as shown by
2563       // shiftright128[] and maskhigh128[]
2564       // 1 <= x <= 33
2565       // kx = 10^(-x) = ten2mk128[ind - 1]
2566       // C* = C1 * 10^(-x)
2567       // the approximation of 10^(-x) was rounded up to 118 bits
2568       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2569       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2570         Cstar.w[1] = P256.w[3];
2571         Cstar.w[0] = P256.w[2];
2572         fstar.w[3] = 0;
2573         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
2574         fstar.w[1] = P256.w[1];
2575         fstar.w[0] = P256.w[0];
2576       } else {  // 22 <= ind - 1 <= 33
2577         Cstar.w[1] = 0;
2578         Cstar.w[0] = P256.w[3];
2579         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
2580         fstar.w[2] = P256.w[2];
2581         fstar.w[1] = P256.w[1];
2582         fstar.w[0] = P256.w[0];
2583       }
2584       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2585       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2586       // C* = floor(C*) (logical right shift; C has p decimal digits,
2587       //     correct by Property 1)
2588       // n = C* * 10^(e+x)
2589
2590       // shift right C* by Ex-128 = shiftright128[ind]
2591       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
2592       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2593         Cstar.w[0] =
2594           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2595         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2596       } else {  // 22 <= ind - 1 <= 33
2597         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
2598       }
2599       // determine inexactness of the rounding of C*
2600       // if (0 < f* < 10^(-x)) then
2601       //   the result is exact
2602       // else // if (f* > T*) then
2603       //   the result is inexact
2604       if (ind - 1 <= 2) {
2605         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2606             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2607                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2608           // set the inexact flag
2609           *pfpsf |= INEXACT_EXCEPTION;
2610         }       // else the result is exact
2611       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
2612         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2613             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2614                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2615           // set the inexact flag
2616           *pfpsf |= INEXACT_EXCEPTION;
2617         }       // else the result is exact
2618       } else {  // if 22 <= ind <= 33
2619         if (fstar.w[3] || fstar.w[2]
2620             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2621             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2622                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2623           // set the inexact flag
2624           *pfpsf |= INEXACT_EXCEPTION;
2625         }       // else the result is exact
2626       }
2627
2628       res = Cstar.w[0]; // the result is positive
2629     } else if (exp == 0) {
2630       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2631       // res = C (exact)
2632       res = C1.w[0];
2633     } else {
2634       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2635       // res = C * 10^exp (exact) - must fit in 64 bits
2636       res = C1.w[0] * ten2k64[exp];
2637     }
2638   }
2639 }
2640
2641 BID_RETURN (res);
2642 }
2643
2644 /*****************************************************************************
2645  *  BID128_to_uint64_rninta
2646  ****************************************************************************/
2647
2648 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
2649                                           bid128_to_uint64_rninta, x)
2650
2651      UINT64 res;
2652      UINT64 x_sign;
2653      UINT64 x_exp;
2654      int exp;                   // unbiased exponent
2655   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2656      UINT64 tmp64;
2657      BID_UI64DOUBLE tmp1;
2658      unsigned int x_nr_bits;
2659      int q, ind, shift;
2660      UINT128 C1, C;
2661      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
2662      UINT256 P256;
2663
2664   // unpack x
2665 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
2666 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
2667 C1.w[1] = x.w[1] & MASK_COEFF;
2668 C1.w[0] = x.w[0];
2669
2670   // check for NaN or Infinity
2671 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2672     // x is special
2673 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
2674   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
2675     // set invalid flag
2676     *pfpsf |= INVALID_EXCEPTION;
2677     // return Integer Indefinite
2678     res = 0x8000000000000000ull;
2679   } else {      // x is QNaN
2680     // set invalid flag
2681     *pfpsf |= INVALID_EXCEPTION;
2682     // return Integer Indefinite
2683     res = 0x8000000000000000ull;
2684   }
2685   BID_RETURN (res);
2686 } else {        // x is not a NaN, so it must be infinity
2687   if (!x_sign) {        // x is +inf
2688     // set invalid flag
2689     *pfpsf |= INVALID_EXCEPTION;
2690     // return Integer Indefinite
2691     res = 0x8000000000000000ull;
2692   } else {      // x is -inf
2693     // set invalid flag
2694     *pfpsf |= INVALID_EXCEPTION;
2695     // return Integer Indefinite
2696     res = 0x8000000000000000ull;
2697   }
2698   BID_RETURN (res);
2699 }
2700 }
2701   // check for non-canonical values (after the check for special values)
2702 if ((C1.w[1] > 0x0001ed09bead87c0ull)
2703     || (C1.w[1] == 0x0001ed09bead87c0ull
2704         && (C1.w[0] > 0x378d8e63ffffffffull))
2705     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2706   res = 0x0000000000000000ull;
2707   BID_RETURN (res);
2708 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2709   // x is 0
2710   res = 0x0000000000000000ull;
2711   BID_RETURN (res);
2712 } else {        // x is not special and is not zero
2713
2714   // q = nr. of decimal digits in x
2715   //  determine first the nr. of bits in x
2716   if (C1.w[1] == 0) {
2717     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
2718       // split the 64-bit value in two 32-bit halves to avoid rounding errors
2719       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
2720         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
2721         x_nr_bits =
2722           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2723       } else {  // x < 2^32
2724         tmp1.d = (double) (C1.w[0]);    // exact conversion
2725         x_nr_bits =
2726           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2727       }
2728     } else {    // if x < 2^53
2729       tmp1.d = (double) C1.w[0];        // exact conversion
2730       x_nr_bits =
2731         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2732     }
2733   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2734     tmp1.d = (double) C1.w[1];  // exact conversion
2735     x_nr_bits =
2736       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2737   }
2738   q = nr_digits[x_nr_bits - 1].digits;
2739   if (q == 0) {
2740     q = nr_digits[x_nr_bits - 1].digits1;
2741     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2742         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2743             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2744       q++;
2745   }
2746   exp = (x_exp >> 49) - 6176;
2747   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
2748     // set invalid flag
2749     *pfpsf |= INVALID_EXCEPTION;
2750     // return Integer Indefinite
2751     res = 0x8000000000000000ull;
2752     BID_RETURN (res);
2753   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
2754     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2755     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
2756     // the cases that do not fit are identified here; the ones that fit
2757     // fall through and will be handled with other cases further,
2758     // under '1 <= q + exp <= 20'
2759     if (x_sign) {       // if n < 0 and q + exp = 20
2760       // if n <= -1/2 then n cannot be converted to uint64 with RN
2761       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1/2
2762       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x05, 1<=q<=34
2763       // <=> C * 10^(21-q) >= 0x05, 1<=q<=34
2764       if (q == 21) {
2765         // C >= 5
2766         if (C1.w[1] != 0 || C1.w[0] >= 0x05ull) {
2767           // set invalid flag
2768           *pfpsf |= INVALID_EXCEPTION;
2769           // return Integer Indefinite
2770           res = 0x8000000000000000ull;
2771           BID_RETURN (res);
2772         }
2773         // else cases that can be rounded to 64-bit unsigned int fall through
2774         // to '1 <= q + exp <= 20'
2775       } else {
2776         // if 1 <= q <= 20
2777         //   C * 10^(21-q) >= 5 is true because C >= 1 and 10^(21-q) >= 10
2778         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2779         //  C >= 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
2780         // set invalid flag
2781         *pfpsf |= INVALID_EXCEPTION;
2782         // return Integer Indefinite
2783         res = 0x8000000000000000ull;
2784         BID_RETURN (res);
2785       }
2786     } else {    // if n > 0 and q + exp = 20
2787       // if n >= 2^64 - 1/2 then n is too large
2788       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
2789       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
2790       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
2791       // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
2792       if (q == 1) {
2793         // C * 10^20 >= 0x9fffffffffffffffb
2794         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
2795         if (C.w[1] > 0x09 || (C.w[1] == 0x09
2796                               && C.w[0] >= 0xfffffffffffffffbull)) {
2797           // set invalid flag
2798           *pfpsf |= INVALID_EXCEPTION;
2799           // return Integer Indefinite
2800           res = 0x8000000000000000ull;
2801           BID_RETURN (res);
2802         }
2803         // else cases that can be rounded to a 64-bit int fall through
2804         // to '1 <= q + exp <= 20'
2805       } else if (q <= 19) {
2806         // C * 10^(21-q) >= 0x9fffffffffffffffb
2807         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
2808         if (C.w[1] > 0x09 || (C.w[1] == 0x09
2809                               && C.w[0] >= 0xfffffffffffffffbull)) {
2810           // set invalid flag
2811           *pfpsf |= INVALID_EXCEPTION;
2812           // return Integer Indefinite
2813           res = 0x8000000000000000ull;
2814           BID_RETURN (res);
2815         }
2816         // else cases that can be rounded to a 64-bit int fall through
2817         // to '1 <= q + exp <= 20'
2818       } else if (q == 20) {
2819         // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
2820         C.w[0] = C1.w[0] + C1.w[0];
2821         C.w[1] = C1.w[1] + C1.w[1];
2822         if (C.w[0] < C1.w[0])
2823           C.w[1]++;
2824         if (C.w[1] > 0x01 || (C.w[1] == 0x01
2825                               && C.w[0] >= 0xffffffffffffffffull)) {
2826           // set invalid flag
2827           *pfpsf |= INVALID_EXCEPTION;
2828           // return Integer Indefinite
2829           res = 0x8000000000000000ull;
2830           BID_RETURN (res);
2831         }
2832         // else cases that can be rounded to a 64-bit int fall through
2833         // to '1 <= q + exp <= 20'
2834       } else if (q == 21) {
2835         // C >= 0x9fffffffffffffffb
2836         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
2837                                && C1.w[0] >= 0xfffffffffffffffbull)) {
2838           // set invalid flag
2839           *pfpsf |= INVALID_EXCEPTION;
2840           // return Integer Indefinite
2841           res = 0x8000000000000000ull;
2842           BID_RETURN (res);
2843         }
2844         // else cases that can be rounded to a 64-bit int fall through
2845         // to '1 <= q + exp <= 20'
2846       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2847         // C  >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
2848         C.w[1] = 0x09;
2849         C.w[0] = 0xfffffffffffffffbull;
2850         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
2851         if (C1.w[1] > C.w[1]
2852             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2853           // set invalid flag
2854           *pfpsf |= INVALID_EXCEPTION;
2855           // return Integer Indefinite
2856           res = 0x8000000000000000ull;
2857           BID_RETURN (res);
2858         }
2859         // else cases that can be rounded to a 64-bit int fall through
2860         // to '1 <= q + exp <= 20'
2861       }
2862     }
2863   }
2864   // n is not too large to be converted to int64 if -1/2 < n < 2^64 - 1/2
2865   // Note: some of the cases tested for above fall through to this point
2866   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
2867     // return 0
2868     res = 0x0000000000000000ull;
2869     BID_RETURN (res);
2870   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
2871     // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
2872     //   res = 0
2873     // else if x > 0
2874     //   res = +1
2875     // else // if x < 0
2876     //   invalid exc
2877     ind = q - 1;
2878     if (ind <= 18) {    // 0 <= ind <= 18
2879       if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
2880         res = 0x0000000000000000ull;    // return 0
2881       } else if (!x_sign) {     // n > 0
2882         res = 0x00000001;       // return +1
2883       } else {
2884         // set invalid flag
2885         *pfpsf |= INVALID_EXCEPTION;
2886         // return Integer Indefinite
2887         res = 0x8000000000000000ull;
2888         BID_RETURN (res);
2889       }
2890     } else {    // 19 <= ind <= 33
2891       if ((C1.w[1] < midpoint128[ind - 19].w[1])
2892           || ((C1.w[1] == midpoint128[ind - 19].w[1])
2893               && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
2894         res = 0x0000000000000000ull;    // return 0
2895       } else if (!x_sign) {     // n > 0
2896         res = 0x00000001;       // return +1
2897       } else {
2898         // set invalid flag
2899         *pfpsf |= INVALID_EXCEPTION;
2900         // return Integer Indefinite
2901         res = 0x8000000000000000ull;
2902         BID_RETURN (res);
2903       }
2904     }
2905   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
2906     // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
2907     // to nearest to a 64-bit unsigned signed integer
2908     if (x_sign) {       // x <= -1
2909       // set invalid flag
2910       *pfpsf |= INVALID_EXCEPTION;
2911       // return Integer Indefinite
2912       res = 0x8000000000000000ull;
2913       BID_RETURN (res);
2914     }
2915     // 1 <= x < 2^64-1/2 so x can be rounded
2916     // to nearest to a 64-bit unsigned integer
2917     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
2918       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
2919       // chop off ind digits from the lower part of C1
2920       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2921       tmp64 = C1.w[0];
2922       if (ind <= 19) {
2923         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2924       } else {
2925         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2926         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2927       }
2928       if (C1.w[0] < tmp64)
2929         C1.w[1]++;
2930       // calculate C* and f*
2931       // C* is actually floor(C*) in this case
2932       // C* and f* need shifting and masking, as shown by
2933       // shiftright128[] and maskhigh128[]
2934       // 1 <= x <= 33
2935       // kx = 10^(-x) = ten2mk128[ind - 1]
2936       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2937       // the approximation of 10^(-x) was rounded up to 118 bits
2938       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2939       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2940         Cstar.w[1] = P256.w[3];
2941         Cstar.w[0] = P256.w[2];
2942       } else {  // 22 <= ind - 1 <= 33
2943         Cstar.w[1] = 0;
2944         Cstar.w[0] = P256.w[3];
2945       }
2946       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2947       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2948       // if (0 < f* < 10^(-x)) then the result is a midpoint
2949       //   if floor(C*) is even then C* = floor(C*) - logical right
2950       //       shift; C* has p decimal digits, correct by Prop. 1)
2951       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
2952       //       shift; C* has p decimal digits, correct by Pr. 1)
2953       // else
2954       //   C* = floor(C*) (logical right shift; C has p decimal digits,
2955       //       correct by Property 1)
2956       // n = C* * 10^(e+x)
2957
2958       // shift right C* by Ex-128 = shiftright128[ind]
2959       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
2960       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2961         Cstar.w[0] =
2962           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2963         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2964       } else {  // 22 <= ind - 1 <= 33
2965         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
2966       }
2967
2968       // if the result was a midpoint it was rounded away from zero
2969       res = Cstar.w[0]; // the result is positive
2970     } else if (exp == 0) {
2971       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2972       // res = C (exact)
2973       res = C1.w[0];
2974     } else {
2975       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2976       // res = C * 10^exp (exact) - must fit in 64 bits
2977       res = C1.w[0] * ten2k64[exp];
2978     }
2979   }
2980 }
2981
2982 BID_RETURN (res);
2983 }
2984
2985 /*****************************************************************************
2986  *  BID128_to_uint64_xrninta
2987  ****************************************************************************/
2988
2989 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
2990                                           bid128_to_uint64_xrninta, x)
2991
2992      UINT64 res;
2993      UINT64 x_sign;
2994      UINT64 x_exp;
2995      int exp;                   // unbiased exponent
2996   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2997      UINT64 tmp64, tmp64A;
2998      BID_UI64DOUBLE tmp1;
2999      unsigned int x_nr_bits;
3000      int q, ind, shift;
3001      UINT128 C1, C;
3002      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
3003      UINT256 fstar;
3004      UINT256 P256;
3005
3006   // unpack x
3007 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
3008 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
3009 C1.w[1] = x.w[1] & MASK_COEFF;
3010 C1.w[0] = x.w[0];
3011
3012   // check for NaN or Infinity
3013 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
3014     // x is special
3015 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
3016   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
3017     // set invalid flag
3018     *pfpsf |= INVALID_EXCEPTION;
3019     // return Integer Indefinite
3020     res = 0x8000000000000000ull;
3021   } else {      // x is QNaN
3022     // set invalid flag
3023     *pfpsf |= INVALID_EXCEPTION;
3024     // return Integer Indefinite
3025     res = 0x8000000000000000ull;
3026   }
3027   BID_RETURN (res);
3028 } else {        // x is not a NaN, so it must be infinity
3029   if (!x_sign) {        // x is +inf
3030     // set invalid flag
3031     *pfpsf |= INVALID_EXCEPTION;
3032     // return Integer Indefinite
3033     res = 0x8000000000000000ull;
3034   } else {      // x is -inf
3035     // set invalid flag
3036     *pfpsf |= INVALID_EXCEPTION;
3037     // return Integer Indefinite
3038     res = 0x8000000000000000ull;
3039   }
3040   BID_RETURN (res);
3041 }
3042 }
3043   // check for non-canonical values (after the check for special values)
3044 if ((C1.w[1] > 0x0001ed09bead87c0ull)
3045     || (C1.w[1] == 0x0001ed09bead87c0ull
3046         && (C1.w[0] > 0x378d8e63ffffffffull))
3047     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
3048   res = 0x0000000000000000ull;
3049   BID_RETURN (res);
3050 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
3051   // x is 0
3052   res = 0x0000000000000000ull;
3053   BID_RETURN (res);
3054 } else {        // x is not special and is not zero
3055
3056   // q = nr. of decimal digits in x
3057   //  determine first the nr. of bits in x
3058   if (C1.w[1] == 0) {
3059     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
3060       // split the 64-bit value in two 32-bit halves to avoid rounding errors
3061       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
3062         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
3063         x_nr_bits =
3064           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3065       } else {  // x < 2^32
3066         tmp1.d = (double) (C1.w[0]);    // exact conversion
3067         x_nr_bits =
3068           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3069       }
3070     } else {    // if x < 2^53
3071       tmp1.d = (double) C1.w[0];        // exact conversion
3072       x_nr_bits =
3073         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3074     }
3075   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
3076     tmp1.d = (double) C1.w[1];  // exact conversion
3077     x_nr_bits =
3078       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3079   }
3080   q = nr_digits[x_nr_bits - 1].digits;
3081   if (q == 0) {
3082     q = nr_digits[x_nr_bits - 1].digits1;
3083     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
3084         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
3085             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
3086       q++;
3087   }
3088   exp = (x_exp >> 49) - 6176;
3089
3090   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
3091     // set invalid flag
3092     *pfpsf |= INVALID_EXCEPTION;
3093     // return Integer Indefinite
3094     res = 0x8000000000000000ull;
3095     BID_RETURN (res);
3096   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
3097     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
3098     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
3099     // the cases that do not fit are identified here; the ones that fit
3100     // fall through and will be handled with other cases further,
3101     // under '1 <= q + exp <= 20'
3102     if (x_sign) {       // if n < 0 and q + exp = 20
3103       // if n <= -1/2 then n cannot be converted to uint64 with RN
3104       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1/2
3105       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x05, 1<=q<=34
3106       // <=> C * 10^(21-q) >= 0x05, 1<=q<=34
3107       if (q == 21) {
3108         // C >= 5
3109         if (C1.w[1] != 0 || C1.w[0] >= 0x05ull) {
3110           // set invalid flag
3111           *pfpsf |= INVALID_EXCEPTION;
3112           // return Integer Indefinite
3113           res = 0x8000000000000000ull;
3114           BID_RETURN (res);
3115         }
3116         // else cases that can be rounded to 64-bit unsigned int fall through
3117         // to '1 <= q + exp <= 20'
3118       } else {
3119         // if 1 <= q <= 20
3120         //   C * 10^(21-q) >= 5 is true because C >= 1 and 10^(21-q) >= 10
3121         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
3122         //  C >= 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
3123         // set invalid flag
3124         *pfpsf |= INVALID_EXCEPTION;
3125         // return Integer Indefinite
3126         res = 0x8000000000000000ull;
3127         BID_RETURN (res);
3128       }
3129     } else {    // if n > 0 and q + exp = 20
3130       // if n >= 2^64 - 1/2 then n is too large
3131       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
3132       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
3133       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
3134       // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
3135       if (q == 1) {
3136         // C * 10^20 >= 0x9fffffffffffffffb
3137         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
3138         if (C.w[1] > 0x09 || (C.w[1] == 0x09
3139                               && C.w[0] >= 0xfffffffffffffffbull)) {
3140           // set invalid flag
3141           *pfpsf |= INVALID_EXCEPTION;
3142           // return Integer Indefinite
3143           res = 0x8000000000000000ull;
3144           BID_RETURN (res);
3145         }
3146         // else cases that can be rounded to a 64-bit int fall through
3147         // to '1 <= q + exp <= 20'
3148       } else if (q <= 19) {
3149         // C * 10^(21-q) >= 0x9fffffffffffffffb
3150         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
3151         if (C.w[1] > 0x09 || (C.w[1] == 0x09
3152                               && C.w[0] >= 0xfffffffffffffffbull)) {
3153           // set invalid flag
3154           *pfpsf |= INVALID_EXCEPTION;
3155           // return Integer Indefinite
3156           res = 0x8000000000000000ull;
3157           BID_RETURN (res);
3158         }
3159         // else cases that can be rounded to a 64-bit int fall through
3160         // to '1 <= q + exp <= 20'
3161       } else if (q == 20) {
3162         // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
3163         C.w[0] = C1.w[0] + C1.w[0];
3164         C.w[1] = C1.w[1] + C1.w[1];
3165         if (C.w[0] < C1.w[0])
3166           C.w[1]++;
3167         if (C.w[1] > 0x01 || (C.w[1] == 0x01
3168                               && C.w[0] >= 0xffffffffffffffffull)) {
3169           // set invalid flag
3170           *pfpsf |= INVALID_EXCEPTION;
3171           // return Integer Indefinite
3172           res = 0x8000000000000000ull;
3173           BID_RETURN (res);
3174         }
3175         // else cases that can be rounded to a 64-bit int fall through
3176         // to '1 <= q + exp <= 20'
3177       } else if (q == 21) {
3178         // C >= 0x9fffffffffffffffb
3179         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
3180                                && C1.w[0] >= 0xfffffffffffffffbull)) {
3181           // set invalid flag
3182           *pfpsf |= INVALID_EXCEPTION;
3183           // return Integer Indefinite
3184           res = 0x8000000000000000ull;
3185           BID_RETURN (res);
3186         }
3187         // else cases that can be rounded to a 64-bit int fall through
3188         // to '1 <= q + exp <= 20'
3189       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
3190         // C  >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
3191         C.w[1] = 0x09;
3192         C.w[0] = 0xfffffffffffffffbull;
3193         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
3194         if (C1.w[1] > C.w[1]
3195             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3196           // set invalid flag
3197           *pfpsf |= INVALID_EXCEPTION;
3198           // return Integer Indefinite
3199           res = 0x8000000000000000ull;
3200           BID_RETURN (res);
3201         }
3202         // else cases that can be rounded to a 64-bit int fall through
3203         // to '1 <= q + exp <= 20'
3204       }
3205     }
3206   }
3207   // n is not too large to be converted to int64 if -1/2 < n < 2^64 - 1/2
3208   // Note: some of the cases tested for above fall through to this point
3209   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
3210     // set inexact flag
3211     *pfpsf |= INEXACT_EXCEPTION;
3212     // return 0
3213     res = 0x0000000000000000ull;
3214     BID_RETURN (res);
3215   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
3216     // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3217     //   res = 0
3218     // else if x > 0
3219     //   res = +1
3220     // else // if x < 0
3221     //   invalid exc
3222     ind = q - 1;
3223     if (ind <= 18) {    // 0 <= ind <= 18
3224       if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
3225         res = 0x0000000000000000ull;    // return 0
3226       } else if (!x_sign) {     // n > 0
3227         res = 0x00000001;       // return +1
3228       } else {
3229         res = 0x8000000000000000ull;
3230         // set invalid flag
3231         *pfpsf |= INVALID_EXCEPTION;
3232         // return Integer Indefinite
3233         res = 0x8000000000000000ull;
3234         BID_RETURN (res);
3235       }
3236     } else {    // 19 <= ind <= 33
3237       if ((C1.w[1] < midpoint128[ind - 19].w[1])
3238           || ((C1.w[1] == midpoint128[ind - 19].w[1])
3239               && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
3240         res = 0x0000000000000000ull;    // return 0
3241       } else if (!x_sign) {     // n > 0
3242         res = 0x00000001;       // return +1
3243       } else {
3244         res = 0x8000000000000000ull;
3245         *pfpsf |= INVALID_EXCEPTION;
3246         // return Integer Indefinite
3247         res = 0x8000000000000000ull;
3248         BID_RETURN (res);
3249       }
3250     }
3251     // set inexact flag
3252     *pfpsf |= INEXACT_EXCEPTION;
3253   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
3254     // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
3255     // to nearest to a 64-bit unsigned signed integer
3256     if (x_sign) {       // x <= -1
3257       // set invalid flag
3258       *pfpsf |= INVALID_EXCEPTION;
3259       // return Integer Indefinite
3260       res = 0x8000000000000000ull;
3261       BID_RETURN (res);
3262     }
3263     // 1 <= x < 2^64-1/2 so x can be rounded
3264     // to nearest to a 64-bit unsigned integer
3265     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
3266       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
3267       // chop off ind digits from the lower part of C1
3268       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3269       tmp64 = C1.w[0];
3270       if (ind <= 19) {
3271         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
3272       } else {
3273         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
3274         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
3275       }
3276       if (C1.w[0] < tmp64)
3277         C1.w[1]++;
3278       // calculate C* and f*
3279       // C* is actually floor(C*) in this case
3280       // C* and f* need shifting and masking, as shown by
3281       // shiftright128[] and maskhigh128[]
3282       // 1 <= x <= 33
3283       // kx = 10^(-x) = ten2mk128[ind - 1]
3284       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3285       // the approximation of 10^(-x) was rounded up to 118 bits
3286       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
3287       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
3288         Cstar.w[1] = P256.w[3];
3289         Cstar.w[0] = P256.w[2];
3290         fstar.w[3] = 0;
3291         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
3292         fstar.w[1] = P256.w[1];
3293         fstar.w[0] = P256.w[0];
3294       } else {  // 22 <= ind - 1 <= 33
3295         Cstar.w[1] = 0;
3296         Cstar.w[0] = P256.w[3];
3297         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
3298         fstar.w[2] = P256.w[2];
3299         fstar.w[1] = P256.w[1];
3300         fstar.w[0] = P256.w[0];
3301       }
3302       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3303       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3304       // if (0 < f* < 10^(-x)) then the result is a midpoint
3305       //   if floor(C*) is even then C* = floor(C*) - logical right
3306       //       shift; C* has p decimal digits, correct by Prop. 1)
3307       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
3308       //       shift; C* has p decimal digits, correct by Pr. 1)
3309       // else
3310       //   C* = floor(C*) (logical right shift; C has p decimal digits,
3311       //       correct by Property 1)
3312       // n = C* * 10^(e+x)
3313
3314       // shift right C* by Ex-128 = shiftright128[ind]
3315       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
3316       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
3317         Cstar.w[0] =
3318           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
3319         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3320       } else {  // 22 <= ind - 1 <= 33
3321         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
3322       }
3323       // determine inexactness of the rounding of C*
3324       // if (0 < f* - 1/2 < 10^(-x)) then
3325       //   the result is exact
3326       // else // if (f* - 1/2 > T*) then
3327       //   the result is inexact
3328       if (ind - 1 <= 2) {
3329         if (fstar.w[1] > 0x8000000000000000ull ||
3330             (fstar.w[1] == 0x8000000000000000ull
3331              && fstar.w[0] > 0x0ull)) {
3332           // f* > 1/2 and the result may be exact
3333           tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
3334           if (tmp64 > ten2mk128trunc[ind - 1].w[1]
3335               || (tmp64 == ten2mk128trunc[ind - 1].w[1]
3336                   && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
3337             // set the inexact flag
3338             *pfpsf |= INEXACT_EXCEPTION;
3339           }     // else the result is exact
3340         } else {        // the result is inexact; f2* <= 1/2
3341           // set the inexact flag
3342           *pfpsf |= INEXACT_EXCEPTION;
3343         }
3344       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
3345         if (fstar.w[3] > 0x0 ||
3346             (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
3347             (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
3348              (fstar.w[1] || fstar.w[0]))) {
3349           // f2* > 1/2 and the result may be exact
3350           // Calculate f2* - 1/2
3351           tmp64 = fstar.w[2] - onehalf128[ind - 1];
3352           tmp64A = fstar.w[3];
3353           if (tmp64 > fstar.w[2])
3354             tmp64A--;
3355           if (tmp64A || tmp64
3356               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
3357               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
3358                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
3359             // set the inexact flag
3360             *pfpsf |= INEXACT_EXCEPTION;
3361           }     // else the result is exact
3362         } else {        // the result is inexact; f2* <= 1/2
3363           // set the inexact flag
3364           *pfpsf |= INEXACT_EXCEPTION;
3365         }
3366       } else {  // if 22 <= ind <= 33
3367         if (fstar.w[3] > onehalf128[ind - 1] ||
3368             (fstar.w[3] == onehalf128[ind - 1] &&
3369              (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
3370           // f2* > 1/2 and the result may be exact
3371           // Calculate f2* - 1/2
3372           tmp64 = fstar.w[3] - onehalf128[ind - 1];
3373           if (tmp64 || fstar.w[2]
3374               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
3375               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
3376                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
3377             // set the inexact flag
3378             *pfpsf |= INEXACT_EXCEPTION;
3379           }     // else the result is exact
3380         } else {        // the result is inexact; f2* <= 1/2
3381           // set the inexact flag
3382           *pfpsf |= INEXACT_EXCEPTION;
3383         }
3384       }
3385
3386       // if the result was a midpoint it was rounded away from zero
3387       res = Cstar.w[0]; // the result is positive
3388     } else if (exp == 0) {
3389       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
3390       // res = C (exact)
3391       res = C1.w[0];
3392     } else {
3393       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
3394       // res = C * 10^exp (exact) - must fit in 64 bits
3395       res = C1.w[0] * ten2k64[exp];
3396     }
3397   }
3398 }
3399
3400 BID_RETURN (res);
3401 }