libgcc/config/libbid/bid64_to_int32.c

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