libgcc/config/libbid/bid128_to_int64.c

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