libgcc/config/libbid/bid64_to_uint64.c

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