libgcc/config/libbid/bid_round.c

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