libgcc/config/libbid/bid128_to_int64.c

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