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