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