libgcc/config/libbid/bid64_to_int64.c

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