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