Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

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