libgcc/config/libbid/bid128_to_int32.c

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