libgcc/config/libbid/bid128_to_uint64.c

   1 /* Copyright (C) 2007  Free Software Foundation, Inc.
   2
   3 This file is part of GCC.
   4
   5 GCC is free software; you can redistribute it and/or modify it under
   6 the terms of the GNU General Public License as published by the Free
   7 Software Foundation; either version 2, or (at your option) any later
   8 version.
   9
  10 In addition to the permissions in the GNU General Public License, the
  11 Free Software Foundation gives you unlimited permission to link the
  12 compiled version of this file into combinations with other programs,
  13 and to distribute those combinations without any restriction coming
  14 from the use of this file.  (The General Public License restrictions
  15 do apply in other respects; for example, they cover modification of
  16 the file, and distribution when not linked into a combine
  17 executable.)
  18
  19 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  20 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  21 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  22 for more details.
  23
  24 You should have received a copy of the GNU General Public License
  25 along with GCC; see the file COPYING.  If not, write to the Free
  26 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
  27 02110-1301, USA.  */
  28
  29 #include "bid_internal.h"
  30
  31 /*****************************************************************************
  32  *  BID128_to_uint64_rnint
  33  ****************************************************************************/
  34
  35 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
  36                                           bid128_to_uint64_rnint, x)
  37
  38      UINT64 res;
  39      UINT64 x_sign;
  40      UINT64 x_exp;
  41      int exp;                   // unbiased exponent
  42   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  43      UINT64 tmp64;
  44      BID_UI64DOUBLE tmp1;
  45      unsigned int x_nr_bits;
  46      int q, ind, shift;
  47      UINT128 C1, C;
  48      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  49      UINT256 fstar;
  50      UINT256 P256;
  51
  52   // unpack x
  53 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  54 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
  55 C1.w[1] = x.w[1] & MASK_COEFF;
  56 C1.w[0] = x.w[0];
  57
  58   // check for NaN or Infinity
  59 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  60     // x is special
  61 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
  62   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
  63     // set invalid flag
  64     *pfpsf |= INVALID_EXCEPTION;
  65     // return Integer Indefinite
  66     res = 0x8000000000000000ull;
  67   } else {      // x is QNaN
  68     // set invalid flag
  69     *pfpsf |= INVALID_EXCEPTION;
  70     // return Integer Indefinite
  71     res = 0x8000000000000000ull;
  72   }
  73   BID_RETURN (res);
  74 } else {        // x is not a NaN, so it must be infinity
  75   if (!x_sign) {        // x is +inf
  76     // set invalid flag
  77     *pfpsf |= INVALID_EXCEPTION;
  78     // return Integer Indefinite
  79     res = 0x8000000000000000ull;
  80   } else {      // x is -inf
  81     // set invalid flag
  82     *pfpsf |= INVALID_EXCEPTION;
  83     // return Integer Indefinite
  84     res = 0x8000000000000000ull;
  85   }
  86   BID_RETURN (res);
  87 }
  88 }
  89   // check for non-canonical values (after the check for special values)
  90 if ((C1.w[1] > 0x0001ed09bead87c0ull)
  91     || (C1.w[1] == 0x0001ed09bead87c0ull
  92         && (C1.w[0] > 0x378d8e63ffffffffull))
  93     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  94   res = 0x0000000000000000ull;
  95   BID_RETURN (res);
  96 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  97   // x is 0
  98   res = 0x0000000000000000ull;
  99   BID_RETURN (res);
 100 } else {        // x is not special and is not zero
 101
 102   // q = nr. of decimal digits in x
 103   //  determine first the nr. of bits in x
 104   if (C1.w[1] == 0) {
 105     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
 106       // split the 64-bit value in two 32-bit halves to avoid rounding errors
 107       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
 108         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
 109         x_nr_bits =
 110           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 111       } else {  // x < 2^32
 112         tmp1.d = (double) (C1.w[0]);    // exact conversion
 113         x_nr_bits =
 114           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 115       }
 116     } else {    // if x < 2^53
 117       tmp1.d = (double) C1.w[0];        // exact conversion
 118       x_nr_bits =
 119         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 120     }
 121   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
 122     tmp1.d = (double) C1.w[1];  // exact conversion
 123     x_nr_bits =
 124       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 125   }
 126   q = nr_digits[x_nr_bits - 1].digits;
 127   if (q == 0) {
 128     q = nr_digits[x_nr_bits - 1].digits1;
 129     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
 130         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
 131             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
 132       q++;
 133   }
 134   exp = (x_exp >> 49) - 6176;
 135
 136   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
 137     // set invalid flag
 138     *pfpsf |= INVALID_EXCEPTION;
 139     // return Integer Indefinite
 140     res = 0x8000000000000000ull;
 141     BID_RETURN (res);
 142   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
 143     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
 144     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
 145     // the cases that do not fit are identified here; the ones that fit
 146     // fall through and will be handled with other cases further,
 147     // under '1 <= q + exp <= 20'
 148     if (x_sign) {       // if n < 0 and q + exp = 20
 149       // if n < -1/2 then n cannot be converted to uint64 with RN
 150       // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 1/2
 151       // <=> 0.c(0)c(1)...c(q-1) * 10^21 > 0x05, 1<=q<=34
 152       // <=> C * 10^(21-q) > 0x05, 1<=q<=34
 153       if (q == 21) {
 154         // C > 5
 155         if (C1.w[1] != 0 || C1.w[0] > 0x05ull) {
 156           // set invalid flag
 157           *pfpsf |= INVALID_EXCEPTION;
 158           // return Integer Indefinite
 159           res = 0x8000000000000000ull;
 160           BID_RETURN (res);
 161         }
 162         // else cases that can be rounded to 64-bit unsigned int fall through
 163         // to '1 <= q + exp <= 20'
 164       } else {
 165         // if 1 <= q <= 20
 166         //   C * 10^(21-q) > 5 is true because C >= 1 and 10^(21-q) >= 10
 167         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
 168         //   C > 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
 169         // set invalid flag
 170         *pfpsf |= INVALID_EXCEPTION;
 171         // return Integer Indefinite
 172         res = 0x8000000000000000ull;
 173         BID_RETURN (res);
 174       }
 175     } else {    // if n > 0 and q + exp = 20
 176       // if n >= 2^64 - 1/2 then n is too large
 177       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
 178       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
 179       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
 180       // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
 181       if (q == 1) {
 182         // C * 10^20 >= 0x9fffffffffffffffb
 183         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
 184         if (C.w[1] > 0x09 || (C.w[1] == 0x09
 185                               && C.w[0] >= 0xfffffffffffffffbull)) {
 186           // set invalid flag
 187           *pfpsf |= INVALID_EXCEPTION;
 188           // return Integer Indefinite
 189           res = 0x8000000000000000ull;
 190           BID_RETURN (res);
 191         }
 192         // else cases that can be rounded to a 64-bit int fall through
 193         // to '1 <= q + exp <= 20'
 194       } else if (q <= 19) {
 195         // C * 10^(21-q) >= 0x9fffffffffffffffb
 196         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
 197         if (C.w[1] > 0x09 || (C.w[1] == 0x09
 198                               && C.w[0] >= 0xfffffffffffffffbull)) {
 199           // set invalid flag
 200           *pfpsf |= INVALID_EXCEPTION;
 201           // return Integer Indefinite
 202           res = 0x8000000000000000ull;
 203           BID_RETURN (res);
 204         }
 205         // else cases that can be rounded to a 64-bit int fall through
 206         // to '1 <= q + exp <= 20'
 207       } else if (q == 20) {
 208         // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
 209         C.w[0] = C1.w[0] + C1.w[0];
 210         C.w[1] = C1.w[1] + C1.w[1];
 211         if (C.w[0] < C1.w[0])
 212           C.w[1]++;
 213         if (C.w[1] > 0x01 || (C.w[1] == 0x01
 214                               && C.w[0] >= 0xffffffffffffffffull)) {
 215           // set invalid flag
 216           *pfpsf |= INVALID_EXCEPTION;
 217           // return Integer Indefinite
 218           res = 0x8000000000000000ull;
 219           BID_RETURN (res);
 220         }
 221         // else cases that can be rounded to a 64-bit int fall through
 222         // to '1 <= q + exp <= 20'
 223       } else if (q == 21) {
 224         // C >= 0x9fffffffffffffffb
 225         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
 226                                && C1.w[0] >= 0xfffffffffffffffbull)) {
 227           // set invalid flag
 228           *pfpsf |= INVALID_EXCEPTION;
 229           // return Integer Indefinite
 230           res = 0x8000000000000000ull;
 231           BID_RETURN (res);
 232         }
 233         // else cases that can be rounded to a 64-bit int fall through
 234         // to '1 <= q + exp <= 20'
 235       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
 236         // C  >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
 237         C.w[1] = 0x09;
 238         C.w[0] = 0xfffffffffffffffbull;
 239         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
 240         if (C1.w[1] > C.w[1]
 241             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
 242           // set invalid flag
 243           *pfpsf |= INVALID_EXCEPTION;
 244           // return Integer Indefinite
 245           res = 0x8000000000000000ull;
 246           BID_RETURN (res);
 247         }
 248         // else cases that can be rounded to a 64-bit int fall through
 249         // to '1 <= q + exp <= 20'
 250       }
 251     }
 252   }
 253   // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
 254   // Note: some of the cases tested for above fall through to this point
 255   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
 256     // return 0
 257     res = 0x0000000000000000ull;
 258     BID_RETURN (res);
 259   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
 260     // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
 261     //   res = 0
 262     // else if x > 0
 263     //   res = +1
 264     // else // if x < 0
 265     //   invalid exc
 266     ind = q - 1;
 267     if (ind <= 18) {    // 0 <= ind <= 18
 268       if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
 269         res = 0x0000000000000000ull;    // return 0
 270       } else if (!x_sign) {     // n > 0
 271         res = 0x00000001;       // return +1
 272       } else {
 273         res = 0x8000000000000000ull;
 274         *pfpsf |= INVALID_EXCEPTION;
 275       }
 276     } else {    // 19 <= ind <= 33
 277       if ((C1.w[1] < midpoint128[ind - 19].w[1])
 278           || ((C1.w[1] == midpoint128[ind - 19].w[1])
 279               && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
 280         res = 0x0000000000000000ull;    // return 0
 281       } else if (!x_sign) {     // n > 0
 282         res = 0x00000001;       // return +1
 283       } else {
 284         res = 0x8000000000000000ull;
 285         *pfpsf |= INVALID_EXCEPTION;
 286       }
 287     }
 288   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
 289     // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
 290     // to nearest to a 64-bit unsigned signed integer
 291     if (x_sign) {       // x <= -1
 292       // set invalid flag
 293       *pfpsf |= INVALID_EXCEPTION;
 294       // return Integer Indefinite
 295       res = 0x8000000000000000ull;
 296       BID_RETURN (res);
 297     }
 298     // 1 <= x < 2^64-1/2 so x can be rounded
 299     // to nearest to a 64-bit unsigned integer
 300     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
 301       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
 302       // chop off ind digits from the lower part of C1
 303       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
 304       tmp64 = C1.w[0];
 305       if (ind <= 19) {
 306         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
 307       } else {
 308         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
 309         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
 310       }
 311       if (C1.w[0] < tmp64)
 312         C1.w[1]++;
 313       // calculate C* and f*
 314       // C* is actually floor(C*) in this case
 315       // C* and f* need shifting and masking, as shown by
 316       // shiftright128[] and maskhigh128[]
 317       // 1 <= x <= 33
 318       // kx = 10^(-x) = ten2mk128[ind - 1]
 319       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 320       // the approximation of 10^(-x) was rounded up to 118 bits
 321       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
 322       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
 323         Cstar.w[1] = P256.w[3];
 324         Cstar.w[0] = P256.w[2];
 325         fstar.w[3] = 0;
 326         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
 327         fstar.w[1] = P256.w[1];
 328         fstar.w[0] = P256.w[0];
 329       } else {  // 22 <= ind - 1 <= 33
 330         Cstar.w[1] = 0;
 331         Cstar.w[0] = P256.w[3];
 332         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
 333         fstar.w[2] = P256.w[2];
 334         fstar.w[1] = P256.w[1];
 335         fstar.w[0] = P256.w[0];
 336       }
 337       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
 338       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
 339       // if (0 < f* < 10^(-x)) then the result is a midpoint
 340       //   if floor(C*) is even then C* = floor(C*) - logical right
 341       //       shift; C* has p decimal digits, correct by Prop. 1)
 342       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
 343       //       shift; C* has p decimal digits, correct by Pr. 1)
 344       // else
 345       //   C* = floor(C*) (logical right shift; C has p decimal digits,
 346       //       correct by Property 1)
 347       // n = C* * 10^(e+x)
 348
 349       // shift right C* by Ex-128 = shiftright128[ind]
 350       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
 351       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
 352         Cstar.w[0] =
 353           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
 354         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
 355       } else {  // 22 <= ind - 1 <= 33
 356         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
 357       }
 358       // if the result was a midpoint it was rounded away from zero, so
 359       // it will need a correction
 360       // check for midpoints
 361       if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
 362           && (fstar.w[1] || fstar.w[0])
 363           && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
 364               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
 365                   && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
 366         // the result is a midpoint; round to nearest
 367         if (Cstar.w[0] & 0x01) {        // Cstar.w[0] is odd; MP in [EVEN, ODD]
 368           // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
 369           Cstar.w[0]--; // Cstar.w[0] is now even
 370         }       // else MP in [ODD, EVEN]
 371       }
 372       res = Cstar.w[0]; // the result is positive
 373     } else if (exp == 0) {
 374       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
 375       // res = C (exact)
 376       res = C1.w[0];
 377     } else {
 378       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
 379       // res = C * 10^exp (exact) - must fit in 64 bits
 380       res = C1.w[0] * ten2k64[exp];
 381     }
 382   }
 383 }
 384
 385 BID_RETURN (res);
 386 }
 387
 388 /*****************************************************************************
 389  *  BID128_to_uint64_xrnint
 390  ****************************************************************************/
 391
 392 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
 393                                           bid128_to_uint64_xrnint, x)
 394
 395      UINT64 res;
 396      UINT64 x_sign;
 397      UINT64 x_exp;
 398      int exp;                   // unbiased exponent
 399   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
 400      UINT64 tmp64, tmp64A;
 401      BID_UI64DOUBLE tmp1;
 402      unsigned int x_nr_bits;
 403      int q, ind, shift;
 404      UINT128 C1, C;
 405      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
 406      UINT256 fstar;
 407      UINT256 P256;
 408
 409   // unpack x
 410 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
 411 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
 412 C1.w[1] = x.w[1] & MASK_COEFF;
 413 C1.w[0] = x.w[0];
 414
 415   // check for NaN or Infinity
 416 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
 417     // x is special
 418 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
 419   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
 420     // set invalid flag
 421     *pfpsf |= INVALID_EXCEPTION;
 422     // return Integer Indefinite
 423     res = 0x8000000000000000ull;
 424   } else {      // x is QNaN
 425     // set invalid flag
 426     *pfpsf |= INVALID_EXCEPTION;
 427     // return Integer Indefinite
 428     res = 0x8000000000000000ull;
 429   }
 430   BID_RETURN (res);
 431 } else {        // x is not a NaN, so it must be infinity
 432   if (!x_sign) {        // x is +inf
 433     // set invalid flag
 434     *pfpsf |= INVALID_EXCEPTION;
 435     // return Integer Indefinite
 436     res = 0x8000000000000000ull;
 437   } else {      // x is -inf
 438     // set invalid flag
 439     *pfpsf |= INVALID_EXCEPTION;
 440     // return Integer Indefinite
 441     res = 0x8000000000000000ull;
 442   }
 443   BID_RETURN (res);
 444 }
 445 }
 446   // check for non-canonical values (after the check for special values)
 447 if ((C1.w[1] > 0x0001ed09bead87c0ull)
 448     || (C1.w[1] == 0x0001ed09bead87c0ull
 449         && (C1.w[0] > 0x378d8e63ffffffffull))
 450     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
 451   res = 0x0000000000000000ull;
 452   BID_RETURN (res);
 453 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
 454   // x is 0
 455   res = 0x0000000000000000ull;
 456   BID_RETURN (res);
 457 } else {        // x is not special and is not zero
 458
 459   // q = nr. of decimal digits in x
 460   //  determine first the nr. of bits in x
 461   if (C1.w[1] == 0) {
 462     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
 463       // split the 64-bit value in two 32-bit halves to avoid rounding errors
 464       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
 465         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
 466         x_nr_bits =
 467           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 468       } else {  // x < 2^32
 469         tmp1.d = (double) (C1.w[0]);    // exact conversion
 470         x_nr_bits =
 471           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 472       }
 473     } else {    // if x < 2^53
 474       tmp1.d = (double) C1.w[0];        // exact conversion
 475       x_nr_bits =
 476         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 477     }
 478   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
 479     tmp1.d = (double) C1.w[1];  // exact conversion
 480     x_nr_bits =
 481       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 482   }
 483   q = nr_digits[x_nr_bits - 1].digits;
 484   if (q == 0) {
 485     q = nr_digits[x_nr_bits - 1].digits1;
 486     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
 487         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
 488             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
 489       q++;
 490   }
 491   exp = (x_exp >> 49) - 6176;
 492
 493   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
 494     // set invalid flag
 495     *pfpsf |= INVALID_EXCEPTION;
 496     // return Integer Indefinite
 497     res = 0x8000000000000000ull;
 498     BID_RETURN (res);
 499   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
 500     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
 501     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
 502     // the cases that do not fit are identified here; the ones that fit
 503     // fall through and will be handled with other cases further,
 504     // under '1 <= q + exp <= 20'
 505     if (x_sign) {       // if n < 0 and q + exp = 20
 506       // if n < -1/2 then n cannot be converted to uint64 with RN
 507       // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 1/2
 508       // <=> 0.c(0)c(1)...c(q-1) * 10^21 > 0x05, 1<=q<=34
 509       // <=> C * 10^(21-q) > 0x05, 1<=q<=34
 510       if (q == 21) {
 511         // C > 5
 512         if (C1.w[1] != 0 || C1.w[0] > 0x05ull) {
 513           // set invalid flag
 514           *pfpsf |= INVALID_EXCEPTION;
 515           // return Integer Indefinite
 516           res = 0x8000000000000000ull;
 517           BID_RETURN (res);
 518         }
 519         // else cases that can be rounded to 64-bit unsigned int fall through
 520         // to '1 <= q + exp <= 20'
 521       } else {
 522         // if 1 <= q <= 20
 523         //   C * 10^(21-q) > 5 is true because C >= 1 and 10^(21-q) >= 10
 524         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
 525         //   C > 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
 526         // set invalid flag
 527         *pfpsf |= INVALID_EXCEPTION;
 528         // return Integer Indefinite
 529         res = 0x8000000000000000ull;
 530         BID_RETURN (res);
 531       }
 532     } else {    // if n > 0 and q + exp = 20
 533       // if n >= 2^64 - 1/2 then n is too large
 534       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
 535       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
 536       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
 537       // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
 538       if (q == 1) {
 539         // C * 10^20 >= 0x9fffffffffffffffb
 540         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
 541         if (C.w[1] > 0x09 || (C.w[1] == 0x09
 542                               && C.w[0] >= 0xfffffffffffffffbull)) {
 543           // set invalid flag
 544           *pfpsf |= INVALID_EXCEPTION;
 545           // return Integer Indefinite
 546           res = 0x8000000000000000ull;
 547           BID_RETURN (res);
 548         }
 549         // else cases that can be rounded to a 64-bit int fall through
 550         // to '1 <= q + exp <= 20'
 551       } else if (q <= 19) {
 552         // C * 10^(21-q) >= 0x9fffffffffffffffb
 553         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
 554         if (C.w[1] > 0x09 || (C.w[1] == 0x09
 555                               && C.w[0] >= 0xfffffffffffffffbull)) {
 556           // set invalid flag
 557           *pfpsf |= INVALID_EXCEPTION;
 558           // return Integer Indefinite
 559           res = 0x8000000000000000ull;
 560           BID_RETURN (res);
 561         }
 562         // else cases that can be rounded to a 64-bit int fall through
 563         // to '1 <= q + exp <= 20'
 564       } else if (q == 20) {
 565         // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
 566         C.w[0] = C1.w[0] + C1.w[0];
 567         C.w[1] = C1.w[1] + C1.w[1];
 568         if (C.w[0] < C1.w[0])
 569           C.w[1]++;
 570         if (C.w[1] > 0x01 || (C.w[1] == 0x01
 571                               && C.w[0] >= 0xffffffffffffffffull)) {
 572           // set invalid flag
 573           *pfpsf |= INVALID_EXCEPTION;
 574           // return Integer Indefinite
 575           res = 0x8000000000000000ull;
 576           BID_RETURN (res);
 577         }
 578         // else cases that can be rounded to a 64-bit int fall through
 579         // to '1 <= q + exp <= 20'
 580       } else if (q == 21) {
 581         // C >= 0x9fffffffffffffffb
 582         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
 583                                && C1.w[0] >= 0xfffffffffffffffbull)) {
 584           // set invalid flag
 585           *pfpsf |= INVALID_EXCEPTION;
 586           // return Integer Indefinite
 587           res = 0x8000000000000000ull;
 588           BID_RETURN (res);
 589         }
 590         // else cases that can be rounded to a 64-bit int fall through
 591         // to '1 <= q + exp <= 20'
 592       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
 593         // C  >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
 594         C.w[1] = 0x09;
 595         C.w[0] = 0xfffffffffffffffbull;
 596         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
 597         if (C1.w[1] > C.w[1]
 598             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
 599           // set invalid flag
 600           *pfpsf |= INVALID_EXCEPTION;
 601           // return Integer Indefinite
 602           res = 0x8000000000000000ull;
 603           BID_RETURN (res);
 604         }
 605         // else cases that can be rounded to a 64-bit int fall through
 606         // to '1 <= q + exp <= 20'
 607       }
 608     }
 609   }
 610   // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
 611   // Note: some of the cases tested for above fall through to this point
 612   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
 613     // set inexact flag
 614     *pfpsf |= INEXACT_EXCEPTION;
 615     // return 0
 616     res = 0x0000000000000000ull;
 617     BID_RETURN (res);
 618   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
 619     // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
 620     //   res = 0
 621     // else if x > 0
 622     //   res = +1
 623     // else // if x < 0
 624     //   invalid exc
 625     ind = q - 1;
 626     if (ind <= 18) {    // 0 <= ind <= 18
 627       if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
 628         res = 0x0000000000000000ull;    // return 0
 629       } else if (!x_sign) {     // n > 0
 630         res = 0x00000001;       // return +1
 631       } else {
 632         res = 0x8000000000000000ull;
 633         *pfpsf |= INVALID_EXCEPTION;
 634         BID_RETURN (res);
 635       }
 636     } else {    // 19 <= ind <= 33
 637       if ((C1.w[1] < midpoint128[ind - 19].w[1])
 638           || ((C1.w[1] == midpoint128[ind - 19].w[1])
 639               && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
 640         res = 0x0000000000000000ull;    // return 0
 641       } else if (!x_sign) {     // n > 0
 642         res = 0x00000001;       // return +1
 643       } else {
 644         res = 0x8000000000000000ull;
 645         *pfpsf |= INVALID_EXCEPTION;
 646         BID_RETURN (res);
 647       }
 648     }
 649     // set inexact flag
 650     *pfpsf |= INEXACT_EXCEPTION;
 651   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
 652     // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
 653     // to nearest to a 64-bit unsigned signed integer
 654     if (x_sign) {       // x <= -1
 655       // set invalid flag
 656       *pfpsf |= INVALID_EXCEPTION;
 657       // return Integer Indefinite
 658       res = 0x8000000000000000ull;
 659       BID_RETURN (res);
 660     }
 661     // 1 <= x < 2^64-1/2 so x can be rounded
 662     // to nearest to a 64-bit unsigned integer
 663     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
 664       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
 665       // chop off ind digits from the lower part of C1
 666       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
 667       tmp64 = C1.w[0];
 668       if (ind <= 19) {
 669         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
 670       } else {
 671         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
 672         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
 673       }
 674       if (C1.w[0] < tmp64)
 675         C1.w[1]++;
 676       // calculate C* and f*
 677       // C* is actually floor(C*) in this case
 678       // C* and f* need shifting and masking, as shown by
 679       // shiftright128[] and maskhigh128[]
 680       // 1 <= x <= 33
 681       // kx = 10^(-x) = ten2mk128[ind - 1]
 682       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
 683       // the approximation of 10^(-x) was rounded up to 118 bits
 684       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
 685       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
 686         Cstar.w[1] = P256.w[3];
 687         Cstar.w[0] = P256.w[2];
 688         fstar.w[3] = 0;
 689         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
 690         fstar.w[1] = P256.w[1];
 691         fstar.w[0] = P256.w[0];
 692       } else {  // 22 <= ind - 1 <= 33
 693         Cstar.w[1] = 0;
 694         Cstar.w[0] = P256.w[3];
 695         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
 696         fstar.w[2] = P256.w[2];
 697         fstar.w[1] = P256.w[1];
 698         fstar.w[0] = P256.w[0];
 699       }
 700       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
 701       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
 702       // if (0 < f* < 10^(-x)) then the result is a midpoint
 703       //   if floor(C*) is even then C* = floor(C*) - logical right
 704       //       shift; C* has p decimal digits, correct by Prop. 1)
 705       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
 706       //       shift; C* has p decimal digits, correct by Pr. 1)
 707       // else
 708       //   C* = floor(C*) (logical right shift; C has p decimal digits,
 709       //       correct by Property 1)
 710       // n = C* * 10^(e+x)
 711
 712       // shift right C* by Ex-128 = shiftright128[ind]
 713       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
 714       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
 715         Cstar.w[0] =
 716           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
 717         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
 718       } else {  // 22 <= ind - 1 <= 33
 719         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
 720       }
 721       // determine inexactness of the rounding of C*
 722       // if (0 < f* - 1/2 < 10^(-x)) then
 723       //   the result is exact
 724       // else // if (f* - 1/2 > T*) then
 725       //   the result is inexact
 726       if (ind - 1 <= 2) {
 727         if (fstar.w[1] > 0x8000000000000000ull ||
 728             (fstar.w[1] == 0x8000000000000000ull
 729              && fstar.w[0] > 0x0ull)) {
 730           // f* > 1/2 and the result may be exact
 731           tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
 732           if (tmp64 > ten2mk128trunc[ind - 1].w[1]
 733               || (tmp64 == ten2mk128trunc[ind - 1].w[1]
 734                   && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
 735             // set the inexact flag
 736             *pfpsf |= INEXACT_EXCEPTION;
 737           }     // else the result is exact
 738         } else {        // the result is inexact; f2* <= 1/2
 739           // set the inexact flag
 740           *pfpsf |= INEXACT_EXCEPTION;
 741         }
 742       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
 743         if (fstar.w[3] > 0x0 ||
 744             (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
 745             (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
 746              (fstar.w[1] || fstar.w[0]))) {
 747           // f2* > 1/2 and the result may be exact
 748           // Calculate f2* - 1/2
 749           tmp64 = fstar.w[2] - onehalf128[ind - 1];
 750           tmp64A = fstar.w[3];
 751           if (tmp64 > fstar.w[2])
 752             tmp64A--;
 753           if (tmp64A || tmp64
 754               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
 755               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
 756                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
 757             // set the inexact flag
 758             *pfpsf |= INEXACT_EXCEPTION;
 759           }     // else the result is exact
 760         } else {        // the result is inexact; f2* <= 1/2
 761           // set the inexact flag
 762           *pfpsf |= INEXACT_EXCEPTION;
 763         }
 764       } else {  // if 22 <= ind <= 33
 765         if (fstar.w[3] > onehalf128[ind - 1] ||
 766             (fstar.w[3] == onehalf128[ind - 1] &&
 767              (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
 768           // f2* > 1/2 and the result may be exact
 769           // Calculate f2* - 1/2
 770           tmp64 = fstar.w[3] - onehalf128[ind - 1];
 771           if (tmp64 || fstar.w[2]
 772               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
 773               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
 774                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
 775             // set the inexact flag
 776             *pfpsf |= INEXACT_EXCEPTION;
 777           }     // else the result is exact
 778         } else {        // the result is inexact; f2* <= 1/2
 779           // set the inexact flag
 780           *pfpsf |= INEXACT_EXCEPTION;
 781         }
 782       }
 783
 784       // if the result was a midpoint it was rounded away from zero, so
 785       // it will need a correction
 786       // check for midpoints
 787       if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
 788           && (fstar.w[1] || fstar.w[0])
 789           && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
 790               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
 791                   && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
 792         // the result is a midpoint; round to nearest
 793         if (Cstar.w[0] & 0x01) {        // Cstar.w[0] is odd; MP in [EVEN, ODD]
 794           // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
 795           Cstar.w[0]--; // Cstar.w[0] is now even
 796         }       // else MP in [ODD, EVEN]
 797       }
 798       res = Cstar.w[0]; // the result is positive
 799     } else if (exp == 0) {
 800       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
 801       // res = C (exact)
 802       res = C1.w[0];
 803     } else {
 804       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
 805       // res = C * 10^exp (exact) - must fit in 64 bits
 806       res = C1.w[0] * ten2k64[exp];
 807     }
 808   }
 809 }
 810
 811 BID_RETURN (res);
 812 }
 813
 814 /*****************************************************************************
 815  *  BID128_to_uint64_floor
 816  ****************************************************************************/
 817
 818 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
 819                                           bid128_to_uint64_floor, x)
 820
 821      UINT64 res;
 822      UINT64 x_sign;
 823      UINT64 x_exp;
 824      int exp;                   // unbiased exponent
 825   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
 826      BID_UI64DOUBLE tmp1;
 827      unsigned int x_nr_bits;
 828      int q, ind, shift;
 829      UINT128 C1, C;
 830      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
 831      UINT256 P256;
 832
 833   // unpack x
 834 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
 835 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
 836 C1.w[1] = x.w[1] & MASK_COEFF;
 837 C1.w[0] = x.w[0];
 838
 839   // check for NaN or Infinity
 840 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
 841     // x is special
 842 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
 843   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
 844     // set invalid flag
 845     *pfpsf |= INVALID_EXCEPTION;
 846     // return Integer Indefinite
 847     res = 0x8000000000000000ull;
 848   } else {      // x is QNaN
 849     // set invalid flag
 850     *pfpsf |= INVALID_EXCEPTION;
 851     // return Integer Indefinite
 852     res = 0x8000000000000000ull;
 853   }
 854   BID_RETURN (res);
 855 } else {        // x is not a NaN, so it must be infinity
 856   if (!x_sign) {        // x is +inf
 857     // set invalid flag
 858     *pfpsf |= INVALID_EXCEPTION;
 859     // return Integer Indefinite
 860     res = 0x8000000000000000ull;
 861   } else {      // x is -inf
 862     // set invalid flag
 863     *pfpsf |= INVALID_EXCEPTION;
 864     // return Integer Indefinite
 865     res = 0x8000000000000000ull;
 866   }
 867   BID_RETURN (res);
 868 }
 869 }
 870   // check for non-canonical values (after the check for special values)
 871 if ((C1.w[1] > 0x0001ed09bead87c0ull)
 872     || (C1.w[1] == 0x0001ed09bead87c0ull
 873         && (C1.w[0] > 0x378d8e63ffffffffull))
 874     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
 875   res = 0x0000000000000000ull;
 876   BID_RETURN (res);
 877 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
 878   // x is 0
 879   res = 0x0000000000000000ull;
 880   BID_RETURN (res);
 881 } else {        // x is not special and is not zero
 882
 883   // if n < 0 then n cannot be converted to uint64 with RM
 884   if (x_sign) { // if n < 0 and q + exp = 20
 885     // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 0
 886     // set invalid flag
 887     *pfpsf |= INVALID_EXCEPTION;
 888     // return Integer Indefinite
 889     res = 0x8000000000000000ull;
 890     BID_RETURN (res);
 891   }
 892   // q = nr. of decimal digits in x
 893   //  determine first the nr. of bits in x
 894   if (C1.w[1] == 0) {
 895     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
 896       // split the 64-bit value in two 32-bit halves to avoid rounding errors
 897       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
 898         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
 899         x_nr_bits =
 900           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 901       } else {  // x < 2^32
 902         tmp1.d = (double) (C1.w[0]);    // exact conversion
 903         x_nr_bits =
 904           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 905       }
 906     } else {    // if x < 2^53
 907       tmp1.d = (double) C1.w[0];        // exact conversion
 908       x_nr_bits =
 909         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 910     }
 911   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
 912     tmp1.d = (double) C1.w[1];  // exact conversion
 913     x_nr_bits =
 914       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 915   }
 916   q = nr_digits[x_nr_bits - 1].digits;
 917   if (q == 0) {
 918     q = nr_digits[x_nr_bits - 1].digits1;
 919     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
 920         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
 921             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
 922       q++;
 923   }
 924   exp = (x_exp >> 49) - 6176;
 925   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
 926     // set invalid flag
 927     *pfpsf |= INVALID_EXCEPTION;
 928     // return Integer Indefinite
 929     res = 0x8000000000000000ull;
 930     BID_RETURN (res);
 931   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
 932     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
 933     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
 934     // the cases that do not fit are identified here; the ones that fit
 935     // fall through and will be handled with other cases further,
 936     // under '1 <= q + exp <= 20'
 937     // if n > 0 and q + exp = 20
 938     // if n >= 2^64 then n is too large
 939     // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
 940     // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
 941     // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
 942     // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
 943     if (q == 1) {
 944       // C * 10^20 >= 0xa0000000000000000
 945       __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);    // 10^20 * C
 946       if (C.w[1] >= 0x0a) {
 947         // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
 948         // set invalid flag
 949         *pfpsf |= INVALID_EXCEPTION;
 950         // return Integer Indefinite
 951         res = 0x8000000000000000ull;
 952         BID_RETURN (res);
 953       }
 954       // else cases that can be rounded to a 64-bit int fall through
 955       // to '1 <= q + exp <= 20'
 956     } else if (q <= 19) {
 957       // C * 10^(21-q) >= 0xa0000000000000000
 958       __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
 959       if (C.w[1] >= 0x0a) {
 960         // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
 961         // set invalid flag
 962         *pfpsf |= INVALID_EXCEPTION;
 963         // return Integer Indefinite
 964         res = 0x8000000000000000ull;
 965         BID_RETURN (res);
 966       }
 967       // else cases that can be rounded to a 64-bit int fall through
 968       // to '1 <= q + exp <= 20'
 969     } else if (q == 20) {
 970       // C >= 0x10000000000000000
 971       if (C1.w[1] >= 0x01) {
 972         // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
 973         // set invalid flag
 974         *pfpsf |= INVALID_EXCEPTION;
 975         // return Integer Indefinite
 976         res = 0x8000000000000000ull;
 977         BID_RETURN (res);
 978       }
 979       // else cases that can be rounded to a 64-bit int fall through
 980       // to '1 <= q + exp <= 20'
 981     } else if (q == 21) {
 982       // C >= 0xa0000000000000000
 983       if (C1.w[1] >= 0x0a) {
 984         // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
 985         // set invalid flag
 986         *pfpsf |= INVALID_EXCEPTION;
 987         // return Integer Indefinite
 988         res = 0x8000000000000000ull;
 989         BID_RETURN (res);
 990       }
 991       // else cases that can be rounded to a 64-bit int fall through
 992       // to '1 <= q + exp <= 20'
 993     } else {    // if 22 <= q <= 34 => 1 <= q - 21 <= 13
 994       // C  >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
 995       C.w[1] = 0x0a;
 996       C.w[0] = 0x0000000000000000ull;
 997       __mul_128x64_to_128 (C, ten2k64[q - 21], C);
 998       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
 999         // set invalid flag
1000         *pfpsf |= INVALID_EXCEPTION;
1001         // return Integer Indefinite
1002         res = 0x8000000000000000ull;
1003         BID_RETURN (res);
1004       }
1005       // else cases that can be rounded to a 64-bit int fall through
1006       // to '1 <= q + exp <= 20'
1007     }
1008   }
1009   // n is not too large to be converted to int64 if 0 <= n < 2^64
1010   // Note: some of the cases tested for above fall through to this point
1011   if ((q + exp) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
1012     // return 0
1013     res = 0x0000000000000000ull;
1014     BID_RETURN (res);
1015   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1016     // 1 <= x < 2^64 so x can be rounded
1017     // down to a 64-bit unsigned signed integer
1018     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1019       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1020       // chop off ind digits from the lower part of C1
1021       // C1 fits in 127 bits
1022       // calculate C* and f*
1023       // C* is actually floor(C*) in this case
1024       // C* and f* need shifting and masking, as shown by
1025       // shiftright128[] and maskhigh128[]
1026       // 1 <= x <= 33
1027       // kx = 10^(-x) = ten2mk128[ind - 1]
1028       // C* = C1 * 10^(-x)
1029       // the approximation of 10^(-x) was rounded up to 118 bits
1030       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1031       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1032         Cstar.w[1] = P256.w[3];
1033         Cstar.w[0] = P256.w[2];
1034       } else {  // 22 <= ind - 1 <= 33
1035         Cstar.w[1] = 0;
1036         Cstar.w[0] = P256.w[3];
1037       }
1038       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1039       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1040       // C* = floor(C*) (logical right shift; C has p decimal digits,
1041       //     correct by Property 1)
1042       // n = C* * 10^(e+x)
1043
1044       // shift right C* by Ex-128 = shiftright128[ind]
1045       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
1046       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1047         Cstar.w[0] =
1048           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1049         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1050       } else {  // 22 <= ind - 1 <= 33
1051         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
1052       }
1053       res = Cstar.w[0]; // the result is positive
1054     } else if (exp == 0) {
1055       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
1056       // res = C (exact)
1057       res = C1.w[0];
1058     } else {
1059       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
1060       // res = C * 10^exp (exact) - must fit in 64 bits
1061       res = C1.w[0] * ten2k64[exp];
1062     }
1063   }
1064 }
1065
1066 BID_RETURN (res);
1067 }
1068
1069 /*****************************************************************************
1070  *  BID128_to_uint64_xfloor
1071  ****************************************************************************/
1072
1073 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
1074                                           bid128_to_uint64_xfloor, x)
1075
1076      UINT64 res;
1077      UINT64 x_sign;
1078      UINT64 x_exp;
1079      int exp;                   // unbiased exponent
1080   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1081      BID_UI64DOUBLE tmp1;
1082      unsigned int x_nr_bits;
1083      int q, ind, shift;
1084      UINT128 C1, C;
1085      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
1086      UINT256 fstar;
1087      UINT256 P256;
1088
1089   // unpack x
1090 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1091 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
1092 C1.w[1] = x.w[1] & MASK_COEFF;
1093 C1.w[0] = x.w[0];
1094
1095   // check for NaN or Infinity
1096 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1097     // x is special
1098 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1099   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1100     // set invalid flag
1101     *pfpsf |= INVALID_EXCEPTION;
1102     // return Integer Indefinite
1103     res = 0x8000000000000000ull;
1104   } else {      // x is QNaN
1105     // set invalid flag
1106     *pfpsf |= INVALID_EXCEPTION;
1107     // return Integer Indefinite
1108     res = 0x8000000000000000ull;
1109   }
1110   BID_RETURN (res);
1111 } else {        // x is not a NaN, so it must be infinity
1112   if (!x_sign) {        // x is +inf
1113     // set invalid flag
1114     *pfpsf |= INVALID_EXCEPTION;
1115     // return Integer Indefinite
1116     res = 0x8000000000000000ull;
1117   } else {      // x is -inf
1118     // set invalid flag
1119     *pfpsf |= INVALID_EXCEPTION;
1120     // return Integer Indefinite
1121     res = 0x8000000000000000ull;
1122   }
1123   BID_RETURN (res);
1124 }
1125 }
1126   // check for non-canonical values (after the check for special values)
1127 if ((C1.w[1] > 0x0001ed09bead87c0ull)
1128     || (C1.w[1] == 0x0001ed09bead87c0ull
1129         && (C1.w[0] > 0x378d8e63ffffffffull))
1130     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1131   res = 0x0000000000000000ull;
1132   BID_RETURN (res);
1133 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1134   // x is 0
1135   res = 0x0000000000000000ull;
1136   BID_RETURN (res);
1137 } else {        // x is not special and is not zero
1138
1139   // if n < 0 then n cannot be converted to uint64 with RM
1140   if (x_sign) { // if n < 0 and q + exp = 20
1141     // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 0
1142     // set invalid flag
1143     *pfpsf |= INVALID_EXCEPTION;
1144     // return Integer Indefinite
1145     res = 0x8000000000000000ull;
1146     BID_RETURN (res);
1147   }
1148   // q = nr. of decimal digits in x
1149   //  determine first the nr. of bits in x
1150   if (C1.w[1] == 0) {
1151     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
1152       // split the 64-bit value in two 32-bit halves to avoid rounding errors
1153       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
1154         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
1155         x_nr_bits =
1156           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1157       } else {  // x < 2^32
1158         tmp1.d = (double) (C1.w[0]);    // exact conversion
1159         x_nr_bits =
1160           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1161       }
1162     } else {    // if x < 2^53
1163       tmp1.d = (double) C1.w[0];        // exact conversion
1164       x_nr_bits =
1165         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1166     }
1167   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1168     tmp1.d = (double) C1.w[1];  // exact conversion
1169     x_nr_bits =
1170       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1171   }
1172   q = nr_digits[x_nr_bits - 1].digits;
1173   if (q == 0) {
1174     q = nr_digits[x_nr_bits - 1].digits1;
1175     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1176         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1177             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1178       q++;
1179   }
1180   exp = (x_exp >> 49) - 6176;
1181   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1182     // set invalid flag
1183     *pfpsf |= INVALID_EXCEPTION;
1184     // return Integer Indefinite
1185     res = 0x8000000000000000ull;
1186     BID_RETURN (res);
1187   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1188     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1189     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1190     // the cases that do not fit are identified here; the ones that fit
1191     // fall through and will be handled with other cases further,
1192     // under '1 <= q + exp <= 20'
1193     // if n > 0 and q + exp = 20
1194     // if n >= 2^64 then n is too large
1195     // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
1196     // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
1197     // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
1198     // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
1199     if (q == 1) {
1200       // C * 10^20 >= 0xa0000000000000000
1201       __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);    // 10^20 * C
1202       if (C.w[1] >= 0x0a) {
1203         // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
1204         // set invalid flag
1205         *pfpsf |= INVALID_EXCEPTION;
1206         // return Integer Indefinite
1207         res = 0x8000000000000000ull;
1208         BID_RETURN (res);
1209       }
1210       // else cases that can be rounded to a 64-bit int fall through
1211       // to '1 <= q + exp <= 20'
1212     } else if (q <= 19) {
1213       // C * 10^(21-q) >= 0xa0000000000000000
1214       __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
1215       if (C.w[1] >= 0x0a) {
1216         // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
1217         // set invalid flag
1218         *pfpsf |= INVALID_EXCEPTION;
1219         // return Integer Indefinite
1220         res = 0x8000000000000000ull;
1221         BID_RETURN (res);
1222       }
1223       // else cases that can be rounded to a 64-bit int fall through
1224       // to '1 <= q + exp <= 20'
1225     } else if (q == 20) {
1226       // C >= 0x10000000000000000
1227       if (C1.w[1] >= 0x01) {
1228         // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
1229         // set invalid flag
1230         *pfpsf |= INVALID_EXCEPTION;
1231         // return Integer Indefinite
1232         res = 0x8000000000000000ull;
1233         BID_RETURN (res);
1234       }
1235       // else cases that can be rounded to a 64-bit int fall through
1236       // to '1 <= q + exp <= 20'
1237     } else if (q == 21) {
1238       // C >= 0xa0000000000000000
1239       if (C1.w[1] >= 0x0a) {
1240         // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
1241         // set invalid flag
1242         *pfpsf |= INVALID_EXCEPTION;
1243         // return Integer Indefinite
1244         res = 0x8000000000000000ull;
1245         BID_RETURN (res);
1246       }
1247       // else cases that can be rounded to a 64-bit int fall through
1248       // to '1 <= q + exp <= 20'
1249     } else {    // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1250       // C  >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
1251       C.w[1] = 0x0a;
1252       C.w[0] = 0x0000000000000000ull;
1253       __mul_128x64_to_128 (C, ten2k64[q - 21], C);
1254       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
1255         // set invalid flag
1256         *pfpsf |= INVALID_EXCEPTION;
1257         // return Integer Indefinite
1258         res = 0x8000000000000000ull;
1259         BID_RETURN (res);
1260       }
1261       // else cases that can be rounded to a 64-bit int fall through
1262       // to '1 <= q + exp <= 20'
1263     }
1264   }
1265   // n is not too large to be converted to int64 if 0 <= n < 2^64
1266   // Note: some of the cases tested for above fall through to this point
1267   if ((q + exp) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
1268     // set inexact flag
1269     *pfpsf |= INEXACT_EXCEPTION;
1270     // return 0
1271     res = 0x0000000000000000ull;
1272     BID_RETURN (res);
1273   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1274     // 1 <= x < 2^64 so x can be rounded
1275     // down to a 64-bit unsigned signed integer
1276     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1277       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1278       // chop off ind digits from the lower part of C1
1279       // C1 fits in 127 bits
1280       // calculate C* and f*
1281       // C* is actually floor(C*) in this case
1282       // C* and f* need shifting and masking, as shown by
1283       // shiftright128[] and maskhigh128[]
1284       // 1 <= x <= 33
1285       // kx = 10^(-x) = ten2mk128[ind - 1]
1286       // C* = C1 * 10^(-x)
1287       // the approximation of 10^(-x) was rounded up to 118 bits
1288       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1289       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1290         Cstar.w[1] = P256.w[3];
1291         Cstar.w[0] = P256.w[2];
1292         fstar.w[3] = 0;
1293         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1294         fstar.w[1] = P256.w[1];
1295         fstar.w[0] = P256.w[0];
1296       } else {  // 22 <= ind - 1 <= 33
1297         Cstar.w[1] = 0;
1298         Cstar.w[0] = P256.w[3];
1299         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1300         fstar.w[2] = P256.w[2];
1301         fstar.w[1] = P256.w[1];
1302         fstar.w[0] = P256.w[0];
1303       }
1304       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1305       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1306       // C* = floor(C*) (logical right shift; C has p decimal digits,
1307       //     correct by Property 1)
1308       // n = C* * 10^(e+x)
1309
1310       // shift right C* by Ex-128 = shiftright128[ind]
1311       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
1312       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1313         Cstar.w[0] =
1314           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1315         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1316       } else {  // 22 <= ind - 1 <= 33
1317         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
1318       }
1319       // determine inexactness of the rounding of C*
1320       // if (0 < f* < 10^(-x)) then
1321       //   the result is exact
1322       // else // if (f* > T*) then
1323       //   the result is inexact
1324       if (ind - 1 <= 2) {
1325         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1] ||
1326             (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] &&
1327              fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1328           // set the inexact flag
1329           *pfpsf |= INEXACT_EXCEPTION;
1330         }       // else the result is exact
1331       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
1332         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1333             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1334                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1335           // set the inexact flag
1336           *pfpsf |= INEXACT_EXCEPTION;
1337         }       // else the result is exact
1338       } else {  // if 22 <= ind <= 33
1339         if (fstar.w[3] || fstar.w[2]
1340             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1341             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1342                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1343           // set the inexact flag
1344           *pfpsf |= INEXACT_EXCEPTION;
1345         }       // else the result is exact
1346       }
1347
1348       res = Cstar.w[0]; // the result is positive
1349     } else if (exp == 0) {
1350       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
1351       // res = C (exact)
1352       res = C1.w[0];
1353     } else {
1354       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
1355       // res = C * 10^exp (exact) - must fit in 64 bits
1356       res = C1.w[0] * ten2k64[exp];
1357     }
1358   }
1359 }
1360
1361 BID_RETURN (res);
1362 }
1363
1364 /*****************************************************************************
1365  *  BID128_to_uint64_ceil
1366  ****************************************************************************/
1367
1368 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64, bid128_to_uint64_ceil,
1369                                           x)
1370
1371      UINT64 res;
1372      UINT64 x_sign;
1373      UINT64 x_exp;
1374      int exp;                   // unbiased exponent
1375   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1376      BID_UI64DOUBLE tmp1;
1377      unsigned int x_nr_bits;
1378      int q, ind, shift;
1379      UINT128 C1, C;
1380      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
1381      UINT256 fstar;
1382      UINT256 P256;
1383
1384   // unpack x
1385 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1386 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
1387 C1.w[1] = x.w[1] & MASK_COEFF;
1388 C1.w[0] = x.w[0];
1389
1390   // check for NaN or Infinity
1391 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1392     // x is special
1393 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1394   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1395     // set invalid flag
1396     *pfpsf |= INVALID_EXCEPTION;
1397     // return Integer Indefinite
1398     res = 0x8000000000000000ull;
1399   } else {      // x is QNaN
1400     // set invalid flag
1401     *pfpsf |= INVALID_EXCEPTION;
1402     // return Integer Indefinite
1403     res = 0x8000000000000000ull;
1404   }
1405   BID_RETURN (res);
1406 } else {        // x is not a NaN, so it must be infinity
1407   if (!x_sign) {        // x is +inf
1408     // set invalid flag
1409     *pfpsf |= INVALID_EXCEPTION;
1410     // return Integer Indefinite
1411     res = 0x8000000000000000ull;
1412   } else {      // x is -inf
1413     // set invalid flag
1414     *pfpsf |= INVALID_EXCEPTION;
1415     // return Integer Indefinite
1416     res = 0x8000000000000000ull;
1417   }
1418   BID_RETURN (res);
1419 }
1420 }
1421   // check for non-canonical values (after the check for special values)
1422 if ((C1.w[1] > 0x0001ed09bead87c0ull)
1423     || (C1.w[1] == 0x0001ed09bead87c0ull
1424         && (C1.w[0] > 0x378d8e63ffffffffull))
1425     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1426   res = 0x0000000000000000ull;
1427   BID_RETURN (res);
1428 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1429   // x is 0
1430   res = 0x0000000000000000ull;
1431   BID_RETURN (res);
1432 } else {        // x is not special and is not zero
1433
1434   // q = nr. of decimal digits in x
1435   //  determine first the nr. of bits in x
1436   if (C1.w[1] == 0) {
1437     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
1438       // split the 64-bit value in two 32-bit halves to avoid rounding errors
1439       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
1440         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
1441         x_nr_bits =
1442           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1443       } else {  // x < 2^32
1444         tmp1.d = (double) (C1.w[0]);    // exact conversion
1445         x_nr_bits =
1446           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1447       }
1448     } else {    // if x < 2^53
1449       tmp1.d = (double) C1.w[0];        // exact conversion
1450       x_nr_bits =
1451         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1452     }
1453   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1454     tmp1.d = (double) C1.w[1];  // exact conversion
1455     x_nr_bits =
1456       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1457   }
1458   q = nr_digits[x_nr_bits - 1].digits;
1459   if (q == 0) {
1460     q = nr_digits[x_nr_bits - 1].digits1;
1461     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1462         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1463             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1464       q++;
1465   }
1466   exp = (x_exp >> 49) - 6176;
1467   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1468     // set invalid flag
1469     *pfpsf |= INVALID_EXCEPTION;
1470     // return Integer Indefinite
1471     res = 0x8000000000000000ull;
1472     BID_RETURN (res);
1473   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1474     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1475     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1476     // the cases that do not fit are identified here; the ones that fit
1477     // fall through and will be handled with other cases further,
1478     // under '1 <= q + exp <= 20'
1479     if (x_sign) {       // if n < 0 and q + exp = 20
1480       // if n <= -1 then n cannot be converted to uint64 with RZ
1481       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
1482       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
1483       // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
1484       if (q == 21) {
1485         // C >= a
1486         if (C1.w[1] != 0 || C1.w[0] >= 0x0aull) {
1487           // set invalid flag
1488           *pfpsf |= INVALID_EXCEPTION;
1489           // return Integer Indefinite
1490           res = 0x8000000000000000ull;
1491           BID_RETURN (res);
1492         }
1493         // else cases that can be rounded to 64-bit unsigned int fall through
1494         // to '1 <= q + exp <= 20'
1495       } else {
1496         // if 1 <= q <= 20
1497         //   C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
1498         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1499         //  C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
1500         // set invalid flag
1501         *pfpsf |= INVALID_EXCEPTION;
1502         // return Integer Indefinite
1503         res = 0x8000000000000000ull;
1504         BID_RETURN (res);
1505       }
1506     } else {    // if n > 0 and q + exp = 20
1507       // if n > 2^64 - 1 then n is too large
1508       // <=> c(0)c(1)...c(19).c(20)...c(q-1) > 2^64 - 1
1509       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 > 2^64 - 1
1510       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 > 10 * (2^64 - 1)
1511       // <=> C * 10^(21-q) > 0x9fffffffffffffff6, 1<=q<=34
1512       if (q == 1) {
1513         // C * 10^20 > 0x9fffffffffffffff6
1514         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
1515         if (C.w[1] > 0x09 || (C.w[1] == 0x09
1516                               && C.w[0] > 0xfffffffffffffff6ull)) {
1517           // set invalid flag
1518           *pfpsf |= INVALID_EXCEPTION;
1519           // return Integer Indefinite
1520           res = 0x8000000000000000ull;
1521           BID_RETURN (res);
1522         }
1523         // else cases that can be rounded to a 64-bit int fall through
1524         // to '1 <= q + exp <= 20'
1525       } else if (q <= 19) {
1526         // C * 10^(21-q) > 0x9fffffffffffffff6
1527         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
1528         if (C.w[1] > 0x09 || (C.w[1] == 0x09
1529                               && C.w[0] > 0xfffffffffffffff6ull)) {
1530           // set invalid flag
1531           *pfpsf |= INVALID_EXCEPTION;
1532           // return Integer Indefinite
1533           res = 0x8000000000000000ull;
1534           BID_RETURN (res);
1535         }
1536         // else cases that can be rounded to a 64-bit int fall through
1537         // to '1 <= q + exp <= 20'
1538       } else if (q == 20) {
1539         // C > 0xffffffffffffffff
1540         if (C1.w[1]) {
1541           // set invalid flag
1542           *pfpsf |= INVALID_EXCEPTION;
1543           // return Integer Indefinite
1544           res = 0x8000000000000000ull;
1545           BID_RETURN (res);
1546         }
1547         // else cases that can be rounded to a 64-bit int fall through
1548         // to '1 <= q + exp <= 20'
1549       } else if (q == 21) {
1550         // C > 0x9fffffffffffffff6
1551         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
1552                                && C1.w[0] > 0xfffffffffffffff6ull)) {
1553           // set invalid flag
1554           *pfpsf |= INVALID_EXCEPTION;
1555           // return Integer Indefinite
1556           res = 0x8000000000000000ull;
1557           BID_RETURN (res);
1558         }
1559         // else cases that can be rounded to a 64-bit int fall through
1560         // to '1 <= q + exp <= 20'
1561       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1562         // C  > 10^(q-21) * 0x9fffffffffffffff6 max 44 bits x 68 bits
1563         C.w[1] = 0x09;
1564         C.w[0] = 0xfffffffffffffff6ull;
1565         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
1566         if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
1567           // set invalid flag
1568           *pfpsf |= INVALID_EXCEPTION;
1569           // return Integer Indefinite
1570           res = 0x8000000000000000ull;
1571           BID_RETURN (res);
1572         }
1573         // else cases that can be rounded to a 64-bit int fall through
1574         // to '1 <= q + exp <= 20'
1575       }
1576     }
1577   }
1578   // n is not too large to be converted to int64 if -1 < n <= 2^64 - 1
1579   // Note: some of the cases tested for above fall through to this point
1580   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1581     // return 0 or 1
1582     if (x_sign)
1583       res = 0x0000000000000000ull;
1584     else
1585       res = 0x0000000000000001ull;
1586     BID_RETURN (res);
1587   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1588     // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
1589     // to zero to a 64-bit unsigned signed integer
1590     if (x_sign) {       // x <= -1
1591       // set invalid flag
1592       *pfpsf |= INVALID_EXCEPTION;
1593       // return Integer Indefinite
1594       res = 0x8000000000000000ull;
1595       BID_RETURN (res);
1596     }
1597     // 1 <= x <= 2^64 - 1 so x can be rounded
1598     // to zero to a 64-bit unsigned integer
1599     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1600       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1601       // chop off ind digits from the lower part of C1
1602       // C1 fits in 127 bits
1603       // calculate C* and f*
1604       // C* is actually floor(C*) in this case
1605       // C* and f* need shifting and masking, as shown by
1606       // shiftright128[] and maskhigh128[]
1607       // 1 <= x <= 33
1608       // kx = 10^(-x) = ten2mk128[ind - 1]
1609       // C* = C1 * 10^(-x)
1610       // the approximation of 10^(-x) was rounded up to 118 bits
1611       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1612       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1613         Cstar.w[1] = P256.w[3];
1614         Cstar.w[0] = P256.w[2];
1615         fstar.w[3] = 0;
1616         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1617         fstar.w[1] = P256.w[1];
1618         fstar.w[0] = P256.w[0];
1619       } else {  // 22 <= ind - 1 <= 33
1620         Cstar.w[1] = 0;
1621         Cstar.w[0] = P256.w[3];
1622         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1623         fstar.w[2] = P256.w[2];
1624         fstar.w[1] = P256.w[1];
1625         fstar.w[0] = P256.w[0];
1626       }
1627       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1628       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1629       // C* = floor(C*) (logical right shift; C has p decimal digits,
1630       //     correct by Property 1)
1631       // n = C* * 10^(e+x)
1632
1633       // shift right C* by Ex-128 = shiftright128[ind]
1634       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
1635       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1636         Cstar.w[0] =
1637           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1638         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1639       } else {  // 22 <= ind - 1 <= 33
1640         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
1641       }
1642       // if the result is positive and inexact, need to add 1 to it
1643
1644       // determine inexactness of the rounding of C*
1645       // if (0 < f* < 10^(-x)) then
1646       //   the result is exact
1647       // else // if (f* > T*) then
1648       //   the result is inexact
1649       if (ind - 1 <= 2) {
1650         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1651             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1652                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1653           if (!x_sign) {        // positive and inexact
1654             Cstar.w[0]++;
1655             if (Cstar.w[0] == 0x0)
1656               Cstar.w[1]++;
1657           }
1658         }       // else the result is exact
1659       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
1660         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1661             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1662                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1663           if (!x_sign) {        // positive and inexact
1664             Cstar.w[0]++;
1665             if (Cstar.w[0] == 0x0)
1666               Cstar.w[1]++;
1667           }
1668         }       // else the result is exact
1669       } else {  // if 22 <= ind <= 33
1670         if (fstar.w[3] || fstar.w[2]
1671             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1672             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1673                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1674           if (!x_sign) {        // positive and inexact
1675             Cstar.w[0]++;
1676             if (Cstar.w[0] == 0x0)
1677               Cstar.w[1]++;
1678           }
1679         }       // else the result is exact
1680       }
1681
1682       res = Cstar.w[0]; // the result is positive
1683     } else if (exp == 0) {
1684       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
1685       // res = C (exact)
1686       res = C1.w[0];
1687     } else {
1688       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
1689       // res = C * 10^exp (exact) - must fit in 64 bits
1690       res = C1.w[0] * ten2k64[exp];
1691     }
1692   }
1693 }
1694
1695 BID_RETURN (res);
1696 }
1697
1698 /*****************************************************************************
1699  *  BID128_to_uint64_xceil
1700  ****************************************************************************/
1701
1702 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
1703                                           bid128_to_uint64_xceil, x)
1704
1705      UINT64 res;
1706      UINT64 x_sign;
1707      UINT64 x_exp;
1708      int exp;                   // unbiased exponent
1709   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1710      BID_UI64DOUBLE tmp1;
1711      unsigned int x_nr_bits;
1712      int q, ind, shift;
1713      UINT128 C1, C;
1714      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
1715      UINT256 fstar;
1716      UINT256 P256;
1717
1718   // unpack x
1719 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1720 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
1721 C1.w[1] = x.w[1] & MASK_COEFF;
1722 C1.w[0] = x.w[0];
1723
1724   // check for NaN or Infinity
1725 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1726     // x is special
1727 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1728   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1729     // set invalid flag
1730     *pfpsf |= INVALID_EXCEPTION;
1731     // return Integer Indefinite
1732     res = 0x8000000000000000ull;
1733   } else {      // x is QNaN
1734     // set invalid flag
1735     *pfpsf |= INVALID_EXCEPTION;
1736     // return Integer Indefinite
1737     res = 0x8000000000000000ull;
1738   }
1739   BID_RETURN (res);
1740 } else {        // x is not a NaN, so it must be infinity
1741   if (!x_sign) {        // x is +inf
1742     // set invalid flag
1743     *pfpsf |= INVALID_EXCEPTION;
1744     // return Integer Indefinite
1745     res = 0x8000000000000000ull;
1746   } else {      // x is -inf
1747     // set invalid flag
1748     *pfpsf |= INVALID_EXCEPTION;
1749     // return Integer Indefinite
1750     res = 0x8000000000000000ull;
1751   }
1752   BID_RETURN (res);
1753 }
1754 }
1755   // check for non-canonical values (after the check for special values)
1756 if ((C1.w[1] > 0x0001ed09bead87c0ull)
1757     || (C1.w[1] == 0x0001ed09bead87c0ull
1758         && (C1.w[0] > 0x378d8e63ffffffffull))
1759     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1760   res = 0x0000000000000000ull;
1761   BID_RETURN (res);
1762 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1763   // x is 0
1764   res = 0x0000000000000000ull;
1765   BID_RETURN (res);
1766 } else {        // x is not special and is not zero
1767
1768   // q = nr. of decimal digits in x
1769   //  determine first the nr. of bits in x
1770   if (C1.w[1] == 0) {
1771     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
1772       // split the 64-bit value in two 32-bit halves to avoid rounding errors
1773       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
1774         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
1775         x_nr_bits =
1776           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1777       } else {  // x < 2^32
1778         tmp1.d = (double) (C1.w[0]);    // exact conversion
1779         x_nr_bits =
1780           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1781       }
1782     } else {    // if x < 2^53
1783       tmp1.d = (double) C1.w[0];        // exact conversion
1784       x_nr_bits =
1785         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1786     }
1787   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1788     tmp1.d = (double) C1.w[1];  // exact conversion
1789     x_nr_bits =
1790       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1791   }
1792   q = nr_digits[x_nr_bits - 1].digits;
1793   if (q == 0) {
1794     q = nr_digits[x_nr_bits - 1].digits1;
1795     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1796         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1797             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1798       q++;
1799   }
1800   exp = (x_exp >> 49) - 6176;
1801   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1802     // set invalid flag
1803     *pfpsf |= INVALID_EXCEPTION;
1804     // return Integer Indefinite
1805     res = 0x8000000000000000ull;
1806     BID_RETURN (res);
1807   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1808     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1809     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1810     // the cases that do not fit are identified here; the ones that fit
1811     // fall through and will be handled with other cases further,
1812     // under '1 <= q + exp <= 20'
1813     if (x_sign) {       // if n < 0 and q + exp = 20
1814       // if n <= -1 then n cannot be converted to uint64 with RZ
1815       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
1816       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
1817       // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
1818       if (q == 21) {
1819         // C >= a
1820         if (C1.w[1] != 0 || C1.w[0] >= 0x0aull) {
1821           // set invalid flag
1822           *pfpsf |= INVALID_EXCEPTION;
1823           // return Integer Indefinite
1824           res = 0x8000000000000000ull;
1825           BID_RETURN (res);
1826         }
1827         // else cases that can be rounded to 64-bit unsigned int fall through
1828         // to '1 <= q + exp <= 20'
1829       } else {
1830         // if 1 <= q <= 20
1831         //   C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
1832         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1833         //  C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
1834         // set invalid flag
1835         *pfpsf |= INVALID_EXCEPTION;
1836         // return Integer Indefinite
1837         res = 0x8000000000000000ull;
1838         BID_RETURN (res);
1839       }
1840     } else {    // if n > 0 and q + exp = 20
1841       // if n > 2^64 - 1 then n is too large
1842       // <=> c(0)c(1)...c(19).c(20)...c(q-1) > 2^64 - 1
1843       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 > 2^64 - 1
1844       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 > 10 * (2^64 - 1)
1845       // <=> C * 10^(21-q) > 0x9fffffffffffffff6, 1<=q<=34
1846       if (q == 1) {
1847         // C * 10^20 > 0x9fffffffffffffff6
1848         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
1849         if (C.w[1] > 0x09 || (C.w[1] == 0x09
1850                               && C.w[0] > 0xfffffffffffffff6ull)) {
1851           // set invalid flag
1852           *pfpsf |= INVALID_EXCEPTION;
1853           // return Integer Indefinite
1854           res = 0x8000000000000000ull;
1855           BID_RETURN (res);
1856         }
1857         // else cases that can be rounded to a 64-bit int fall through
1858         // to '1 <= q + exp <= 20'
1859       } else if (q <= 19) {
1860         // C * 10^(21-q) > 0x9fffffffffffffff6
1861         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
1862         if (C.w[1] > 0x09 || (C.w[1] == 0x09
1863                               && C.w[0] > 0xfffffffffffffff6ull)) {
1864           // set invalid flag
1865           *pfpsf |= INVALID_EXCEPTION;
1866           // return Integer Indefinite
1867           res = 0x8000000000000000ull;
1868           BID_RETURN (res);
1869         }
1870         // else cases that can be rounded to a 64-bit int fall through
1871         // to '1 <= q + exp <= 20'
1872       } else if (q == 20) {
1873         // C > 0xffffffffffffffff
1874         if (C1.w[1]) {
1875           // set invalid flag
1876           *pfpsf |= INVALID_EXCEPTION;
1877           // return Integer Indefinite
1878           res = 0x8000000000000000ull;
1879           BID_RETURN (res);
1880         }
1881         // else cases that can be rounded to a 64-bit int fall through
1882         // to '1 <= q + exp <= 20'
1883       } else if (q == 21) {
1884         // C > 0x9fffffffffffffff6
1885         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
1886                                && C1.w[0] > 0xfffffffffffffff6ull)) {
1887           // set invalid flag
1888           *pfpsf |= INVALID_EXCEPTION;
1889           // return Integer Indefinite
1890           res = 0x8000000000000000ull;
1891           BID_RETURN (res);
1892         }
1893         // else cases that can be rounded to a 64-bit int fall through
1894         // to '1 <= q + exp <= 20'
1895       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1896         // C  > 10^(q-21) * 0x9fffffffffffffff6 max 44 bits x 68 bits
1897         C.w[1] = 0x09;
1898         C.w[0] = 0xfffffffffffffff6ull;
1899         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
1900         if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
1901           // set invalid flag
1902           *pfpsf |= INVALID_EXCEPTION;
1903           // return Integer Indefinite
1904           res = 0x8000000000000000ull;
1905           BID_RETURN (res);
1906         }
1907         // else cases that can be rounded to a 64-bit int fall through
1908         // to '1 <= q + exp <= 20'
1909       }
1910     }
1911   }
1912   // n is not too large to be converted to int64 if -1 < n <= 2^64 - 1
1913   // Note: some of the cases tested for above fall through to this point
1914   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1915     // set inexact flag
1916     *pfpsf |= INEXACT_EXCEPTION;
1917     // return 0 or 1
1918     if (x_sign)
1919       res = 0x0000000000000000ull;
1920     else
1921       res = 0x0000000000000001ull;
1922     BID_RETURN (res);
1923   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1924     // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
1925     // to zero to a 64-bit unsigned signed integer
1926     if (x_sign) {       // x <= -1
1927       // set invalid flag
1928       *pfpsf |= INVALID_EXCEPTION;
1929       // return Integer Indefinite
1930       res = 0x8000000000000000ull;
1931       BID_RETURN (res);
1932     }
1933     // 1 <= x <= 2^64 - 1 so x can be rounded
1934     // to zero to a 64-bit unsigned integer
1935     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1936       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1937       // chop off ind digits from the lower part of C1
1938       // C1 fits in 127 bits
1939       // calculate C* and f*
1940       // C* is actually floor(C*) in this case
1941       // C* and f* need shifting and masking, as shown by
1942       // shiftright128[] and maskhigh128[]
1943       // 1 <= x <= 33
1944       // kx = 10^(-x) = ten2mk128[ind - 1]
1945       // C* = C1 * 10^(-x)
1946       // the approximation of 10^(-x) was rounded up to 118 bits
1947       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1948       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1949         Cstar.w[1] = P256.w[3];
1950         Cstar.w[0] = P256.w[2];
1951         fstar.w[3] = 0;
1952         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1953         fstar.w[1] = P256.w[1];
1954         fstar.w[0] = P256.w[0];
1955       } else {  // 22 <= ind - 1 <= 33
1956         Cstar.w[1] = 0;
1957         Cstar.w[0] = P256.w[3];
1958         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1959         fstar.w[2] = P256.w[2];
1960         fstar.w[1] = P256.w[1];
1961         fstar.w[0] = P256.w[0];
1962       }
1963       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1964       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1965       // C* = floor(C*) (logical right shift; C has p decimal digits,
1966       //     correct by Property 1)
1967       // n = C* * 10^(e+x)
1968
1969       // shift right C* by Ex-128 = shiftright128[ind]
1970       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
1971       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1972         Cstar.w[0] =
1973           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1974         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1975       } else {  // 22 <= ind - 1 <= 33
1976         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
1977       }
1978       // if the result is positive and inexact, need to add 1 to it
1979
1980       // determine inexactness of the rounding of C*
1981       // if (0 < f* < 10^(-x)) then
1982       //   the result is exact
1983       // else // if (f* > T*) then
1984       //   the result is inexact
1985       if (ind - 1 <= 2) {
1986         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1987             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1988                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1989           if (!x_sign) {        // positive and inexact
1990             Cstar.w[0]++;
1991             if (Cstar.w[0] == 0x0)
1992               Cstar.w[1]++;
1993           }
1994           // set the inexact flag
1995           *pfpsf |= INEXACT_EXCEPTION;
1996         }       // else the result is exact
1997       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
1998         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1999             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2000                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2001           if (!x_sign) {        // positive and inexact
2002             Cstar.w[0]++;
2003             if (Cstar.w[0] == 0x0)
2004               Cstar.w[1]++;
2005           }
2006           // set the inexact flag
2007           *pfpsf |= INEXACT_EXCEPTION;
2008         }       // else the result is exact
2009       } else {  // if 22 <= ind <= 33
2010         if (fstar.w[3] || fstar.w[2]
2011             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2012             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2013                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2014           if (!x_sign) {        // positive and inexact
2015             Cstar.w[0]++;
2016             if (Cstar.w[0] == 0x0)
2017               Cstar.w[1]++;
2018           }
2019           // set the inexact flag
2020           *pfpsf |= INEXACT_EXCEPTION;
2021         }       // else the result is exact
2022       }
2023
2024       res = Cstar.w[0]; // the result is positive
2025     } else if (exp == 0) {
2026       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2027       // res = C (exact)
2028       res = C1.w[0];
2029     } else {
2030       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2031       // res = C * 10^exp (exact) - must fit in 64 bits
2032       res = C1.w[0] * ten2k64[exp];
2033     }
2034   }
2035 }
2036
2037 BID_RETURN (res);
2038 }
2039
2040 /*****************************************************************************
2041  *  BID128_to_uint64_int
2042  ****************************************************************************/
2043
2044 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64, bid128_to_uint64_int,
2045                                           x)
2046
2047      UINT64 res;
2048      UINT64 x_sign;
2049      UINT64 x_exp;
2050      int exp;                   // unbiased exponent
2051   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2052      BID_UI64DOUBLE tmp1;
2053      unsigned int x_nr_bits;
2054      int q, ind, shift;
2055      UINT128 C1, C;
2056      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
2057      UINT256 P256;
2058
2059   // unpack x
2060 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
2061 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
2062 C1.w[1] = x.w[1] & MASK_COEFF;
2063 C1.w[0] = x.w[0];
2064
2065   // check for NaN or Infinity
2066 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2067     // x is special
2068 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
2069   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
2070     // set invalid flag
2071     *pfpsf |= INVALID_EXCEPTION;
2072     // return Integer Indefinite
2073     res = 0x8000000000000000ull;
2074   } else {      // x is QNaN
2075     // set invalid flag
2076     *pfpsf |= INVALID_EXCEPTION;
2077     // return Integer Indefinite
2078     res = 0x8000000000000000ull;
2079   }
2080   BID_RETURN (res);
2081 } else {        // x is not a NaN, so it must be infinity
2082   if (!x_sign) {        // x is +inf
2083     // set invalid flag
2084     *pfpsf |= INVALID_EXCEPTION;
2085     // return Integer Indefinite
2086     res = 0x8000000000000000ull;
2087   } else {      // x is -inf
2088     // set invalid flag
2089     *pfpsf |= INVALID_EXCEPTION;
2090     // return Integer Indefinite
2091     res = 0x8000000000000000ull;
2092   }
2093   BID_RETURN (res);
2094 }
2095 }
2096   // check for non-canonical values (after the check for special values)
2097 if ((C1.w[1] > 0x0001ed09bead87c0ull)
2098     || (C1.w[1] == 0x0001ed09bead87c0ull
2099         && (C1.w[0] > 0x378d8e63ffffffffull))
2100     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2101   res = 0x0000000000000000ull;
2102   BID_RETURN (res);
2103 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2104   // x is 0
2105   res = 0x0000000000000000ull;
2106   BID_RETURN (res);
2107 } else {        // x is not special and is not zero
2108
2109   // q = nr. of decimal digits in x
2110   //  determine first the nr. of bits in x
2111   if (C1.w[1] == 0) {
2112     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
2113       // split the 64-bit value in two 32-bit halves to avoid rounding errors
2114       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
2115         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
2116         x_nr_bits =
2117           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2118       } else {  // x < 2^32
2119         tmp1.d = (double) (C1.w[0]);    // exact conversion
2120         x_nr_bits =
2121           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2122       }
2123     } else {    // if x < 2^53
2124       tmp1.d = (double) C1.w[0];        // exact conversion
2125       x_nr_bits =
2126         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2127     }
2128   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2129     tmp1.d = (double) C1.w[1];  // exact conversion
2130     x_nr_bits =
2131       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2132   }
2133   q = nr_digits[x_nr_bits - 1].digits;
2134   if (q == 0) {
2135     q = nr_digits[x_nr_bits - 1].digits1;
2136     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2137         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2138             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2139       q++;
2140   }
2141   exp = (x_exp >> 49) - 6176;
2142
2143   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
2144     // set invalid flag
2145     *pfpsf |= INVALID_EXCEPTION;
2146     // return Integer Indefinite
2147     res = 0x8000000000000000ull;
2148     BID_RETURN (res);
2149   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
2150     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2151     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
2152     // the cases that do not fit are identified here; the ones that fit
2153     // fall through and will be handled with other cases further,
2154     // under '1 <= q + exp <= 20'
2155     if (x_sign) {       // if n < 0 and q + exp = 20
2156       // if n <= -1 then n cannot be converted to uint64 with RZ
2157       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
2158       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
2159       // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
2160       if (q == 21) {
2161         // C >= a
2162         if (C1.w[1] != 0 || C1.w[0] >= 0x0aull) {
2163           // set invalid flag
2164           *pfpsf |= INVALID_EXCEPTION;
2165           // return Integer Indefinite
2166           res = 0x8000000000000000ull;
2167           BID_RETURN (res);
2168         }
2169         // else cases that can be rounded to 64-bit unsigned int fall through
2170         // to '1 <= q + exp <= 20'
2171       } else {
2172         // if 1 <= q <= 20
2173         //   C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
2174         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2175         //  C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
2176         // set invalid flag
2177         *pfpsf |= INVALID_EXCEPTION;
2178         // return Integer Indefinite
2179         res = 0x8000000000000000ull;
2180         BID_RETURN (res);
2181       }
2182     } else {    // if n > 0 and q + exp = 20
2183       // if n >= 2^64 then n is too large
2184       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
2185       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
2186       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
2187       // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
2188       if (q == 1) {
2189         // C * 10^20 >= 0xa0000000000000000
2190         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
2191         if (C.w[1] >= 0x0a) {
2192           // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2193           // set invalid flag
2194           *pfpsf |= INVALID_EXCEPTION;
2195           // return Integer Indefinite
2196           res = 0x8000000000000000ull;
2197           BID_RETURN (res);
2198         }
2199         // else cases that can be rounded to a 64-bit int fall through
2200         // to '1 <= q + exp <= 20'
2201       } else if (q <= 19) {
2202         // C * 10^(21-q) >= 0xa0000000000000000
2203         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
2204         if (C.w[1] >= 0x0a) {
2205           // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2206           // set invalid flag
2207           *pfpsf |= INVALID_EXCEPTION;
2208           // return Integer Indefinite
2209           res = 0x8000000000000000ull;
2210           BID_RETURN (res);
2211         }
2212         // else cases that can be rounded to a 64-bit int fall through
2213         // to '1 <= q + exp <= 20'
2214       } else if (q == 20) {
2215         // C >= 0x10000000000000000
2216         if (C1.w[1] >= 0x01) {
2217           // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
2218           // set invalid flag
2219           *pfpsf |= INVALID_EXCEPTION;
2220           // return Integer Indefinite
2221           res = 0x8000000000000000ull;
2222           BID_RETURN (res);
2223         }
2224         // else cases that can be rounded to a 64-bit int fall through
2225         // to '1 <= q + exp <= 20'
2226       } else if (q == 21) {
2227         // C >= 0xa0000000000000000
2228         if (C1.w[1] >= 0x0a) {
2229           // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
2230           // set invalid flag
2231           *pfpsf |= INVALID_EXCEPTION;
2232           // return Integer Indefinite
2233           res = 0x8000000000000000ull;
2234           BID_RETURN (res);
2235         }
2236         // else cases that can be rounded to a 64-bit int fall through
2237         // to '1 <= q + exp <= 20'
2238       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2239         // C  >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
2240         C.w[1] = 0x0a;
2241         C.w[0] = 0x0000000000000000ull;
2242         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
2243         if (C1.w[1] > C.w[1]
2244             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2245           // set invalid flag
2246           *pfpsf |= INVALID_EXCEPTION;
2247           // return Integer Indefinite
2248           res = 0x8000000000000000ull;
2249           BID_RETURN (res);
2250         }
2251         // else cases that can be rounded to a 64-bit int fall through
2252         // to '1 <= q + exp <= 20'
2253       }
2254     }
2255   }
2256   // n is not too large to be converted to int64 if -1 < n < 2^64
2257   // Note: some of the cases tested for above fall through to this point
2258   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
2259     // return 0
2260     res = 0x0000000000000000ull;
2261     BID_RETURN (res);
2262   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
2263     // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
2264     // to zero to a 64-bit unsigned signed integer
2265     if (x_sign) {       // x <= -1
2266       // set invalid flag
2267       *pfpsf |= INVALID_EXCEPTION;
2268       // return Integer Indefinite
2269       res = 0x8000000000000000ull;
2270       BID_RETURN (res);
2271     }
2272     // 1 <= x < 2^64 so x can be rounded
2273     // to zero to a 64-bit unsigned integer
2274     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
2275       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
2276       // chop off ind digits from the lower part of C1
2277       // C1 fits in 127 bits
2278       // calculate C* and f*
2279       // C* is actually floor(C*) in this case
2280       // C* and f* need shifting and masking, as shown by
2281       // shiftright128[] and maskhigh128[]
2282       // 1 <= x <= 33
2283       // kx = 10^(-x) = ten2mk128[ind - 1]
2284       // C* = C1 * 10^(-x)
2285       // the approximation of 10^(-x) was rounded up to 118 bits
2286       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2287       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2288         Cstar.w[1] = P256.w[3];
2289         Cstar.w[0] = P256.w[2];
2290       } else {  // 22 <= ind - 1 <= 33
2291         Cstar.w[1] = 0;
2292         Cstar.w[0] = P256.w[3];
2293       }
2294       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2295       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2296       // C* = floor(C*) (logical right shift; C has p decimal digits,
2297       //     correct by Property 1)
2298       // n = C* * 10^(e+x)
2299
2300       // shift right C* by Ex-128 = shiftright128[ind]
2301       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
2302       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2303         Cstar.w[0] =
2304           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2305         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2306       } else {  // 22 <= ind - 1 <= 33
2307         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
2308       }
2309       res = Cstar.w[0]; // the result is positive
2310     } else if (exp == 0) {
2311       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2312       // res = C (exact)
2313       res = C1.w[0];
2314     } else {
2315       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2316       // res = C * 10^exp (exact) - must fit in 64 bits
2317       res = C1.w[0] * ten2k64[exp];
2318     }
2319   }
2320 }
2321
2322 BID_RETURN (res);
2323 }
2324
2325 /*****************************************************************************
2326  *  BID128_to_uint64_xint
2327  ****************************************************************************/
2328
2329 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64, bid128_to_uint64_xint,
2330                                           x)
2331
2332      UINT64 res;
2333      UINT64 x_sign;
2334      UINT64 x_exp;
2335      int exp;                   // unbiased exponent
2336   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2337      BID_UI64DOUBLE tmp1;
2338      unsigned int x_nr_bits;
2339      int q, ind, shift;
2340      UINT128 C1, C;
2341      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
2342      UINT256 fstar;
2343      UINT256 P256;
2344
2345   // unpack x
2346 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
2347 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
2348 C1.w[1] = x.w[1] & MASK_COEFF;
2349 C1.w[0] = x.w[0];
2350
2351   // check for NaN or Infinity
2352 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2353     // x is special
2354 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
2355   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
2356     // set invalid flag
2357     *pfpsf |= INVALID_EXCEPTION;
2358     // return Integer Indefinite
2359     res = 0x8000000000000000ull;
2360   } else {      // x is QNaN
2361     // set invalid flag
2362     *pfpsf |= INVALID_EXCEPTION;
2363     // return Integer Indefinite
2364     res = 0x8000000000000000ull;
2365   }
2366   BID_RETURN (res);
2367 } else {        // x is not a NaN, so it must be infinity
2368   if (!x_sign) {        // x is +inf
2369     // set invalid flag
2370     *pfpsf |= INVALID_EXCEPTION;
2371     // return Integer Indefinite
2372     res = 0x8000000000000000ull;
2373   } else {      // x is -inf
2374     // set invalid flag
2375     *pfpsf |= INVALID_EXCEPTION;
2376     // return Integer Indefinite
2377     res = 0x8000000000000000ull;
2378   }
2379   BID_RETURN (res);
2380 }
2381 }
2382   // check for non-canonical values (after the check for special values)
2383 if ((C1.w[1] > 0x0001ed09bead87c0ull)
2384     || (C1.w[1] == 0x0001ed09bead87c0ull
2385         && (C1.w[0] > 0x378d8e63ffffffffull))
2386     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2387   res = 0x0000000000000000ull;
2388   BID_RETURN (res);
2389 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2390   // x is 0
2391   res = 0x0000000000000000ull;
2392   BID_RETURN (res);
2393 } else {        // x is not special and is not zero
2394
2395   // q = nr. of decimal digits in x
2396   //  determine first the nr. of bits in x
2397   if (C1.w[1] == 0) {
2398     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
2399       // split the 64-bit value in two 32-bit halves to avoid rounding errors
2400       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
2401         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
2402         x_nr_bits =
2403           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2404       } else {  // x < 2^32
2405         tmp1.d = (double) (C1.w[0]);    // exact conversion
2406         x_nr_bits =
2407           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2408       }
2409     } else {    // if x < 2^53
2410       tmp1.d = (double) C1.w[0];        // exact conversion
2411       x_nr_bits =
2412         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2413     }
2414   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2415     tmp1.d = (double) C1.w[1];  // exact conversion
2416     x_nr_bits =
2417       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2418   }
2419   q = nr_digits[x_nr_bits - 1].digits;
2420   if (q == 0) {
2421     q = nr_digits[x_nr_bits - 1].digits1;
2422     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2423         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2424             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2425       q++;
2426   }
2427   exp = (x_exp >> 49) - 6176;
2428   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
2429     // set invalid flag
2430     *pfpsf |= INVALID_EXCEPTION;
2431     // return Integer Indefinite
2432     res = 0x8000000000000000ull;
2433     BID_RETURN (res);
2434   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
2435     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2436     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
2437     // the cases that do not fit are identified here; the ones that fit
2438     // fall through and will be handled with other cases further,
2439     // under '1 <= q + exp <= 20'
2440     if (x_sign) {       // if n < 0 and q + exp = 20
2441       // if n <= -1 then n cannot be converted to uint64 with RZ
2442       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
2443       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
2444       // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
2445       if (q == 21) {
2446         // C >= a
2447         if (C1.w[1] != 0 || C1.w[0] >= 0x0aull) {
2448           // set invalid flag
2449           *pfpsf |= INVALID_EXCEPTION;
2450           // return Integer Indefinite
2451           res = 0x8000000000000000ull;
2452           BID_RETURN (res);
2453         }
2454         // else cases that can be rounded to 64-bit unsigned int fall through
2455         // to '1 <= q + exp <= 20'
2456       } else {
2457         // if 1 <= q <= 20
2458         //   C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
2459         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2460         //  C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
2461         // set invalid flag
2462         *pfpsf |= INVALID_EXCEPTION;
2463         // return Integer Indefinite
2464         res = 0x8000000000000000ull;
2465         BID_RETURN (res);
2466       }
2467     } else {    // if n > 0 and q + exp = 20
2468       // if n >= 2^64 then n is too large
2469       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
2470       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
2471       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
2472       // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
2473       if (q == 1) {
2474         // C * 10^20 >= 0xa0000000000000000
2475         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
2476         if (C.w[1] >= 0x0a) {
2477           // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2478           // set invalid flag
2479           *pfpsf |= INVALID_EXCEPTION;
2480           // return Integer Indefinite
2481           res = 0x8000000000000000ull;
2482           BID_RETURN (res);
2483         }
2484         // else cases that can be rounded to a 64-bit int fall through
2485         // to '1 <= q + exp <= 20'
2486       } else if (q <= 19) {
2487         // C * 10^(21-q) >= 0xa0000000000000000
2488         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
2489         if (C.w[1] >= 0x0a) {
2490           // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2491           // set invalid flag
2492           *pfpsf |= INVALID_EXCEPTION;
2493           // return Integer Indefinite
2494           res = 0x8000000000000000ull;
2495           BID_RETURN (res);
2496         }
2497         // else cases that can be rounded to a 64-bit int fall through
2498         // to '1 <= q + exp <= 20'
2499       } else if (q == 20) {
2500         // C >= 0x10000000000000000
2501         if (C1.w[1] >= 0x01) {
2502           // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
2503           // set invalid flag
2504           *pfpsf |= INVALID_EXCEPTION;
2505           // return Integer Indefinite
2506           res = 0x8000000000000000ull;
2507           BID_RETURN (res);
2508         }
2509         // else cases that can be rounded to a 64-bit int fall through
2510         // to '1 <= q + exp <= 20'
2511       } else if (q == 21) {
2512         // C >= 0xa0000000000000000
2513         if (C1.w[1] >= 0x0a) {
2514           // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
2515           // set invalid flag
2516           *pfpsf |= INVALID_EXCEPTION;
2517           // return Integer Indefinite
2518           res = 0x8000000000000000ull;
2519           BID_RETURN (res);
2520         }
2521         // else cases that can be rounded to a 64-bit int fall through
2522         // to '1 <= q + exp <= 20'
2523       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2524         // C  >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
2525         C.w[1] = 0x0a;
2526         C.w[0] = 0x0000000000000000ull;
2527         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
2528         if (C1.w[1] > C.w[1]
2529             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2530           // set invalid flag
2531           *pfpsf |= INVALID_EXCEPTION;
2532           // return Integer Indefinite
2533           res = 0x8000000000000000ull;
2534           BID_RETURN (res);
2535         }
2536         // else cases that can be rounded to a 64-bit int fall through
2537         // to '1 <= q + exp <= 20'
2538       }
2539     }
2540   }
2541   // n is not too large to be converted to int64 if -1 < n < 2^64
2542   // Note: some of the cases tested for above fall through to this point
2543   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
2544     // set inexact flag
2545     *pfpsf |= INEXACT_EXCEPTION;
2546     // return 0
2547     res = 0x0000000000000000ull;
2548     BID_RETURN (res);
2549   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
2550     // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
2551     // to zero to a 64-bit unsigned signed integer
2552     if (x_sign) {       // x <= -1
2553       // set invalid flag
2554       *pfpsf |= INVALID_EXCEPTION;
2555       // return Integer Indefinite
2556       res = 0x8000000000000000ull;
2557       BID_RETURN (res);
2558     }
2559     // 1 <= x < 2^64 so x can be rounded
2560     // to zero to a 64-bit unsigned integer
2561     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
2562       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
2563       // chop off ind digits from the lower part of C1
2564       // C1 fits in 127 bits
2565       // calculate C* and f*
2566       // C* is actually floor(C*) in this case
2567       // C* and f* need shifting and masking, as shown by
2568       // shiftright128[] and maskhigh128[]
2569       // 1 <= x <= 33
2570       // kx = 10^(-x) = ten2mk128[ind - 1]
2571       // C* = C1 * 10^(-x)
2572       // the approximation of 10^(-x) was rounded up to 118 bits
2573       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2574       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2575         Cstar.w[1] = P256.w[3];
2576         Cstar.w[0] = P256.w[2];
2577         fstar.w[3] = 0;
2578         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
2579         fstar.w[1] = P256.w[1];
2580         fstar.w[0] = P256.w[0];
2581       } else {  // 22 <= ind - 1 <= 33
2582         Cstar.w[1] = 0;
2583         Cstar.w[0] = P256.w[3];
2584         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
2585         fstar.w[2] = P256.w[2];
2586         fstar.w[1] = P256.w[1];
2587         fstar.w[0] = P256.w[0];
2588       }
2589       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2590       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2591       // C* = floor(C*) (logical right shift; C has p decimal digits,
2592       //     correct by Property 1)
2593       // n = C* * 10^(e+x)
2594
2595       // shift right C* by Ex-128 = shiftright128[ind]
2596       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
2597       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2598         Cstar.w[0] =
2599           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2600         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2601       } else {  // 22 <= ind - 1 <= 33
2602         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
2603       }
2604       // determine inexactness of the rounding of C*
2605       // if (0 < f* < 10^(-x)) then
2606       //   the result is exact
2607       // else // if (f* > T*) then
2608       //   the result is inexact
2609       if (ind - 1 <= 2) {
2610         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2611             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2612                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2613           // set the inexact flag
2614           *pfpsf |= INEXACT_EXCEPTION;
2615         }       // else the result is exact
2616       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
2617         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2618             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2619                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2620           // set the inexact flag
2621           *pfpsf |= INEXACT_EXCEPTION;
2622         }       // else the result is exact
2623       } else {  // if 22 <= ind <= 33
2624         if (fstar.w[3] || fstar.w[2]
2625             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2626             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2627                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2628           // set the inexact flag
2629           *pfpsf |= INEXACT_EXCEPTION;
2630         }       // else the result is exact
2631       }
2632
2633       res = Cstar.w[0]; // the result is positive
2634     } else if (exp == 0) {
2635       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2636       // res = C (exact)
2637       res = C1.w[0];
2638     } else {
2639       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2640       // res = C * 10^exp (exact) - must fit in 64 bits
2641       res = C1.w[0] * ten2k64[exp];
2642     }
2643   }
2644 }
2645
2646 BID_RETURN (res);
2647 }
2648
2649 /*****************************************************************************
2650  *  BID128_to_uint64_rninta
2651  ****************************************************************************/
2652
2653 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
2654                                           bid128_to_uint64_rninta, x)
2655
2656      UINT64 res;
2657      UINT64 x_sign;
2658      UINT64 x_exp;
2659      int exp;                   // unbiased exponent
2660   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2661      UINT64 tmp64;
2662      BID_UI64DOUBLE tmp1;
2663      unsigned int x_nr_bits;
2664      int q, ind, shift;
2665      UINT128 C1, C;
2666      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
2667      UINT256 P256;
2668
2669   // unpack x
2670 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
2671 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
2672 C1.w[1] = x.w[1] & MASK_COEFF;
2673 C1.w[0] = x.w[0];
2674
2675   // check for NaN or Infinity
2676 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2677     // x is special
2678 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
2679   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
2680     // set invalid flag
2681     *pfpsf |= INVALID_EXCEPTION;
2682     // return Integer Indefinite
2683     res = 0x8000000000000000ull;
2684   } else {      // x is QNaN
2685     // set invalid flag
2686     *pfpsf |= INVALID_EXCEPTION;
2687     // return Integer Indefinite
2688     res = 0x8000000000000000ull;
2689   }
2690   BID_RETURN (res);
2691 } else {        // x is not a NaN, so it must be infinity
2692   if (!x_sign) {        // x is +inf
2693     // set invalid flag
2694     *pfpsf |= INVALID_EXCEPTION;
2695     // return Integer Indefinite
2696     res = 0x8000000000000000ull;
2697   } else {      // x is -inf
2698     // set invalid flag
2699     *pfpsf |= INVALID_EXCEPTION;
2700     // return Integer Indefinite
2701     res = 0x8000000000000000ull;
2702   }
2703   BID_RETURN (res);
2704 }
2705 }
2706   // check for non-canonical values (after the check for special values)
2707 if ((C1.w[1] > 0x0001ed09bead87c0ull)
2708     || (C1.w[1] == 0x0001ed09bead87c0ull
2709         && (C1.w[0] > 0x378d8e63ffffffffull))
2710     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2711   res = 0x0000000000000000ull;
2712   BID_RETURN (res);
2713 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2714   // x is 0
2715   res = 0x0000000000000000ull;
2716   BID_RETURN (res);
2717 } else {        // x is not special and is not zero
2718
2719   // q = nr. of decimal digits in x
2720   //  determine first the nr. of bits in x
2721   if (C1.w[1] == 0) {
2722     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
2723       // split the 64-bit value in two 32-bit halves to avoid rounding errors
2724       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
2725         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
2726         x_nr_bits =
2727           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2728       } else {  // x < 2^32
2729         tmp1.d = (double) (C1.w[0]);    // exact conversion
2730         x_nr_bits =
2731           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2732       }
2733     } else {    // if x < 2^53
2734       tmp1.d = (double) C1.w[0];        // exact conversion
2735       x_nr_bits =
2736         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2737     }
2738   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2739     tmp1.d = (double) C1.w[1];  // exact conversion
2740     x_nr_bits =
2741       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2742   }
2743   q = nr_digits[x_nr_bits - 1].digits;
2744   if (q == 0) {
2745     q = nr_digits[x_nr_bits - 1].digits1;
2746     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2747         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2748             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2749       q++;
2750   }
2751   exp = (x_exp >> 49) - 6176;
2752   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
2753     // set invalid flag
2754     *pfpsf |= INVALID_EXCEPTION;
2755     // return Integer Indefinite
2756     res = 0x8000000000000000ull;
2757     BID_RETURN (res);
2758   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
2759     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2760     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
2761     // the cases that do not fit are identified here; the ones that fit
2762     // fall through and will be handled with other cases further,
2763     // under '1 <= q + exp <= 20'
2764     if (x_sign) {       // if n < 0 and q + exp = 20
2765       // if n <= -1/2 then n cannot be converted to uint64 with RN
2766       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1/2
2767       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x05, 1<=q<=34
2768       // <=> C * 10^(21-q) >= 0x05, 1<=q<=34
2769       if (q == 21) {
2770         // C >= 5
2771         if (C1.w[1] != 0 || C1.w[0] >= 0x05ull) {
2772           // set invalid flag
2773           *pfpsf |= INVALID_EXCEPTION;
2774           // return Integer Indefinite
2775           res = 0x8000000000000000ull;
2776           BID_RETURN (res);
2777         }
2778         // else cases that can be rounded to 64-bit unsigned int fall through
2779         // to '1 <= q + exp <= 20'
2780       } else {
2781         // if 1 <= q <= 20
2782         //   C * 10^(21-q) >= 5 is true because C >= 1 and 10^(21-q) >= 10
2783         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2784         //  C >= 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
2785         // set invalid flag
2786         *pfpsf |= INVALID_EXCEPTION;
2787         // return Integer Indefinite
2788         res = 0x8000000000000000ull;
2789         BID_RETURN (res);
2790       }
2791     } else {    // if n > 0 and q + exp = 20
2792       // if n >= 2^64 - 1/2 then n is too large
2793       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
2794       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
2795       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
2796       // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
2797       if (q == 1) {
2798         // C * 10^20 >= 0x9fffffffffffffffb
2799         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
2800         if (C.w[1] > 0x09 || (C.w[1] == 0x09
2801                               && C.w[0] >= 0xfffffffffffffffbull)) {
2802           // set invalid flag
2803           *pfpsf |= INVALID_EXCEPTION;
2804           // return Integer Indefinite
2805           res = 0x8000000000000000ull;
2806           BID_RETURN (res);
2807         }
2808         // else cases that can be rounded to a 64-bit int fall through
2809         // to '1 <= q + exp <= 20'
2810       } else if (q <= 19) {
2811         // C * 10^(21-q) >= 0x9fffffffffffffffb
2812         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
2813         if (C.w[1] > 0x09 || (C.w[1] == 0x09
2814                               && C.w[0] >= 0xfffffffffffffffbull)) {
2815           // set invalid flag
2816           *pfpsf |= INVALID_EXCEPTION;
2817           // return Integer Indefinite
2818           res = 0x8000000000000000ull;
2819           BID_RETURN (res);
2820         }
2821         // else cases that can be rounded to a 64-bit int fall through
2822         // to '1 <= q + exp <= 20'
2823       } else if (q == 20) {
2824         // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
2825         C.w[0] = C1.w[0] + C1.w[0];
2826         C.w[1] = C1.w[1] + C1.w[1];
2827         if (C.w[0] < C1.w[0])
2828           C.w[1]++;
2829         if (C.w[1] > 0x01 || (C.w[1] == 0x01
2830                               && C.w[0] >= 0xffffffffffffffffull)) {
2831           // set invalid flag
2832           *pfpsf |= INVALID_EXCEPTION;
2833           // return Integer Indefinite
2834           res = 0x8000000000000000ull;
2835           BID_RETURN (res);
2836         }
2837         // else cases that can be rounded to a 64-bit int fall through
2838         // to '1 <= q + exp <= 20'
2839       } else if (q == 21) {
2840         // C >= 0x9fffffffffffffffb
2841         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
2842                                && C1.w[0] >= 0xfffffffffffffffbull)) {
2843           // set invalid flag
2844           *pfpsf |= INVALID_EXCEPTION;
2845           // return Integer Indefinite
2846           res = 0x8000000000000000ull;
2847           BID_RETURN (res);
2848         }
2849         // else cases that can be rounded to a 64-bit int fall through
2850         // to '1 <= q + exp <= 20'
2851       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2852         // C  >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
2853         C.w[1] = 0x09;
2854         C.w[0] = 0xfffffffffffffffbull;
2855         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
2856         if (C1.w[1] > C.w[1]
2857             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2858           // set invalid flag
2859           *pfpsf |= INVALID_EXCEPTION;
2860           // return Integer Indefinite
2861           res = 0x8000000000000000ull;
2862           BID_RETURN (res);
2863         }
2864         // else cases that can be rounded to a 64-bit int fall through
2865         // to '1 <= q + exp <= 20'
2866       }
2867     }
2868   }
2869   // n is not too large to be converted to int64 if -1/2 < n < 2^64 - 1/2
2870   // Note: some of the cases tested for above fall through to this point
2871   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
2872     // return 0
2873     res = 0x0000000000000000ull;
2874     BID_RETURN (res);
2875   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
2876     // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
2877     //   res = 0
2878     // else if x > 0
2879     //   res = +1
2880     // else // if x < 0
2881     //   invalid exc
2882     ind = q - 1;
2883     if (ind <= 18) {    // 0 <= ind <= 18
2884       if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
2885         res = 0x0000000000000000ull;    // return 0
2886       } else if (!x_sign) {     // n > 0
2887         res = 0x00000001;       // return +1
2888       } else {
2889         // set invalid flag
2890         *pfpsf |= INVALID_EXCEPTION;
2891         // return Integer Indefinite
2892         res = 0x8000000000000000ull;
2893         BID_RETURN (res);
2894       }
2895     } else {    // 19 <= ind <= 33
2896       if ((C1.w[1] < midpoint128[ind - 19].w[1])
2897           || ((C1.w[1] == midpoint128[ind - 19].w[1])
2898               && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
2899         res = 0x0000000000000000ull;    // return 0
2900       } else if (!x_sign) {     // n > 0
2901         res = 0x00000001;       // return +1
2902       } else {
2903         // set invalid flag
2904         *pfpsf |= INVALID_EXCEPTION;
2905         // return Integer Indefinite
2906         res = 0x8000000000000000ull;
2907         BID_RETURN (res);
2908       }
2909     }
2910   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
2911     // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
2912     // to nearest to a 64-bit unsigned signed integer
2913     if (x_sign) {       // x <= -1
2914       // set invalid flag
2915       *pfpsf |= INVALID_EXCEPTION;
2916       // return Integer Indefinite
2917       res = 0x8000000000000000ull;
2918       BID_RETURN (res);
2919     }
2920     // 1 <= x < 2^64-1/2 so x can be rounded
2921     // to nearest to a 64-bit unsigned integer
2922     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
2923       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
2924       // chop off ind digits from the lower part of C1
2925       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2926       tmp64 = C1.w[0];
2927       if (ind <= 19) {
2928         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2929       } else {
2930         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2931         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2932       }
2933       if (C1.w[0] < tmp64)
2934         C1.w[1]++;
2935       // calculate C* and f*
2936       // C* is actually floor(C*) in this case
2937       // C* and f* need shifting and masking, as shown by
2938       // shiftright128[] and maskhigh128[]
2939       // 1 <= x <= 33
2940       // kx = 10^(-x) = ten2mk128[ind - 1]
2941       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2942       // the approximation of 10^(-x) was rounded up to 118 bits
2943       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2944       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2945         Cstar.w[1] = P256.w[3];
2946         Cstar.w[0] = P256.w[2];
2947       } else {  // 22 <= ind - 1 <= 33
2948         Cstar.w[1] = 0;
2949         Cstar.w[0] = P256.w[3];
2950       }
2951       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2952       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2953       // if (0 < f* < 10^(-x)) then the result is a midpoint
2954       //   if floor(C*) is even then C* = floor(C*) - logical right
2955       //       shift; C* has p decimal digits, correct by Prop. 1)
2956       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
2957       //       shift; C* has p decimal digits, correct by Pr. 1)
2958       // else
2959       //   C* = floor(C*) (logical right shift; C has p decimal digits,
2960       //       correct by Property 1)
2961       // n = C* * 10^(e+x)
2962
2963       // shift right C* by Ex-128 = shiftright128[ind]
2964       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
2965       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2966         Cstar.w[0] =
2967           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2968         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2969       } else {  // 22 <= ind - 1 <= 33
2970         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
2971       }
2972
2973       // if the result was a midpoint it was rounded away from zero
2974       res = Cstar.w[0]; // the result is positive
2975     } else if (exp == 0) {
2976       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2977       // res = C (exact)
2978       res = C1.w[0];
2979     } else {
2980       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2981       // res = C * 10^exp (exact) - must fit in 64 bits
2982       res = C1.w[0] * ten2k64[exp];
2983     }
2984   }
2985 }
2986
2987 BID_RETURN (res);
2988 }
2989
2990 /*****************************************************************************
2991  *  BID128_to_uint64_xrninta
2992  ****************************************************************************/
2993
2994 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
2995                                           bid128_to_uint64_xrninta, x)
2996
2997      UINT64 res;
2998      UINT64 x_sign;
2999      UINT64 x_exp;
3000      int exp;                   // unbiased exponent
3001   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
3002      UINT64 tmp64, tmp64A;
3003      BID_UI64DOUBLE tmp1;
3004      unsigned int x_nr_bits;
3005      int q, ind, shift;
3006      UINT128 C1, C;
3007      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
3008      UINT256 fstar;
3009      UINT256 P256;
3010
3011   // unpack x
3012 x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
3013 x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
3014 C1.w[1] = x.w[1] & MASK_COEFF;
3015 C1.w[0] = x.w[0];
3016
3017   // check for NaN or Infinity
3018 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
3019     // x is special
3020 if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
3021   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
3022     // set invalid flag
3023     *pfpsf |= INVALID_EXCEPTION;
3024     // return Integer Indefinite
3025     res = 0x8000000000000000ull;
3026   } else {      // x is QNaN
3027     // set invalid flag
3028     *pfpsf |= INVALID_EXCEPTION;
3029     // return Integer Indefinite
3030     res = 0x8000000000000000ull;
3031   }
3032   BID_RETURN (res);
3033 } else {        // x is not a NaN, so it must be infinity
3034   if (!x_sign) {        // x is +inf
3035     // set invalid flag
3036     *pfpsf |= INVALID_EXCEPTION;
3037     // return Integer Indefinite
3038     res = 0x8000000000000000ull;
3039   } else {      // x is -inf
3040     // set invalid flag
3041     *pfpsf |= INVALID_EXCEPTION;
3042     // return Integer Indefinite
3043     res = 0x8000000000000000ull;
3044   }
3045   BID_RETURN (res);
3046 }
3047 }
3048   // check for non-canonical values (after the check for special values)
3049 if ((C1.w[1] > 0x0001ed09bead87c0ull)
3050     || (C1.w[1] == 0x0001ed09bead87c0ull
3051         && (C1.w[0] > 0x378d8e63ffffffffull))
3052     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
3053   res = 0x0000000000000000ull;
3054   BID_RETURN (res);
3055 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
3056   // x is 0
3057   res = 0x0000000000000000ull;
3058   BID_RETURN (res);
3059 } else {        // x is not special and is not zero
3060
3061   // q = nr. of decimal digits in x
3062   //  determine first the nr. of bits in x
3063   if (C1.w[1] == 0) {
3064     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
3065       // split the 64-bit value in two 32-bit halves to avoid rounding errors
3066       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
3067         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
3068         x_nr_bits =
3069           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3070       } else {  // x < 2^32
3071         tmp1.d = (double) (C1.w[0]);    // exact conversion
3072         x_nr_bits =
3073           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3074       }
3075     } else {    // if x < 2^53
3076       tmp1.d = (double) C1.w[0];        // exact conversion
3077       x_nr_bits =
3078         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3079     }
3080   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
3081     tmp1.d = (double) C1.w[1];  // exact conversion
3082     x_nr_bits =
3083       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3084   }
3085   q = nr_digits[x_nr_bits - 1].digits;
3086   if (q == 0) {
3087     q = nr_digits[x_nr_bits - 1].digits1;
3088     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
3089         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
3090             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
3091       q++;
3092   }
3093   exp = (x_exp >> 49) - 6176;
3094
3095   if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
3096     // set invalid flag
3097     *pfpsf |= INVALID_EXCEPTION;
3098     // return Integer Indefinite
3099     res = 0x8000000000000000ull;
3100     BID_RETURN (res);
3101   } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
3102     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
3103     // so x rounded to an integer may or may not fit in an unsigned 64-bit int
3104     // the cases that do not fit are identified here; the ones that fit
3105     // fall through and will be handled with other cases further,
3106     // under '1 <= q + exp <= 20'
3107     if (x_sign) {       // if n < 0 and q + exp = 20
3108       // if n <= -1/2 then n cannot be converted to uint64 with RN
3109       // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1/2
3110       // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x05, 1<=q<=34
3111       // <=> C * 10^(21-q) >= 0x05, 1<=q<=34
3112       if (q == 21) {
3113         // C >= 5
3114         if (C1.w[1] != 0 || C1.w[0] >= 0x05ull) {
3115           // set invalid flag
3116           *pfpsf |= INVALID_EXCEPTION;
3117           // return Integer Indefinite
3118           res = 0x8000000000000000ull;
3119           BID_RETURN (res);
3120         }
3121         // else cases that can be rounded to 64-bit unsigned int fall through
3122         // to '1 <= q + exp <= 20'
3123       } else {
3124         // if 1 <= q <= 20
3125         //   C * 10^(21-q) >= 5 is true because C >= 1 and 10^(21-q) >= 10
3126         // if 22 <= q <= 34 => 1 <= q - 21 <= 13
3127         //  C >= 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
3128         // set invalid flag
3129         *pfpsf |= INVALID_EXCEPTION;
3130         // return Integer Indefinite
3131         res = 0x8000000000000000ull;
3132         BID_RETURN (res);
3133       }
3134     } else {    // if n > 0 and q + exp = 20
3135       // if n >= 2^64 - 1/2 then n is too large
3136       // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
3137       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
3138       // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
3139       // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
3140       if (q == 1) {
3141         // C * 10^20 >= 0x9fffffffffffffffb
3142         __mul_128x64_to_128 (C, C1.w[0], ten2k128[0]);  // 10^20 * C
3143         if (C.w[1] > 0x09 || (C.w[1] == 0x09
3144                               && C.w[0] >= 0xfffffffffffffffbull)) {
3145           // set invalid flag
3146           *pfpsf |= INVALID_EXCEPTION;
3147           // return Integer Indefinite
3148           res = 0x8000000000000000ull;
3149           BID_RETURN (res);
3150         }
3151         // else cases that can be rounded to a 64-bit int fall through
3152         // to '1 <= q + exp <= 20'
3153       } else if (q <= 19) {
3154         // C * 10^(21-q) >= 0x9fffffffffffffffb
3155         __mul_64x64_to_128MACH (C, C1.w[0], ten2k64[21 - q]);
3156         if (C.w[1] > 0x09 || (C.w[1] == 0x09
3157                               && C.w[0] >= 0xfffffffffffffffbull)) {
3158           // set invalid flag
3159           *pfpsf |= INVALID_EXCEPTION;
3160           // return Integer Indefinite
3161           res = 0x8000000000000000ull;
3162           BID_RETURN (res);
3163         }
3164         // else cases that can be rounded to a 64-bit int fall through
3165         // to '1 <= q + exp <= 20'
3166       } else if (q == 20) {
3167         // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
3168         C.w[0] = C1.w[0] + C1.w[0];
3169         C.w[1] = C1.w[1] + C1.w[1];
3170         if (C.w[0] < C1.w[0])
3171           C.w[1]++;
3172         if (C.w[1] > 0x01 || (C.w[1] == 0x01
3173                               && C.w[0] >= 0xffffffffffffffffull)) {
3174           // set invalid flag
3175           *pfpsf |= INVALID_EXCEPTION;
3176           // return Integer Indefinite
3177           res = 0x8000000000000000ull;
3178           BID_RETURN (res);
3179         }
3180         // else cases that can be rounded to a 64-bit int fall through
3181         // to '1 <= q + exp <= 20'
3182       } else if (q == 21) {
3183         // C >= 0x9fffffffffffffffb
3184         if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
3185                                && C1.w[0] >= 0xfffffffffffffffbull)) {
3186           // set invalid flag
3187           *pfpsf |= INVALID_EXCEPTION;
3188           // return Integer Indefinite
3189           res = 0x8000000000000000ull;
3190           BID_RETURN (res);
3191         }
3192         // else cases that can be rounded to a 64-bit int fall through
3193         // to '1 <= q + exp <= 20'
3194       } else {  // if 22 <= q <= 34 => 1 <= q - 21 <= 13
3195         // C  >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
3196         C.w[1] = 0x09;
3197         C.w[0] = 0xfffffffffffffffbull;
3198         __mul_128x64_to_128 (C, ten2k64[q - 21], C);
3199         if (C1.w[1] > C.w[1]
3200             || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3201           // set invalid flag
3202           *pfpsf |= INVALID_EXCEPTION;
3203           // return Integer Indefinite
3204           res = 0x8000000000000000ull;
3205           BID_RETURN (res);
3206         }
3207         // else cases that can be rounded to a 64-bit int fall through
3208         // to '1 <= q + exp <= 20'
3209       }
3210     }
3211   }
3212   // n is not too large to be converted to int64 if -1/2 < n < 2^64 - 1/2
3213   // Note: some of the cases tested for above fall through to this point
3214   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
3215     // set inexact flag
3216     *pfpsf |= INEXACT_EXCEPTION;
3217     // return 0
3218     res = 0x0000000000000000ull;
3219     BID_RETURN (res);
3220   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
3221     // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3222     //   res = 0
3223     // else if x > 0
3224     //   res = +1
3225     // else // if x < 0
3226     //   invalid exc
3227     ind = q - 1;
3228     if (ind <= 18) {    // 0 <= ind <= 18
3229       if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
3230         res = 0x0000000000000000ull;    // return 0
3231       } else if (!x_sign) {     // n > 0
3232         res = 0x00000001;       // return +1
3233       } else {
3234         res = 0x8000000000000000ull;
3235         // set invalid flag
3236         *pfpsf |= INVALID_EXCEPTION;
3237         // return Integer Indefinite
3238         res = 0x8000000000000000ull;
3239         BID_RETURN (res);
3240       }
3241     } else {    // 19 <= ind <= 33
3242       if ((C1.w[1] < midpoint128[ind - 19].w[1])
3243           || ((C1.w[1] == midpoint128[ind - 19].w[1])
3244               && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
3245         res = 0x0000000000000000ull;    // return 0
3246       } else if (!x_sign) {     // n > 0
3247         res = 0x00000001;       // return +1
3248       } else {
3249         res = 0x8000000000000000ull;
3250         *pfpsf |= INVALID_EXCEPTION;
3251         // return Integer Indefinite
3252         res = 0x8000000000000000ull;
3253         BID_RETURN (res);
3254       }
3255     }
3256     // set inexact flag
3257     *pfpsf |= INEXACT_EXCEPTION;
3258   } else {      // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
3259     // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
3260     // to nearest to a 64-bit unsigned signed integer
3261     if (x_sign) {       // x <= -1
3262       // set invalid flag
3263       *pfpsf |= INVALID_EXCEPTION;
3264       // return Integer Indefinite
3265       res = 0x8000000000000000ull;
3266       BID_RETURN (res);
3267     }
3268     // 1 <= x < 2^64-1/2 so x can be rounded
3269     // to nearest to a 64-bit unsigned integer
3270     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
3271       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
3272       // chop off ind digits from the lower part of C1
3273       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3274       tmp64 = C1.w[0];
3275       if (ind <= 19) {
3276         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
3277       } else {
3278         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
3279         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
3280       }
3281       if (C1.w[0] < tmp64)
3282         C1.w[1]++;
3283       // calculate C* and f*
3284       // C* is actually floor(C*) in this case
3285       // C* and f* need shifting and masking, as shown by
3286       // shiftright128[] and maskhigh128[]
3287       // 1 <= x <= 33
3288       // kx = 10^(-x) = ten2mk128[ind - 1]
3289       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3290       // the approximation of 10^(-x) was rounded up to 118 bits
3291       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
3292       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
3293         Cstar.w[1] = P256.w[3];
3294         Cstar.w[0] = P256.w[2];
3295         fstar.w[3] = 0;
3296         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
3297         fstar.w[1] = P256.w[1];
3298         fstar.w[0] = P256.w[0];
3299       } else {  // 22 <= ind - 1 <= 33
3300         Cstar.w[1] = 0;
3301         Cstar.w[0] = P256.w[3];
3302         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
3303         fstar.w[2] = P256.w[2];
3304         fstar.w[1] = P256.w[1];
3305         fstar.w[0] = P256.w[0];
3306       }
3307       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3308       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3309       // if (0 < f* < 10^(-x)) then the result is a midpoint
3310       //   if floor(C*) is even then C* = floor(C*) - logical right
3311       //       shift; C* has p decimal digits, correct by Prop. 1)
3312       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
3313       //       shift; C* has p decimal digits, correct by Pr. 1)
3314       // else
3315       //   C* = floor(C*) (logical right shift; C has p decimal digits,
3316       //       correct by Property 1)
3317       // n = C* * 10^(e+x)
3318
3319       // shift right C* by Ex-128 = shiftright128[ind]
3320       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
3321       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
3322         Cstar.w[0] =
3323           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
3324         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3325       } else {  // 22 <= ind - 1 <= 33
3326         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
3327       }
3328       // determine inexactness of the rounding of C*
3329       // if (0 < f* - 1/2 < 10^(-x)) then
3330       //   the result is exact
3331       // else // if (f* - 1/2 > T*) then
3332       //   the result is inexact
3333       if (ind - 1 <= 2) {
3334         if (fstar.w[1] > 0x8000000000000000ull ||
3335             (fstar.w[1] == 0x8000000000000000ull
3336              && fstar.w[0] > 0x0ull)) {
3337           // f* > 1/2 and the result may be exact
3338           tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
3339           if (tmp64 > ten2mk128trunc[ind - 1].w[1]
3340               || (tmp64 == ten2mk128trunc[ind - 1].w[1]
3341                   && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
3342             // set the inexact flag
3343             *pfpsf |= INEXACT_EXCEPTION;
3344           }     // else the result is exact
3345         } else {        // the result is inexact; f2* <= 1/2
3346           // set the inexact flag
3347           *pfpsf |= INEXACT_EXCEPTION;
3348         }
3349       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
3350         if (fstar.w[3] > 0x0 ||
3351             (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
3352             (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
3353              (fstar.w[1] || fstar.w[0]))) {
3354           // f2* > 1/2 and the result may be exact
3355           // Calculate f2* - 1/2
3356           tmp64 = fstar.w[2] - onehalf128[ind - 1];
3357           tmp64A = fstar.w[3];
3358           if (tmp64 > fstar.w[2])
3359             tmp64A--;
3360           if (tmp64A || tmp64
3361               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
3362               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
3363                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
3364             // set the inexact flag
3365             *pfpsf |= INEXACT_EXCEPTION;
3366           }     // else the result is exact
3367         } else {        // the result is inexact; f2* <= 1/2
3368           // set the inexact flag
3369           *pfpsf |= INEXACT_EXCEPTION;
3370         }
3371       } else {  // if 22 <= ind <= 33
3372         if (fstar.w[3] > onehalf128[ind - 1] ||
3373             (fstar.w[3] == onehalf128[ind - 1] &&
3374              (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
3375           // f2* > 1/2 and the result may be exact
3376           // Calculate f2* - 1/2
3377           tmp64 = fstar.w[3] - onehalf128[ind - 1];
3378           if (tmp64 || fstar.w[2]
3379               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
3380               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
3381                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
3382             // set the inexact flag
3383             *pfpsf |= INEXACT_EXCEPTION;
3384           }     // else the result is exact
3385         } else {        // the result is inexact; f2* <= 1/2
3386           // set the inexact flag
3387           *pfpsf |= INEXACT_EXCEPTION;
3388         }
3389       }
3390
3391       // if the result was a midpoint it was rounded away from zero
3392       res = Cstar.w[0]; // the result is positive
3393     } else if (exp == 0) {
3394       // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
3395       // res = C (exact)
3396       res = C1.w[0];
3397     } else {
3398       // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
3399       // res = C * 10^exp (exact) - must fit in 64 bits
3400       res = C1.w[0] * ten2k64[exp];
3401     }
3402   }
3403 }
3404
3405 BID_RETURN (res);
3406 }