Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

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