libgcc/config/libbid/bid64_to_int32.c

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