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. /*****************************************************************************
  25.  *
  26.  *  BID64 encoding:
  27.  * ****************************************
  28.  *  63  62              53 52           0
  29.  * |---|------------------|--------------|
  30.  * | S |  Biased Exp (E)  |  Coeff (c)   |
  31.  * |---|------------------|--------------|
  32.  *
  33.  * bias = 398
  34.  * number = (-1)^s * 10^(E-398) * c
  35.  * coefficient range - 0 to (2^53)-1
  36.  * COEFF_MAX = 2^53-1 = 9007199254740991
  37.  *
  38.  *****************************************************************************
  39.  *
  40.  * BID128 encoding:
  41.  *   1-bit sign
  42.  *   14-bit biased exponent in [0x21, 0x3020] = [33, 12320]
  43.  *         unbiased exponent in [-6176, 6111]; exponent bias = 6176
  44.  *   113-bit unsigned binary integer coefficient (49-bit high + 64-bit low)
  45.  *   Note: 10^34-1 ~ 2^112.945555... < 2^113 => coefficient fits in 113 bits
  46.  *
  47.  *   Note: assume invalid encodings are not passed to this function
  48.  *
  49.  * Round a number C with q decimal digits, represented as a binary integer
  50.  * to q - x digits. Six different routines are provided for different values
  51.  * of q. The maximum value of q used in the library is q = 3 * P - 1 where
  52.  * P = 16 or P = 34 (so q <= 111 decimal digits).
  53.  * The partitioning is based on the following, where Kx is the scaled
  54.  * integer representing the value of 10^(-x) rounded up to a number of bits
  55.  * sufficient to ensure correct rounding:
  56.  *
  57.  * --------------------------------------------------------------------------
  58.  * q    x           max. value of  a            max number      min. number
  59.  *                                              of bits in C    of bits in Kx
  60.  * --------------------------------------------------------------------------
  61.  *
  62.  *                          GROUP 1: 64 bits
  63.  *                          round64_2_18 ()
  64.  *
  65.  * 2    [1,1]       10^1 - 1 < 2^3.33            4              4
  66.  * ...  ...         ...                         ...             ...
  67.  * 18   [1,17]      10^18 - 1 < 2^59.80         60              61
  68.  *
  69.  *                        GROUP 2: 128 bits
  70.  *                        round128_19_38 ()
  71.  *
  72.  * 19   [1,18]      10^19 - 1 < 2^63.11         64              65
  73.  * 20   [1,19]      10^20 - 1 < 2^66.44         67              68
  74.  * ...  ...         ...                         ...             ...
  75.  * 38   [1,37]      10^38 - 1 < 2^126.24        127             128
  76.  *
  77.  *                        GROUP 3: 192 bits
  78.  *                        round192_39_57 ()
  79.  *
  80.  * 39   [1,38]      10^39 - 1 < 2^129.56        130             131
  81.  * ...  ...         ...                         ...             ...
  82.  * 57   [1,56]      10^57 - 1 < 2^189.35        190             191
  83.  *
  84.  *                        GROUP 4: 256 bits
  85.  *                        round256_58_76 ()
  86.  *
  87.  * 58   [1,57]      10^58 - 1 < 2^192.68        193             194
  88.  * ...  ...         ...                         ...             ...
  89.  * 76   [1,75]      10^76 - 1 < 2^252.47        253             254
  90.  *
  91.  *                        GROUP 5: 320 bits
  92.  *                        round320_77_96 ()
  93.  *
  94.  * 77   [1,76]      10^77 - 1 < 2^255.79        256             257
  95.  * 78   [1,77]      10^78 - 1 < 2^259.12        260             261
  96.  * ...  ...         ...                         ...             ...
  97.  * 96   [1,95]      10^96 - 1 < 2^318.91        319             320
  98.  *
  99.  *                        GROUP 6: 384 bits
  100.  *                        round384_97_115 ()
  101.  *
  102.  * 97   [1,96]      10^97 - 1 < 2^322.23        323             324
  103.  * ...  ...         ...                         ...             ...
  104.  * 115  [1,114]     10^115 - 1 < 2^382.03       383             384
  105.  *
  106.  ****************************************************************************/
  107.  
  108. #include "bid_internal.h"
  109.  
  110. void
  111. round64_2_18 (int q,
  112.               int x,
  113.               UINT64 C,
  114.               UINT64 * ptr_Cstar,
  115.               int *incr_exp,
  116.               int *ptr_is_midpoint_lt_even,
  117.               int *ptr_is_midpoint_gt_even,
  118.               int *ptr_is_inexact_lt_midpoint,
  119.               int *ptr_is_inexact_gt_midpoint) {
  120.  
  121.   UINT128 P128;
  122.   UINT128 fstar;
  123.   UINT64 Cstar;
  124.   UINT64 tmp64;
  125.   int shift;
  126.   int ind;
  127.  
  128.   // Note:
  129.   //    In round128_2_18() positive numbers with 2 <= q <= 18 will be
  130.   //    rounded to nearest only for 1 <= x <= 3:
  131.   //     x = 1 or x = 2 when q = 17
  132.   //     x = 2 or x = 3 when q = 18
  133.   // However, for generality and possible uses outside the frame of IEEE 754R
  134.   // this implementation works for 1 <= x <= q - 1
  135.  
  136.   // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even,
  137.   // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are
  138.   // initialized to 0 by the caller
  139.  
  140.   // round a number C with q decimal digits, 2 <= q <= 18
  141.   // to q - x digits, 1 <= x <= 17
  142.   // C = C + 1/2 * 10^x where the result C fits in 64 bits
  143.   // (because the largest value is 999999999999999999 + 50000000000000000 =
  144.   // 0x0e92596fd628ffff, which fits in 60 bits)
  145.   ind = x - 1;  // 0 <= ind <= 16
  146.   C = C + midpoint64[ind];
  147.   // kx ~= 10^(-x), kx = Kx64[ind] * 2^(-Ex), 0 <= ind <= 16
  148.   // P128 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx
  149.   // the approximation kx of 10^(-x) was rounded up to 64 bits
  150.   __mul_64x64_to_128MACH (P128, C, Kx64[ind]);
  151.   // calculate C* = floor (P128) and f*
  152.   // Cstar = P128 >> Ex
  153.   // fstar = low Ex bits of P128
  154.   shift = Ex64m64[ind]; // in [3, 56]
  155.   Cstar = P128.w[1] >> shift;
  156.   fstar.w[1] = P128.w[1] & mask64[ind];
  157.   fstar.w[0] = P128.w[0];
  158.   // the top Ex bits of 10^(-x) are T* = ten2mxtrunc64[ind], e.g.
  159.   // if x=1, T*=ten2mxtrunc64[0]=0xcccccccccccccccc
  160.   // if (0 < f* < 10^(-x)) then the result is a midpoint
  161.   //   if floor(C*) is even then C* = floor(C*) - logical right
  162.   //       shift; C* has q - x decimal digits, correct by Prop. 1)
  163.   //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
  164.   //       shift; C* has q - x decimal digits, correct by Pr. 1)
  165.   // else
  166.   //   C* = floor(C*) (logical right shift; C has q - x decimal digits,
  167.   //       correct by Property 1)
  168.   // in the caling function n = C* * 10^(e+x)
  169.  
  170.   // determine inexactness of the rounding of C*
  171.   // if (0 < f* - 1/2 < 10^(-x)) then
  172.   //   the result is exact
  173.   // else // if (f* - 1/2 > T*) then
  174.   //   the result is inexact
  175.   if (fstar.w[1] > half64[ind] ||
  176.       (fstar.w[1] == half64[ind] && fstar.w[0])) {
  177.     // f* > 1/2 and the result may be exact
  178.     // Calculate f* - 1/2
  179.     tmp64 = fstar.w[1] - half64[ind];
  180.     if (tmp64 || fstar.w[0] > ten2mxtrunc64[ind]) {     // f* - 1/2 > 10^(-x)
  181.       *ptr_is_inexact_lt_midpoint = 1;
  182.     }   // else the result is exact
  183.   } else {      // the result is inexact; f2* <= 1/2
  184.     *ptr_is_inexact_gt_midpoint = 1;
  185.   }
  186.   // check for midpoints (could do this before determining inexactness)
  187.   if (fstar.w[1] == 0 && fstar.w[0] <= ten2mxtrunc64[ind]) {
  188.     // the result is a midpoint
  189.     if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD]
  190.       // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
  191.       Cstar--;  // Cstar is now even
  192.       *ptr_is_midpoint_gt_even = 1;
  193.       *ptr_is_inexact_lt_midpoint = 0;
  194.       *ptr_is_inexact_gt_midpoint = 0;
  195.     } else {    // else MP in [ODD, EVEN]
  196.       *ptr_is_midpoint_lt_even = 1;
  197.       *ptr_is_inexact_lt_midpoint = 0;
  198.       *ptr_is_inexact_gt_midpoint = 0;
  199.     }
  200.   }
  201.   // check for rounding overflow, which occurs if Cstar = 10^(q-x)
  202.   ind = q - x;  // 1 <= ind <= q - 1
  203.   if (Cstar == ten2k64[ind]) {  // if  Cstar = 10^(q-x)
  204.     Cstar = ten2k64[ind - 1];   // Cstar = 10^(q-x-1)
  205.     *incr_exp = 1;
  206.   } else {      // 10^33 <= Cstar <= 10^34 - 1
  207.     *incr_exp = 0;
  208.   }
  209.   *ptr_Cstar = Cstar;
  210. }
  211.  
  212.  
  213. void
  214. round128_19_38 (int q,
  215.                 int x,
  216.                 UINT128 C,
  217.                 UINT128 * ptr_Cstar,
  218.                 int *incr_exp,
  219.                 int *ptr_is_midpoint_lt_even,
  220.                 int *ptr_is_midpoint_gt_even,
  221.                 int *ptr_is_inexact_lt_midpoint,
  222.                 int *ptr_is_inexact_gt_midpoint) {
  223.  
  224.   UINT256 P256;
  225.   UINT256 fstar;
  226.   UINT128 Cstar;
  227.   UINT64 tmp64;
  228.   int shift;
  229.   int ind;
  230.  
  231.   // Note:
  232.   //    In round128_19_38() positive numbers with 19 <= q <= 38 will be
  233.   //    rounded to nearest only for 1 <= x <= 23:
  234.   //     x = 3 or x = 4 when q = 19
  235.   //     x = 4 or x = 5 when q = 20
  236.   //     ...
  237.   //     x = 18 or x = 19 when q = 34
  238.   //     x = 1 or x = 2 or x = 19 or x = 20 when q = 35
  239.   //     x = 2 or x = 3 or x = 20 or x = 21 when q = 36
  240.   //     x = 3 or x = 4 or x = 21 or x = 22 when q = 37
  241.   //     x = 4 or x = 5 or x = 22 or x = 23 when q = 38
  242.   // However, for generality and possible uses outside the frame of IEEE 754R
  243.   // this implementation works for 1 <= x <= q - 1
  244.  
  245.   // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even,
  246.   // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are
  247.   // initialized to 0 by the caller
  248.  
  249.   // round a number C with q decimal digits, 19 <= q <= 38
  250.   // to q - x digits, 1 <= x <= 37
  251.   // C = C + 1/2 * 10^x where the result C fits in 128 bits
  252.   // (because the largest value is 99999999999999999999999999999999999999 +
  253.   // 5000000000000000000000000000000000000 =
  254.   // 0x4efe43b0c573e7e68a043d8fffffffff, which fits is 127 bits)
  255.  
  256.   ind = x - 1;  // 0 <= ind <= 36
  257.   if (ind <= 18) {      // if 0 <= ind <= 18
  258.     tmp64 = C.w[0];
  259.     C.w[0] = C.w[0] + midpoint64[ind];
  260.     if (C.w[0] < tmp64)
  261.       C.w[1]++;
  262.   } else {      // if 19 <= ind <= 37
  263.     tmp64 = C.w[0];
  264.     C.w[0] = C.w[0] + midpoint128[ind - 19].w[0];
  265.     if (C.w[0] < tmp64) {
  266.       C.w[1]++;
  267.     }
  268.     C.w[1] = C.w[1] + midpoint128[ind - 19].w[1];
  269.   }
  270.   // kx ~= 10^(-x), kx = Kx128[ind] * 2^(-Ex), 0 <= ind <= 36
  271.   // P256 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx
  272.   // the approximation kx of 10^(-x) was rounded up to 128 bits
  273.   __mul_128x128_to_256 (P256, C, Kx128[ind]);
  274.   // calculate C* = floor (P256) and f*
  275.   // Cstar = P256 >> Ex
  276.   // fstar = low Ex bits of P256
  277.   shift = Ex128m128[ind];       // in [2, 63] but have to consider two cases
  278.   if (ind <= 18) {      // if 0 <= ind <= 18
  279.     Cstar.w[0] = (P256.w[2] >> shift) | (P256.w[3] << (64 - shift));
  280.     Cstar.w[1] = (P256.w[3] >> shift);
  281.     fstar.w[0] = P256.w[0];
  282.     fstar.w[1] = P256.w[1];
  283.     fstar.w[2] = P256.w[2] & mask128[ind];
  284.     fstar.w[3] = 0x0ULL;
  285.   } else {      // if 19 <= ind <= 37
  286.     Cstar.w[0] = P256.w[3] >> shift;
  287.     Cstar.w[1] = 0x0ULL;
  288.     fstar.w[0] = P256.w[0];
  289.     fstar.w[1] = P256.w[1];
  290.     fstar.w[2] = P256.w[2];
  291.     fstar.w[3] = P256.w[3] & mask128[ind];
  292.   }
  293.   // the top Ex bits of 10^(-x) are T* = ten2mxtrunc64[ind], e.g.
  294.   // if x=1, T*=ten2mxtrunc128[0]=0xcccccccccccccccccccccccccccccccc
  295.   // if (0 < f* < 10^(-x)) then the result is a midpoint
  296.   //   if floor(C*) is even then C* = floor(C*) - logical right
  297.   //       shift; C* has q - x decimal digits, correct by Prop. 1)
  298.   //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
  299.   //       shift; C* has q - x decimal digits, correct by Pr. 1)
  300.   // else
  301.   //   C* = floor(C*) (logical right shift; C has q - x decimal digits,
  302.   //       correct by Property 1)
  303.   // in the caling function n = C* * 10^(e+x)
  304.  
  305.   // determine inexactness of the rounding of C*
  306.   // if (0 < f* - 1/2 < 10^(-x)) then
  307.   //   the result is exact
  308.   // else // if (f* - 1/2 > T*) then
  309.   //   the result is inexact
  310.   if (ind <= 18) {      // if 0 <= ind <= 18
  311.     if (fstar.w[2] > half128[ind] ||
  312.         (fstar.w[2] == half128[ind] && (fstar.w[1] || fstar.w[0]))) {
  313.       // f* > 1/2 and the result may be exact
  314.       // Calculate f* - 1/2
  315.       tmp64 = fstar.w[2] - half128[ind];
  316.       if (tmp64 || fstar.w[1] > ten2mxtrunc128[ind].w[1] || (fstar.w[1] == ten2mxtrunc128[ind].w[1] && fstar.w[0] > ten2mxtrunc128[ind].w[0])) {        // f* - 1/2 > 10^(-x)
  317.         *ptr_is_inexact_lt_midpoint = 1;
  318.       } // else the result is exact
  319.     } else {    // the result is inexact; f2* <= 1/2
  320.       *ptr_is_inexact_gt_midpoint = 1;
  321.     }
  322.   } else {      // if 19 <= ind <= 37
  323.     if (fstar.w[3] > half128[ind] || (fstar.w[3] == half128[ind] &&
  324.                                       (fstar.w[2] || fstar.w[1]
  325.                                        || fstar.w[0]))) {
  326.       // f* > 1/2 and the result may be exact
  327.       // Calculate f* - 1/2
  328.       tmp64 = fstar.w[3] - half128[ind];
  329.       if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mxtrunc128[ind].w[1] || (fstar.w[1] == ten2mxtrunc128[ind].w[1] && fstar.w[0] > ten2mxtrunc128[ind].w[0])) {  // f* - 1/2 > 10^(-x)
  330.         *ptr_is_inexact_lt_midpoint = 1;
  331.       } // else the result is exact
  332.     } else {    // the result is inexact; f2* <= 1/2
  333.       *ptr_is_inexact_gt_midpoint = 1;
  334.     }
  335.   }
  336.   // check for midpoints (could do this before determining inexactness)
  337.   if (fstar.w[3] == 0 && fstar.w[2] == 0 &&
  338.       (fstar.w[1] < ten2mxtrunc128[ind].w[1] ||
  339.        (fstar.w[1] == ten2mxtrunc128[ind].w[1] &&
  340.         fstar.w[0] <= ten2mxtrunc128[ind].w[0]))) {
  341.     // the result is a midpoint
  342.     if (Cstar.w[0] & 0x01) {    // Cstar is odd; MP in [EVEN, ODD]
  343.       // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
  344.       Cstar.w[0]--;     // Cstar is now even
  345.       if (Cstar.w[0] == 0xffffffffffffffffULL) {
  346.         Cstar.w[1]--;
  347.       }
  348.       *ptr_is_midpoint_gt_even = 1;
  349.       *ptr_is_inexact_lt_midpoint = 0;
  350.       *ptr_is_inexact_gt_midpoint = 0;
  351.     } else {    // else MP in [ODD, EVEN]
  352.       *ptr_is_midpoint_lt_even = 1;
  353.       *ptr_is_inexact_lt_midpoint = 0;
  354.       *ptr_is_inexact_gt_midpoint = 0;
  355.     }
  356.   }
  357.   // check for rounding overflow, which occurs if Cstar = 10^(q-x)
  358.   ind = q - x;  // 1 <= ind <= q - 1
  359.   if (ind <= 19) {
  360.     if (Cstar.w[1] == 0x0ULL && Cstar.w[0] == ten2k64[ind]) {
  361.       // if  Cstar = 10^(q-x)
  362.       Cstar.w[0] = ten2k64[ind - 1];    // Cstar = 10^(q-x-1)
  363.       *incr_exp = 1;
  364.     } else {
  365.       *incr_exp = 0;
  366.     }
  367.   } else if (ind == 20) {
  368.     // if ind = 20
  369.     if (Cstar.w[1] == ten2k128[0].w[1]
  370.         && Cstar.w[0] == ten2k128[0].w[0]) {
  371.       // if  Cstar = 10^(q-x)
  372.       Cstar.w[0] = ten2k64[19]; // Cstar = 10^(q-x-1)
  373.       Cstar.w[1] = 0x0ULL;
  374.       *incr_exp = 1;
  375.     } else {
  376.       *incr_exp = 0;
  377.     }
  378.   } else {      // if 21 <= ind <= 37
  379.     if (Cstar.w[1] == ten2k128[ind - 20].w[1] &&
  380.         Cstar.w[0] == ten2k128[ind - 20].w[0]) {
  381.       // if  Cstar = 10^(q-x)
  382.       Cstar.w[0] = ten2k128[ind - 21].w[0];     // Cstar = 10^(q-x-1)
  383.       Cstar.w[1] = ten2k128[ind - 21].w[1];
  384.       *incr_exp = 1;
  385.     } else {
  386.       *incr_exp = 0;
  387.     }
  388.   }
  389.   ptr_Cstar->w[1] = Cstar.w[1];
  390.   ptr_Cstar->w[0] = Cstar.w[0];
  391. }
  392.  
  393.  
  394. void
  395. round192_39_57 (int q,
  396.                 int x,
  397.                 UINT192 C,
  398.                 UINT192 * ptr_Cstar,
  399.                 int *incr_exp,
  400.                 int *ptr_is_midpoint_lt_even,
  401.                 int *ptr_is_midpoint_gt_even,
  402.                 int *ptr_is_inexact_lt_midpoint,
  403.                 int *ptr_is_inexact_gt_midpoint) {
  404.  
  405.   UINT384 P384;
  406.   UINT384 fstar;
  407.   UINT192 Cstar;
  408.   UINT64 tmp64;
  409.   int shift;
  410.   int ind;
  411.  
  412.   // Note:
  413.   //    In round192_39_57() positive numbers with 39 <= q <= 57 will be
  414.   //    rounded to nearest only for 5 <= x <= 42:
  415.   //     x = 23 or x = 24 or x = 5 or x = 6 when q = 39
  416.   //     x = 24 or x = 25 or x = 6 or x = 7 when q = 40
  417.   //     ...
  418.   //     x = 41 or x = 42 or x = 23 or x = 24 when q = 57
  419.   // However, for generality and possible uses outside the frame of IEEE 754R
  420.   // this implementation works for 1 <= x <= q - 1
  421.  
  422.   // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even,
  423.   // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are
  424.   // initialized to 0 by the caller
  425.  
  426.   // round a number C with q decimal digits, 39 <= q <= 57
  427.   // to q - x digits, 1 <= x <= 56
  428.   // C = C + 1/2 * 10^x where the result C fits in 192 bits
  429.   // (because the largest value is
  430.   // 999999999999999999999999999999999999999999999999999999999 +
  431.   //  50000000000000000000000000000000000000000000000000000000 =
  432.   // 0x2ad282f212a1da846afdaf18c034ff09da7fffffffffffff, which fits in 190 bits)
  433.   ind = x - 1;  // 0 <= ind <= 55
  434.   if (ind <= 18) {      // if 0 <= ind <= 18
  435.     tmp64 = C.w[0];
  436.     C.w[0] = C.w[0] + midpoint64[ind];
  437.     if (C.w[0] < tmp64) {
  438.       C.w[1]++;
  439.       if (C.w[1] == 0x0) {
  440.         C.w[2]++;
  441.       }
  442.     }
  443.   } else if (ind <= 37) {       // if 19 <= ind <= 37
  444.     tmp64 = C.w[0];
  445.     C.w[0] = C.w[0] + midpoint128[ind - 19].w[0];
  446.     if (C.w[0] < tmp64) {
  447.       C.w[1]++;
  448.       if (C.w[1] == 0x0) {
  449.         C.w[2]++;
  450.       }
  451.     }
  452.     tmp64 = C.w[1];
  453.     C.w[1] = C.w[1] + midpoint128[ind - 19].w[1];
  454.     if (C.w[1] < tmp64) {
  455.       C.w[2]++;
  456.     }
  457.   } else {      // if 38 <= ind <= 57 (actually ind <= 55)
  458.     tmp64 = C.w[0];
  459.     C.w[0] = C.w[0] + midpoint192[ind - 38].w[0];
  460.     if (C.w[0] < tmp64) {
  461.       C.w[1]++;
  462.       if (C.w[1] == 0x0ull) {
  463.         C.w[2]++;
  464.       }
  465.     }
  466.     tmp64 = C.w[1];
  467.     C.w[1] = C.w[1] + midpoint192[ind - 38].w[1];
  468.     if (C.w[1] < tmp64) {
  469.       C.w[2]++;
  470.     }
  471.     C.w[2] = C.w[2] + midpoint192[ind - 38].w[2];
  472.   }
  473.   // kx ~= 10^(-x), kx = Kx192[ind] * 2^(-Ex), 0 <= ind <= 55
  474.   // P384 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx
  475.   // the approximation kx of 10^(-x) was rounded up to 192 bits
  476.   __mul_192x192_to_384 (P384, C, Kx192[ind]);
  477.   // calculate C* = floor (P384) and f*
  478.   // Cstar = P384 >> Ex
  479.   // fstar = low Ex bits of P384
  480.   shift = Ex192m192[ind];       // in [1, 63] but have to consider three cases
  481.   if (ind <= 18) {      // if 0 <= ind <= 18
  482.     Cstar.w[2] = (P384.w[5] >> shift);
  483.     Cstar.w[1] = (P384.w[5] << (64 - shift)) | (P384.w[4] >> shift);
  484.     Cstar.w[0] = (P384.w[4] << (64 - shift)) | (P384.w[3] >> shift);
  485.     fstar.w[5] = 0x0ULL;
  486.     fstar.w[4] = 0x0ULL;
  487.     fstar.w[3] = P384.w[3] & mask192[ind];
  488.     fstar.w[2] = P384.w[2];
  489.     fstar.w[1] = P384.w[1];
  490.     fstar.w[0] = P384.w[0];
  491.   } else if (ind <= 37) {       // if 19 <= ind <= 37
  492.     Cstar.w[2] = 0x0ULL;
  493.     Cstar.w[1] = P384.w[5] >> shift;
  494.     Cstar.w[0] = (P384.w[5] << (64 - shift)) | (P384.w[4] >> shift);
  495.     fstar.w[5] = 0x0ULL;
  496.     fstar.w[4] = P384.w[4] & mask192[ind];
  497.     fstar.w[3] = P384.w[3];
  498.     fstar.w[2] = P384.w[2];
  499.     fstar.w[1] = P384.w[1];
  500.     fstar.w[0] = P384.w[0];
  501.   } else {      // if 38 <= ind <= 57
  502.     Cstar.w[2] = 0x0ULL;
  503.     Cstar.w[1] = 0x0ULL;
  504.     Cstar.w[0] = P384.w[5] >> shift;
  505.     fstar.w[5] = P384.w[5] & mask192[ind];
  506.     fstar.w[4] = P384.w[4];
  507.     fstar.w[3] = P384.w[3];
  508.     fstar.w[2] = P384.w[2];
  509.     fstar.w[1] = P384.w[1];
  510.     fstar.w[0] = P384.w[0];
  511.   }
  512.  
  513.   // the top Ex bits of 10^(-x) are T* = ten2mxtrunc192[ind], e.g. if x=1,
  514.   // T*=ten2mxtrunc192[0]=0xcccccccccccccccccccccccccccccccccccccccccccccccc
  515.   // if (0 < f* < 10^(-x)) then the result is a midpoint
  516.   //   if floor(C*) is even then C* = floor(C*) - logical right
  517.   //       shift; C* has q - x decimal digits, correct by Prop. 1)
  518.   //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
  519.   //       shift; C* has q - x decimal digits, correct by Pr. 1)
  520.   // else
  521.   //   C* = floor(C*) (logical right shift; C has q - x decimal digits,
  522.   //       correct by Property 1)
  523.   // in the caling function n = C* * 10^(e+x)
  524.  
  525.   // determine inexactness of the rounding of C*
  526.   // if (0 < f* - 1/2 < 10^(-x)) then
  527.   //   the result is exact
  528.   // else // if (f* - 1/2 > T*) then
  529.   //   the result is inexact
  530.   if (ind <= 18) {      // if 0 <= ind <= 18
  531.     if (fstar.w[3] > half192[ind] || (fstar.w[3] == half192[ind] &&
  532.                                       (fstar.w[2] || fstar.w[1]
  533.                                        || fstar.w[0]))) {
  534.       // f* > 1/2 and the result may be exact
  535.       // Calculate f* - 1/2
  536.       tmp64 = fstar.w[3] - half192[ind];
  537.       if (tmp64 || fstar.w[2] > ten2mxtrunc192[ind].w[2] || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] > ten2mxtrunc192[ind].w[1]) || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] == ten2mxtrunc192[ind].w[1] && fstar.w[0] > ten2mxtrunc192[ind].w[0])) { // f* - 1/2 > 10^(-x)
  538.         *ptr_is_inexact_lt_midpoint = 1;
  539.       } // else the result is exact
  540.     } else {    // the result is inexact; f2* <= 1/2
  541.       *ptr_is_inexact_gt_midpoint = 1;
  542.     }
  543.   } else if (ind <= 37) {       // if 19 <= ind <= 37
  544.     if (fstar.w[4] > half192[ind] || (fstar.w[4] == half192[ind] &&
  545.                                       (fstar.w[3] || fstar.w[2]
  546.                                        || fstar.w[1] || fstar.w[0]))) {
  547.       // f* > 1/2 and the result may be exact
  548.       // Calculate f* - 1/2
  549.       tmp64 = fstar.w[4] - half192[ind];
  550.       if (tmp64 || fstar.w[3] || fstar.w[2] > ten2mxtrunc192[ind].w[2] || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] > ten2mxtrunc192[ind].w[1]) || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] == ten2mxtrunc192[ind].w[1] && fstar.w[0] > ten2mxtrunc192[ind].w[0])) {   // f* - 1/2 > 10^(-x)
  551.         *ptr_is_inexact_lt_midpoint = 1;
  552.       } // else the result is exact
  553.     } else {    // the result is inexact; f2* <= 1/2
  554.       *ptr_is_inexact_gt_midpoint = 1;
  555.     }
  556.   } else {      // if 38 <= ind <= 55
  557.     if (fstar.w[5] > half192[ind] || (fstar.w[5] == half192[ind] &&
  558.                                       (fstar.w[4] || fstar.w[3]
  559.                                        || fstar.w[2] || fstar.w[1]
  560.                                        || fstar.w[0]))) {
  561.       // f* > 1/2 and the result may be exact
  562.       // Calculate f* - 1/2
  563.       tmp64 = fstar.w[5] - half192[ind];
  564.       if (tmp64 || fstar.w[4] || fstar.w[3] || fstar.w[2] > ten2mxtrunc192[ind].w[2] || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] > ten2mxtrunc192[ind].w[1]) || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] == ten2mxtrunc192[ind].w[1] && fstar.w[0] > ten2mxtrunc192[ind].w[0])) {     // f* - 1/2 > 10^(-x)
  565.         *ptr_is_inexact_lt_midpoint = 1;
  566.       } // else the result is exact
  567.     } else {    // the result is inexact; f2* <= 1/2
  568.       *ptr_is_inexact_gt_midpoint = 1;
  569.     }
  570.   }
  571.   // check for midpoints (could do this before determining inexactness)
  572.   if (fstar.w[5] == 0 && fstar.w[4] == 0 && fstar.w[3] == 0 &&
  573.       (fstar.w[2] < ten2mxtrunc192[ind].w[2] ||
  574.        (fstar.w[2] == ten2mxtrunc192[ind].w[2] &&
  575.         fstar.w[1] < ten2mxtrunc192[ind].w[1]) ||
  576.        (fstar.w[2] == ten2mxtrunc192[ind].w[2] &&
  577.         fstar.w[1] == ten2mxtrunc192[ind].w[1] &&
  578.         fstar.w[0] <= ten2mxtrunc192[ind].w[0]))) {
  579.     // the result is a midpoint
  580.     if (Cstar.w[0] & 0x01) {    // Cstar is odd; MP in [EVEN, ODD]
  581.       // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
  582.       Cstar.w[0]--;     // Cstar is now even
  583.       if (Cstar.w[0] == 0xffffffffffffffffULL) {
  584.         Cstar.w[1]--;
  585.         if (Cstar.w[1] == 0xffffffffffffffffULL) {
  586.           Cstar.w[2]--;
  587.         }
  588.       }
  589.       *ptr_is_midpoint_gt_even = 1;
  590.       *ptr_is_inexact_lt_midpoint = 0;
  591.       *ptr_is_inexact_gt_midpoint = 0;
  592.     } else {    // else MP in [ODD, EVEN]
  593.       *ptr_is_midpoint_lt_even = 1;
  594.       *ptr_is_inexact_lt_midpoint = 0;
  595.       *ptr_is_inexact_gt_midpoint = 0;
  596.     }
  597.   }
  598.   // check for rounding overflow, which occurs if Cstar = 10^(q-x)
  599.   ind = q - x;  // 1 <= ind <= q - 1
  600.   if (ind <= 19) {
  601.     if (Cstar.w[2] == 0x0ULL && Cstar.w[1] == 0x0ULL &&
  602.         Cstar.w[0] == ten2k64[ind]) {
  603.       // if  Cstar = 10^(q-x)
  604.       Cstar.w[0] = ten2k64[ind - 1];    // Cstar = 10^(q-x-1)
  605.       *incr_exp = 1;
  606.     } else {
  607.       *incr_exp = 0;
  608.     }
  609.   } else if (ind == 20) {
  610.     // if ind = 20
  611.     if (Cstar.w[2] == 0x0ULL && Cstar.w[1] == ten2k128[0].w[1] &&
  612.         Cstar.w[0] == ten2k128[0].w[0]) {
  613.       // if  Cstar = 10^(q-x)
  614.       Cstar.w[0] = ten2k64[19]; // Cstar = 10^(q-x-1)
  615.       Cstar.w[1] = 0x0ULL;
  616.       *incr_exp = 1;
  617.     } else {
  618.       *incr_exp = 0;
  619.     }
  620.   } else if (ind <= 38) {       // if 21 <= ind <= 38
  621.     if (Cstar.w[2] == 0x0ULL && Cstar.w[1] == ten2k128[ind - 20].w[1] &&
  622.         Cstar.w[0] == ten2k128[ind - 20].w[0]) {
  623.       // if  Cstar = 10^(q-x)
  624.       Cstar.w[0] = ten2k128[ind - 21].w[0];     // Cstar = 10^(q-x-1)
  625.       Cstar.w[1] = ten2k128[ind - 21].w[1];
  626.       *incr_exp = 1;
  627.     } else {
  628.       *incr_exp = 0;
  629.     }
  630.   } else if (ind == 39) {
  631.     if (Cstar.w[2] == ten2k256[0].w[2] && Cstar.w[1] == ten2k256[0].w[1]
  632.         && Cstar.w[0] == ten2k256[0].w[0]) {
  633.       // if  Cstar = 10^(q-x)
  634.       Cstar.w[0] = ten2k128[18].w[0];   // Cstar = 10^(q-x-1)
  635.       Cstar.w[1] = ten2k128[18].w[1];
  636.       Cstar.w[2] = 0x0ULL;
  637.       *incr_exp = 1;
  638.     } else {
  639.       *incr_exp = 0;
  640.     }
  641.   } else {      // if 40 <= ind <= 56
  642.     if (Cstar.w[2] == ten2k256[ind - 39].w[2] &&
  643.         Cstar.w[1] == ten2k256[ind - 39].w[1] &&
  644.         Cstar.w[0] == ten2k256[ind - 39].w[0]) {
  645.       // if  Cstar = 10^(q-x)
  646.       Cstar.w[0] = ten2k256[ind - 40].w[0];     // Cstar = 10^(q-x-1)
  647.       Cstar.w[1] = ten2k256[ind - 40].w[1];
  648.       Cstar.w[2] = ten2k256[ind - 40].w[2];
  649.       *incr_exp = 1;
  650.     } else {
  651.       *incr_exp = 0;
  652.     }
  653.   }
  654.   ptr_Cstar->w[2] = Cstar.w[2];
  655.   ptr_Cstar->w[1] = Cstar.w[1];
  656.   ptr_Cstar->w[0] = Cstar.w[0];
  657. }
  658.  
  659.  
  660. void
  661. round256_58_76 (int q,
  662.                 int x,
  663.                 UINT256 C,
  664.                 UINT256 * ptr_Cstar,
  665.                 int *incr_exp,
  666.                 int *ptr_is_midpoint_lt_even,
  667.                 int *ptr_is_midpoint_gt_even,
  668.                 int *ptr_is_inexact_lt_midpoint,
  669.                 int *ptr_is_inexact_gt_midpoint) {
  670.  
  671.   UINT512 P512;
  672.   UINT512 fstar;
  673.   UINT256 Cstar;
  674.   UINT64 tmp64;
  675.   int shift;
  676.   int ind;
  677.  
  678.   // Note:
  679.   //    In round256_58_76() positive numbers with 58 <= q <= 76 will be
  680.   //    rounded to nearest only for 24 <= x <= 61:
  681.   //     x = 42 or x = 43 or x = 24 or x = 25 when q = 58
  682.   //     x = 43 or x = 44 or x = 25 or x = 26 when q = 59
  683.   //     ...
  684.   //     x = 60 or x = 61 or x = 42 or x = 43 when q = 76
  685.   // However, for generality and possible uses outside the frame of IEEE 754R
  686.   // this implementation works for 1 <= x <= q - 1
  687.  
  688.   // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even,
  689.   // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are
  690.   // initialized to 0 by the caller
  691.  
  692.   // round a number C with q decimal digits, 58 <= q <= 76
  693.   // to q - x digits, 1 <= x <= 75
  694.   // C = C + 1/2 * 10^x where the result C fits in 256 bits
  695.   // (because the largest value is 9999999999999999999999999999999999999999
  696.   //     999999999999999999999999999999999999 + 500000000000000000000000000
  697.   //     000000000000000000000000000000000000000000000000 =
  698.   //     0x1736ca15d27a56cae15cf0e7b403d1f2bd6ebb0a50dc83ffffffffffffffffff,
  699.   // which fits in 253 bits)
  700.   ind = x - 1;  // 0 <= ind <= 74
  701.   if (ind <= 18) {      // if 0 <= ind <= 18
  702.     tmp64 = C.w[0];
  703.     C.w[0] = C.w[0] + midpoint64[ind];
  704.     if (C.w[0] < tmp64) {
  705.       C.w[1]++;
  706.       if (C.w[1] == 0x0) {
  707.         C.w[2]++;
  708.         if (C.w[2] == 0x0) {
  709.           C.w[3]++;
  710.         }
  711.       }
  712.     }
  713.   } else if (ind <= 37) {       // if 19 <= ind <= 37
  714.     tmp64 = C.w[0];
  715.     C.w[0] = C.w[0] + midpoint128[ind - 19].w[0];
  716.     if (C.w[0] < tmp64) {
  717.       C.w[1]++;
  718.       if (C.w[1] == 0x0) {
  719.         C.w[2]++;
  720.         if (C.w[2] == 0x0) {
  721.           C.w[3]++;
  722.         }
  723.       }
  724.     }
  725.     tmp64 = C.w[1];
  726.     C.w[1] = C.w[1] + midpoint128[ind - 19].w[1];
  727.     if (C.w[1] < tmp64) {
  728.       C.w[2]++;
  729.       if (C.w[2] == 0x0) {
  730.         C.w[3]++;
  731.       }
  732.     }
  733.   } else if (ind <= 57) {       // if 38 <= ind <= 57
  734.     tmp64 = C.w[0];
  735.     C.w[0] = C.w[0] + midpoint192[ind - 38].w[0];
  736.     if (C.w[0] < tmp64) {
  737.       C.w[1]++;
  738.       if (C.w[1] == 0x0ull) {
  739.         C.w[2]++;
  740.         if (C.w[2] == 0x0) {
  741.           C.w[3]++;
  742.         }
  743.       }
  744.     }
  745.     tmp64 = C.w[1];
  746.     C.w[1] = C.w[1] + midpoint192[ind - 38].w[1];
  747.     if (C.w[1] < tmp64) {
  748.       C.w[2]++;
  749.       if (C.w[2] == 0x0) {
  750.         C.w[3]++;
  751.       }
  752.     }
  753.     tmp64 = C.w[2];
  754.     C.w[2] = C.w[2] + midpoint192[ind - 38].w[2];
  755.     if (C.w[2] < tmp64) {
  756.       C.w[3]++;
  757.     }
  758.   } else {      // if 58 <= ind <= 76 (actually 58 <= ind <= 74)
  759.     tmp64 = C.w[0];
  760.     C.w[0] = C.w[0] + midpoint256[ind - 58].w[0];
  761.     if (C.w[0] < tmp64) {
  762.       C.w[1]++;
  763.       if (C.w[1] == 0x0ull) {
  764.         C.w[2]++;
  765.         if (C.w[2] == 0x0) {
  766.           C.w[3]++;
  767.         }
  768.       }
  769.     }
  770.     tmp64 = C.w[1];
  771.     C.w[1] = C.w[1] + midpoint256[ind - 58].w[1];
  772.     if (C.w[1] < tmp64) {
  773.       C.w[2]++;
  774.       if (C.w[2] == 0x0) {
  775.         C.w[3]++;
  776.       }
  777.     }
  778.     tmp64 = C.w[2];
  779.     C.w[2] = C.w[2] + midpoint256[ind - 58].w[2];
  780.     if (C.w[2] < tmp64) {
  781.       C.w[3]++;
  782.     }
  783.     C.w[3] = C.w[3] + midpoint256[ind - 58].w[3];
  784.   }
  785.   // kx ~= 10^(-x), kx = Kx256[ind] * 2^(-Ex), 0 <= ind <= 74
  786.   // P512 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx
  787.   // the approximation kx of 10^(-x) was rounded up to 192 bits
  788.   __mul_256x256_to_512 (P512, C, Kx256[ind]);
  789.   // calculate C* = floor (P512) and f*
  790.   // Cstar = P512 >> Ex
  791.   // fstar = low Ex bits of P512
  792.   shift = Ex256m256[ind];       // in [0, 63] but have to consider four cases
  793.   if (ind <= 18) {      // if 0 <= ind <= 18
  794.     Cstar.w[3] = (P512.w[7] >> shift);
  795.     Cstar.w[2] = (P512.w[7] << (64 - shift)) | (P512.w[6] >> shift);
  796.     Cstar.w[1] = (P512.w[6] << (64 - shift)) | (P512.w[5] >> shift);
  797.     Cstar.w[0] = (P512.w[5] << (64 - shift)) | (P512.w[4] >> shift);
  798.     fstar.w[7] = 0x0ULL;
  799.     fstar.w[6] = 0x0ULL;
  800.     fstar.w[5] = 0x0ULL;
  801.     fstar.w[4] = P512.w[4] & mask256[ind];
  802.     fstar.w[3] = P512.w[3];
  803.     fstar.w[2] = P512.w[2];
  804.     fstar.w[1] = P512.w[1];
  805.     fstar.w[0] = P512.w[0];
  806.   } else if (ind <= 37) {       // if 19 <= ind <= 37
  807.     Cstar.w[3] = 0x0ULL;
  808.     Cstar.w[2] = P512.w[7] >> shift;
  809.     Cstar.w[1] = (P512.w[7] << (64 - shift)) | (P512.w[6] >> shift);
  810.     Cstar.w[0] = (P512.w[6] << (64 - shift)) | (P512.w[5] >> shift);
  811.     fstar.w[7] = 0x0ULL;
  812.     fstar.w[6] = 0x0ULL;
  813.     fstar.w[5] = P512.w[5] & mask256[ind];
  814.     fstar.w[4] = P512.w[4];
  815.     fstar.w[3] = P512.w[3];
  816.     fstar.w[2] = P512.w[2];
  817.     fstar.w[1] = P512.w[1];
  818.     fstar.w[0] = P512.w[0];
  819.   } else if (ind <= 56) {       // if 38 <= ind <= 56
  820.     Cstar.w[3] = 0x0ULL;
  821.     Cstar.w[2] = 0x0ULL;
  822.     Cstar.w[1] = P512.w[7] >> shift;
  823.     Cstar.w[0] = (P512.w[7] << (64 - shift)) | (P512.w[6] >> shift);
  824.     fstar.w[7] = 0x0ULL;
  825.     fstar.w[6] = P512.w[6] & mask256[ind];
  826.     fstar.w[5] = P512.w[5];
  827.     fstar.w[4] = P512.w[4];
  828.     fstar.w[3] = P512.w[3];
  829.     fstar.w[2] = P512.w[2];
  830.     fstar.w[1] = P512.w[1];
  831.     fstar.w[0] = P512.w[0];
  832.   } else if (ind == 57) {
  833.     Cstar.w[3] = 0x0ULL;
  834.     Cstar.w[2] = 0x0ULL;
  835.     Cstar.w[1] = 0x0ULL;
  836.     Cstar.w[0] = P512.w[7];
  837.     fstar.w[7] = 0x0ULL;
  838.     fstar.w[6] = P512.w[6];
  839.     fstar.w[5] = P512.w[5];
  840.     fstar.w[4] = P512.w[4];
  841.     fstar.w[3] = P512.w[3];
  842.     fstar.w[2] = P512.w[2];
  843.     fstar.w[1] = P512.w[1];
  844.     fstar.w[0] = P512.w[0];
  845.   } else {      // if 58 <= ind <= 74
  846.     Cstar.w[3] = 0x0ULL;
  847.     Cstar.w[2] = 0x0ULL;
  848.     Cstar.w[1] = 0x0ULL;
  849.     Cstar.w[0] = P512.w[7] >> shift;
  850.     fstar.w[7] = P512.w[7] & mask256[ind];
  851.     fstar.w[6] = P512.w[6];
  852.     fstar.w[5] = P512.w[5];
  853.     fstar.w[4] = P512.w[4];
  854.     fstar.w[3] = P512.w[3];
  855.     fstar.w[2] = P512.w[2];
  856.     fstar.w[1] = P512.w[1];
  857.     fstar.w[0] = P512.w[0];
  858.   }
  859.  
  860.   // the top Ex bits of 10^(-x) are T* = ten2mxtrunc256[ind], e.g. if x=1,
  861.   // T*=ten2mxtrunc256[0]=
  862.   //     0xcccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  863.   // if (0 < f* < 10^(-x)) then the result is a midpoint
  864.   //   if floor(C*) is even then C* = floor(C*) - logical right
  865.   //       shift; C* has q - x decimal digits, correct by Prop. 1)
  866.   //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
  867.   //       shift; C* has q - x decimal digits, correct by Pr. 1)
  868.   // else
  869.   //   C* = floor(C*) (logical right shift; C has q - x decimal digits,
  870.   //       correct by Property 1)
  871.   // in the caling function n = C* * 10^(e+x)
  872.  
  873.   // determine inexactness of the rounding of C*
  874.   // if (0 < f* - 1/2 < 10^(-x)) then
  875.   //   the result is exact
  876.   // else // if (f* - 1/2 > T*) then
  877.   //   the result is inexact
  878.   if (ind <= 18) {      // if 0 <= ind <= 18
  879.     if (fstar.w[4] > half256[ind] || (fstar.w[4] == half256[ind] &&
  880.                                       (fstar.w[3] || fstar.w[2]
  881.                                        || fstar.w[1] || fstar.w[0]))) {
  882.       // f* > 1/2 and the result may be exact
  883.       // Calculate f* - 1/2
  884.       tmp64 = fstar.w[4] - half256[ind];
  885.       if (tmp64 || fstar.w[3] > ten2mxtrunc256[ind].w[2] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) {        // f* - 1/2 > 10^(-x)
  886.         *ptr_is_inexact_lt_midpoint = 1;
  887.       } // else the result is exact
  888.     } else {    // the result is inexact; f2* <= 1/2
  889.       *ptr_is_inexact_gt_midpoint = 1;
  890.     }
  891.   } else if (ind <= 37) {       // if 19 <= ind <= 37
  892.     if (fstar.w[5] > half256[ind] || (fstar.w[5] == half256[ind] &&
  893.                                       (fstar.w[4] || fstar.w[3]
  894.                                        || fstar.w[2] || fstar.w[1]
  895.                                        || fstar.w[0]))) {
  896.       // f* > 1/2 and the result may be exact
  897.       // Calculate f* - 1/2
  898.       tmp64 = fstar.w[5] - half256[ind];
  899.       if (tmp64 || fstar.w[4] || fstar.w[3] > ten2mxtrunc256[ind].w[3] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) {  // f* - 1/2 > 10^(-x)
  900.         *ptr_is_inexact_lt_midpoint = 1;
  901.       } // else the result is exact
  902.     } else {    // the result is inexact; f2* <= 1/2
  903.       *ptr_is_inexact_gt_midpoint = 1;
  904.     }
  905.   } else if (ind <= 57) {       // if 38 <= ind <= 57
  906.     if (fstar.w[6] > half256[ind] || (fstar.w[6] == half256[ind] &&
  907.                                       (fstar.w[5] || fstar.w[4]
  908.                                        || fstar.w[3] || fstar.w[2]
  909.                                        || fstar.w[1] || fstar.w[0]))) {
  910.       // f* > 1/2 and the result may be exact
  911.       // Calculate f* - 1/2
  912.       tmp64 = fstar.w[6] - half256[ind];
  913.       if (tmp64 || fstar.w[5] || fstar.w[4] || fstar.w[3] > ten2mxtrunc256[ind].w[3] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) {    // f* - 1/2 > 10^(-x)
  914.         *ptr_is_inexact_lt_midpoint = 1;
  915.       } // else the result is exact
  916.     } else {    // the result is inexact; f2* <= 1/2
  917.       *ptr_is_inexact_gt_midpoint = 1;
  918.     }
  919.   } else {      // if 58 <= ind <= 74
  920.     if (fstar.w[7] > half256[ind] || (fstar.w[7] == half256[ind] &&
  921.                                       (fstar.w[6] || fstar.w[5]
  922.                                        || fstar.w[4] || fstar.w[3]
  923.                                        || fstar.w[2] || fstar.w[1]
  924.                                        || fstar.w[0]))) {
  925.       // f* > 1/2 and the result may be exact
  926.       // Calculate f* - 1/2
  927.       tmp64 = fstar.w[7] - half256[ind];
  928.       if (tmp64 || fstar.w[6] || fstar.w[5] || fstar.w[4] || fstar.w[3] > ten2mxtrunc256[ind].w[3] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) {      // f* - 1/2 > 10^(-x)
  929.         *ptr_is_inexact_lt_midpoint = 1;
  930.       } // else the result is exact
  931.     } else {    // the result is inexact; f2* <= 1/2
  932.       *ptr_is_inexact_gt_midpoint = 1;
  933.     }
  934.   }
  935.   // check for midpoints (could do this before determining inexactness)
  936.   if (fstar.w[7] == 0 && fstar.w[6] == 0 &&
  937.       fstar.w[5] == 0 && fstar.w[4] == 0 &&
  938.       (fstar.w[3] < ten2mxtrunc256[ind].w[3] ||
  939.        (fstar.w[3] == ten2mxtrunc256[ind].w[3] &&
  940.         fstar.w[2] < ten2mxtrunc256[ind].w[2]) ||
  941.        (fstar.w[3] == ten2mxtrunc256[ind].w[3] &&
  942.         fstar.w[2] == ten2mxtrunc256[ind].w[2] &&
  943.         fstar.w[1] < ten2mxtrunc256[ind].w[1]) ||
  944.        (fstar.w[3] == ten2mxtrunc256[ind].w[3] &&
  945.         fstar.w[2] == ten2mxtrunc256[ind].w[2] &&
  946.         fstar.w[1] == ten2mxtrunc256[ind].w[1] &&
  947.         fstar.w[0] <= ten2mxtrunc256[ind].w[0]))) {
  948.     // the result is a midpoint
  949.     if (Cstar.w[0] & 0x01) {    // Cstar is odd; MP in [EVEN, ODD]
  950.       // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
  951.       Cstar.w[0]--;     // Cstar is now even
  952.       if (Cstar.w[0] == 0xffffffffffffffffULL) {
  953.         Cstar.w[1]--;
  954.         if (Cstar.w[1] == 0xffffffffffffffffULL) {
  955.           Cstar.w[2]--;
  956.           if (Cstar.w[2] == 0xffffffffffffffffULL) {
  957.             Cstar.w[3]--;
  958.           }
  959.         }
  960.       }
  961.       *ptr_is_midpoint_gt_even = 1;
  962.       *ptr_is_inexact_lt_midpoint = 0;
  963.       *ptr_is_inexact_gt_midpoint = 0;
  964.     } else {    // else MP in [ODD, EVEN]
  965.       *ptr_is_midpoint_lt_even = 1;
  966.       *ptr_is_inexact_lt_midpoint = 0;
  967.       *ptr_is_inexact_gt_midpoint = 0;
  968.     }
  969.   }
  970.   // check for rounding overflow, which occurs if Cstar = 10^(q-x)
  971.   ind = q - x;  // 1 <= ind <= q - 1
  972.   if (ind <= 19) {
  973.     if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == 0x0ULL &&
  974.         Cstar.w[1] == 0x0ULL && Cstar.w[0] == ten2k64[ind]) {
  975.       // if  Cstar = 10^(q-x)
  976.       Cstar.w[0] = ten2k64[ind - 1];    // Cstar = 10^(q-x-1)
  977.       *incr_exp = 1;
  978.     } else {
  979.       *incr_exp = 0;
  980.     }
  981.   } else if (ind == 20) {
  982.     // if ind = 20
  983.     if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == 0x0ULL &&
  984.         Cstar.w[1] == ten2k128[0].w[1]
  985.         && Cstar.w[0] == ten2k128[0].w[0]) {
  986.       // if  Cstar = 10^(q-x)
  987.       Cstar.w[0] = ten2k64[19]; // Cstar = 10^(q-x-1)
  988.       Cstar.w[1] = 0x0ULL;
  989.       *incr_exp = 1;
  990.     } else {
  991.       *incr_exp = 0;
  992.     }
  993.   } else if (ind <= 38) {       // if 21 <= ind <= 38
  994.     if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == 0x0ULL &&
  995.         Cstar.w[1] == ten2k128[ind - 20].w[1] &&
  996.         Cstar.w[0] == ten2k128[ind - 20].w[0]) {
  997.       // if  Cstar = 10^(q-x)
  998.       Cstar.w[0] = ten2k128[ind - 21].w[0];     // Cstar = 10^(q-x-1)
  999.       Cstar.w[1] = ten2k128[ind - 21].w[1];
  1000.       *incr_exp = 1;
  1001.     } else {
  1002.       *incr_exp = 0;
  1003.     }
  1004.   } else if (ind == 39) {
  1005.     if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == ten2k256[0].w[2] &&
  1006.         Cstar.w[1] == ten2k256[0].w[1]
  1007.         && Cstar.w[0] == ten2k256[0].w[0]) {
  1008.       // if  Cstar = 10^(q-x)
  1009.       Cstar.w[0] = ten2k128[18].w[0];   // Cstar = 10^(q-x-1)
  1010.       Cstar.w[1] = ten2k128[18].w[1];
  1011.       Cstar.w[2] = 0x0ULL;
  1012.       *incr_exp = 1;
  1013.     } else {
  1014.       *incr_exp = 0;
  1015.     }
  1016.   } else if (ind <= 57) {       // if 40 <= ind <= 57
  1017.     if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == ten2k256[ind - 39].w[2] &&
  1018.         Cstar.w[1] == ten2k256[ind - 39].w[1] &&
  1019.         Cstar.w[0] == ten2k256[ind - 39].w[0]) {
  1020.       // if  Cstar = 10^(q-x)
  1021.       Cstar.w[0] = ten2k256[ind - 40].w[0];     // Cstar = 10^(q-x-1)
  1022.       Cstar.w[1] = ten2k256[ind - 40].w[1];
  1023.       Cstar.w[2] = ten2k256[ind - 40].w[2];
  1024.       *incr_exp = 1;
  1025.     } else {
  1026.       *incr_exp = 0;
  1027.     }
  1028.     // else if (ind == 58) is not needed becauae we do not have ten2k192[] yet
  1029.   } else {      // if 58 <= ind <= 77 (actually 58 <= ind <= 74)
  1030.     if (Cstar.w[3] == ten2k256[ind - 39].w[3] &&
  1031.         Cstar.w[2] == ten2k256[ind - 39].w[2] &&
  1032.         Cstar.w[1] == ten2k256[ind - 39].w[1] &&
  1033.         Cstar.w[0] == ten2k256[ind - 39].w[0]) {
  1034.       // if  Cstar = 10^(q-x)
  1035.       Cstar.w[0] = ten2k256[ind - 40].w[0];     // Cstar = 10^(q-x-1)
  1036.       Cstar.w[1] = ten2k256[ind - 40].w[1];
  1037.       Cstar.w[2] = ten2k256[ind - 40].w[2];
  1038.       Cstar.w[3] = ten2k256[ind - 40].w[3];
  1039.       *incr_exp = 1;
  1040.     } else {
  1041.       *incr_exp = 0;
  1042.     }
  1043.   }
  1044.   ptr_Cstar->w[3] = Cstar.w[3];
  1045.   ptr_Cstar->w[2] = Cstar.w[2];
  1046.   ptr_Cstar->w[1] = Cstar.w[1];
  1047.   ptr_Cstar->w[0] = Cstar.w[0];
  1048.  
  1049. }
  1050.