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_int64_rnint
  28.  ****************************************************************************/
  29.  
  30. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (SINT64, bid128_to_int64_rnint,
  31.                                           x)
  32.  
  33.      SINT64 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.   if ((q + exp) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
  131.     // set invalid flag
  132.     *pfpsf |= INVALID_EXCEPTION;
  133.     // return Integer Indefinite
  134.     res = 0x8000000000000000ull;
  135.     BID_RETURN (res);
  136.   } else if ((q + exp) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
  137.     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
  138.     // so x rounded to an integer may or may not fit in a signed 64-bit int
  139.     // the cases that do not fit are identified here; the ones that fit
  140.     // fall through and will be handled with other cases further,
  141.     // under '1 <= q + exp <= 19'
  142.     if (x_sign) {       // if n < 0 and q + exp = 19
  143.       // if n < -2^63 - 1/2 then n is too large
  144.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63+1/2
  145.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 5*(2^64+1), 1<=q<=34
  146.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 0x50000000000000005, 1<=q<=34
  147.       C.w[1] = 0x0000000000000005ull;
  148.       C.w[0] = 0000000000000005ull;
  149.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  150.         // 10^(20-q) is 64-bit, and so is C1
  151.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  152.       } else if (q == 20) {
  153.         ;       // C1 * 10^0 = C1
  154.       } else {  // if 21 <= q <= 34
  155.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  156.       }
  157.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
  158.         // set invalid flag
  159.         *pfpsf |= INVALID_EXCEPTION;
  160.         // return Integer Indefinite
  161.         res = 0x8000000000000000ull;
  162.         BID_RETURN (res);
  163.       }
  164.       // else cases that can be rounded to a 64-bit int fall through
  165.       // to '1 <= q + exp <= 19'
  166.     } else {    // if n > 0 and q + exp = 19
  167.       // if n >= 2^63 - 1/2 then n is too large
  168.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63-1/2
  169.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64-1), 1<=q<=34
  170.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x4fffffffffffffffb, 1<=q<=34
  171.       C.w[1] = 0x0000000000000004ull;
  172.       C.w[0] = 0xfffffffffffffffbull;
  173.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  174.         // 10^(20-q) is 64-bit, and so is C1
  175.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  176.       } else if (q == 20) {
  177.         ;       // C1 * 10^0 = C1
  178.       } else {  // if 21 <= q <= 34
  179.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  180.       }
  181.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  182.         // set invalid flag
  183.         *pfpsf |= INVALID_EXCEPTION;
  184.         // return Integer Indefinite
  185.         res = 0x8000000000000000ull;
  186.         BID_RETURN (res);
  187.       }
  188.       // else cases that can be rounded to a 64-bit int fall through
  189.       // to '1 <= q + exp <= 19'
  190.     }
  191.   }
  192.   // n is not too large to be converted to int64: -2^63-1/2 <= n < 2^63-1/2
  193.   // Note: some of the cases tested for above fall through to this point
  194.   // Restore C1 which may have been modified above
  195.   C1.w[1] = x.w[1] & MASK_COEFF;
  196.   C1.w[0] = x.w[0];
  197.   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
  198.     // return 0
  199.     res = 0x0000000000000000ull;
  200.     BID_RETURN (res);
  201.   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
  202.     // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
  203.     //   res = 0
  204.     // else
  205.     //   res = +/-1
  206.     ind = q - 1;
  207.     if (ind <= 18) {    // 0 <= ind <= 18
  208.       if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
  209.         res = 0x0000000000000000ull;    // return 0
  210.       } else if (x_sign) {      // n < 0
  211.         res = 0xffffffffffffffffull;    // return -1
  212.       } else {  // n > 0
  213.         res = 0x0000000000000001ull;    // return +1
  214.       }
  215.     } else {    // 19 <= ind <= 33
  216.       if ((C1.w[1] < midpoint128[ind - 19].w[1])
  217.           || ((C1.w[1] == midpoint128[ind - 19].w[1])
  218.               && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
  219.         res = 0x0000000000000000ull;    // return 0
  220.       } else if (x_sign) {      // n < 0
  221.         res = 0xffffffffffffffffull;    // return -1
  222.       } else {  // n > 0
  223.         res = 0x0000000000000001ull;    // return +1
  224.       }
  225.     }
  226.   } else {      // if (1 <= q + exp <= 19, 1 <= q <= 34, -33 <= exp <= 18)
  227.     // -2^63-1/2 <= x <= -1 or 1 <= x < 2^63-1/2 so x can be rounded
  228.     // to nearest to a 64-bit signed integer
  229.     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 19
  230.       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
  231.       // chop off ind digits from the lower part of C1
  232.       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  233.       tmp64 = C1.w[0];
  234.       if (ind <= 19) {
  235.         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  236.       } else {
  237.         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  238.         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  239.       }
  240.       if (C1.w[0] < tmp64)
  241.         C1.w[1]++;
  242.       // calculate C* and f*
  243.       // C* is actually floor(C*) in this case
  244.       // C* and f* need shifting and masking, as shown by
  245.       // shiftright128[] and maskhigh128[]
  246.       // 1 <= x <= 33
  247.       // kx = 10^(-x) = ten2mk128[ind - 1]
  248.       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  249.       // the approximation of 10^(-x) was rounded up to 118 bits
  250.       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  251.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  252.         Cstar.w[1] = P256.w[3];
  253.         Cstar.w[0] = P256.w[2];
  254.         fstar.w[3] = 0;
  255.         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  256.         fstar.w[1] = P256.w[1];
  257.         fstar.w[0] = P256.w[0];
  258.       } else {  // 22 <= ind - 1 <= 33
  259.         Cstar.w[1] = 0;
  260.         Cstar.w[0] = P256.w[3];
  261.         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  262.         fstar.w[2] = P256.w[2];
  263.         fstar.w[1] = P256.w[1];
  264.         fstar.w[0] = P256.w[0];
  265.       }
  266.       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  267.       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  268.       // if (0 < f* < 10^(-x)) then the result is a midpoint
  269.       //   if floor(C*) is even then C* = floor(C*) - logical right
  270.       //       shift; C* has p decimal digits, correct by Prop. 1)
  271.       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
  272.       //       shift; C* has p decimal digits, correct by Pr. 1)
  273.       // else
  274.       //   C* = floor(C*) (logical right shift; C has p decimal digits,
  275.       //       correct by Property 1)
  276.       // n = C* * 10^(e+x)
  277.  
  278.       // shift right C* by Ex-128 = shiftright128[ind]
  279.       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
  280.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  281.         Cstar.w[0] =
  282.           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  283.         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  284.       } else {  // 22 <= ind - 1 <= 33
  285.         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
  286.       }
  287.       // if the result was a midpoint it was rounded away from zero, so
  288.       // it will need a correction
  289.       // check for midpoints
  290.       if ((fstar.w[3] == 0) && (fstar.w[2] == 0) &&
  291.           (fstar.w[1] || fstar.w[0]) &&
  292.           (fstar.w[1] < ten2mk128trunc[ind - 1].w[1] ||
  293.            (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] &&
  294.             fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
  295.         // the result is a midpoint; round to nearest
  296.         if (Cstar.w[0] & 0x01) {        // Cstar.w[0] is odd; MP in [EVEN, ODD]
  297.           // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
  298.           Cstar.w[0]--; // Cstar.w[0] is now even
  299.         }       // else MP in [ODD, EVEN]
  300.       }
  301.       if (x_sign)
  302.         res = -Cstar.w[0];
  303.       else
  304.         res = Cstar.w[0];
  305.     } else if (exp == 0) {
  306.       // 1 <= q <= 19
  307.       // res = +/-C (exact)
  308.       if (x_sign)
  309.         res = -C1.w[0];
  310.       else
  311.         res = C1.w[0];
  312.     } else {    // if (exp>0) => 1 <= exp <= 18, 1 <= q < 18, 2 <= q + exp <= 19
  313.       // res = +/-C * 10^exp (exact) where this fits in 64-bit integer
  314.       if (x_sign)
  315.         res = -C1.w[0] * ten2k64[exp];
  316.       else
  317.         res = C1.w[0] * ten2k64[exp];
  318.     }
  319.   }
  320. }
  321.  
  322. BID_RETURN (res);
  323. }
  324.  
  325. /*****************************************************************************
  326.  *  BID128_to_int64_xrnint
  327.  ****************************************************************************/
  328.  
  329. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (SINT64,
  330.                                           bid128_to_int64_xrnint, x)
  331.  
  332.      SINT64 res;
  333.      UINT64 x_sign;
  334.      UINT64 x_exp;
  335.      int exp;                   // unbiased exponent
  336.   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  337.      UINT64 tmp64, tmp64A;
  338.      BID_UI64DOUBLE tmp1;
  339.      unsigned int x_nr_bits;
  340.      int q, ind, shift;
  341.      UINT128 C1, C;
  342.      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  343.      UINT256 fstar;
  344.      UINT256 P256;
  345.  
  346.   // unpack x
  347. x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  348. x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
  349. C1.w[1] = x.w[1] & MASK_COEFF;
  350. C1.w[0] = x.w[0];
  351.  
  352.   // check for NaN or Infinity
  353. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  354.     // x is special
  355. if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
  356.   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
  357.     // set invalid flag
  358.     *pfpsf |= INVALID_EXCEPTION;
  359.     // return Integer Indefinite
  360.     res = 0x8000000000000000ull;
  361.   } else {      // x is QNaN
  362.     // set invalid flag
  363.     *pfpsf |= INVALID_EXCEPTION;
  364.     // return Integer Indefinite
  365.     res = 0x8000000000000000ull;
  366.   }
  367.   BID_RETURN (res);
  368. } else {        // x is not a NaN, so it must be infinity
  369.   if (!x_sign) {        // x is +inf
  370.     // set invalid flag
  371.     *pfpsf |= INVALID_EXCEPTION;
  372.     // return Integer Indefinite
  373.     res = 0x8000000000000000ull;
  374.   } else {      // x is -inf
  375.     // set invalid flag
  376.     *pfpsf |= INVALID_EXCEPTION;
  377.     // return Integer Indefinite
  378.     res = 0x8000000000000000ull;
  379.   }
  380.   BID_RETURN (res);
  381. }
  382. }
  383.   // check for non-canonical values (after the check for special values)
  384. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  385.     || (C1.w[1] == 0x0001ed09bead87c0ull
  386.         && (C1.w[0] > 0x378d8e63ffffffffull))
  387.     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  388.   res = 0x0000000000000000ull;
  389.   BID_RETURN (res);
  390. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  391.   // x is 0
  392.   res = 0x0000000000000000ull;
  393.   BID_RETURN (res);
  394. } else {        // x is not special and is not zero
  395.  
  396.   // q = nr. of decimal digits in x
  397.   //  determine first the nr. of bits in x
  398.   if (C1.w[1] == 0) {
  399.     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
  400.       // split the 64-bit value in two 32-bit halves to avoid rounding errors
  401.       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
  402.         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
  403.         x_nr_bits =
  404.           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  405.       } else {  // x < 2^32
  406.         tmp1.d = (double) (C1.w[0]);    // exact conversion
  407.         x_nr_bits =
  408.           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  409.       }
  410.     } else {    // if x < 2^53
  411.       tmp1.d = (double) C1.w[0];        // exact conversion
  412.       x_nr_bits =
  413.         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  414.     }
  415.   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  416.     tmp1.d = (double) C1.w[1];  // exact conversion
  417.     x_nr_bits =
  418.       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  419.   }
  420.   q = nr_digits[x_nr_bits - 1].digits;
  421.   if (q == 0) {
  422.     q = nr_digits[x_nr_bits - 1].digits1;
  423.     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  424.         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  425.             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  426.       q++;
  427.   }
  428.   exp = (x_exp >> 49) - 6176;
  429.   if ((q + exp) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
  430.     // set invalid flag
  431.     *pfpsf |= INVALID_EXCEPTION;
  432.     // return Integer Indefinite
  433.     res = 0x8000000000000000ull;
  434.     BID_RETURN (res);
  435.   } else if ((q + exp) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
  436.     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
  437.     // so x rounded to an integer may or may not fit in a signed 64-bit int
  438.     // the cases that do not fit are identified here; the ones that fit
  439.     // fall through and will be handled with other cases further,
  440.     // under '1 <= q + exp <= 19'
  441.     if (x_sign) {       // if n < 0 and q + exp = 19
  442.       // if n < -2^63 - 1/2 then n is too large
  443.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63+1/2
  444.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 5*(2^64+1), 1<=q<=34
  445.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 0x50000000000000005, 1<=q<=34
  446.       C.w[1] = 0x0000000000000005ull;
  447.       C.w[0] = 0000000000000005ull;
  448.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  449.         // 10^(20-q) is 64-bit, and so is C1
  450.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  451.       } else if (q == 20) {
  452.         ;       // C1 * 10^0 = C1
  453.       } else {  // if 21 <= q <= 34
  454.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  455.       }
  456.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
  457.         // set invalid flag
  458.         *pfpsf |= INVALID_EXCEPTION;
  459.         // return Integer Indefinite
  460.         res = 0x8000000000000000ull;
  461.         BID_RETURN (res);
  462.       }
  463.       // else cases that can be rounded to a 64-bit int fall through
  464.       // to '1 <= q + exp <= 19'
  465.     } else {    // if n > 0 and q + exp = 19
  466.       // if n >= 2^63 - 1/2 then n is too large
  467.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63-1/2
  468.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64-1), 1<=q<=34
  469.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x4fffffffffffffffb, 1<=q<=34
  470.       C.w[1] = 0x0000000000000004ull;
  471.       C.w[0] = 0xfffffffffffffffbull;
  472.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  473.         // 10^(20-q) is 64-bit, and so is C1
  474.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  475.       } else if (q == 20) {
  476.         ;       // C1 * 10^0 = C1
  477.       } else {  // if 21 <= q <= 34
  478.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  479.       }
  480.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  481.         // set invalid flag
  482.         *pfpsf |= INVALID_EXCEPTION;
  483.         // return Integer Indefinite
  484.         res = 0x8000000000000000ull;
  485.         BID_RETURN (res);
  486.       }
  487.       // else cases that can be rounded to a 64-bit int fall through
  488.       // to '1 <= q + exp <= 19'
  489.     }
  490.   }
  491.   // n is not too large to be converted to int64: -2^63-1/2 <= n < 2^63-1/2
  492.   // Note: some of the cases tested for above fall through to this point
  493.   // Restore C1 which may have been modified above
  494.   C1.w[1] = x.w[1] & MASK_COEFF;
  495.   C1.w[0] = x.w[0];
  496.   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
  497.     // set inexact flag
  498.     *pfpsf |= INEXACT_EXCEPTION;
  499.     // return 0
  500.     res = 0x0000000000000000ull;
  501.     BID_RETURN (res);
  502.   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
  503.     // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
  504.     //   res = 0
  505.     // else
  506.     //   res = +/-1
  507.     ind = q - 1;
  508.     if (ind <= 18) {    // 0 <= ind <= 18
  509.       if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
  510.         res = 0x0000000000000000ull;    // return 0
  511.       } else if (x_sign) {      // n < 0
  512.         res = 0xffffffffffffffffull;    // return -1
  513.       } else {  // n > 0
  514.         res = 0x0000000000000001ull;    // return +1
  515.       }
  516.     } else {    // 19 <= ind <= 33
  517.       if ((C1.w[1] < midpoint128[ind - 19].w[1])
  518.           || ((C1.w[1] == midpoint128[ind - 19].w[1])
  519.               && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
  520.         res = 0x0000000000000000ull;    // return 0
  521.       } else if (x_sign) {      // n < 0
  522.         res = 0xffffffffffffffffull;    // return -1
  523.       } else {  // n > 0
  524.         res = 0x0000000000000001ull;    // return +1
  525.       }
  526.     }
  527.     // set inexact flag
  528.     *pfpsf |= INEXACT_EXCEPTION;
  529.   } else {      // if (1 <= q + exp <= 19, 1 <= q <= 34, -33 <= exp <= 18)
  530.     // -2^63-1/2 <= x <= -1 or 1 <= x < 2^63-1/2 so x can be rounded
  531.     // to nearest to a 64-bit signed integer
  532.     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 19
  533.       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
  534.       // chop off ind digits from the lower part of C1
  535.       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  536.       tmp64 = C1.w[0];
  537.       if (ind <= 19) {
  538.         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  539.       } else {
  540.         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  541.         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  542.       }
  543.       if (C1.w[0] < tmp64)
  544.         C1.w[1]++;
  545.       // calculate C* and f*
  546.       // C* is actually floor(C*) in this case
  547.       // C* and f* need shifting and masking, as shown by
  548.       // shiftright128[] and maskhigh128[]
  549.       // 1 <= x <= 33
  550.       // kx = 10^(-x) = ten2mk128[ind - 1]
  551.       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  552.       // the approximation of 10^(-x) was rounded up to 118 bits
  553.       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  554.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  555.         Cstar.w[1] = P256.w[3];
  556.         Cstar.w[0] = P256.w[2];
  557.         fstar.w[3] = 0;
  558.         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  559.         fstar.w[1] = P256.w[1];
  560.         fstar.w[0] = P256.w[0];
  561.       } else {  // 22 <= ind - 1 <= 33
  562.         Cstar.w[1] = 0;
  563.         Cstar.w[0] = P256.w[3];
  564.         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  565.         fstar.w[2] = P256.w[2];
  566.         fstar.w[1] = P256.w[1];
  567.         fstar.w[0] = P256.w[0];
  568.       }
  569.       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  570.       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  571.       // if (0 < f* < 10^(-x)) then the result is a midpoint
  572.       //   if floor(C*) is even then C* = floor(C*) - logical right
  573.       //       shift; C* has p decimal digits, correct by Prop. 1)
  574.       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
  575.       //       shift; C* has p decimal digits, correct by Pr. 1)
  576.       // else
  577.       //   C* = floor(C*) (logical right shift; C has p decimal digits,
  578.       //       correct by Property 1)
  579.       // n = C* * 10^(e+x)
  580.  
  581.       // shift right C* by Ex-128 = shiftright128[ind]
  582.       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
  583.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  584.         Cstar.w[0] =
  585.           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  586.         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  587.       } else {  // 22 <= ind - 1 <= 33
  588.         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
  589.       }
  590.       // determine inexactness of the rounding of C*
  591.       // if (0 < f* - 1/2 < 10^(-x)) then
  592.       //   the result is exact
  593.       // else // if (f* - 1/2 > T*) then
  594.       //   the result is inexact
  595.       if (ind - 1 <= 2) {
  596.         if (fstar.w[1] > 0x8000000000000000ull ||
  597.             (fstar.w[1] == 0x8000000000000000ull
  598.              && fstar.w[0] > 0x0ull)) {
  599.           // f* > 1/2 and the result may be exact
  600.           tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
  601.           if (tmp64 > ten2mk128trunc[ind - 1].w[1]
  602.               || (tmp64 == ten2mk128trunc[ind - 1].w[1]
  603.                   && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
  604.             // set the inexact flag
  605.             *pfpsf |= INEXACT_EXCEPTION;
  606.           }     // else the result is exact
  607.         } else {        // the result is inexact; f2* <= 1/2
  608.           // set the inexact flag
  609.           *pfpsf |= INEXACT_EXCEPTION;
  610.         }
  611.       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
  612.         if (fstar.w[3] > 0x0 ||
  613.             (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
  614.             (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
  615.              (fstar.w[1] || fstar.w[0]))) {
  616.           // f2* > 1/2 and the result may be exact
  617.           // Calculate f2* - 1/2
  618.           tmp64 = fstar.w[2] - onehalf128[ind - 1];
  619.           tmp64A = fstar.w[3];
  620.           if (tmp64 > fstar.w[2])
  621.             tmp64A--;
  622.           if (tmp64A || tmp64
  623.               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  624.               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  625.                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  626.             // set the inexact flag
  627.             *pfpsf |= INEXACT_EXCEPTION;
  628.           }     // else the result is exact
  629.         } else {        // the result is inexact; f2* <= 1/2
  630.           // set the inexact flag
  631.           *pfpsf |= INEXACT_EXCEPTION;
  632.         }
  633.       } else {  // if 22 <= ind <= 33
  634.         if (fstar.w[3] > onehalf128[ind - 1] ||
  635.             (fstar.w[3] == onehalf128[ind - 1] &&
  636.              (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
  637.           // f2* > 1/2 and the result may be exact
  638.           // Calculate f2* - 1/2
  639.           tmp64 = fstar.w[3] - onehalf128[ind - 1];
  640.           if (tmp64 || fstar.w[2]
  641.               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  642.               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  643.                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  644.             // set the inexact flag
  645.             *pfpsf |= INEXACT_EXCEPTION;
  646.           }     // else the result is exact
  647.         } else {        // the result is inexact; f2* <= 1/2
  648.           // set the inexact flag
  649.           *pfpsf |= INEXACT_EXCEPTION;
  650.         }
  651.       }
  652.  
  653.       // if the result was a midpoint it was rounded away from zero, so
  654.       // it will need a correction
  655.       // check for midpoints
  656.       if ((fstar.w[3] == 0) && (fstar.w[2] == 0) &&
  657.           (fstar.w[1] || fstar.w[0]) &&
  658.           (fstar.w[1] < ten2mk128trunc[ind - 1].w[1] ||
  659.            (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] &&
  660.             fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
  661.         // the result is a midpoint; round to nearest
  662.         if (Cstar.w[0] & 0x01) {        // Cstar.w[0] is odd; MP in [EVEN, ODD]
  663.           // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
  664.           Cstar.w[0]--; // Cstar.w[0] is now even
  665.         }       // else MP in [ODD, EVEN]
  666.       }
  667.       if (x_sign)
  668.         res = -Cstar.w[0];
  669.       else
  670.         res = Cstar.w[0];
  671.     } else if (exp == 0) {
  672.       // 1 <= q <= 19
  673.       // res = +/-C (exact)
  674.       if (x_sign)
  675.         res = -C1.w[0];
  676.       else
  677.         res = C1.w[0];
  678.     } else {    // if (exp>0) => 1 <= exp <= 18, 1 <= q < 18, 2 <= q + exp <= 19
  679.       // res = +/-C * 10^exp (exact) where this fits in 64-bit integer
  680.       if (x_sign)
  681.         res = -C1.w[0] * ten2k64[exp];
  682.       else
  683.         res = C1.w[0] * ten2k64[exp];
  684.     }
  685.   }
  686. }
  687.  
  688. BID_RETURN (res);
  689. }
  690.  
  691. /*****************************************************************************
  692.  *  BID128_to_int64_floor
  693.  ****************************************************************************/
  694.  
  695. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (SINT64, bid128_to_int64_floor,
  696.                                           x)
  697.  
  698.      SINT64 res;
  699.      UINT64 x_sign;
  700.      UINT64 x_exp;
  701.      int exp;                   // unbiased exponent
  702.   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  703.      BID_UI64DOUBLE tmp1;
  704.      unsigned int x_nr_bits;
  705.      int q, ind, shift;
  706.      UINT128 C1, C;
  707.      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  708.      UINT256 fstar;
  709.      UINT256 P256;
  710.  
  711.   // unpack x
  712. x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  713. x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
  714. C1.w[1] = x.w[1] & MASK_COEFF;
  715. C1.w[0] = x.w[0];
  716.  
  717.   // check for NaN or Infinity
  718. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  719.     // x is special
  720. if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
  721.   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
  722.     // set invalid flag
  723.     *pfpsf |= INVALID_EXCEPTION;
  724.     // return Integer Indefinite
  725.     res = 0x8000000000000000ull;
  726.   } else {      // x is QNaN
  727.     // set invalid flag
  728.     *pfpsf |= INVALID_EXCEPTION;
  729.     // return Integer Indefinite
  730.     res = 0x8000000000000000ull;
  731.   }
  732.   BID_RETURN (res);
  733. } else {        // x is not a NaN, so it must be infinity
  734.   if (!x_sign) {        // x is +inf
  735.     // set invalid flag
  736.     *pfpsf |= INVALID_EXCEPTION;
  737.     // return Integer Indefinite
  738.     res = 0x8000000000000000ull;
  739.   } else {      // x is -inf
  740.     // set invalid flag
  741.     *pfpsf |= INVALID_EXCEPTION;
  742.     // return Integer Indefinite
  743.     res = 0x8000000000000000ull;
  744.   }
  745.   BID_RETURN (res);
  746. }
  747. }
  748.   // check for non-canonical values (after the check for special values)
  749. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  750.     || (C1.w[1] == 0x0001ed09bead87c0ull
  751.         && (C1.w[0] > 0x378d8e63ffffffffull))
  752.     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  753.   res = 0x0000000000000000ull;
  754.   BID_RETURN (res);
  755. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  756.   // x is 0
  757.   res = 0x0000000000000000ull;
  758.   BID_RETURN (res);
  759. } else {        // x is not special and is not zero
  760.  
  761.   // q = nr. of decimal digits in x
  762.   //  determine first the nr. of bits in x
  763.   if (C1.w[1] == 0) {
  764.     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
  765.       // split the 64-bit value in two 32-bit halves to avoid rounding errors
  766.       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
  767.         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
  768.         x_nr_bits =
  769.           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  770.       } else {  // x < 2^32
  771.         tmp1.d = (double) (C1.w[0]);    // exact conversion
  772.         x_nr_bits =
  773.           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  774.       }
  775.     } else {    // if x < 2^53
  776.       tmp1.d = (double) C1.w[0];        // exact conversion
  777.       x_nr_bits =
  778.         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  779.     }
  780.   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  781.     tmp1.d = (double) C1.w[1];  // exact conversion
  782.     x_nr_bits =
  783.       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  784.   }
  785.   q = nr_digits[x_nr_bits - 1].digits;
  786.   if (q == 0) {
  787.     q = nr_digits[x_nr_bits - 1].digits1;
  788.     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  789.         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  790.             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  791.       q++;
  792.   }
  793.   exp = (x_exp >> 49) - 6176;
  794.  
  795.   if ((q + exp) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
  796.     // set invalid flag
  797.     *pfpsf |= INVALID_EXCEPTION;
  798.     // return Integer Indefinite
  799.     res = 0x8000000000000000ull;
  800.     BID_RETURN (res);
  801.   } else if ((q + exp) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
  802.     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
  803.     // so x rounded to an integer may or may not fit in a signed 64-bit int
  804.     // the cases that do not fit are identified here; the ones that fit
  805.     // fall through and will be handled with other cases further,
  806.     // under '1 <= q + exp <= 19'
  807.     if (x_sign) {       // if n < 0 and q + exp = 19
  808.       // if n < -2^63 then n is too large
  809.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63
  810.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 10*2^63, 1<=q<=34
  811.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 0x50000000000000000, 1<=q<=34
  812.       C.w[1] = 0x0000000000000005ull;
  813.       C.w[0] = 0x0000000000000000ull;
  814.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  815.         // 10^(20-q) is 64-bit, and so is C1
  816.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  817.       } else if (q == 20) {
  818.         ;       // C1 * 10^0 = C1
  819.       } else {  // if 21 <= q <= 34
  820.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  821.       }
  822.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
  823.         // set invalid flag
  824.         *pfpsf |= INVALID_EXCEPTION;
  825.         // return Integer Indefinite
  826.         res = 0x8000000000000000ull;
  827.         BID_RETURN (res);
  828.       }
  829.       // else cases that can be rounded to a 64-bit int fall through
  830.       // to '1 <= q + exp <= 19'
  831.     } else {    // if n > 0 and q + exp = 19
  832.       // if n >= 2^63 then n is too large
  833.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63
  834.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*2^64, 1<=q<=34
  835.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000000, 1<=q<=34
  836.       C.w[1] = 0x0000000000000005ull;
  837.       C.w[0] = 0x0000000000000000ull;
  838.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  839.         // 10^(20-q) is 64-bit, and so is C1
  840.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  841.       } else if (q == 20) {
  842.         ;       // C1 * 10^0 = C1
  843.       } else {  // if 21 <= q <= 34
  844.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  845.       }
  846.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  847.         // set invalid flag
  848.         *pfpsf |= INVALID_EXCEPTION;
  849.         // return Integer Indefinite
  850.         res = 0x8000000000000000ull;
  851.         BID_RETURN (res);
  852.       }
  853.       // else cases that can be rounded to a 64-bit int fall through
  854.       // to '1 <= q + exp <= 19'
  855.     }
  856.   }
  857.   // n is not too large to be converted to int64: -2^63-1 < n < 2^63
  858.   // Note: some of the cases tested for above fall through to this point
  859.   // Restore C1 which may have been modified above
  860.   C1.w[1] = x.w[1] & MASK_COEFF;
  861.   C1.w[0] = x.w[0];
  862.   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
  863.     // return -1 or 0
  864.     if (x_sign)
  865.       res = 0xffffffffffffffffull;
  866.     else
  867.       res = 0x0000000000000000ull;
  868.     BID_RETURN (res);
  869.   } else {      // if (1 <= q + exp <= 19, 1 <= q <= 34, -33 <= exp <= 18)
  870.     // -2^63 <= x <= -1 or 1 <= x < 2^63 so x can be rounded
  871.     // toward zero to a 64-bit signed integer
  872.     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 19
  873.       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
  874.       // chop off ind digits from the lower part of C1
  875.       // C1 fits in 127 bits
  876.       // calculate C* and f*
  877.       // C* is actually floor(C*) in this case
  878.       // C* and f* need shifting and masking, as shown by
  879.       // shiftright128[] and maskhigh128[]
  880.       // 1 <= x <= 33
  881.       // kx = 10^(-x) = ten2mk128[ind - 1]
  882.       // C* = C1 * 10^(-x)
  883.       // the approximation of 10^(-x) was rounded up to 118 bits
  884.       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  885.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  886.         Cstar.w[1] = P256.w[3];
  887.         Cstar.w[0] = P256.w[2];
  888.         fstar.w[3] = 0;
  889.         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  890.         fstar.w[1] = P256.w[1];
  891.         fstar.w[0] = P256.w[0];
  892.       } else {  // 22 <= ind - 1 <= 33
  893.         Cstar.w[1] = 0;
  894.         Cstar.w[0] = P256.w[3];
  895.         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  896.         fstar.w[2] = P256.w[2];
  897.         fstar.w[1] = P256.w[1];
  898.         fstar.w[0] = P256.w[0];
  899.       }
  900.       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  901.       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  902.       // C* = floor(C*) (logical right shift; C has p decimal digits,
  903.       //     correct by Property 1)
  904.       // n = C* * 10^(e+x)
  905.  
  906.       // shift right C* by Ex-128 = shiftright128[ind]
  907.       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
  908.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  909.         Cstar.w[0] =
  910.           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  911.         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  912.       } else {  // 22 <= ind - 1 <= 33
  913.         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
  914.       }
  915.       // if the result is negative and inexact, need to add 1 to it
  916.  
  917.       // determine inexactness of the rounding of C*
  918.       // if (0 < f* < 10^(-x)) then
  919.       //   the result is exact
  920.       // else // if (f* > T*) then
  921.       //   the result is inexact
  922.       if (ind - 1 <= 2) {
  923.         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1] ||
  924.             (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] &&
  925.              fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  926.           if (x_sign) { // positive and inexact
  927.             Cstar.w[0]++;
  928.             if (Cstar.w[0] == 0x0)
  929.               Cstar.w[1]++;
  930.           }
  931.         }       // else the result is exact
  932.       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
  933.         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  934.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  935.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  936.           if (x_sign) { // positive and inexact
  937.             Cstar.w[0]++;
  938.             if (Cstar.w[0] == 0x0)
  939.               Cstar.w[1]++;
  940.           }
  941.         }       // else the result is exact
  942.       } else {  // if 22 <= ind <= 33
  943.         if (fstar.w[3] || fstar.w[2]
  944.             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  945.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  946.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  947.           if (x_sign) { // positive and inexact
  948.             Cstar.w[0]++;
  949.             if (Cstar.w[0] == 0x0)
  950.               Cstar.w[1]++;
  951.           }
  952.         }       // else the result is exact
  953.       }
  954.  
  955.       if (x_sign)
  956.         res = -Cstar.w[0];
  957.       else
  958.         res = Cstar.w[0];
  959.     } else if (exp == 0) {
  960.       // 1 <= q <= 19
  961.       // res = +/-C (exact)
  962.       if (x_sign)
  963.         res = -C1.w[0];
  964.       else
  965.         res = C1.w[0];
  966.     } else {    // if (exp>0) => 1 <= exp <= 18, 1 <= q < 18, 2 <= q + exp <= 19
  967.       // res = +/-C * 10^exp (exact) where this fits in 64-bit integer
  968.       if (x_sign)
  969.         res = -C1.w[0] * ten2k64[exp];
  970.       else
  971.         res = C1.w[0] * ten2k64[exp];
  972.     }
  973.   }
  974. }
  975.  
  976. BID_RETURN (res);
  977. }
  978.  
  979. /*****************************************************************************
  980.  *  BID128_to_int64_xfloor
  981.  ****************************************************************************/
  982.  
  983. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (SINT64,
  984.                                           bid128_to_int64_xfloor, x)
  985.  
  986.      SINT64 res;
  987.      UINT64 x_sign;
  988.      UINT64 x_exp;
  989.      int exp;                   // unbiased exponent
  990.   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  991.      BID_UI64DOUBLE tmp1;
  992.      unsigned int x_nr_bits;
  993.      int q, ind, shift;
  994.      UINT128 C1, C;
  995.      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  996.      UINT256 fstar;
  997.      UINT256 P256;
  998.  
  999.   // unpack x
  1000. x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  1001. x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
  1002. C1.w[1] = x.w[1] & MASK_COEFF;
  1003. C1.w[0] = x.w[0];
  1004.  
  1005.   // check for NaN or Infinity
  1006. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  1007.     // x is special
  1008. if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
  1009.   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
  1010.     // set invalid flag
  1011.     *pfpsf |= INVALID_EXCEPTION;
  1012.     // return Integer Indefinite
  1013.     res = 0x8000000000000000ull;
  1014.   } else {      // x is QNaN
  1015.     // set invalid flag
  1016.     *pfpsf |= INVALID_EXCEPTION;
  1017.     // return Integer Indefinite
  1018.     res = 0x8000000000000000ull;
  1019.   }
  1020.   BID_RETURN (res);
  1021. } else {        // x is not a NaN, so it must be infinity
  1022.   if (!x_sign) {        // x is +inf
  1023.     // set invalid flag
  1024.     *pfpsf |= INVALID_EXCEPTION;
  1025.     // return Integer Indefinite
  1026.     res = 0x8000000000000000ull;
  1027.   } else {      // x is -inf
  1028.     // set invalid flag
  1029.     *pfpsf |= INVALID_EXCEPTION;
  1030.     // return Integer Indefinite
  1031.     res = 0x8000000000000000ull;
  1032.   }
  1033.   BID_RETURN (res);
  1034. }
  1035. }
  1036.   // check for non-canonical values (after the check for special values)
  1037. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  1038.     || (C1.w[1] == 0x0001ed09bead87c0ull
  1039.         && (C1.w[0] > 0x378d8e63ffffffffull))
  1040.     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  1041.   res = 0x0000000000000000ull;
  1042.   BID_RETURN (res);
  1043. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  1044.   // x is 0
  1045.   res = 0x0000000000000000ull;
  1046.   BID_RETURN (res);
  1047. } else {        // x is not special and is not zero
  1048.  
  1049.   // q = nr. of decimal digits in x
  1050.   //  determine first the nr. of bits in x
  1051.   if (C1.w[1] == 0) {
  1052.     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
  1053.       // split the 64-bit value in two 32-bit halves to avoid rounding errors
  1054.       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
  1055.         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
  1056.         x_nr_bits =
  1057.           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1058.       } else {  // x < 2^32
  1059.         tmp1.d = (double) (C1.w[0]);    // exact conversion
  1060.         x_nr_bits =
  1061.           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1062.       }
  1063.     } else {    // if x < 2^53
  1064.       tmp1.d = (double) C1.w[0];        // exact conversion
  1065.       x_nr_bits =
  1066.         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1067.     }
  1068.   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  1069.     tmp1.d = (double) C1.w[1];  // exact conversion
  1070.     x_nr_bits =
  1071.       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1072.   }
  1073.   q = nr_digits[x_nr_bits - 1].digits;
  1074.   if (q == 0) {
  1075.     q = nr_digits[x_nr_bits - 1].digits1;
  1076.     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  1077.         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  1078.             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  1079.       q++;
  1080.   }
  1081.   exp = (x_exp >> 49) - 6176;
  1082.   if ((q + exp) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
  1083.     // set invalid flag
  1084.     *pfpsf |= INVALID_EXCEPTION;
  1085.     // return Integer Indefinite
  1086.     res = 0x8000000000000000ull;
  1087.     BID_RETURN (res);
  1088.   } else if ((q + exp) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
  1089.     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
  1090.     // so x rounded to an integer may or may not fit in a signed 64-bit int
  1091.     // the cases that do not fit are identified here; the ones that fit
  1092.     // fall through and will be handled with other cases further,
  1093.     // under '1 <= q + exp <= 19'
  1094.     if (x_sign) {       // if n < 0 and q + exp = 19
  1095.       // if n < -2^63 then n is too large
  1096.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63
  1097.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 10*2^63, 1<=q<=34
  1098.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 0x50000000000000000, 1<=q<=34
  1099.       C.w[1] = 0x0000000000000005ull;
  1100.       C.w[0] = 0x0000000000000000ull;
  1101.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  1102.         // 10^(20-q) is 64-bit, and so is C1
  1103.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  1104.       } else if (q == 20) {
  1105.         ;       // C1 * 10^0 = C1
  1106.       } else {  // if 21 <= q <= 34
  1107.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  1108.       }
  1109.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
  1110.         // set invalid flag
  1111.         *pfpsf |= INVALID_EXCEPTION;
  1112.         // return Integer Indefinite
  1113.         res = 0x8000000000000000ull;
  1114.         BID_RETURN (res);
  1115.       }
  1116.       // else cases that can be rounded to a 64-bit int fall through
  1117.       // to '1 <= q + exp <= 19'
  1118.     } else {    // if n > 0 and q + exp = 19
  1119.       // if n >= 2^63 then n is too large
  1120.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63
  1121.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*2^64, 1<=q<=34
  1122.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000000, 1<=q<=34
  1123.       C.w[1] = 0x0000000000000005ull;
  1124.       C.w[0] = 0x0000000000000000ull;
  1125.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  1126.         // 10^(20-q) is 64-bit, and so is C1
  1127.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  1128.       } else if (q == 20) {
  1129.         ;       // C1 * 10^0 = C1
  1130.       } else {  // if 21 <= q <= 34
  1131.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  1132.       }
  1133.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  1134.         // set invalid flag
  1135.         *pfpsf |= INVALID_EXCEPTION;
  1136.         // return Integer Indefinite
  1137.         res = 0x8000000000000000ull;
  1138.         BID_RETURN (res);
  1139.       }
  1140.       // else cases that can be rounded to a 64-bit int fall through
  1141.       // to '1 <= q + exp <= 19'
  1142.     }
  1143.   }
  1144.   // n is not too large to be converted to int64: -2^63-1 < n < 2^63
  1145.   // Note: some of the cases tested for above fall through to this point
  1146.   // Restore C1 which may have been modified above
  1147.   C1.w[1] = x.w[1] & MASK_COEFF;
  1148.   C1.w[0] = x.w[0];
  1149.   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
  1150.     // set inexact flag
  1151.     *pfpsf |= INEXACT_EXCEPTION;
  1152.     // return -1 or 0
  1153.     if (x_sign)
  1154.       res = 0xffffffffffffffffull;
  1155.     else
  1156.       res = 0x0000000000000000ull;
  1157.     BID_RETURN (res);
  1158.   } else {      // if (1 <= q + exp <= 19, 1 <= q <= 34, -33 <= exp <= 18)
  1159.     // -2^63 <= x <= -1 or 1 <= x < 2^63 so x can be rounded
  1160.     // toward zero to a 64-bit signed integer
  1161.     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 19
  1162.       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
  1163.       // chop off ind digits from the lower part of C1
  1164.       // C1 fits in 127 bits
  1165.       // calculate C* and f*
  1166.       // C* is actually floor(C*) in this case
  1167.       // C* and f* need shifting and masking, as shown by
  1168.       // shiftright128[] and maskhigh128[]
  1169.       // 1 <= x <= 33
  1170.       // kx = 10^(-x) = ten2mk128[ind - 1]
  1171.       // C* = C1 * 10^(-x)
  1172.       // the approximation of 10^(-x) was rounded up to 118 bits
  1173.       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  1174.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  1175.         Cstar.w[1] = P256.w[3];
  1176.         Cstar.w[0] = P256.w[2];
  1177.         fstar.w[3] = 0;
  1178.         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  1179.         fstar.w[1] = P256.w[1];
  1180.         fstar.w[0] = P256.w[0];
  1181.       } else {  // 22 <= ind - 1 <= 33
  1182.         Cstar.w[1] = 0;
  1183.         Cstar.w[0] = P256.w[3];
  1184.         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  1185.         fstar.w[2] = P256.w[2];
  1186.         fstar.w[1] = P256.w[1];
  1187.         fstar.w[0] = P256.w[0];
  1188.       }
  1189.       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  1190.       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  1191.       // C* = floor(C*) (logical right shift; C has p decimal digits,
  1192.       //     correct by Property 1)
  1193.       // n = C* * 10^(e+x)
  1194.  
  1195.       // shift right C* by Ex-128 = shiftright128[ind]
  1196.       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
  1197.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  1198.         Cstar.w[0] =
  1199.           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  1200.         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  1201.       } else {  // 22 <= ind - 1 <= 33
  1202.         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
  1203.       }
  1204.       // if the result is negative and inexact, need to add 1 to it
  1205.  
  1206.       // determine inexactness of the rounding of C*
  1207.       // if (0 < f* < 10^(-x)) then
  1208.       //   the result is exact
  1209.       // else // if (f* > T*) then
  1210.       //   the result is inexact
  1211.       if (ind - 1 <= 2) {
  1212.         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1] ||
  1213.             (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] &&
  1214.              fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1215.           if (x_sign) { // positive and inexact
  1216.             Cstar.w[0]++;
  1217.             if (Cstar.w[0] == 0x0)
  1218.               Cstar.w[1]++;
  1219.           }
  1220.           // set the inexact flag
  1221.           *pfpsf |= INEXACT_EXCEPTION;
  1222.         }       // else the result is exact
  1223.       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
  1224.         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1225.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1226.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1227.           if (x_sign) { // positive and inexact
  1228.             Cstar.w[0]++;
  1229.             if (Cstar.w[0] == 0x0)
  1230.               Cstar.w[1]++;
  1231.           }
  1232.           // set the inexact flag
  1233.           *pfpsf |= INEXACT_EXCEPTION;
  1234.         }       // else the result is exact
  1235.       } else {  // if 22 <= ind <= 33
  1236.         if (fstar.w[3] || fstar.w[2]
  1237.             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1238.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1239.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1240.           if (x_sign) { // positive and inexact
  1241.             Cstar.w[0]++;
  1242.             if (Cstar.w[0] == 0x0)
  1243.               Cstar.w[1]++;
  1244.           }
  1245.           // set the inexact flag
  1246.           *pfpsf |= INEXACT_EXCEPTION;
  1247.         }       // else the result is exact
  1248.       }
  1249.  
  1250.       if (x_sign)
  1251.         res = -Cstar.w[0];
  1252.       else
  1253.         res = Cstar.w[0];
  1254.     } else if (exp == 0) {
  1255.       // 1 <= q <= 19
  1256.       // res = +/-C (exact)
  1257.       if (x_sign)
  1258.         res = -C1.w[0];
  1259.       else
  1260.         res = C1.w[0];
  1261.     } else {    // if (exp>0) => 1 <= exp <= 18, 1 <= q < 18, 2 <= q + exp <= 19
  1262.       // res = +/-C * 10^exp (exact) where this fits in 64-bit integer
  1263.       if (x_sign)
  1264.         res = -C1.w[0] * ten2k64[exp];
  1265.       else
  1266.         res = C1.w[0] * ten2k64[exp];
  1267.     }
  1268.   }
  1269. }
  1270.  
  1271. BID_RETURN (res);
  1272. }
  1273.  
  1274. /*****************************************************************************
  1275.  *  BID128_to_int64_ceil
  1276.  ****************************************************************************/
  1277.  
  1278. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (SINT64, bid128_to_int64_ceil,
  1279.                                           x)
  1280.  
  1281.      SINT64 res;
  1282.      UINT64 x_sign;
  1283.      UINT64 x_exp;
  1284.      int exp;                   // unbiased exponent
  1285.   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  1286.      BID_UI64DOUBLE tmp1;
  1287.      unsigned int x_nr_bits;
  1288.      int q, ind, shift;
  1289.      UINT128 C1, C;
  1290.      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  1291.      UINT256 fstar;
  1292.      UINT256 P256;
  1293.  
  1294.   // unpack x
  1295. x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  1296. x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
  1297. C1.w[1] = x.w[1] & MASK_COEFF;
  1298. C1.w[0] = x.w[0];
  1299.  
  1300.   // check for NaN or Infinity
  1301. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  1302.     // x is special
  1303. if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
  1304.   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
  1305.     // set invalid flag
  1306.     *pfpsf |= INVALID_EXCEPTION;
  1307.     // return Integer Indefinite
  1308.     res = 0x8000000000000000ull;
  1309.   } else {      // x is QNaN
  1310.     // set invalid flag
  1311.     *pfpsf |= INVALID_EXCEPTION;
  1312.     // return Integer Indefinite
  1313.     res = 0x8000000000000000ull;
  1314.   }
  1315.   BID_RETURN (res);
  1316. } else {        // x is not a NaN, so it must be infinity
  1317.   if (!x_sign) {        // x is +inf
  1318.     // set invalid flag
  1319.     *pfpsf |= INVALID_EXCEPTION;
  1320.     // return Integer Indefinite
  1321.     res = 0x8000000000000000ull;
  1322.   } else {      // x is -inf
  1323.     // set invalid flag
  1324.     *pfpsf |= INVALID_EXCEPTION;
  1325.     // return Integer Indefinite
  1326.     res = 0x8000000000000000ull;
  1327.   }
  1328.   BID_RETURN (res);
  1329. }
  1330. }
  1331.   // check for non-canonical values (after the check for special values)
  1332. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  1333.     || (C1.w[1] == 0x0001ed09bead87c0ull
  1334.         && (C1.w[0] > 0x378d8e63ffffffffull))
  1335.     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  1336.   res = 0x0000000000000000ull;
  1337.   BID_RETURN (res);
  1338. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  1339.   // x is 0
  1340.   res = 0x0000000000000000ull;
  1341.   BID_RETURN (res);
  1342. } else {        // x is not special and is not zero
  1343.  
  1344.   // q = nr. of decimal digits in x
  1345.   //  determine first the nr. of bits in x
  1346.   if (C1.w[1] == 0) {
  1347.     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
  1348.       // split the 64-bit value in two 32-bit halves to avoid rounding errors
  1349.       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
  1350.         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
  1351.         x_nr_bits =
  1352.           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1353.       } else {  // x < 2^32
  1354.         tmp1.d = (double) (C1.w[0]);    // exact conversion
  1355.         x_nr_bits =
  1356.           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1357.       }
  1358.     } else {    // if x < 2^53
  1359.       tmp1.d = (double) C1.w[0];        // exact conversion
  1360.       x_nr_bits =
  1361.         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1362.     }
  1363.   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  1364.     tmp1.d = (double) C1.w[1];  // exact conversion
  1365.     x_nr_bits =
  1366.       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1367.   }
  1368.   q = nr_digits[x_nr_bits - 1].digits;
  1369.   if (q == 0) {
  1370.     q = nr_digits[x_nr_bits - 1].digits1;
  1371.     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  1372.         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  1373.             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  1374.       q++;
  1375.   }
  1376.   exp = (x_exp >> 49) - 6176;
  1377.   if ((q + exp) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
  1378.     // set invalid flag
  1379.     *pfpsf |= INVALID_EXCEPTION;
  1380.     // return Integer Indefinite
  1381.     res = 0x8000000000000000ull;
  1382.     BID_RETURN (res);
  1383.   } else if ((q + exp) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
  1384.     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
  1385.     // so x rounded to an integer may or may not fit in a signed 64-bit int
  1386.     // the cases that do not fit are identified here; the ones that fit
  1387.     // fall through and will be handled with other cases further,
  1388.     // under '1 <= q + exp <= 19'
  1389.     if (x_sign) {       // if n < 0 and q + exp = 19
  1390.       // if n <= -2^63 - 1 then n is too large
  1391.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1
  1392.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 5*(2^64+2), 1<=q<=34
  1393.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 0x5000000000000000a, 1<=q<=34
  1394.       C.w[1] = 0x0000000000000005ull;
  1395.       C.w[0] = 0x000000000000000aull;
  1396.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  1397.         // 10^(20-q) is 64-bit, and so is C1
  1398.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  1399.       } else if (q == 20) {
  1400.         ;       // C1 * 10^0 = C1
  1401.       } else {  // if 21 <= q <= 34
  1402.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  1403.       }
  1404.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  1405.         // set invalid flag
  1406.         *pfpsf |= INVALID_EXCEPTION;
  1407.         // return Integer Indefinite
  1408.         res = 0x8000000000000000ull;
  1409.         BID_RETURN (res);
  1410.       }
  1411.       // else cases that can be rounded to a 64-bit int fall through
  1412.       // to '1 <= q + exp <= 19'
  1413.     } else {    // if n > 0 and q + exp = 19
  1414.       // if n > 2^63 - 1 then n is too large
  1415.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63 - 1
  1416.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 > 10*(2^63-1), 1<=q<=34
  1417.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 > 0x4fffffffffffffff6, 1<=q<=34
  1418.       C.w[1] = 0x0000000000000004ull;
  1419.       C.w[0] = 0xfffffffffffffff6ull;
  1420.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  1421.         // 10^(20-q) is 64-bit, and so is C1
  1422.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  1423.       } else if (q == 20) {
  1424.         ;       // C1 * 10^0 = C1
  1425.       } else {  // if 21 <= q <= 34
  1426.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  1427.       }
  1428.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
  1429.         // set invalid flag
  1430.         *pfpsf |= INVALID_EXCEPTION;
  1431.         // return Integer Indefinite
  1432.         res = 0x8000000000000000ull;
  1433.         BID_RETURN (res);
  1434.       }
  1435.       // else cases that can be rounded to a 64-bit int fall through
  1436.       // to '1 <= q + exp <= 19'
  1437.     }
  1438.   }
  1439.   // n is not too large to be converted to int64: -2^63-1 < n <= 2^63 - 1
  1440.   // Note: some of the cases tested for above fall through to this point
  1441.   // Restore C1 which may have been modified above
  1442.   C1.w[1] = x.w[1] & MASK_COEFF;
  1443.   C1.w[0] = x.w[0];
  1444.   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
  1445.     // return 0 or 1
  1446.     if (x_sign)
  1447.       res = 0x0000000000000000ull;
  1448.     else
  1449.       res = 0x0000000000000001ull;
  1450.     BID_RETURN (res);
  1451.   } else {      // if (1 <= q + exp <= 19, 1 <= q <= 34, -33 <= exp <= 18)
  1452.     // -2^63-1 < x <= -1 or 1 <= x <= 2^63 - 1 so x can be rounded
  1453.     // up to a 64-bit signed integer
  1454.     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 19
  1455.       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
  1456.       // chop off ind digits from the lower part of C1
  1457.       // C1 fits in 127 bits
  1458.       // calculate C* and f*
  1459.       // C* is actually floor(C*) in this case
  1460.       // C* and f* need shifting and masking, as shown by
  1461.       // shiftright128[] and maskhigh128[]
  1462.       // 1 <= x <= 33
  1463.       // kx = 10^(-x) = ten2mk128[ind - 1]
  1464.       // C* = C1 * 10^(-x)
  1465.       // the approximation of 10^(-x) was rounded up to 118 bits
  1466.       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  1467.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  1468.         Cstar.w[1] = P256.w[3];
  1469.         Cstar.w[0] = P256.w[2];
  1470.         fstar.w[3] = 0;
  1471.         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  1472.         fstar.w[1] = P256.w[1];
  1473.         fstar.w[0] = P256.w[0];
  1474.       } else {  // 22 <= ind - 1 <= 33
  1475.         Cstar.w[1] = 0;
  1476.         Cstar.w[0] = P256.w[3];
  1477.         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  1478.         fstar.w[2] = P256.w[2];
  1479.         fstar.w[1] = P256.w[1];
  1480.         fstar.w[0] = P256.w[0];
  1481.       }
  1482.       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  1483.       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  1484.       // C* = floor(C*) (logical right shift; C has p decimal digits,
  1485.       //     correct by Property 1)
  1486.       // n = C* * 10^(e+x)
  1487.  
  1488.       // shift right C* by Ex-128 = shiftright128[ind]
  1489.       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
  1490.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  1491.         Cstar.w[0] =
  1492.           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  1493.         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  1494.       } else {  // 22 <= ind - 1 <= 33
  1495.         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
  1496.       }
  1497.       // if the result is positive and inexact, need to add 1 to it
  1498.  
  1499.       // determine inexactness of the rounding of C*
  1500.       // if (0 < f* < 10^(-x)) then
  1501.       //   the result is exact
  1502.       // else // if (f* > T*) then
  1503.       //   the result is inexact
  1504.       if (ind - 1 <= 2) {
  1505.         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1506.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1507.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1508.           if (!x_sign) {        // positive and inexact
  1509.             Cstar.w[0]++;
  1510.             if (Cstar.w[0] == 0x0)
  1511.               Cstar.w[1]++;
  1512.           }
  1513.         }       // else the result is exact
  1514.       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
  1515.         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1516.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1517.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1518.           if (!x_sign) {        // positive and inexact
  1519.             Cstar.w[0]++;
  1520.             if (Cstar.w[0] == 0x0)
  1521.               Cstar.w[1]++;
  1522.           }
  1523.         }       // else the result is exact
  1524.       } else {  // if 22 <= ind <= 33
  1525.         if (fstar.w[3] || fstar.w[2]
  1526.             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1527.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1528.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1529.           if (!x_sign) {        // positive and inexact
  1530.             Cstar.w[0]++;
  1531.             if (Cstar.w[0] == 0x0)
  1532.               Cstar.w[1]++;
  1533.           }
  1534.         }       // else the result is exact
  1535.       }
  1536.       if (x_sign)
  1537.         res = -Cstar.w[0];
  1538.       else
  1539.         res = Cstar.w[0];
  1540.     } else if (exp == 0) {
  1541.       // 1 <= q <= 19
  1542.       // res = +/-C (exact)
  1543.       if (x_sign)
  1544.         res = -C1.w[0];
  1545.       else
  1546.         res = C1.w[0];
  1547.     } else {    // if (exp>0) => 1 <= exp <= 18, 1 <= q < 18, 2 <= q + exp <= 19
  1548.       // res = +/-C * 10^exp (exact) where this fits in 64-bit integer
  1549.       if (x_sign)
  1550.         res = -C1.w[0] * ten2k64[exp];
  1551.       else
  1552.         res = C1.w[0] * ten2k64[exp];
  1553.     }
  1554.   }
  1555. }
  1556.  
  1557. BID_RETURN (res);
  1558. }
  1559.  
  1560. /*****************************************************************************
  1561.  *  BID128_to_int64_xceil
  1562.  ****************************************************************************/
  1563.  
  1564. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (SINT64, bid128_to_int64_xceil,
  1565.                                           x)
  1566.  
  1567.      SINT64 res;
  1568.      UINT64 x_sign;
  1569.      UINT64 x_exp;
  1570.      int exp;                   // unbiased exponent
  1571.   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  1572.      BID_UI64DOUBLE tmp1;
  1573.      unsigned int x_nr_bits;
  1574.      int q, ind, shift;
  1575.      UINT128 C1, C;
  1576.      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  1577.      UINT256 fstar;
  1578.      UINT256 P256;
  1579.  
  1580.   // unpack x
  1581. x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  1582. x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
  1583. C1.w[1] = x.w[1] & MASK_COEFF;
  1584. C1.w[0] = x.w[0];
  1585.  
  1586.   // check for NaN or Infinity
  1587. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  1588.     // x is special
  1589. if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
  1590.   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
  1591.     // set invalid flag
  1592.     *pfpsf |= INVALID_EXCEPTION;
  1593.     // return Integer Indefinite
  1594.     res = 0x8000000000000000ull;
  1595.   } else {      // x is QNaN
  1596.     // set invalid flag
  1597.     *pfpsf |= INVALID_EXCEPTION;
  1598.     // return Integer Indefinite
  1599.     res = 0x8000000000000000ull;
  1600.   }
  1601.   BID_RETURN (res);
  1602. } else {        // x is not a NaN, so it must be infinity
  1603.   if (!x_sign) {        // x is +inf
  1604.     // set invalid flag
  1605.     *pfpsf |= INVALID_EXCEPTION;
  1606.     // return Integer Indefinite
  1607.     res = 0x8000000000000000ull;
  1608.   } else {      // x is -inf
  1609.     // set invalid flag
  1610.     *pfpsf |= INVALID_EXCEPTION;
  1611.     // return Integer Indefinite
  1612.     res = 0x8000000000000000ull;
  1613.   }
  1614.   BID_RETURN (res);
  1615. }
  1616. }
  1617.   // check for non-canonical values (after the check for special values)
  1618. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  1619.     || (C1.w[1] == 0x0001ed09bead87c0ull
  1620.         && (C1.w[0] > 0x378d8e63ffffffffull))
  1621.     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  1622.   res = 0x0000000000000000ull;
  1623.   BID_RETURN (res);
  1624. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  1625.   // x is 0
  1626.   res = 0x0000000000000000ull;
  1627.   BID_RETURN (res);
  1628. } else {        // x is not special and is not zero
  1629.  
  1630.   // q = nr. of decimal digits in x
  1631.   //  determine first the nr. of bits in x
  1632.   if (C1.w[1] == 0) {
  1633.     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
  1634.       // split the 64-bit value in two 32-bit halves to avoid rounding errors
  1635.       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
  1636.         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
  1637.         x_nr_bits =
  1638.           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1639.       } else {  // x < 2^32
  1640.         tmp1.d = (double) (C1.w[0]);    // exact conversion
  1641.         x_nr_bits =
  1642.           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1643.       }
  1644.     } else {    // if x < 2^53
  1645.       tmp1.d = (double) C1.w[0];        // exact conversion
  1646.       x_nr_bits =
  1647.         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1648.     }
  1649.   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  1650.     tmp1.d = (double) C1.w[1];  // exact conversion
  1651.     x_nr_bits =
  1652.       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1653.   }
  1654.   q = nr_digits[x_nr_bits - 1].digits;
  1655.   if (q == 0) {
  1656.     q = nr_digits[x_nr_bits - 1].digits1;
  1657.     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  1658.         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  1659.             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  1660.       q++;
  1661.   }
  1662.   exp = (x_exp >> 49) - 6176;
  1663.   if ((q + exp) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
  1664.     // set invalid flag
  1665.     *pfpsf |= INVALID_EXCEPTION;
  1666.     // return Integer Indefinite
  1667.     res = 0x8000000000000000ull;
  1668.     BID_RETURN (res);
  1669.   } else if ((q + exp) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
  1670.     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
  1671.     // so x rounded to an integer may or may not fit in a signed 64-bit int
  1672.     // the cases that do not fit are identified here; the ones that fit
  1673.     // fall through and will be handled with other cases further,
  1674.     // under '1 <= q + exp <= 19'
  1675.     if (x_sign) {       // if n < 0 and q + exp = 19
  1676.       // if n <= -2^63 - 1 then n is too large
  1677.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1
  1678.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 5*(2^64+2), 1<=q<=34
  1679.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 0x5000000000000000a, 1<=q<=34
  1680.       C.w[1] = 0x0000000000000005ull;
  1681.       C.w[0] = 0x000000000000000aull;
  1682.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  1683.         // 10^(20-q) is 64-bit, and so is C1
  1684.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  1685.       } else if (q == 20) {
  1686.         ;       // C1 * 10^0 = C1
  1687.       } else {  // if 21 <= q <= 34
  1688.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  1689.       }
  1690.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  1691.         // set invalid flag
  1692.         *pfpsf |= INVALID_EXCEPTION;
  1693.         // return Integer Indefinite
  1694.         res = 0x8000000000000000ull;
  1695.         BID_RETURN (res);
  1696.       }
  1697.       // else cases that can be rounded to a 64-bit int fall through
  1698.       // to '1 <= q + exp <= 19'
  1699.     } else {    // if n > 0 and q + exp = 19
  1700.       // if n > 2^63 - 1 then n is too large
  1701.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63 - 1
  1702.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 > 10*(2^63-1), 1<=q<=34
  1703.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 > 0x4fffffffffffffff6, 1<=q<=34
  1704.       C.w[1] = 0x0000000000000004ull;
  1705.       C.w[0] = 0xfffffffffffffff6ull;
  1706.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  1707.         // 10^(20-q) is 64-bit, and so is C1
  1708.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  1709.       } else if (q == 20) {
  1710.         ;       // C1 * 10^0 = C1
  1711.       } else {  // if 21 <= q <= 34
  1712.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  1713.       }
  1714.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
  1715.         // set invalid flag
  1716.         *pfpsf |= INVALID_EXCEPTION;
  1717.         // return Integer Indefinite
  1718.         res = 0x8000000000000000ull;
  1719.         BID_RETURN (res);
  1720.       }
  1721.       // else cases that can be rounded to a 64-bit int fall through
  1722.       // to '1 <= q + exp <= 19'
  1723.     }
  1724.   }
  1725.   // n is not too large to be converted to int64: -2^63-1 < n <= 2^63 - 1
  1726.   // Note: some of the cases tested for above fall through to this point
  1727.   // Restore C1 which may have been modified above
  1728.   C1.w[1] = x.w[1] & MASK_COEFF;
  1729.   C1.w[0] = x.w[0];
  1730.   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
  1731.     // set inexact flag
  1732.     *pfpsf |= INEXACT_EXCEPTION;
  1733.     // return 0 or 1
  1734.     if (x_sign)
  1735.       res = 0x0000000000000000ull;
  1736.     else
  1737.       res = 0x0000000000000001ull;
  1738.     BID_RETURN (res);
  1739.   } else {      // if (1 <= q + exp <= 19, 1 <= q <= 34, -33 <= exp <= 18)
  1740.     // -2^63-1 < x <= -1 or 1 <= x <= 2^63 - 1 so x can be rounded
  1741.     // up to a 64-bit signed integer
  1742.     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 19
  1743.       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
  1744.       // chop off ind digits from the lower part of C1
  1745.       // C1 fits in 127 bits
  1746.       // calculate C* and f*
  1747.       // C* is actually floor(C*) in this case
  1748.       // C* and f* need shifting and masking, as shown by
  1749.       // shiftright128[] and maskhigh128[]
  1750.       // 1 <= x <= 33
  1751.       // kx = 10^(-x) = ten2mk128[ind - 1]
  1752.       // C* = C1 * 10^(-x)
  1753.       // the approximation of 10^(-x) was rounded up to 118 bits
  1754.       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  1755.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  1756.         Cstar.w[1] = P256.w[3];
  1757.         Cstar.w[0] = P256.w[2];
  1758.         fstar.w[3] = 0;
  1759.         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  1760.         fstar.w[1] = P256.w[1];
  1761.         fstar.w[0] = P256.w[0];
  1762.       } else {  // 22 <= ind - 1 <= 33
  1763.         Cstar.w[1] = 0;
  1764.         Cstar.w[0] = P256.w[3];
  1765.         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  1766.         fstar.w[2] = P256.w[2];
  1767.         fstar.w[1] = P256.w[1];
  1768.         fstar.w[0] = P256.w[0];
  1769.       }
  1770.       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  1771.       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  1772.       // C* = floor(C*) (logical right shift; C has p decimal digits,
  1773.       //     correct by Property 1)
  1774.       // n = C* * 10^(e+x)
  1775.  
  1776.       // shift right C* by Ex-128 = shiftright128[ind]
  1777.       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
  1778.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  1779.         Cstar.w[0] =
  1780.           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  1781.         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  1782.       } else {  // 22 <= ind - 1 <= 33
  1783.         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
  1784.       }
  1785.       // if the result is positive and inexact, need to add 1 to it
  1786.  
  1787.       // determine inexactness of the rounding of C*
  1788.       // if (0 < f* < 10^(-x)) then
  1789.       //   the result is exact
  1790.       // else // if (f* > T*) then
  1791.       //   the result is inexact
  1792.       if (ind - 1 <= 2) {
  1793.         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1794.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1795.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1796.           if (!x_sign) {        // positive and inexact
  1797.             Cstar.w[0]++;
  1798.             if (Cstar.w[0] == 0x0)
  1799.               Cstar.w[1]++;
  1800.           }
  1801.           // set the inexact flag
  1802.           *pfpsf |= INEXACT_EXCEPTION;
  1803.         }       // else the result is exact
  1804.       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
  1805.         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1806.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1807.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1808.           if (!x_sign) {        // positive and inexact
  1809.             Cstar.w[0]++;
  1810.             if (Cstar.w[0] == 0x0)
  1811.               Cstar.w[1]++;
  1812.           }
  1813.           // set the inexact flag
  1814.           *pfpsf |= INEXACT_EXCEPTION;
  1815.         }       // else the result is exact
  1816.       } else {  // if 22 <= ind <= 33
  1817.         if (fstar.w[3] || fstar.w[2]
  1818.             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1819.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1820.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1821.           if (!x_sign) {        // positive and inexact
  1822.             Cstar.w[0]++;
  1823.             if (Cstar.w[0] == 0x0)
  1824.               Cstar.w[1]++;
  1825.           }
  1826.           // set the inexact flag
  1827.           *pfpsf |= INEXACT_EXCEPTION;
  1828.         }       // else the result is exact
  1829.       }
  1830.  
  1831.       if (x_sign)
  1832.         res = -Cstar.w[0];
  1833.       else
  1834.         res = Cstar.w[0];
  1835.     } else if (exp == 0) {
  1836.       // 1 <= q <= 19
  1837.       // res = +/-C (exact)
  1838.       if (x_sign)
  1839.         res = -C1.w[0];
  1840.       else
  1841.         res = C1.w[0];
  1842.     } else {    // if (exp>0) => 1 <= exp <= 18, 1 <= q < 18, 2 <= q + exp <= 19
  1843.       // res = +/-C * 10^exp (exact) where this fits in 64-bit integer
  1844.       if (x_sign)
  1845.         res = -C1.w[0] * ten2k64[exp];
  1846.       else
  1847.         res = C1.w[0] * ten2k64[exp];
  1848.     }
  1849.   }
  1850. }
  1851.  
  1852. BID_RETURN (res);
  1853. }
  1854.  
  1855. /*****************************************************************************
  1856.  *  BID128_to_int64_int
  1857.  ****************************************************************************/
  1858.  
  1859. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (SINT64, bid128_to_int64_int,
  1860.                                           x)
  1861.  
  1862.      SINT64 res;
  1863.      UINT64 x_sign;
  1864.      UINT64 x_exp;
  1865.      int exp;                   // unbiased exponent
  1866.   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  1867.      BID_UI64DOUBLE tmp1;
  1868.      unsigned int x_nr_bits;
  1869.      int q, ind, shift;
  1870.      UINT128 C1, C;
  1871.      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  1872.      UINT256 P256;
  1873.  
  1874.   // unpack x
  1875. x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  1876. x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
  1877. C1.w[1] = x.w[1] & MASK_COEFF;
  1878. C1.w[0] = x.w[0];
  1879.  
  1880.   // check for NaN or Infinity
  1881. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  1882.     // x is special
  1883. if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
  1884.   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
  1885.     // set invalid flag
  1886.     *pfpsf |= INVALID_EXCEPTION;
  1887.     // return Integer Indefinite
  1888.     res = 0x8000000000000000ull;
  1889.   } else {      // x is QNaN
  1890.     // set invalid flag
  1891.     *pfpsf |= INVALID_EXCEPTION;
  1892.     // return Integer Indefinite
  1893.     res = 0x8000000000000000ull;
  1894.   }
  1895.   BID_RETURN (res);
  1896. } else {        // x is not a NaN, so it must be infinity
  1897.   if (!x_sign) {        // x is +inf
  1898.     // set invalid flag
  1899.     *pfpsf |= INVALID_EXCEPTION;
  1900.     // return Integer Indefinite
  1901.     res = 0x8000000000000000ull;
  1902.   } else {      // x is -inf
  1903.     // set invalid flag
  1904.     *pfpsf |= INVALID_EXCEPTION;
  1905.     // return Integer Indefinite
  1906.     res = 0x8000000000000000ull;
  1907.   }
  1908.   BID_RETURN (res);
  1909. }
  1910. }
  1911.   // check for non-canonical values (after the check for special values)
  1912. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  1913.     || (C1.w[1] == 0x0001ed09bead87c0ull
  1914.         && (C1.w[0] > 0x378d8e63ffffffffull))
  1915.     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  1916.   res = 0x0000000000000000ull;
  1917.   BID_RETURN (res);
  1918. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  1919.   // x is 0
  1920.   res = 0x0000000000000000ull;
  1921.   BID_RETURN (res);
  1922. } else {        // x is not special and is not zero
  1923.  
  1924.   // q = nr. of decimal digits in x
  1925.   //  determine first the nr. of bits in x
  1926.   if (C1.w[1] == 0) {
  1927.     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
  1928.       // split the 64-bit value in two 32-bit halves to avoid rounding errors
  1929.       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
  1930.         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
  1931.         x_nr_bits =
  1932.           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1933.       } else {  // x < 2^32
  1934.         tmp1.d = (double) (C1.w[0]);    // exact conversion
  1935.         x_nr_bits =
  1936.           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1937.       }
  1938.     } else {    // if x < 2^53
  1939.       tmp1.d = (double) C1.w[0];        // exact conversion
  1940.       x_nr_bits =
  1941.         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1942.     }
  1943.   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  1944.     tmp1.d = (double) C1.w[1];  // exact conversion
  1945.     x_nr_bits =
  1946.       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1947.   }
  1948.   q = nr_digits[x_nr_bits - 1].digits;
  1949.   if (q == 0) {
  1950.     q = nr_digits[x_nr_bits - 1].digits1;
  1951.     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  1952.         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  1953.             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  1954.       q++;
  1955.   }
  1956.   exp = (x_exp >> 49) - 6176;
  1957.   if ((q + exp) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
  1958.     // set invalid flag
  1959.     *pfpsf |= INVALID_EXCEPTION;
  1960.     // return Integer Indefinite
  1961.     res = 0x8000000000000000ull;
  1962.     BID_RETURN (res);
  1963.   } else if ((q + exp) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
  1964.     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
  1965.     // so x rounded to an integer may or may not fit in a signed 64-bit int
  1966.     // the cases that do not fit are identified here; the ones that fit
  1967.     // fall through and will be handled with other cases further,
  1968.     // under '1 <= q + exp <= 19'
  1969.     if (x_sign) {       // if n < 0 and q + exp = 19
  1970.       // if n <= -2^63 - 1 then n is too large
  1971.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1
  1972.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64+2), 1<=q<=34
  1973.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 0x5000000000000000a, 1<=q<=34
  1974.       C.w[1] = 0x0000000000000005ull;
  1975.       C.w[0] = 0x000000000000000aull;
  1976.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  1977.         // 10^(20-q) is 64-bit, and so is C1
  1978.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  1979.       } else if (q == 20) {
  1980.         ;       // C1 * 10^0 = C1
  1981.       } else {  // if 21 <= q <= 34
  1982.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  1983.       }
  1984.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  1985.         // set invalid flag
  1986.         *pfpsf |= INVALID_EXCEPTION;
  1987.         // return Integer Indefinite
  1988.         res = 0x8000000000000000ull;
  1989.         BID_RETURN (res);
  1990.       }
  1991.       // else cases that can be rounded to a 64-bit int fall through
  1992.       // to '1 <= q + exp <= 19'
  1993.     } else {    // if n > 0 and q + exp = 19
  1994.       // if n >= 2^63 then n is too large
  1995.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63
  1996.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*2^64, 1<=q<=34
  1997.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000000, 1<=q<=34
  1998.       C.w[1] = 0x0000000000000005ull;
  1999.       C.w[0] = 0x0000000000000000ull;
  2000.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  2001.         // 10^(20-q) is 64-bit, and so is C1
  2002.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  2003.       } else if (q == 20) {
  2004.         ;       // C1 * 10^0 = C1
  2005.       } else {  // if 21 <= q <= 34
  2006.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  2007.       }
  2008.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2009.         // set invalid flag
  2010.         *pfpsf |= INVALID_EXCEPTION;
  2011.         // return Integer Indefinite
  2012.         res = 0x8000000000000000ull;
  2013.         BID_RETURN (res);
  2014.       }
  2015.       // else cases that can be rounded to a 64-bit int fall through
  2016.       // to '1 <= q + exp <= 19'
  2017.     }
  2018.   }
  2019.   // n is not too large to be converted to int64: -2^63-1 < n < 2^63
  2020.   // Note: some of the cases tested for above fall through to this point
  2021.   // Restore C1 which may have been modified above
  2022.   C1.w[1] = x.w[1] & MASK_COEFF;
  2023.   C1.w[0] = x.w[0];
  2024.   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
  2025.     // return 0
  2026.     res = 0x0000000000000000ull;
  2027.     BID_RETURN (res);
  2028.   } else {      // if (1 <= q + exp <= 19, 1 <= q <= 34, -33 <= exp <= 18)
  2029.     // -2^63-1 < x <= -1 or 1 <= x < 2^63 so x can be rounded
  2030.     // toward zero to a 64-bit signed integer
  2031.     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 19
  2032.       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
  2033.       // chop off ind digits from the lower part of C1
  2034.       // C1 fits in 127 bits
  2035.       // calculate C* and f*
  2036.       // C* is actually floor(C*) in this case
  2037.       // C* and f* need shifting and masking, as shown by
  2038.       // shiftright128[] and maskhigh128[]
  2039.       // 1 <= x <= 33
  2040.       // kx = 10^(-x) = ten2mk128[ind - 1]
  2041.       // C* = C1 * 10^(-x)
  2042.       // the approximation of 10^(-x) was rounded up to 118 bits
  2043.       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  2044.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  2045.         Cstar.w[1] = P256.w[3];
  2046.         Cstar.w[0] = P256.w[2];
  2047.       } else {  // 22 <= ind - 1 <= 33
  2048.         Cstar.w[1] = 0;
  2049.         Cstar.w[0] = P256.w[3];
  2050.       }
  2051.       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  2052.       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  2053.       // C* = floor(C*) (logical right shift; C has p decimal digits,
  2054.       //     correct by Property 1)
  2055.       // n = C* * 10^(e+x)
  2056.  
  2057.       // shift right C* by Ex-128 = shiftright128[ind]
  2058.       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
  2059.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  2060.         Cstar.w[0] =
  2061.           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  2062.         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  2063.       } else {  // 22 <= ind - 1 <= 33
  2064.         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
  2065.       }
  2066.       if (x_sign)
  2067.         res = -Cstar.w[0];
  2068.       else
  2069.         res = Cstar.w[0];
  2070.     } else if (exp == 0) {
  2071.       // 1 <= q <= 19
  2072.       // res = +/-C (exact)
  2073.       if (x_sign)
  2074.         res = -C1.w[0];
  2075.       else
  2076.         res = C1.w[0];
  2077.     } else {    // if (exp>0) => 1 <= exp <= 18, 1 <= q < 18, 2 <= q + exp <= 19
  2078.       // res = +/-C * 10^exp (exact) where this fits in 64-bit integer
  2079.       if (x_sign)
  2080.         res = -C1.w[0] * ten2k64[exp];
  2081.       else
  2082.         res = C1.w[0] * ten2k64[exp];
  2083.     }
  2084.   }
  2085. }
  2086.  
  2087. BID_RETURN (res);
  2088. }
  2089.  
  2090. /*****************************************************************************
  2091.  *  BID128_to_xint64_xint
  2092.  ****************************************************************************/
  2093.  
  2094. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (SINT64, bid128_to_int64_xint,
  2095.                                           x)
  2096.  
  2097.      SINT64 res;
  2098.      UINT64 x_sign;
  2099.      UINT64 x_exp;
  2100.      int exp;                   // unbiased exponent
  2101.   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  2102.      BID_UI64DOUBLE tmp1;
  2103.      unsigned int x_nr_bits;
  2104.      int q, ind, shift;
  2105.      UINT128 C1, C;
  2106.      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  2107.      UINT256 fstar;
  2108.      UINT256 P256;
  2109.  
  2110.   // unpack x
  2111. x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  2112. x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
  2113. C1.w[1] = x.w[1] & MASK_COEFF;
  2114. C1.w[0] = x.w[0];
  2115.  
  2116.   // check for NaN or Infinity
  2117. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  2118.     // x is special
  2119. if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
  2120.   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
  2121.     // set invalid flag
  2122.     *pfpsf |= INVALID_EXCEPTION;
  2123.     // return Integer Indefinite
  2124.     res = 0x8000000000000000ull;
  2125.   } else {      // x is QNaN
  2126.     // set invalid flag
  2127.     *pfpsf |= INVALID_EXCEPTION;
  2128.     // return Integer Indefinite
  2129.     res = 0x8000000000000000ull;
  2130.   }
  2131.   BID_RETURN (res);
  2132. } else {        // x is not a NaN, so it must be infinity
  2133.   if (!x_sign) {        // x is +inf
  2134.     // set invalid flag
  2135.     *pfpsf |= INVALID_EXCEPTION;
  2136.     // return Integer Indefinite
  2137.     res = 0x8000000000000000ull;
  2138.   } else {      // x is -inf
  2139.     // set invalid flag
  2140.     *pfpsf |= INVALID_EXCEPTION;
  2141.     // return Integer Indefinite
  2142.     res = 0x8000000000000000ull;
  2143.   }
  2144.   BID_RETURN (res);
  2145. }
  2146. }
  2147.   // check for non-canonical values (after the check for special values)
  2148. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  2149.     || (C1.w[1] == 0x0001ed09bead87c0ull
  2150.         && (C1.w[0] > 0x378d8e63ffffffffull))
  2151.     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  2152.   res = 0x0000000000000000ull;
  2153.   BID_RETURN (res);
  2154. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  2155.   // x is 0
  2156.   res = 0x0000000000000000ull;
  2157.   BID_RETURN (res);
  2158. } else {        // x is not special and is not zero
  2159.  
  2160.   // q = nr. of decimal digits in x
  2161.   //  determine first the nr. of bits in x
  2162.   if (C1.w[1] == 0) {
  2163.     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
  2164.       // split the 64-bit value in two 32-bit halves to avoid rounding errors
  2165.       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
  2166.         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
  2167.         x_nr_bits =
  2168.           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2169.       } else {  // x < 2^32
  2170.         tmp1.d = (double) (C1.w[0]);    // exact conversion
  2171.         x_nr_bits =
  2172.           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2173.       }
  2174.     } else {    // if x < 2^53
  2175.       tmp1.d = (double) C1.w[0];        // exact conversion
  2176.       x_nr_bits =
  2177.         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2178.     }
  2179.   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  2180.     tmp1.d = (double) C1.w[1];  // exact conversion
  2181.     x_nr_bits =
  2182.       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2183.   }
  2184.   q = nr_digits[x_nr_bits - 1].digits;
  2185.   if (q == 0) {
  2186.     q = nr_digits[x_nr_bits - 1].digits1;
  2187.     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  2188.         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  2189.             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  2190.       q++;
  2191.   }
  2192.   exp = (x_exp >> 49) - 6176;
  2193.   if ((q + exp) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
  2194.     // set invalid flag
  2195.     *pfpsf |= INVALID_EXCEPTION;
  2196.     // return Integer Indefinite
  2197.     res = 0x8000000000000000ull;
  2198.     BID_RETURN (res);
  2199.   } else if ((q + exp) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
  2200.     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
  2201.     // so x rounded to an integer may or may not fit in a signed 64-bit int
  2202.     // the cases that do not fit are identified here; the ones that fit
  2203.     // fall through and will be handled with other cases further,
  2204.     // under '1 <= q + exp <= 19'
  2205.     if (x_sign) {       // if n < 0 and q + exp = 19
  2206.       // if n <= -2^63 - 1 then n is too large
  2207.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1
  2208.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64+2), 1<=q<=34
  2209.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 0x5000000000000000a, 1<=q<=34
  2210.       C.w[1] = 0x0000000000000005ull;
  2211.       C.w[0] = 0x000000000000000aull;
  2212.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  2213.         // 10^(20-q) is 64-bit, and so is C1
  2214.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  2215.       } else if (q == 20) {
  2216.         ;       // C1 * 10^0 = C1
  2217.       } else {  // if 21 <= q <= 34
  2218.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  2219.       }
  2220.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2221.         // set invalid flag
  2222.         *pfpsf |= INVALID_EXCEPTION;
  2223.         // return Integer Indefinite
  2224.         res = 0x8000000000000000ull;
  2225.         BID_RETURN (res);
  2226.       }
  2227.       // else cases that can be rounded to a 64-bit int fall through
  2228.       // to '1 <= q + exp <= 19'
  2229.     } else {    // if n > 0 and q + exp = 19
  2230.       // if n >= 2^63 then n is too large
  2231.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63
  2232.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*2^64, 1<=q<=34
  2233.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000000, 1<=q<=34
  2234.       C.w[1] = 0x0000000000000005ull;
  2235.       C.w[0] = 0x0000000000000000ull;
  2236.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  2237.         // 10^(20-q) is 64-bit, and so is C1
  2238.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  2239.       } else if (q == 20) {
  2240.         ;       // C1 * 10^0 = C1
  2241.       } else {  // if 21 <= q <= 34
  2242.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  2243.       }
  2244.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2245.         // set invalid flag
  2246.         *pfpsf |= INVALID_EXCEPTION;
  2247.         // return Integer Indefinite
  2248.         res = 0x8000000000000000ull;
  2249.         BID_RETURN (res);
  2250.       }
  2251.       // else cases that can be rounded to a 64-bit int fall through
  2252.       // to '1 <= q + exp <= 19'
  2253.     }
  2254.   }
  2255.   // n is not too large to be converted to int64: -2^63-1 < n < 2^63
  2256.   // Note: some of the cases tested for above fall through to this point
  2257.   // Restore C1 which may have been modified above
  2258.   C1.w[1] = x.w[1] & MASK_COEFF;
  2259.   C1.w[0] = x.w[0];
  2260.   if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
  2261.     // set inexact flag
  2262.     *pfpsf |= INEXACT_EXCEPTION;
  2263.     // return 0
  2264.     res = 0x0000000000000000ull;
  2265.     BID_RETURN (res);
  2266.   } else {      // if (1 <= q + exp <= 19, 1 <= q <= 34, -33 <= exp <= 18)
  2267.     // -2^63-1 < x <= -1 or 1 <= x < 2^63 so x can be rounded
  2268.     // toward zero to a 64-bit signed integer
  2269.     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 19
  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.         fstar.w[3] = 0;
  2286.         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  2287.         fstar.w[1] = P256.w[1];
  2288.         fstar.w[0] = P256.w[0];
  2289.       } else {  // 22 <= ind - 1 <= 33
  2290.         Cstar.w[1] = 0;
  2291.         Cstar.w[0] = P256.w[3];
  2292.         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  2293.         fstar.w[2] = P256.w[2];
  2294.         fstar.w[1] = P256.w[1];
  2295.         fstar.w[0] = P256.w[0];
  2296.       }
  2297.       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  2298.       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  2299.       // C* = floor(C*) (logical right shift; C has p decimal digits,
  2300.       //     correct by Property 1)
  2301.       // n = C* * 10^(e+x)
  2302.  
  2303.       // shift right C* by Ex-128 = shiftright128[ind]
  2304.       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
  2305.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  2306.         Cstar.w[0] =
  2307.           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  2308.         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  2309.       } else {  // 22 <= ind - 1 <= 33
  2310.         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
  2311.       }
  2312.       // determine inexactness of the rounding of C*
  2313.       // if (0 < f* < 10^(-x)) then
  2314.       //   the result is exact
  2315.       // else // if (f* > T*) then
  2316.       //   the result is inexact
  2317.       if (ind - 1 <= 2) {
  2318.         if (fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2319.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2320.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2321.           // set the inexact flag
  2322.           *pfpsf |= INEXACT_EXCEPTION;
  2323.         }       // else the result is exact
  2324.       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
  2325.         if (fstar.w[2] || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2326.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2327.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2328.           // set the inexact flag
  2329.           *pfpsf |= INEXACT_EXCEPTION;
  2330.         }
  2331.       } else {  // if 22 <= ind <= 33
  2332.         if (fstar.w[3] || fstar.w[2]
  2333.             || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2334.             || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2335.                 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2336.           // set the inexact flag
  2337.           *pfpsf |= INEXACT_EXCEPTION;
  2338.         }
  2339.       }
  2340.  
  2341.       if (x_sign)
  2342.         res = -Cstar.w[0];
  2343.       else
  2344.         res = Cstar.w[0];
  2345.     } else if (exp == 0) {
  2346.       // 1 <= q <= 19
  2347.       // res = +/-C (exact)
  2348.       if (x_sign)
  2349.         res = -C1.w[0];
  2350.       else
  2351.         res = C1.w[0];
  2352.     } else {    // if (exp>0) => 1 <= exp <= 18, 1 <= q < 18, 2 <= q + exp <= 19
  2353.       // res = +/-C * 10^exp (exact) where this fits in 64-bit integer
  2354.       if (x_sign)
  2355.         res = -C1.w[0] * ten2k64[exp];
  2356.       else
  2357.         res = C1.w[0] * ten2k64[exp];
  2358.     }
  2359.   }
  2360. }
  2361.  
  2362. BID_RETURN (res);
  2363. }
  2364.  
  2365. /*****************************************************************************
  2366.  *  BID128_to_int64_rninta
  2367.  ****************************************************************************/
  2368.  
  2369. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (SINT64,
  2370.                                           bid128_to_int64_rninta, x)
  2371.  
  2372.      SINT64 res;
  2373.      UINT64 x_sign;
  2374.      UINT64 x_exp;
  2375.      int exp;                   // unbiased exponent
  2376.   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  2377.      UINT64 tmp64;
  2378.      BID_UI64DOUBLE tmp1;
  2379.      unsigned int x_nr_bits;
  2380.      int q, ind, shift;
  2381.      UINT128 C1, C;
  2382.      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  2383.      UINT256 P256;
  2384.  
  2385.   // unpack x
  2386. x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
  2387. x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
  2388. C1.w[1] = x.w[1] & MASK_COEFF;
  2389. C1.w[0] = x.w[0];
  2390.  
  2391.   // check for NaN or Infinity
  2392. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  2393.     // x is special
  2394. if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
  2395.   if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
  2396.     // set invalid flag
  2397.     *pfpsf |= INVALID_EXCEPTION;
  2398.     // return Integer Indefinite
  2399.     res = 0x8000000000000000ull;
  2400.   } else {      // x is QNaN
  2401.     // set invalid flag
  2402.     *pfpsf |= INVALID_EXCEPTION;
  2403.     // return Integer Indefinite
  2404.     res = 0x8000000000000000ull;
  2405.   }
  2406.   BID_RETURN (res);
  2407. } else {        // x is not a NaN, so it must be infinity
  2408.   if (!x_sign) {        // x is +inf
  2409.     // set invalid flag
  2410.     *pfpsf |= INVALID_EXCEPTION;
  2411.     // return Integer Indefinite
  2412.     res = 0x8000000000000000ull;
  2413.   } else {      // x is -inf
  2414.     // set invalid flag
  2415.     *pfpsf |= INVALID_EXCEPTION;
  2416.     // return Integer Indefinite
  2417.     res = 0x8000000000000000ull;
  2418.   }
  2419.   BID_RETURN (res);
  2420. }
  2421. }
  2422.   // check for non-canonical values (after the check for special values)
  2423. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  2424.     || (C1.w[1] == 0x0001ed09bead87c0ull
  2425.         && (C1.w[0] > 0x378d8e63ffffffffull))
  2426.     || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  2427.   res = 0x0000000000000000ull;
  2428.   BID_RETURN (res);
  2429. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  2430.   // x is 0
  2431.   res = 0x0000000000000000ull;
  2432.   BID_RETURN (res);
  2433. } else {        // x is not special and is not zero
  2434.  
  2435.   // q = nr. of decimal digits in x
  2436.   //  determine first the nr. of bits in x
  2437.   if (C1.w[1] == 0) {
  2438.     if (C1.w[0] >= 0x0020000000000000ull) {     // x >= 2^53
  2439.       // split the 64-bit value in two 32-bit halves to avoid rounding errors
  2440.       if (C1.w[0] >= 0x0000000100000000ull) {   // x >= 2^32
  2441.         tmp1.d = (double) (C1.w[0] >> 32);      // exact conversion
  2442.         x_nr_bits =
  2443.           33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2444.       } else {  // x < 2^32
  2445.         tmp1.d = (double) (C1.w[0]);    // exact conversion
  2446.         x_nr_bits =
  2447.           1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2448.       }
  2449.     } else {    // if x < 2^53
  2450.       tmp1.d = (double) C1.w[0];        // exact conversion
  2451.       x_nr_bits =
  2452.         1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2453.     }
  2454.   } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  2455.     tmp1.d = (double) C1.w[1];  // exact conversion
  2456.     x_nr_bits =
  2457.       65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2458.   }
  2459.   q = nr_digits[x_nr_bits - 1].digits;
  2460.   if (q == 0) {
  2461.     q = nr_digits[x_nr_bits - 1].digits1;
  2462.     if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  2463.         || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  2464.             && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  2465.       q++;
  2466.   }
  2467.   exp = (x_exp >> 49) - 6176;
  2468.   if ((q + exp) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
  2469.     // set invalid flag
  2470.     *pfpsf |= INVALID_EXCEPTION;
  2471.     // return Integer Indefinite
  2472.     res = 0x8000000000000000ull;
  2473.     BID_RETURN (res);
  2474.   } else if ((q + exp) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
  2475.     // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
  2476.     // so x rounded to an integer may or may not fit in a signed 64-bit int
  2477.     // the cases that do not fit are identified here; the ones that fit
  2478.     // fall through and will be handled with other cases further,
  2479.     // under '1 <= q + exp <= 19'
  2480.     if (x_sign) {       // if n < 0 and q + exp = 19
  2481.       // if n <= -2^63 - 1/2 then n is too large
  2482.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1/2
  2483.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64+1), 1<=q<=34
  2484.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000005, 1<=q<=34
  2485.       C.w[1] = 0x0000000000000005ull;
  2486.       C.w[0] = 0000000000000005ull;
  2487.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  2488.         // 10^(20-q) is 64-bit, and so is C1
  2489.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  2490.       } else if (q == 20) {
  2491.         ;       // C1 * 10^0 = C1
  2492.       } else {  // if 21 <= q <= 34
  2493.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  2494.       }
  2495.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2496.         // set invalid flag
  2497.         *pfpsf |= INVALID_EXCEPTION;
  2498.         // return Integer Indefinite
  2499.         res = 0x8000000000000000ull;
  2500.         BID_RETURN (res);
  2501.       }
  2502.       // else cases that can be rounded to a 64-bit int fall through
  2503.       // to '1 <= q + exp <= 19'
  2504.     } else {    // if n > 0 and q + exp = 19
  2505.       // if n >= 2^63 - 1/2 then n is too large
  2506.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63-1/2
  2507.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64-1), 1<=q<=34
  2508.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x4fffffffffffffffb, 1<=q<=34
  2509.       C.w[1] = 0x0000000000000004ull;
  2510.       C.w[0] = 0xfffffffffffffffbull;
  2511.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  2512.         // 10^(20-q) is 64-bit, and so is C1
  2513.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  2514.       } else if (q == 20) {
  2515.         ;       // C1 * 10^0 = C1
  2516.       } else {  // if 21 <= q <= 34
  2517.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  2518.       }
  2519.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2520.         // set invalid flag
  2521.         *pfpsf |= INVALID_EXCEPTION;
  2522.         // return Integer Indefinite
  2523.         res = 0x8000000000000000ull;
  2524.         BID_RETURN (res);
  2525.       }
  2526.       // else cases that can be rounded to a 64-bit int fall through
  2527.       // to '1 <= q + exp <= 19'
  2528.     }
  2529.   }
  2530.   // n is not too large to be converted to int64: -2^63-1/2 <= n < 2^63-1/2
  2531.   // Note: some of the cases tested for above fall through to this point
  2532.   // Restore C1 which may have been modified above
  2533.   C1.w[1] = x.w[1] & MASK_COEFF;
  2534.   C1.w[0] = x.w[0];
  2535.   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
  2536.     // return 0
  2537.     res = 0x0000000000000000ull;
  2538.     BID_RETURN (res);
  2539.   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
  2540.     // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
  2541.     //   res = 0
  2542.     // else
  2543.     //   res = +/-1
  2544.     ind = q - 1;
  2545.     if (ind <= 18) {    // 0 <= ind <= 18
  2546.       if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
  2547.         res = 0x0000000000000000ull;    // return 0
  2548.       } else if (x_sign) {      // n < 0
  2549.         res = 0xffffffffffffffffull;    // return -1
  2550.       } else {  // n > 0
  2551.         res = 0x0000000000000001ull;    // return +1
  2552.       }
  2553.     } else {    // 19 <= ind <= 33
  2554.       if ((C1.w[1] < midpoint128[ind - 19].w[1])
  2555.           || ((C1.w[1] == midpoint128[ind - 19].w[1])
  2556.               && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
  2557.         res = 0x0000000000000000ull;    // return 0
  2558.       } else if (x_sign) {      // n < 0
  2559.         res = 0xffffffffffffffffull;    // return -1
  2560.       } else {  // n > 0
  2561.         res = 0x0000000000000001ull;    // return +1
  2562.       }
  2563.     }
  2564.   } else {      // if (1 <= q + exp <= 19, 1 <= q <= 34, -33 <= exp <= 18)
  2565.     // -2^63-1/2 <= x <= -1 or 1 <= x < 2^63-1/2 so x can be rounded
  2566.     // to nearest to a 64-bit signed integer
  2567.     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 19
  2568.       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
  2569.       // chop off ind digits from the lower part of C1
  2570.       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  2571.       tmp64 = C1.w[0];
  2572.       if (ind <= 19) {
  2573.         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  2574.       } else {
  2575.         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  2576.         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  2577.       }
  2578.       if (C1.w[0] < tmp64)
  2579.         C1.w[1]++;
  2580.       // calculate C* and f*
  2581.       // C* is actually floor(C*) in this case
  2582.       // C* and f* need shifting and masking, as shown by
  2583.       // shiftright128[] and maskhigh128[]
  2584.       // 1 <= x <= 33
  2585.       // kx = 10^(-x) = ten2mk128[ind - 1]
  2586.       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  2587.       // the approximation of 10^(-x) was rounded up to 118 bits
  2588.       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  2589.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  2590.         Cstar.w[1] = P256.w[3];
  2591.         Cstar.w[0] = P256.w[2];
  2592.       } else {  // 22 <= ind - 1 <= 33
  2593.         Cstar.w[1] = 0;
  2594.         Cstar.w[0] = P256.w[3];
  2595.       }
  2596.       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  2597.       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  2598.       // if (0 < f* < 10^(-x)) then the result is a midpoint
  2599.       //   if floor(C*) is even then C* = floor(C*) - logical right
  2600.       //       shift; C* has p decimal digits, correct by Prop. 1)
  2601.       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
  2602.       //       shift; C* has p decimal digits, correct by Pr. 1)
  2603.       // else
  2604.       //   C* = floor(C*) (logical right shift; C has p decimal digits,
  2605.       //       correct by Property 1)
  2606.       // n = C* * 10^(e+x)
  2607.  
  2608.       // shift right C* by Ex-128 = shiftright128[ind]
  2609.       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
  2610.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  2611.         Cstar.w[0] =
  2612.           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  2613.         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  2614.       } else {  // 22 <= ind - 1 <= 33
  2615.         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
  2616.       }
  2617.  
  2618.       // if the result was a midpoint it was rounded away from zero
  2619.       if (x_sign)
  2620.         res = -Cstar.w[0];
  2621.       else
  2622.         res = Cstar.w[0];
  2623.     } else if (exp == 0) {
  2624.       // 1 <= q <= 19
  2625.       // res = +/-C (exact)
  2626.       if (x_sign)
  2627.         res = -C1.w[0];
  2628.       else
  2629.         res = C1.w[0];
  2630.     } else {    // if (exp>0) => 1 <= exp <= 18, 1 <= q < 18, 2 <= q + exp <= 19
  2631.       // res = +/-C * 10^exp (exact) where this fits in 64-bit integer
  2632.       if (x_sign)
  2633.         res = -C1.w[0] * ten2k64[exp];
  2634.       else
  2635.         res = C1.w[0] * ten2k64[exp];
  2636.     }
  2637.   }
  2638. }
  2639.  
  2640. BID_RETURN (res);
  2641. }
  2642.  
  2643. /*****************************************************************************
  2644.  *  BID128_to_int64_xrninta
  2645.  ****************************************************************************/
  2646.  
  2647. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (SINT64,
  2648.                                           bid128_to_int64_xrninta, x)
  2649.  
  2650.      SINT64 res;
  2651.      UINT64 x_sign;
  2652.      UINT64 x_exp;
  2653.      int exp;                   // unbiased exponent
  2654.   // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  2655.      UINT64 tmp64, tmp64A;
  2656.      BID_UI64DOUBLE tmp1;
  2657.      unsigned int x_nr_bits;
  2658.      int q, ind, shift;
  2659.      UINT128 C1, C;
  2660.      UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
  2661.      UINT256 fstar;
  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) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
  2748.     // set invalid flag
  2749.     *pfpsf |= INVALID_EXCEPTION;
  2750.     // return Integer Indefinite
  2751.     res = 0x8000000000000000ull;
  2752.     BID_RETURN (res);
  2753.   } else if ((q + exp) == 19) { // x = c(0)c(1)...c(18).c(19)...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 a signed 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 <= 19'
  2759.     if (x_sign) {       // if n < 0 and q + exp = 19
  2760.       // if n <= -2^63 - 1/2 then n is too large
  2761.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1/2
  2762.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64+1), 1<=q<=34
  2763.       // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000005, 1<=q<=34
  2764.       C.w[1] = 0x0000000000000005ull;
  2765.       C.w[0] = 0000000000000005ull;
  2766.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  2767.         // 10^(20-q) is 64-bit, and so is C1
  2768.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  2769.       } else if (q == 20) {
  2770.         ;       // C1 * 10^0 = C1
  2771.       } else {  // if 21 <= q <= 34
  2772.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  2773.       }
  2774.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2775.         // set invalid flag
  2776.         *pfpsf |= INVALID_EXCEPTION;
  2777.         // return Integer Indefinite
  2778.         res = 0x8000000000000000ull;
  2779.         BID_RETURN (res);
  2780.       }
  2781.       // else cases that can be rounded to a 64-bit int fall through
  2782.       // to '1 <= q + exp <= 19'
  2783.     } else {    // if n > 0 and q + exp = 19
  2784.       // if n >= 2^63 - 1/2 then n is too large
  2785.       // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63-1/2
  2786.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64-1), 1<=q<=34
  2787.       // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x4fffffffffffffffb, 1<=q<=34
  2788.       C.w[1] = 0x0000000000000004ull;
  2789.       C.w[0] = 0xfffffffffffffffbull;
  2790.       if (q <= 19) {    // 1 <= q <= 19 => 1 <= 20-q <= 19 =>
  2791.         // 10^(20-q) is 64-bit, and so is C1
  2792.         __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[20 - q]);
  2793.       } else if (q == 20) {
  2794.         ;       // C1 * 10^0 = C1
  2795.       } else {  // if 21 <= q <= 34
  2796.         __mul_128x64_to_128 (C, ten2k64[q - 20], C);    // max 47-bit x 67-bit
  2797.       }
  2798.       if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2799.         // set invalid flag
  2800.         *pfpsf |= INVALID_EXCEPTION;
  2801.         // return Integer Indefinite
  2802.         res = 0x8000000000000000ull;
  2803.         BID_RETURN (res);
  2804.       }
  2805.       // else cases that can be rounded to a 64-bit int fall through
  2806.       // to '1 <= q + exp <= 19'
  2807.     }
  2808.   }
  2809.   // n is not too large to be converted to int64: -2^63-1/2 <= n < 2^63-1/2
  2810.   // Note: some of the cases tested for above fall through to this point
  2811.   // Restore C1 which may have been modified above
  2812.   C1.w[1] = x.w[1] & MASK_COEFF;
  2813.   C1.w[0] = x.w[0];
  2814.   if ((q + exp) < 0) {  // n = +/-0.0...c(0)c(1)...c(q-1)
  2815.     // set inexact flag
  2816.     *pfpsf |= INEXACT_EXCEPTION;
  2817.     // return 0
  2818.     res = 0x0000000000000000ull;
  2819.     BID_RETURN (res);
  2820.   } else if ((q + exp) == 0) {  // n = +/-0.c(0)c(1)...c(q-1)
  2821.     // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
  2822.     //   res = 0
  2823.     // else
  2824.     //   res = +/-1
  2825.     ind = q - 1;
  2826.     if (ind <= 18) {    // 0 <= ind <= 18
  2827.       if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
  2828.         res = 0x0000000000000000ull;    // return 0
  2829.       } else if (x_sign) {      // n < 0
  2830.         res = 0xffffffffffffffffull;    // return -1
  2831.       } else {  // n > 0
  2832.         res = 0x0000000000000001ull;    // return +1
  2833.       }
  2834.     } else {    // 19 <= ind <= 33
  2835.       if ((C1.w[1] < midpoint128[ind - 19].w[1])
  2836.           || ((C1.w[1] == midpoint128[ind - 19].w[1])
  2837.               && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
  2838.         res = 0x0000000000000000ull;    // return 0
  2839.       } else if (x_sign) {      // n < 0
  2840.         res = 0xffffffffffffffffull;    // return -1
  2841.       } else {  // n > 0
  2842.         res = 0x0000000000000001ull;    // return +1
  2843.       }
  2844.     }
  2845.     // set inexact flag
  2846.     *pfpsf |= INEXACT_EXCEPTION;
  2847.   } else {      // if (1 <= q + exp <= 19, 1 <= q <= 34, -33 <= exp <= 18)
  2848.     // -2^63-1/2 <= x <= -1 or 1 <= x < 2^63-1/2 so x can be rounded
  2849.     // to nearest to a 64-bit signed integer
  2850.     if (exp < 0) {      // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 19
  2851.       ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
  2852.       // chop off ind digits from the lower part of C1
  2853.       // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  2854.       tmp64 = C1.w[0];
  2855.       if (ind <= 19) {
  2856.         C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  2857.       } else {
  2858.         C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  2859.         C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  2860.       }
  2861.       if (C1.w[0] < tmp64)
  2862.         C1.w[1]++;
  2863.       // calculate C* and f*
  2864.       // C* is actually floor(C*) in this case
  2865.       // C* and f* need shifting and masking, as shown by
  2866.       // shiftright128[] and maskhigh128[]
  2867.       // 1 <= x <= 33
  2868.       // kx = 10^(-x) = ten2mk128[ind - 1]
  2869.       // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  2870.       // the approximation of 10^(-x) was rounded up to 118 bits
  2871.       __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  2872.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  2873.         Cstar.w[1] = P256.w[3];
  2874.         Cstar.w[0] = P256.w[2];
  2875.         fstar.w[3] = 0;
  2876.         fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  2877.         fstar.w[1] = P256.w[1];
  2878.         fstar.w[0] = P256.w[0];
  2879.       } else {  // 22 <= ind - 1 <= 33
  2880.         Cstar.w[1] = 0;
  2881.         Cstar.w[0] = P256.w[3];
  2882.         fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  2883.         fstar.w[2] = P256.w[2];
  2884.         fstar.w[1] = P256.w[1];
  2885.         fstar.w[0] = P256.w[0];
  2886.       }
  2887.       // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  2888.       // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  2889.       // if (0 < f* < 10^(-x)) then the result is a midpoint
  2890.       //   if floor(C*) is even then C* = floor(C*) - logical right
  2891.       //       shift; C* has p decimal digits, correct by Prop. 1)
  2892.       //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
  2893.       //       shift; C* has p decimal digits, correct by Pr. 1)
  2894.       // else
  2895.       //   C* = floor(C*) (logical right shift; C has p decimal digits,
  2896.       //       correct by Property 1)
  2897.       // n = C* * 10^(e+x)
  2898.  
  2899.       // shift right C* by Ex-128 = shiftright128[ind]
  2900.       shift = shiftright128[ind - 1];   // 0 <= shift <= 102
  2901.       if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
  2902.         Cstar.w[0] =
  2903.           (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  2904.         // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  2905.       } else {  // 22 <= ind - 1 <= 33
  2906.         Cstar.w[0] = (Cstar.w[0] >> (shift - 64));      // 2 <= shift - 64 <= 38
  2907.       }
  2908.       // determine inexactness of the rounding of C*
  2909.       // if (0 < f* - 1/2 < 10^(-x)) then
  2910.       //   the result is exact
  2911.       // else // if (f* - 1/2 > T*) then
  2912.       //   the result is inexact
  2913.       if (ind - 1 <= 2) {
  2914.         if (fstar.w[1] > 0x8000000000000000ull ||
  2915.             (fstar.w[1] == 0x8000000000000000ull
  2916.              && fstar.w[0] > 0x0ull)) {
  2917.           // f* > 1/2 and the result may be exact
  2918.           tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
  2919.           if (tmp64 > ten2mk128trunc[ind - 1].w[1]
  2920.               || (tmp64 == ten2mk128trunc[ind - 1].w[1]
  2921.                   && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
  2922.             // set the inexact flag
  2923.             *pfpsf |= INEXACT_EXCEPTION;
  2924.           }     // else the result is exact
  2925.         } else {        // the result is inexact; f2* <= 1/2
  2926.           // set the inexact flag
  2927.           *pfpsf |= INEXACT_EXCEPTION;
  2928.         }
  2929.       } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
  2930.         if (fstar.w[3] > 0x0 ||
  2931.             (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
  2932.             (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
  2933.              (fstar.w[1] || fstar.w[0]))) {
  2934.           // f2* > 1/2 and the result may be exact
  2935.           // Calculate f2* - 1/2
  2936.           tmp64 = fstar.w[2] - onehalf128[ind - 1];
  2937.           tmp64A = fstar.w[3];
  2938.           if (tmp64 > fstar.w[2])
  2939.             tmp64A--;
  2940.           if (tmp64A || tmp64
  2941.               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2942.               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2943.                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2944.             // set the inexact flag
  2945.             *pfpsf |= INEXACT_EXCEPTION;
  2946.           }     // else the result is exact
  2947.         } else {        // the result is inexact; f2* <= 1/2
  2948.           // set the inexact flag
  2949.           *pfpsf |= INEXACT_EXCEPTION;
  2950.         }
  2951.       } else {  // if 22 <= ind <= 33
  2952.         if (fstar.w[3] > onehalf128[ind - 1] ||
  2953.             (fstar.w[3] == onehalf128[ind - 1] &&
  2954.              (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
  2955.           // f2* > 1/2 and the result may be exact
  2956.           // Calculate f2* - 1/2
  2957.           tmp64 = fstar.w[3] - onehalf128[ind - 1];
  2958.           if (tmp64 || fstar.w[2]
  2959.               || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2960.               || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2961.                   && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2962.             // set the inexact flag
  2963.             *pfpsf |= INEXACT_EXCEPTION;
  2964.           }     // else the result is exact
  2965.         } else {        // the result is inexact; f2* <= 1/2
  2966.           // set the inexact flag
  2967.           *pfpsf |= INEXACT_EXCEPTION;
  2968.         }
  2969.       }
  2970.  
  2971.       // if the result was a midpoint it was rounded away from zero
  2972.       if (x_sign)
  2973.         res = -Cstar.w[0];
  2974.       else
  2975.         res = Cstar.w[0];
  2976.     } else if (exp == 0) {
  2977.       // 1 <= q <= 19
  2978.       // res = +/-C (exact)
  2979.       if (x_sign)
  2980.         res = -C1.w[0];
  2981.       else
  2982.         res = C1.w[0];
  2983.     } else {    // if (exp>0) => 1 <= exp <= 18, 1 <= q < 18, 2 <= q + exp <= 19
  2984.       // res = +/-C * 10^exp (exact) where this fits in 64-bit integer
  2985.       if (x_sign)
  2986.         res = -C1.w[0] * ten2k64[exp];
  2987.       else
  2988.         res = C1.w[0] * ten2k64[exp];
  2989.     }
  2990.   }
  2991. }
  2992.  
  2993. BID_RETURN (res);
  2994. }
  2995.