libgcc/config/libbid/bid128_to_uint32.c

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