libgcc/config/libbid/bid64_to_int64.c

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