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