1 /* Copyright (C) 2007-2018 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
round64_2_18(int q,int x,UINT64 C,UINT64 * ptr_Cstar,int * incr_exp,int * ptr_is_midpoint_lt_even,int * ptr_is_midpoint_gt_even,int * ptr_is_inexact_lt_midpoint,int * ptr_is_inexact_gt_midpoint)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
round128_19_38(int q,int x,UINT128 C,UINT128 * ptr_Cstar,int * incr_exp,int * ptr_is_midpoint_lt_even,int * ptr_is_midpoint_gt_even,int * ptr_is_inexact_lt_midpoint,int * ptr_is_inexact_gt_midpoint)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
round192_39_57(int q,int x,UINT192 C,UINT192 * ptr_Cstar,int * incr_exp,int * ptr_is_midpoint_lt_even,int * ptr_is_midpoint_gt_even,int * ptr_is_inexact_lt_midpoint,int * ptr_is_inexact_gt_midpoint)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
round256_58_76(int q,int x,UINT256 C,UINT256 * ptr_Cstar,int * incr_exp,int * ptr_is_midpoint_lt_even,int * ptr_is_midpoint_gt_even,int * ptr_is_inexact_lt_midpoint,int * ptr_is_inexact_gt_midpoint)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