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