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