1 /* Copyright (C) 2007-2021 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  *  BID128 fma   x * y + z
27  *
28  ****************************************************************************/
29 
30 #include "bid_internal.h"
31 
32 static void
rounding_correction(unsigned int rnd_mode,unsigned int is_inexact_lt_midpoint,unsigned int is_inexact_gt_midpoint,unsigned int is_midpoint_lt_even,unsigned int is_midpoint_gt_even,int unbexp,UINT128 * ptrres,_IDEC_flags * ptrfpsf)33 rounding_correction (unsigned int rnd_mode,
34         	     unsigned int is_inexact_lt_midpoint,
35         	     unsigned int is_inexact_gt_midpoint,
36         	     unsigned int is_midpoint_lt_even,
37         	     unsigned int is_midpoint_gt_even,
38         	     int unbexp,
39         	     UINT128 * ptrres, _IDEC_flags * ptrfpsf) {
40   // unbiased true exponent unbexp may be larger than emax
41 
42   UINT128 res = *ptrres; // expected to have the correct sign and coefficient
43   // (the exponent field is ignored, as unbexp is used instead)
44   UINT64 sign, exp;
45   UINT64 C_hi, C_lo;
46 
47   // general correction from RN to RA, RM, RP, RZ
48   // Note: if the result is negative, then is_inexact_lt_midpoint,
49   // is_inexact_gt_midpoint, is_midpoint_lt_even, and is_midpoint_gt_even
50   // have to be considered as if determined for the absolute value of the
51   // result (so they seem to be reversed)
52 
53   if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
54       is_midpoint_lt_even || is_midpoint_gt_even) {
55     *ptrfpsf |= INEXACT_EXCEPTION;
56   }
57   // apply correction to result calculated with unbounded exponent
58   sign = res.w[1] & MASK_SIGN;
59   exp = (UINT64) (unbexp + 6176) << 49; // valid only if expmin<=unbexp<=expmax
60   C_hi = res.w[1] & MASK_COEFF;
61   C_lo = res.w[0];
62   if ((!sign && ((rnd_mode == ROUNDING_UP && is_inexact_lt_midpoint) ||
63       ((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_UP) &&
64       is_midpoint_gt_even))) ||
65       (sign && ((rnd_mode == ROUNDING_DOWN && is_inexact_lt_midpoint) ||
66       ((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_DOWN) &&
67       is_midpoint_gt_even)))) {
68     // C = C + 1
69     C_lo = C_lo + 1;
70     if (C_lo == 0)
71       C_hi = C_hi + 1;
72     if (C_hi == 0x0001ed09bead87c0ull && C_lo == 0x378d8e6400000000ull) {
73       // C = 10^34 => rounding overflow
74       C_hi = 0x0000314dc6448d93ull;
75       C_lo = 0x38c15b0a00000000ull; // 10^33
76       // exp = exp + EXP_P1;
77       unbexp = unbexp + 1;
78       exp = (UINT64) (unbexp + 6176) << 49;
79     }
80   } else if ((is_midpoint_lt_even || is_inexact_gt_midpoint) &&
81       ((sign && (rnd_mode == ROUNDING_UP || rnd_mode == ROUNDING_TO_ZERO)) ||
82       (!sign && (rnd_mode == ROUNDING_DOWN || rnd_mode == ROUNDING_TO_ZERO)))) {
83     // C = C - 1
84     C_lo = C_lo - 1;
85     if (C_lo == 0xffffffffffffffffull)
86       C_hi--;
87     // check if we crossed into the lower decade
88     if (C_hi == 0x0000314dc6448d93ull && C_lo == 0x38c15b09ffffffffull) {
89       // C = 10^33 - 1
90       if (exp > 0) {
91         C_hi = 0x0001ed09bead87c0ull; // 10^34 - 1
92         C_lo = 0x378d8e63ffffffffull;
93         // exp = exp - EXP_P1;
94         unbexp = unbexp - 1;
95         exp = (UINT64) (unbexp + 6176) << 49;
96       } else { // if exp = 0 the result is tiny & inexact
97         *ptrfpsf |= UNDERFLOW_EXCEPTION;
98       }
99     }
100   } else {
101     ; // the result is already correct
102   }
103   if (unbexp > expmax) { // 6111
104     *ptrfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
105     exp = 0;
106     if (!sign) { // result is positive
107       if (rnd_mode == ROUNDING_UP || rnd_mode == ROUNDING_TIES_AWAY) { // +inf
108         C_hi = 0x7800000000000000ull;
109         C_lo = 0x0000000000000000ull;
110       } else { // res = +MAXFP = (10^34-1) * 10^emax
111         C_hi = 0x5fffed09bead87c0ull;
112         C_lo = 0x378d8e63ffffffffull;
113       }
114     } else { // result is negative
115       if (rnd_mode == ROUNDING_DOWN || rnd_mode == ROUNDING_TIES_AWAY) { // -inf
116         C_hi = 0xf800000000000000ull;
117         C_lo = 0x0000000000000000ull;
118       } else { // res = -MAXFP = -(10^34-1) * 10^emax
119         C_hi = 0xdfffed09bead87c0ull;
120         C_lo = 0x378d8e63ffffffffull;
121       }
122     }
123   }
124   // assemble the result
125   res.w[1] = sign | exp | C_hi;
126   res.w[0] = C_lo;
127   *ptrres = res;
128 }
129 
130 static void
add256(UINT256 x,UINT256 y,UINT256 * pz)131 add256 (UINT256 x, UINT256 y, UINT256 * pz) {
132   // *z = x + yl assume the sum fits in 256 bits
133   UINT256 z;
134   z.w[0] = x.w[0] + y.w[0];
135   if (z.w[0] < x.w[0]) {
136     x.w[1]++;
137     if (x.w[1] == 0x0000000000000000ull) {
138       x.w[2]++;
139       if (x.w[2] == 0x0000000000000000ull) {
140         x.w[3]++;
141       }
142     }
143   }
144   z.w[1] = x.w[1] + y.w[1];
145   if (z.w[1] < x.w[1]) {
146     x.w[2]++;
147     if (x.w[2] == 0x0000000000000000ull) {
148       x.w[3]++;
149     }
150   }
151   z.w[2] = x.w[2] + y.w[2];
152   if (z.w[2] < x.w[2]) {
153     x.w[3]++;
154   }
155   z.w[3] = x.w[3] + y.w[3]; // it was assumed that no carry is possible
156   *pz = z;
157 }
158 
159 static void
sub256(UINT256 x,UINT256 y,UINT256 * pz)160 sub256 (UINT256 x, UINT256 y, UINT256 * pz) {
161   // *z = x - y; assume x >= y
162   UINT256 z;
163   z.w[0] = x.w[0] - y.w[0];
164   if (z.w[0] > x.w[0]) {
165     x.w[1]--;
166     if (x.w[1] == 0xffffffffffffffffull) {
167       x.w[2]--;
168       if (x.w[2] == 0xffffffffffffffffull) {
169         x.w[3]--;
170       }
171     }
172   }
173   z.w[1] = x.w[1] - y.w[1];
174   if (z.w[1] > x.w[1]) {
175     x.w[2]--;
176     if (x.w[2] == 0xffffffffffffffffull) {
177       x.w[3]--;
178     }
179   }
180   z.w[2] = x.w[2] - y.w[2];
181   if (z.w[2] > x.w[2]) {
182     x.w[3]--;
183   }
184   z.w[3] = x.w[3] - y.w[3]; // no borrow possible, because x >= y
185   *pz = z;
186 }
187 
188 
189 static int
nr_digits256(UINT256 R256)190 nr_digits256 (UINT256 R256) {
191   int ind;
192   // determine the number of decimal digits in R256
193   if (R256.w[3] == 0x0 && R256.w[2] == 0x0 && R256.w[1] == 0x0) {
194     // between 1 and 19 digits
195     for (ind = 1; ind <= 19; ind++) {
196       if (R256.w[0] < ten2k64[ind]) {
197         break;
198       }
199     }
200     // ind digits
201   } else if (R256.w[3] == 0x0 && R256.w[2] == 0x0 &&
202              (R256.w[1] < ten2k128[0].w[1] ||
203               (R256.w[1] == ten2k128[0].w[1]
204                && R256.w[0] < ten2k128[0].w[0]))) {
205     // 20 digits
206     ind = 20;
207   } else if (R256.w[3] == 0x0 && R256.w[2] == 0x0) {
208     // between 21 and 38 digits
209     for (ind = 1; ind <= 18; ind++) {
210       if (R256.w[1] < ten2k128[ind].w[1] ||
211           (R256.w[1] == ten2k128[ind].w[1] &&
212            R256.w[0] < ten2k128[ind].w[0])) {
213         break;
214       }
215     }
216     // ind + 20 digits
217     ind = ind + 20;
218   } else if (R256.w[3] == 0x0 &&
219              (R256.w[2] < ten2k256[0].w[2] ||
220               (R256.w[2] == ten2k256[0].w[2] &&
221                R256.w[1] < ten2k256[0].w[1]) ||
222               (R256.w[2] == ten2k256[0].w[2] &&
223                R256.w[1] == ten2k256[0].w[1] &&
224                R256.w[0] < ten2k256[0].w[0]))) {
225     // 39 digits
226     ind = 39;
227   } else {
228     // between 40 and 68 digits
229     for (ind = 1; ind <= 29; ind++) {
230       if (R256.w[3] < ten2k256[ind].w[3] ||
231           (R256.w[3] == ten2k256[ind].w[3] &&
232            R256.w[2] < ten2k256[ind].w[2]) ||
233           (R256.w[3] == ten2k256[ind].w[3] &&
234            R256.w[2] == ten2k256[ind].w[2] &&
235            R256.w[1] < ten2k256[ind].w[1]) ||
236           (R256.w[3] == ten2k256[ind].w[3] &&
237            R256.w[2] == ten2k256[ind].w[2] &&
238            R256.w[1] == ten2k256[ind].w[1] &&
239            R256.w[0] < ten2k256[ind].w[0])) {
240         break;
241       }
242     }
243     // ind + 39 digits
244     ind = ind + 39;
245   }
246   return (ind);
247 }
248 
249 // add/subtract C4 and C3 * 10^scale; this may follow a previous rounding, so
250 // use the rounding information from ptr_is_* to avoid a double rounding error
251 static void
add_and_round(int q3,int q4,int e4,int delta,int p34,UINT64 z_sign,UINT64 p_sign,UINT128 C3,UINT256 C4,int rnd_mode,int * ptr_is_midpoint_lt_even,int * ptr_is_midpoint_gt_even,int * ptr_is_inexact_lt_midpoint,int * ptr_is_inexact_gt_midpoint,_IDEC_flags * ptrfpsf,UINT128 * ptrres)252 add_and_round (int q3,
253                int q4,
254                int e4,
255                int delta,
256                int p34,
257                UINT64 z_sign,
258                UINT64 p_sign,
259                UINT128 C3,
260                UINT256 C4,
261                int rnd_mode,
262                int *ptr_is_midpoint_lt_even,
263                int *ptr_is_midpoint_gt_even,
264                int *ptr_is_inexact_lt_midpoint,
265                int *ptr_is_inexact_gt_midpoint,
266                _IDEC_flags * ptrfpsf, UINT128 * ptrres) {
267 
268   int scale;
269   int x0;
270   int ind;
271   UINT64 R64;
272   UINT128 P128, R128;
273   UINT192 P192, R192;
274   UINT256 R256;
275   int is_midpoint_lt_even = 0;
276   int is_midpoint_gt_even = 0;
277   int is_inexact_lt_midpoint = 0;
278   int is_inexact_gt_midpoint = 0;
279   int is_midpoint_lt_even0 = 0;
280   int is_midpoint_gt_even0 = 0;
281   int is_inexact_lt_midpoint0 = 0;
282   int is_inexact_gt_midpoint0 = 0;
283   int incr_exp = 0;
284   int is_tiny = 0;
285   int lt_half_ulp = 0;
286   int eq_half_ulp = 0;
287   // int gt_half_ulp = 0;
288   UINT128 res = *ptrres;
289 
290   // scale C3 up by 10^(q4-delta-q3), 0 <= q4-delta-q3 <= 2*P34-2 = 66
291   scale = q4 - delta - q3; // 0 <= scale <= 66 (or 0 <= scale <= 68 if this
292   // comes from Cases (2), (3), (4), (5), (6), with 0 <= |delta| <= 1
293 
294   // calculate C3 * 10^scale in R256 (it has at most 67 decimal digits for
295   // Cases (15),(16),(17) and at most 69 for Cases (2),(3),(4),(5),(6))
296   if (scale == 0) {
297     R256.w[3] = 0x0ull;
298     R256.w[2] = 0x0ull;
299     R256.w[1] = C3.w[1];
300     R256.w[0] = C3.w[0];
301   } else if (scale <= 19) { // 10^scale fits in 64 bits
302     P128.w[1] = 0;
303     P128.w[0] = ten2k64[scale];
304     __mul_128x128_to_256 (R256, P128, C3);
305   } else if (scale <= 38) { // 10^scale fits in 128 bits
306     __mul_128x128_to_256 (R256, ten2k128[scale - 20], C3);
307   } else if (scale <= 57) { // 39 <= scale <= 57
308     // 10^scale fits in 192 bits but C3 * 10^scale fits in 223 or 230 bits
309     // (10^67 has 223 bits; 10^69 has 230 bits);
310     // must split the computation:
311     // 10^scale * C3 = 10*38 * 10^(scale-38) * C3 where 10^38 takes 127
312     // bits and so 10^(scale-38) * C3 fits in 128 bits with certainty
313     // Note that 1 <= scale - 38 <= 19 => 10^(scale-38) fits in 64 bits
314     __mul_64x128_to_128 (R128, ten2k64[scale - 38], C3);
315     // now multiply R128 by 10^38
316     __mul_128x128_to_256 (R256, R128, ten2k128[18]);
317   } else { // 58 <= scale <= 66
318     // 10^scale takes between 193 and 220 bits,
319     // and C3 * 10^scale fits in 223 bits (10^67/10^69 has 223/230 bits)
320     // must split the computation:
321     // 10^scale * C3 = 10*38 * 10^(scale-38) * C3 where 10^38 takes 127
322     // bits and so 10^(scale-38) * C3 fits in 128 bits with certainty
323     // Note that 20 <= scale - 38 <= 30 => 10^(scale-38) fits in 128 bits
324     // Calculate first 10^(scale-38) * C3, which fits in 128 bits; because
325     // 10^(scale-38) takes more than 64 bits, C3 will take less than 64
326     __mul_64x128_to_128 (R128, C3.w[0], ten2k128[scale - 58]);
327     // now calculate 10*38 * 10^(scale-38) * C3
328     __mul_128x128_to_256 (R256, R128, ten2k128[18]);
329   }
330   // C3 * 10^scale is now in R256
331 
332   // for Cases (15), (16), (17) C4 > C3 * 10^scale because C4 has at least
333   // one extra digit; for Cases (2), (3), (4), (5), or (6) any order is
334   // possible
335   // add/subtract C4 and C3 * 10^scale; the exponent is e4
336   if (p_sign == z_sign) { // R256 = C4 + R256
337     // calculate R256 = C4 + C3 * 10^scale = C4 + R256 which is exact,
338     // but may require rounding
339     add256 (C4, R256, &R256);
340   } else { // if (p_sign != z_sign) { // R256 = C4 - R256
341     // calculate R256 = C4 - C3 * 10^scale = C4 - R256 or
342     // R256 = C3 * 10^scale - C4 = R256 - C4 which is exact,
343     // but may require rounding
344 
345     // compare first R256 = C3 * 10^scale and C4
346     if (R256.w[3] > C4.w[3] || (R256.w[3] == C4.w[3] && R256.w[2] > C4.w[2]) ||
347         (R256.w[3] == C4.w[3] && R256.w[2] == C4.w[2] && R256.w[1] > C4.w[1]) ||
348         (R256.w[3] == C4.w[3] && R256.w[2] == C4.w[2] && R256.w[1] == C4.w[1] &&
349         R256.w[0] >= C4.w[0])) { // C3 * 10^scale >= C4
350       // calculate R256 = C3 * 10^scale - C4 = R256 - C4, which is exact,
351       // but may require rounding
352       sub256 (R256, C4, &R256);
353       // flip p_sign too, because the result has the sign of z
354       p_sign = z_sign;
355     } else { // if C4 > C3 * 10^scale
356       // calculate R256 = C4 - C3 * 10^scale = C4 - R256, which is exact,
357       // but may require rounding
358       sub256 (C4, R256, &R256);
359     }
360     // if the result is pure zero, the sign depends on the rounding mode
361     // (x*y and z had opposite signs)
362     if (R256.w[3] == 0x0ull && R256.w[2] == 0x0ull &&
363         R256.w[1] == 0x0ull && R256.w[0] == 0x0ull) {
364       if (rnd_mode != ROUNDING_DOWN)
365         p_sign = 0x0000000000000000ull;
366       else
367         p_sign = 0x8000000000000000ull;
368       // the exponent is max (e4, expmin)
369       if (e4 < -6176)
370         e4 = expmin;
371       // assemble result
372       res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49);
373       res.w[0] = 0x0;
374       *ptrres = res;
375       return;
376     }
377   }
378 
379   // determine the number of decimal digits in R256
380   ind = nr_digits256 (R256);
381 
382   // the exact result is (-1)^p_sign * R256 * 10^e4 where q (R256) = ind;
383   // round to the destination precision, with unbounded exponent
384 
385   if (ind <= p34) {
386     // result rounded to the destination precision with unbounded exponent
387     // is exact
388     if (ind + e4 < p34 + expmin) {
389       is_tiny = 1; // applies to all rounding modes
390     }
391     res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | R256.w[1];
392     res.w[0] = R256.w[0];
393     // Note: res is correct only if expmin <= e4 <= expmax
394   } else { // if (ind > p34)
395     // if more than P digits, round to nearest to P digits
396     // round R256 to p34 digits
397     x0 = ind - p34; // 1 <= x0 <= 34 as 35 <= ind <= 68
398     if (ind <= 38) {
399       P128.w[1] = R256.w[1];
400       P128.w[0] = R256.w[0];
401       round128_19_38 (ind, x0, P128, &R128, &incr_exp,
402         	      &is_midpoint_lt_even, &is_midpoint_gt_even,
403         	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
404     } else if (ind <= 57) {
405       P192.w[2] = R256.w[2];
406       P192.w[1] = R256.w[1];
407       P192.w[0] = R256.w[0];
408       round192_39_57 (ind, x0, P192, &R192, &incr_exp,
409         	      &is_midpoint_lt_even, &is_midpoint_gt_even,
410         	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
411       R128.w[1] = R192.w[1];
412       R128.w[0] = R192.w[0];
413     } else { // if (ind <= 68)
414       round256_58_76 (ind, x0, R256, &R256, &incr_exp,
415         	      &is_midpoint_lt_even, &is_midpoint_gt_even,
416         	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
417       R128.w[1] = R256.w[1];
418       R128.w[0] = R256.w[0];
419     }
420     // the rounded result has p34 = 34 digits
421     e4 = e4 + x0 + incr_exp;
422     if (rnd_mode == ROUNDING_TO_NEAREST) {
423       if (e4 < expmin) {
424         is_tiny = 1; // for other rounding modes apply correction
425       }
426     } else {
427       // for RM, RP, RZ, RA apply correction in order to determine tininess
428       // but do not save the result; apply the correction to
429       // (-1)^p_sign * significand * 10^0
430       P128.w[1] = p_sign | 0x3040000000000000ull | R128.w[1];
431       P128.w[0] = R128.w[0];
432       rounding_correction (rnd_mode,
433         		   is_inexact_lt_midpoint,
434         		   is_inexact_gt_midpoint, is_midpoint_lt_even,
435         		   is_midpoint_gt_even, 0, &P128, ptrfpsf);
436       scale = ((P128.w[1] & MASK_EXP) >> 49) - 6176; // -1, 0, or +1
437       // the number of digits in the significand is p34 = 34
438       if (e4 + scale < expmin) {
439         is_tiny = 1;
440       }
441     }
442     ind = p34; // the number of decimal digits in the signifcand of res
443     res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | R128.w[1]; // RN
444     res.w[0] = R128.w[0];
445     // Note: res is correct only if expmin <= e4 <= expmax
446     // set the inexact flag after rounding with bounded exponent, if any
447   }
448   // at this point we have the result rounded with unbounded exponent in
449   // res and we know its tininess:
450   // res = (-1)^p_sign * significand * 10^e4,
451   // where q (significand) = ind <= p34
452   // Note: res is correct only if expmin <= e4 <= expmax
453 
454   // check for overflow if RN
455   if (rnd_mode == ROUNDING_TO_NEAREST && (ind + e4) > (p34 + expmax)) {
456     res.w[1] = p_sign | 0x7800000000000000ull;
457     res.w[0] = 0x0000000000000000ull;
458     *ptrres = res;
459     *ptrfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
460     return; // BID_RETURN (res)
461   } // else not overflow or not RN, so continue
462 
463   // if (e4 >= expmin) we have the result rounded with bounded exponent
464   if (e4 < expmin) {
465     x0 = expmin - e4; // x0 >= 1; the number of digits to chop off of res
466     // where the result rounded [at most] once is
467     //   (-1)^p_sign * significand_res * 10^e4
468 
469     // avoid double rounding error
470     is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
471     is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
472     is_midpoint_lt_even0 = is_midpoint_lt_even;
473     is_midpoint_gt_even0 = is_midpoint_gt_even;
474     is_inexact_lt_midpoint = 0;
475     is_inexact_gt_midpoint = 0;
476     is_midpoint_lt_even = 0;
477     is_midpoint_gt_even = 0;
478 
479     if (x0 > ind) {
480       // nothing is left of res when moving the decimal point left x0 digits
481       is_inexact_lt_midpoint = 1;
482       res.w[1] = p_sign | 0x0000000000000000ull;
483       res.w[0] = 0x0000000000000000ull;
484       e4 = expmin;
485     } else if (x0 == ind) { // 1 <= x0 = ind <= p34 = 34
486       // this is <, =, or > 1/2 ulp
487       // compare the ind-digit value in the significand of res with
488       // 1/2 ulp = 5*10^(ind-1), i.e. determine whether it is
489       // less than, equal to, or greater than 1/2 ulp (significand of res)
490       R128.w[1] = res.w[1] & MASK_COEFF;
491       R128.w[0] = res.w[0];
492       if (ind <= 19) {
493         if (R128.w[0] < midpoint64[ind - 1]) { // < 1/2 ulp
494           lt_half_ulp = 1;
495           is_inexact_lt_midpoint = 1;
496         } else if (R128.w[0] == midpoint64[ind - 1]) { // = 1/2 ulp
497           eq_half_ulp = 1;
498           is_midpoint_gt_even = 1;
499         } else { // > 1/2 ulp
500           // gt_half_ulp = 1;
501           is_inexact_gt_midpoint = 1;
502         }
503       } else { // if (ind <= 38) {
504         if (R128.w[1] < midpoint128[ind - 20].w[1] ||
505             (R128.w[1] == midpoint128[ind - 20].w[1] &&
506             R128.w[0] < midpoint128[ind - 20].w[0])) { // < 1/2 ulp
507           lt_half_ulp = 1;
508           is_inexact_lt_midpoint = 1;
509         } else if (R128.w[1] == midpoint128[ind - 20].w[1] &&
510             R128.w[0] == midpoint128[ind - 20].w[0]) { // = 1/2 ulp
511           eq_half_ulp = 1;
512           is_midpoint_gt_even = 1;
513         } else { // > 1/2 ulp
514           // gt_half_ulp = 1;
515           is_inexact_gt_midpoint = 1;
516         }
517       }
518       if (lt_half_ulp || eq_half_ulp) {
519         // res = +0.0 * 10^expmin
520         res.w[1] = 0x0000000000000000ull;
521         res.w[0] = 0x0000000000000000ull;
522       } else { // if (gt_half_ulp)
523         // res = +1 * 10^expmin
524         res.w[1] = 0x0000000000000000ull;
525         res.w[0] = 0x0000000000000001ull;
526       }
527       res.w[1] = p_sign | res.w[1];
528       e4 = expmin;
529     } else { // if (1 <= x0 <= ind - 1 <= 33)
530       // round the ind-digit result to ind - x0 digits
531 
532       if (ind <= 18) { // 2 <= ind <= 18
533         round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
534         	      &is_midpoint_lt_even, &is_midpoint_gt_even,
535         	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
536         res.w[1] = 0x0;
537         res.w[0] = R64;
538       } else if (ind <= 38) {
539         P128.w[1] = res.w[1] & MASK_COEFF;
540         P128.w[0] = res.w[0];
541         round128_19_38 (ind, x0, P128, &res, &incr_exp,
542         		&is_midpoint_lt_even, &is_midpoint_gt_even,
543         		&is_inexact_lt_midpoint,
544         		&is_inexact_gt_midpoint);
545       }
546       e4 = e4 + x0; // expmin
547       // we want the exponent to be expmin, so if incr_exp = 1 then
548       // multiply the rounded result by 10 - it will still fit in 113 bits
549       if (incr_exp) {
550         // 64 x 128 -> 128
551         P128.w[1] = res.w[1] & MASK_COEFF;
552         P128.w[0] = res.w[0];
553         __mul_64x128_to_128 (res, ten2k64[1], P128);
554       }
555       res.w[1] =
556         p_sign | ((UINT64) (e4 + 6176) << 49) | (res.w[1] & MASK_COEFF);
557       // avoid a double rounding error
558       if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
559           is_midpoint_lt_even) { // double rounding error upward
560         // res = res - 1
561         res.w[0]--;
562         if (res.w[0] == 0xffffffffffffffffull)
563           res.w[1]--;
564         // Note: a double rounding error upward is not possible; for this
565         // the result after the first rounding would have to be 99...95
566         // (35 digits in all), possibly followed by a number of zeros; this
567         // is not possible in Cases (2)-(6) or (15)-(17) which may get here
568         is_midpoint_lt_even = 0;
569         is_inexact_lt_midpoint = 1;
570       } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
571           is_midpoint_gt_even) { // double rounding error downward
572         // res = res + 1
573         res.w[0]++;
574         if (res.w[0] == 0)
575           res.w[1]++;
576         is_midpoint_gt_even = 0;
577         is_inexact_gt_midpoint = 1;
578       } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
579         	 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
580         // if this second rounding was exact the result may still be
581         // inexact because of the first rounding
582         if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
583           is_inexact_gt_midpoint = 1;
584         }
585         if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
586           is_inexact_lt_midpoint = 1;
587         }
588       } else if (is_midpoint_gt_even &&
589         	 (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
590         // pulled up to a midpoint
591         is_inexact_lt_midpoint = 1;
592         is_inexact_gt_midpoint = 0;
593         is_midpoint_lt_even = 0;
594         is_midpoint_gt_even = 0;
595       } else if (is_midpoint_lt_even &&
596         	 (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
597         // pulled down to a midpoint
598         is_inexact_lt_midpoint = 0;
599         is_inexact_gt_midpoint = 1;
600         is_midpoint_lt_even = 0;
601         is_midpoint_gt_even = 0;
602       } else {
603         ;
604       }
605     }
606   }
607   // res contains the correct result
608   // apply correction if not rounding to nearest
609   if (rnd_mode != ROUNDING_TO_NEAREST) {
610     rounding_correction (rnd_mode,
611         		 is_inexact_lt_midpoint, is_inexact_gt_midpoint,
612         		 is_midpoint_lt_even, is_midpoint_gt_even,
613         		 e4, &res, ptrfpsf);
614   }
615   if (is_midpoint_lt_even || is_midpoint_gt_even ||
616       is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
617     // set the inexact flag
618     *ptrfpsf |= INEXACT_EXCEPTION;
619     if (is_tiny)
620       *ptrfpsf |= UNDERFLOW_EXCEPTION;
621   }
622 
623   *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
624   *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
625   *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
626   *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
627   *ptrres = res;
628   return;
629 }
630 
631 
632 #if DECIMAL_CALL_BY_REFERENCE
633 static void
bid128_ext_fma(int * ptr_is_midpoint_lt_even,int * ptr_is_midpoint_gt_even,int * ptr_is_inexact_lt_midpoint,int * ptr_is_inexact_gt_midpoint,UINT128 * pres,UINT128 * px,UINT128 * py,UINT128 * pz _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)634 bid128_ext_fma (int *ptr_is_midpoint_lt_even,
635         	int *ptr_is_midpoint_gt_even,
636         	int *ptr_is_inexact_lt_midpoint,
637         	int *ptr_is_inexact_gt_midpoint, UINT128 * pres,
638         	UINT128 * px, UINT128 * py,
639         	UINT128 *
640         	pz _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
641         	_EXC_INFO_PARAM) {
642   UINT128 x = *px, y = *py, z = *pz;
643 #if !DECIMAL_GLOBAL_ROUNDING
644   unsigned int rnd_mode = *prnd_mode;
645 #endif
646 #else
647 static UINT128
648 bid128_ext_fma (int *ptr_is_midpoint_lt_even,
649         	int *ptr_is_midpoint_gt_even,
650         	int *ptr_is_inexact_lt_midpoint,
651         	int *ptr_is_inexact_gt_midpoint, UINT128 x, UINT128 y,
652         	UINT128 z _RND_MODE_PARAM _EXC_FLAGS_PARAM
653         	_EXC_MASKS_PARAM _EXC_INFO_PARAM) {
654 #endif
655 
656   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
657   UINT64 x_sign, y_sign, z_sign, p_sign, tmp_sign;
658   UINT64 x_exp = 0, y_exp = 0, z_exp = 0, p_exp;
659   int true_p_exp;
660   UINT128 C1, C2, C3;
661   UINT256 C4;
662   int q1 = 0, q2 = 0, q3 = 0, q4;
663   int e1, e2, e3, e4;
664   int scale, ind, delta, x0;
665   int p34 = P34; // used to modify the limit on the number of digits
666   BID_UI64DOUBLE tmp;
667   int x_nr_bits, y_nr_bits, z_nr_bits;
668   unsigned int save_fpsf;
669   int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0;
670   int is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
671   int is_midpoint_lt_even0 = 0, is_midpoint_gt_even0 = 0;
672   int is_inexact_lt_midpoint0 = 0, is_inexact_gt_midpoint0 = 0;
673   int incr_exp = 0;
674   int lsb;
675   int lt_half_ulp = 0;
676   int eq_half_ulp = 0;
677   int gt_half_ulp = 0;
678   int is_tiny = 0;
679   UINT64 R64, tmp64;
680   UINT128 P128, R128;
681   UINT192 P192, R192;
682   UINT256 R256;
683 
684   // the following are based on the table of special cases for fma; the NaN
685   // behavior is similar to that of the IA-64 Architecture fma
686 
687   // identify cases where at least one operand is NaN
688 
689   BID_SWAP128 (x);
690   BID_SWAP128 (y);
691   BID_SWAP128 (z);
692   if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
693     // if x = {0, f, inf, NaN}, y = NaN, z = {0, f, inf, NaN} then res = Q (y)
694     // check first for non-canonical NaN payload
695     if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
696         (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
697          (y.w[0] > 0x38c15b09ffffffffull))) {
698       y.w[1] = y.w[1] & 0xffffc00000000000ull;
699       y.w[0] = 0x0ull;
700     }
701     if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
702       // set invalid flag
703       *pfpsf |= INVALID_EXCEPTION;
704       // return quiet (y)
705       res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
706       res.w[0] = y.w[0];
707     } else { // y is QNaN
708       // return y
709       res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
710       res.w[0] = y.w[0];
711       // if z = SNaN or x = SNaN signal invalid exception
712       if ((z.w[1] & MASK_SNAN) == MASK_SNAN ||
713           (x.w[1] & MASK_SNAN) == MASK_SNAN) {
714         // set invalid flag
715         *pfpsf |= INVALID_EXCEPTION;
716       }
717     }
718     *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
719     *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
720     *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
721     *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
722     BID_SWAP128 (res);
723     BID_RETURN (res)
724   } else if ((z.w[1] & MASK_NAN) == MASK_NAN) { // z is NAN
725     // if x = {0, f, inf, NaN}, y = {0, f, inf}, z = NaN then res = Q (z)
726     // check first for non-canonical NaN payload
727     if (((z.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
728         (((z.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
729          (z.w[0] > 0x38c15b09ffffffffull))) {
730       z.w[1] = z.w[1] & 0xffffc00000000000ull;
731       z.w[0] = 0x0ull;
732     }
733     if ((z.w[1] & MASK_SNAN) == MASK_SNAN) { // z is SNAN
734       // set invalid flag
735       *pfpsf |= INVALID_EXCEPTION;
736       // return quiet (z)
737       res.w[1] = z.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
738       res.w[0] = z.w[0];
739     } else { // z is QNaN
740       // return z
741       res.w[1] = z.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
742       res.w[0] = z.w[0];
743       // if x = SNaN signal invalid exception
744       if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {
745         // set invalid flag
746         *pfpsf |= INVALID_EXCEPTION;
747       }
748     }
749     *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
750     *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
751     *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
752     *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
753     BID_SWAP128 (res);
754     BID_RETURN (res)
755   } else if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
756     // if x = NaN, y = {0, f, inf}, z = {0, f, inf} then res = Q (x)
757     // check first for non-canonical NaN payload
758     if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
759         (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
760          (x.w[0] > 0x38c15b09ffffffffull))) {
761       x.w[1] = x.w[1] & 0xffffc00000000000ull;
762       x.w[0] = 0x0ull;
763     }
764     if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
765       // set invalid flag
766       *pfpsf |= INVALID_EXCEPTION;
767       // return quiet (x)
768       res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
769       res.w[0] = x.w[0];
770     } else { // x is QNaN
771       // return x
772       res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
773       res.w[0] = x.w[0];
774     }
775     *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
776     *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
777     *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
778     *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
779     BID_SWAP128 (res);
780     BID_RETURN (res)
781   }
782   // x, y, z are 0, f, or inf but not NaN => unpack the arguments and check
783   // for non-canonical values
784 
785   x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
786   C1.w[1] = x.w[1] & MASK_COEFF;
787   C1.w[0] = x.w[0];
788   if ((x.w[1] & MASK_ANY_INF) != MASK_INF) { // x != inf
789     // if x is not infinity check for non-canonical values - treated as zero
790     if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
791       // non-canonical
792       x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
793       C1.w[1] = 0; // significand high
794       C1.w[0] = 0; // significand low
795     } else { // G0_G1 != 11
796       x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
797       if (C1.w[1] > 0x0001ed09bead87c0ull ||
798           (C1.w[1] == 0x0001ed09bead87c0ull &&
799            C1.w[0] > 0x378d8e63ffffffffull)) {
800         // x is non-canonical if coefficient is larger than 10^34 -1
801         C1.w[1] = 0;
802         C1.w[0] = 0;
803       } else { // canonical
804         ;
805       }
806     }
807   }
808   y_sign = y.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
809   C2.w[1] = y.w[1] & MASK_COEFF;
810   C2.w[0] = y.w[0];
811   if ((y.w[1] & MASK_ANY_INF) != MASK_INF) { // y != inf
812     // if y is not infinity check for non-canonical values - treated as zero
813     if ((y.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
814       // non-canonical
815       y_exp = (y.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
816       C2.w[1] = 0; // significand high
817       C2.w[0] = 0; // significand low
818     } else { // G0_G1 != 11
819       y_exp = y.w[1] & MASK_EXP; // biased and shifted left 49 bits
820       if (C2.w[1] > 0x0001ed09bead87c0ull ||
821           (C2.w[1] == 0x0001ed09bead87c0ull &&
822            C2.w[0] > 0x378d8e63ffffffffull)) {
823         // y is non-canonical if coefficient is larger than 10^34 -1
824         C2.w[1] = 0;
825         C2.w[0] = 0;
826       } else { // canonical
827         ;
828       }
829     }
830   }
831   z_sign = z.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
832   C3.w[1] = z.w[1] & MASK_COEFF;
833   C3.w[0] = z.w[0];
834   if ((z.w[1] & MASK_ANY_INF) != MASK_INF) { // z != inf
835     // if z is not infinity check for non-canonical values - treated as zero
836     if ((z.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
837       // non-canonical
838       z_exp = (z.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
839       C3.w[1] = 0; // significand high
840       C3.w[0] = 0; // significand low
841     } else { // G0_G1 != 11
842       z_exp = z.w[1] & MASK_EXP; // biased and shifted left 49 bits
843       if (C3.w[1] > 0x0001ed09bead87c0ull ||
844           (C3.w[1] == 0x0001ed09bead87c0ull &&
845            C3.w[0] > 0x378d8e63ffffffffull)) {
846         // z is non-canonical if coefficient is larger than 10^34 -1
847         C3.w[1] = 0;
848         C3.w[0] = 0;
849       } else { // canonical
850         ;
851       }
852     }
853   }
854 
855   p_sign = x_sign ^ y_sign; // sign of the product
856 
857   // identify cases where at least one operand is infinity
858 
859   if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf
860     if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
861       if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf
862         if (p_sign == z_sign) {
863           res.w[1] = z_sign | MASK_INF;
864           res.w[0] = 0x0;
865         } else {
866           // return QNaN Indefinite
867           res.w[1] = 0x7c00000000000000ull;
868           res.w[0] = 0x0000000000000000ull;
869           // set invalid flag
870           *pfpsf |= INVALID_EXCEPTION;
871         }
872       } else { // z = 0 or z = f
873         res.w[1] = p_sign | MASK_INF;
874         res.w[0] = 0x0;
875       }
876     } else if (C2.w[1] != 0 || C2.w[0] != 0) { // y = f
877       if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf
878         if (p_sign == z_sign) {
879           res.w[1] = z_sign | MASK_INF;
880           res.w[0] = 0x0;
881         } else {
882           // return QNaN Indefinite
883           res.w[1] = 0x7c00000000000000ull;
884           res.w[0] = 0x0000000000000000ull;
885           // set invalid flag
886           *pfpsf |= INVALID_EXCEPTION;
887         }
888       } else { // z = 0 or z = f
889         res.w[1] = p_sign | MASK_INF;
890         res.w[0] = 0x0;
891       }
892     } else { // y = 0
893       // return QNaN Indefinite
894       res.w[1] = 0x7c00000000000000ull;
895       res.w[0] = 0x0000000000000000ull;
896       // set invalid flag
897       *pfpsf |= INVALID_EXCEPTION;
898     }
899     *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
900     *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
901     *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
902     *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
903     BID_SWAP128 (res);
904     BID_RETURN (res)
905   } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
906     if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf
907       // x = f, necessarily
908       if ((p_sign != z_sign)
909           || (C1.w[1] == 0x0ull && C1.w[0] == 0x0ull)) {
910         // return QNaN Indefinite
911         res.w[1] = 0x7c00000000000000ull;
912         res.w[0] = 0x0000000000000000ull;
913         // set invalid flag
914         *pfpsf |= INVALID_EXCEPTION;
915       } else {
916         res.w[1] = z_sign | MASK_INF;
917         res.w[0] = 0x0;
918       }
919     } else if (C1.w[1] == 0x0 && C1.w[0] == 0x0) { // x = 0
920       // z = 0, f, inf
921       // return QNaN Indefinite
922       res.w[1] = 0x7c00000000000000ull;
923       res.w[0] = 0x0000000000000000ull;
924       // set invalid flag
925       *pfpsf |= INVALID_EXCEPTION;
926     } else {
927       // x = f and z = 0, f, necessarily
928       res.w[1] = p_sign | MASK_INF;
929       res.w[0] = 0x0;
930     }
931     *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
932     *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
933     *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
934     *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
935     BID_SWAP128 (res);
936     BID_RETURN (res)
937   } else if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf
938     // x = 0, f and y = 0, f, necessarily
939     res.w[1] = z_sign | MASK_INF;
940     res.w[0] = 0x0;
941     *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
942     *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
943     *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
944     *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
945     BID_SWAP128 (res);
946     BID_RETURN (res)
947   }
948 
949   true_p_exp = (x_exp >> 49) - 6176 + (y_exp >> 49) - 6176;
950   if (true_p_exp < -6176)
951     p_exp = 0; // cannot be less than EXP_MIN
952   else
953     p_exp = (UINT64) (true_p_exp + 6176) << 49;
954 
955   if (((C1.w[1] == 0x0 && C1.w[0] == 0x0) || (C2.w[1] == 0x0 && C2.w[0] == 0x0)) && C3.w[1] == 0x0 && C3.w[0] == 0x0) { // (x = 0 or y = 0) and z = 0
956     // the result is 0
957     if (p_exp < z_exp)
958       res.w[1] = p_exp; // preferred exponent
959     else
960       res.w[1] = z_exp; // preferred exponent
961     if (p_sign == z_sign) {
962       res.w[1] |= z_sign;
963       res.w[0] = 0x0;
964     } else { // x * y and z have opposite signs
965       if (rnd_mode == ROUNDING_DOWN) {
966         // res = -0.0
967         res.w[1] |= MASK_SIGN;
968         res.w[0] = 0x0;
969       } else {
970         // res = +0.0
971         // res.w[1] |= 0x0;
972         res.w[0] = 0x0;
973       }
974     }
975     *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
976     *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
977     *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
978     *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
979     BID_SWAP128 (res);
980     BID_RETURN (res)
981   }
982   // from this point on, we may need to know the number of decimal digits
983   // in the significands of x, y, z when x, y, z != 0
984 
985   if (C1.w[1] != 0 || C1.w[0] != 0) { // x = f (non-zero finite)
986     // q1 = nr. of decimal digits in x
987     // determine first the nr. of bits in x
988     if (C1.w[1] == 0) {
989       if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
990         // split the 64-bit value in two 32-bit halves to avoid rounding errors
991         if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
992           tmp.d = (double) (C1.w[0] >> 32); // exact conversion
993           x_nr_bits =
994             33 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
995         } else { // x < 2^32
996           tmp.d = (double) (C1.w[0]); // exact conversion
997           x_nr_bits =
998             1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
999         }
1000       } else { // if x < 2^53
1001         tmp.d = (double) C1.w[0]; // exact conversion
1002         x_nr_bits =
1003           1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
1004       }
1005     } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1006       tmp.d = (double) C1.w[1]; // exact conversion
1007       x_nr_bits =
1008         65 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
1009     }
1010     q1 = nr_digits[x_nr_bits - 1].digits;
1011     if (q1 == 0) {
1012       q1 = nr_digits[x_nr_bits - 1].digits1;
1013       if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1014           (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1015            C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1016         q1++;
1017     }
1018   }
1019 
1020   if (C2.w[1] != 0 || C2.w[0] != 0) { // y = f (non-zero finite)
1021     if (C2.w[1] == 0) {
1022       if (C2.w[0] >= 0x0020000000000000ull) { // y >= 2^53
1023         // split the 64-bit value in two 32-bit halves to avoid rounding errors
1024         if (C2.w[0] >= 0x0000000100000000ull) { // y >= 2^32
1025           tmp.d = (double) (C2.w[0] >> 32); // exact conversion
1026           y_nr_bits =
1027             32 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
1028         } else { // y < 2^32
1029           tmp.d = (double) C2.w[0]; // exact conversion
1030           y_nr_bits =
1031             ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
1032         }
1033       } else { // if y < 2^53
1034         tmp.d = (double) C2.w[0]; // exact conversion
1035         y_nr_bits =
1036           ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
1037       }
1038     } else { // C2.w[1] != 0 => nr. bits = 64 + nr_bits (C2.w[1])
1039       tmp.d = (double) C2.w[1]; // exact conversion
1040       y_nr_bits =
1041         64 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
1042     }
1043 
1044     q2 = nr_digits[y_nr_bits].digits;
1045     if (q2 == 0) {
1046       q2 = nr_digits[y_nr_bits].digits1;
1047       if (C2.w[1] > nr_digits[y_nr_bits].threshold_hi ||
1048           (C2.w[1] == nr_digits[y_nr_bits].threshold_hi &&
1049            C2.w[0] >= nr_digits[y_nr_bits].threshold_lo))
1050         q2++;
1051     }
1052   }
1053 
1054   if (C3.w[1] != 0 || C3.w[0] != 0) { // z = f (non-zero finite)
1055     if (C3.w[1] == 0) {
1056       if (C3.w[0] >= 0x0020000000000000ull) { // z >= 2^53
1057         // split the 64-bit value in two 32-bit halves to avoid rounding errors
1058         if (C3.w[0] >= 0x0000000100000000ull) { // z >= 2^32
1059           tmp.d = (double) (C3.w[0] >> 32); // exact conversion
1060           z_nr_bits =
1061             32 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
1062         } else { // z < 2^32
1063           tmp.d = (double) C3.w[0]; // exact conversion
1064           z_nr_bits =
1065             ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
1066         }
1067       } else { // if z < 2^53
1068         tmp.d = (double) C3.w[0]; // exact conversion
1069         z_nr_bits =
1070           ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
1071       }
1072     } else { // C3.w[1] != 0 => nr. bits = 64 + nr_bits (C3.w[1])
1073       tmp.d = (double) C3.w[1]; // exact conversion
1074       z_nr_bits =
1075         64 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
1076     }
1077 
1078     q3 = nr_digits[z_nr_bits].digits;
1079     if (q3 == 0) {
1080       q3 = nr_digits[z_nr_bits].digits1;
1081       if (C3.w[1] > nr_digits[z_nr_bits].threshold_hi ||
1082           (C3.w[1] == nr_digits[z_nr_bits].threshold_hi &&
1083            C3.w[0] >= nr_digits[z_nr_bits].threshold_lo))
1084         q3++;
1085     }
1086   }
1087 
1088   if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
1089       (C2.w[1] == 0x0 && C2.w[0] == 0x0)) {
1090     // x = 0 or y = 0
1091     // z = f, necessarily; for 0 + z return z, with the preferred exponent
1092     // the result is z, but need to get the preferred exponent
1093     if (z_exp <= p_exp) { // the preferred exponent is z_exp
1094       res.w[1] = z_sign | (z_exp & MASK_EXP) | C3.w[1];
1095       res.w[0] = C3.w[0];
1096     } else { // if (p_exp < z_exp) the preferred exponent is p_exp
1097       // return (C3 * 10^scale) * 10^(z_exp - scale)
1098       // where scale = min (p34-q3, (z_exp-p_exp) >> 49)
1099       scale = p34 - q3;
1100       ind = (z_exp - p_exp) >> 49;
1101       if (ind < scale)
1102         scale = ind;
1103       if (scale == 0) {
1104         res.w[1] = z.w[1]; // & MASK_COEFF, which is redundant
1105         res.w[0] = z.w[0];
1106       } else if (q3 <= 19) { // z fits in 64 bits
1107         if (scale <= 19) { // 10^scale fits in 64 bits
1108           // 64 x 64 C3.w[0] * ten2k64[scale]
1109           __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
1110         } else { // 10^scale fits in 128 bits
1111           // 64 x 128 C3.w[0] * ten2k128[scale - 20]
1112           __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
1113         }
1114       } else { // z fits in 128 bits, but 10^scale must fit in 64 bits
1115         // 64 x 128 ten2k64[scale] * C3
1116         __mul_128x64_to_128 (res, ten2k64[scale], C3);
1117       }
1118       // subtract scale from the exponent
1119       z_exp = z_exp - ((UINT64) scale << 49);
1120       res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
1121     }
1122     *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
1123     *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
1124     *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
1125     *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
1126     BID_SWAP128 (res);
1127     BID_RETURN (res)
1128   } else {
1129     ; // continue with x = f, y = f, z = 0 or x = f, y = f, z = f
1130   }
1131 
1132   e1 = (x_exp >> 49) - 6176; // unbiased exponent of x
1133   e2 = (y_exp >> 49) - 6176; // unbiased exponent of y
1134   e3 = (z_exp >> 49) - 6176; // unbiased exponent of z
1135   e4 = e1 + e2; // unbiased exponent of the exact x * y
1136 
1137   // calculate C1 * C2 and its number of decimal digits, q4
1138 
1139   // the exact product has either q1 + q2 - 1 or q1 + q2 decimal digits
1140   // where 2 <= q1 + q2 <= 68
1141   // calculate C4 = C1 * C2 and determine q
1142   C4.w[3] = C4.w[2] = C4.w[1] = C4.w[0] = 0;
1143   if (q1 + q2 <= 19) { // if 2 <= q1 + q2 <= 19, C4 = C1 * C2 fits in 64 bits
1144     C4.w[0] = C1.w[0] * C2.w[0];
1145     // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
1146     if (C4.w[0] < ten2k64[q1 + q2 - 1])
1147       q4 = q1 + q2 - 1; // q4 in [1, 18]
1148     else
1149       q4 = q1 + q2; // q4 in [2, 19]
1150     // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
1151   } else if (q1 + q2 == 20) { // C4 = C1 * C2 fits in 64 or 128 bits
1152     // q1 <= 19 and q2 <= 19 so both C1 and C2 fit in 64 bits
1153     __mul_64x64_to_128MACH (C4, C1.w[0], C2.w[0]);
1154     // if C4 < 10^(q1+q2-1) = 10^19 then q4 = q1+q2-1 = 19 else q4 = q1+q2 = 20
1155     if (C4.w[1] == 0 && C4.w[0] < ten2k64[19]) { // 19 = q1+q2-1
1156       // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
1157       q4 = 19; // 19 = q1 + q2 - 1
1158     } else {
1159       // if (C4.w[1] == 0)
1160       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
1161       // else
1162       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
1163       q4 = 20; // 20 = q1 + q2
1164     }
1165   } else if (q1 + q2 <= 38) { // 21 <= q1 + q2 <= 38
1166     // C4 = C1 * C2 fits in 64 or 128 bits
1167     // (64 bits possibly, but only when q1 + q2 = 21 and C4 has 20 digits)
1168     // at least one of C1, C2 has at most 19 decimal digits & fits in 64 bits
1169     if (q1 <= 19) {
1170       __mul_128x64_to_128 (C4, C1.w[0], C2);
1171     } else { // q2 <= 19
1172       __mul_128x64_to_128 (C4, C2.w[0], C1);
1173     }
1174     // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
1175     if (C4.w[1] < ten2k128[q1 + q2 - 21].w[1] ||
1176         (C4.w[1] == ten2k128[q1 + q2 - 21].w[1] &&
1177          C4.w[0] < ten2k128[q1 + q2 - 21].w[0])) {
1178       // if (C4.w[1] == 0) // q4 = 20, necessarily
1179       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
1180       // else
1181       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
1182       q4 = q1 + q2 - 1; // q4 in [20, 37]
1183     } else {
1184       // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
1185       q4 = q1 + q2; // q4 in [21, 38]
1186     }
1187   } else if (q1 + q2 == 39) { // C4 = C1 * C2 fits in 128 or 192 bits
1188     // both C1 and C2 fit in 128 bits (actually in 113 bits)
1189     // may replace this by 128x128_to192
1190     __mul_128x128_to_256 (C4, C1, C2); // C4.w[3] is 0
1191     // if C4 < 10^(q1+q2-1) = 10^38 then q4 = q1+q2-1 = 38 else q4 = q1+q2 = 39
1192     if (C4.w[2] == 0 && (C4.w[1] < ten2k128[18].w[1] ||
1193         		 (C4.w[1] == ten2k128[18].w[1]
1194         		  && C4.w[0] < ten2k128[18].w[0]))) {
1195       // 18 = 38 - 20 = q1+q2-1 - 20
1196       // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
1197       q4 = 38; // 38 = q1 + q2 - 1
1198     } else {
1199       // if (C4.w[2] == 0)
1200       // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
1201       // else
1202       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
1203       q4 = 39; // 39 = q1 + q2
1204     }
1205   } else if (q1 + q2 <= 57) { // 40 <= q1 + q2 <= 57
1206     // C4 = C1 * C2 fits in 128 or 192 bits
1207     // (128 bits possibly, but only when q1 + q2 = 40 and C4 has 39 digits)
1208     // both C1 and C2 fit in 128 bits (actually in 113 bits); at most one
1209     // may fit in 64 bits
1210     if (C1.w[1] == 0) { // C1 fits in 64 bits
1211       // __mul_64x128_full (REShi64, RESlo128, A64, B128)
1212       __mul_64x128_full (C4.w[2], C4, C1.w[0], C2);
1213     } else if (C2.w[1] == 0) { // C2 fits in 64 bits
1214       // __mul_64x128_full (REShi64, RESlo128, A64, B128)
1215       __mul_64x128_full (C4.w[2], C4, C2.w[0], C1);
1216     } else { // both C1 and C2 require 128 bits
1217       // may use __mul_128x128_to_192 (C4.w[2], C4.w[0], C2.w[0], C1);
1218       __mul_128x128_to_256 (C4, C1, C2); // C4.w[3] = 0
1219     }
1220     // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
1221     if (C4.w[2] < ten2k256[q1 + q2 - 40].w[2] ||
1222         (C4.w[2] == ten2k256[q1 + q2 - 40].w[2] &&
1223          (C4.w[1] < ten2k256[q1 + q2 - 40].w[1] ||
1224           (C4.w[1] == ten2k256[q1 + q2 - 40].w[1] &&
1225            C4.w[0] < ten2k256[q1 + q2 - 40].w[0])))) {
1226       // if (C4.w[2] == 0) // q4 = 39, necessarily
1227       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
1228       // else
1229       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
1230       q4 = q1 + q2 - 1; // q4 in [39, 56]
1231     } else {
1232       // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
1233       q4 = q1 + q2; // q4 in [40, 57]
1234     }
1235   } else if (q1 + q2 == 58) { // C4 = C1 * C2 fits in 192 or 256 bits
1236     // both C1 and C2 fit in 128 bits (actually in 113 bits); at most one
1237     // may fit in 64 bits
1238     if (C1.w[1] == 0) { // C1 * C2 will fit in 192 bits
1239       __mul_64x128_full (C4.w[2], C4, C1.w[0], C2); // may use 64x128_to_192
1240     } else if (C2.w[1] == 0) { // C1 * C2 will fit in 192 bits
1241       __mul_64x128_full (C4.w[2], C4, C2.w[0], C1); // may use 64x128_to_192
1242     } else { // C1 * C2 will fit in 192 bits or in 256 bits
1243       __mul_128x128_to_256 (C4, C1, C2);
1244     }
1245     // if C4 < 10^(q1+q2-1) = 10^57 then q4 = q1+q2-1 = 57 else q4 = q1+q2 = 58
1246     if (C4.w[3] == 0 && (C4.w[2] < ten2k256[18].w[2] ||
1247         		 (C4.w[2] == ten2k256[18].w[2]
1248         		  && (C4.w[1] < ten2k256[18].w[1]
1249         		      || (C4.w[1] == ten2k256[18].w[1]
1250         			  && C4.w[0] < ten2k256[18].w[0]))))) {
1251       // 18 = 57 - 39 = q1+q2-1 - 39
1252       // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
1253       q4 = 57; // 57 = q1 + q2 - 1
1254     } else {
1255       // if (C4.w[3] == 0)
1256       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
1257       // else
1258       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 256;
1259       q4 = 58; // 58 = q1 + q2
1260     }
1261   } else { // if 59 <= q1 + q2 <= 68
1262     // C4 = C1 * C2 fits in 192 or 256 bits
1263     // (192 bits possibly, but only when q1 + q2 = 59 and C4 has 58 digits)
1264     // both C1 and C2 fit in 128 bits (actually in 113 bits); none fits in
1265     // 64 bits
1266     // may use __mul_128x128_to_192 (C4.w[2], C4.w[0], C2.w[0], C1);
1267     __mul_128x128_to_256 (C4, C1, C2); // C4.w[3] = 0
1268     // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
1269     if (C4.w[3] < ten2k256[q1 + q2 - 40].w[3] ||
1270         (C4.w[3] == ten2k256[q1 + q2 - 40].w[3] &&
1271          (C4.w[2] < ten2k256[q1 + q2 - 40].w[2] ||
1272           (C4.w[2] == ten2k256[q1 + q2 - 40].w[2] &&
1273            (C4.w[1] < ten2k256[q1 + q2 - 40].w[1] ||
1274             (C4.w[1] == ten2k256[q1 + q2 - 40].w[1] &&
1275              C4.w[0] < ten2k256[q1 + q2 - 40].w[0])))))) {
1276       // if (C4.w[3] == 0) // q4 = 58, necessarily
1277       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
1278       // else
1279       //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 256;
1280       q4 = q1 + q2 - 1; // q4 in [58, 67]
1281     } else {
1282       // length of C1 * C2 rounded up to a multiple of 64 bits is len = 256;
1283       q4 = q1 + q2; // q4 in [59, 68]
1284     }
1285   }
1286 
1287   if (C3.w[1] == 0x0 && C3.w[0] == 0x0) { // x = f, y = f, z = 0
1288     save_fpsf = *pfpsf; // sticky bits - caller value must be preserved
1289     *pfpsf = 0;
1290 
1291     if (q4 > p34) {
1292 
1293       // truncate C4 to p34 digits into res
1294       // x = q4-p34, 1 <= x <= 34 because 35 <= q4 <= 68
1295       x0 = q4 - p34;
1296       if (q4 <= 38) {
1297         P128.w[1] = C4.w[1];
1298         P128.w[0] = C4.w[0];
1299         round128_19_38 (q4, x0, P128, &res, &incr_exp,
1300         		&is_midpoint_lt_even, &is_midpoint_gt_even,
1301         		&is_inexact_lt_midpoint,
1302         		&is_inexact_gt_midpoint);
1303       } else if (q4 <= 57) { // 35 <= q4 <= 57
1304         P192.w[2] = C4.w[2];
1305         P192.w[1] = C4.w[1];
1306         P192.w[0] = C4.w[0];
1307         round192_39_57 (q4, x0, P192, &R192, &incr_exp,
1308         		&is_midpoint_lt_even, &is_midpoint_gt_even,
1309         		&is_inexact_lt_midpoint,
1310         		&is_inexact_gt_midpoint);
1311         res.w[0] = R192.w[0];
1312         res.w[1] = R192.w[1];
1313       } else { // if (q4 <= 68)
1314         round256_58_76 (q4, x0, C4, &R256, &incr_exp,
1315         		&is_midpoint_lt_even, &is_midpoint_gt_even,
1316         		&is_inexact_lt_midpoint,
1317         		&is_inexact_gt_midpoint);
1318         res.w[0] = R256.w[0];
1319         res.w[1] = R256.w[1];
1320       }
1321       e4 = e4 + x0;
1322       if (incr_exp) {
1323         e4 = e4 + 1;
1324       }
1325       q4 = p34;
1326       // res is now the coefficient of the result rounded to the destination
1327       // precision, with unbounded exponent; the exponent is e4; q4=digits(res)
1328     } else { // if (q4 <= p34)
1329       // C4 * 10^e4 is the result rounded to the destination precision, with
1330       // unbounded exponent (which is exact)
1331 
1332       if ((q4 + e4 <= p34 + expmax) && (e4 > expmax)) {
1333         // e4 is too large, but can be brought within range by scaling up C4
1334         scale = e4 - expmax; // 1 <= scale < P-q4 <= P-1 => 1 <= scale <= P-2
1335         // res = (C4 * 10^scale) * 10^expmax
1336         if (q4 <= 19) { // C4 fits in 64 bits
1337           if (scale <= 19) { // 10^scale fits in 64 bits
1338             // 64 x 64 C4.w[0] * ten2k64[scale]
1339             __mul_64x64_to_128MACH (res, C4.w[0], ten2k64[scale]);
1340           } else { // 10^scale fits in 128 bits
1341             // 64 x 128 C4.w[0] * ten2k128[scale - 20]
1342             __mul_128x64_to_128 (res, C4.w[0], ten2k128[scale - 20]);
1343           }
1344         } else { // C4 fits in 128 bits, but 10^scale must fit in 64 bits
1345           // 64 x 128 ten2k64[scale] * CC43
1346           __mul_128x64_to_128 (res, ten2k64[scale], C4);
1347         }
1348         e4 = e4 - scale; // expmax
1349         q4 = q4 + scale;
1350       } else {
1351         res.w[1] = C4.w[1];
1352         res.w[0] = C4.w[0];
1353       }
1354       // res is the coefficient of the result rounded to the destination
1355       // precision, with unbounded exponent (it has q4 digits); the exponent
1356       // is e4 (exact result)
1357     }
1358 
1359     // check for overflow
1360     if (q4 + e4 > p34 + expmax) {
1361       if (rnd_mode == ROUNDING_TO_NEAREST) {
1362         res.w[1] = p_sign | 0x7800000000000000ull; // +/-inf
1363         res.w[0] = 0x0000000000000000ull;
1364         *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
1365       } else {
1366         res.w[1] = p_sign | res.w[1];
1367         rounding_correction (rnd_mode,
1368         		     is_inexact_lt_midpoint,
1369         		     is_inexact_gt_midpoint,
1370         		     is_midpoint_lt_even, is_midpoint_gt_even,
1371         		     e4, &res, pfpsf);
1372       }
1373       *pfpsf |= save_fpsf;
1374       *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
1375       *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
1376       *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
1377       *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
1378       BID_SWAP128 (res);
1379       BID_RETURN (res)
1380     }
1381     // check for underflow
1382     if (q4 + e4 < expmin + P34) {
1383       is_tiny = 1; // the result is tiny
1384       if (e4 < expmin) {
1385         // if e4 < expmin, we must truncate more of res
1386         x0 = expmin - e4; // x0 >= 1
1387         is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
1388         is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
1389         is_midpoint_lt_even0 = is_midpoint_lt_even;
1390         is_midpoint_gt_even0 = is_midpoint_gt_even;
1391         is_inexact_lt_midpoint = 0;
1392         is_inexact_gt_midpoint = 0;
1393         is_midpoint_lt_even = 0;
1394         is_midpoint_gt_even = 0;
1395         // the number of decimal digits in res is q4
1396         if (x0 < q4) { // 1 <= x0 <= q4-1 => round res to q4 - x0 digits
1397           if (q4 <= 18) { // 2 <= q4 <= 18, 1 <= x0 <= 17
1398             round64_2_18 (q4, x0, res.w[0], &R64, &incr_exp,
1399         		  &is_midpoint_lt_even, &is_midpoint_gt_even,
1400         		  &is_inexact_lt_midpoint,
1401         		  &is_inexact_gt_midpoint);
1402             if (incr_exp) {
1403               // R64 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 17
1404               R64 = ten2k64[q4 - x0];
1405             }
1406             // res.w[1] = 0; (from above)
1407             res.w[0] = R64;
1408           } else { // if (q4 <= 34)
1409             // 19 <= q4 <= 38
1410             P128.w[1] = res.w[1];
1411             P128.w[0] = res.w[0];
1412             round128_19_38 (q4, x0, P128, &res, &incr_exp,
1413         		    &is_midpoint_lt_even, &is_midpoint_gt_even,
1414         		    &is_inexact_lt_midpoint,
1415         		    &is_inexact_gt_midpoint);
1416             if (incr_exp) {
1417               // increase coefficient by a factor of 10; this will be <= 10^33
1418               // R128 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 37
1419               if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
1420         	// res.w[1] = 0;
1421         	res.w[0] = ten2k64[q4 - x0];
1422               } else { // 20 <= q4 - x0 <= 37
1423         	res.w[0] = ten2k128[q4 - x0 - 20].w[0];
1424         	res.w[1] = ten2k128[q4 - x0 - 20].w[1];
1425               }
1426             }
1427           }
1428           e4 = e4 + x0; // expmin
1429         } else if (x0 == q4) {
1430           // the second rounding is for 0.d(0)d(1)...d(q4-1) * 10^emin
1431           // determine relationship with 1/2 ulp
1432           if (q4 <= 19) {
1433             if (res.w[0] < midpoint64[q4 - 1]) { // < 1/2 ulp
1434               lt_half_ulp = 1;
1435               is_inexact_lt_midpoint = 1;
1436             } else if (res.w[0] == midpoint64[q4 - 1]) { // = 1/2 ulp
1437               eq_half_ulp = 1;
1438               is_midpoint_gt_even = 1;
1439             } else { // > 1/2 ulp
1440               // gt_half_ulp = 1;
1441               is_inexact_gt_midpoint = 1;
1442             }
1443           } else { // if (q4 <= 34)
1444             if (res.w[1] < midpoint128[q4 - 20].w[1] ||
1445                 (res.w[1] == midpoint128[q4 - 20].w[1] &&
1446                 res.w[0] < midpoint128[q4 - 20].w[0])) { // < 1/2 ulp
1447               lt_half_ulp = 1;
1448               is_inexact_lt_midpoint = 1;
1449             } else if (res.w[1] == midpoint128[q4 - 20].w[1] &&
1450                 res.w[0] == midpoint128[q4 - 20].w[0]) { // = 1/2 ulp
1451               eq_half_ulp = 1;
1452               is_midpoint_gt_even = 1;
1453             } else { // > 1/2 ulp
1454               // gt_half_ulp = 1;
1455               is_inexact_gt_midpoint = 1;
1456             }
1457           }
1458           if (lt_half_ulp || eq_half_ulp) {
1459             // res = +0.0 * 10^expmin
1460             res.w[1] = 0x0000000000000000ull;
1461             res.w[0] = 0x0000000000000000ull;
1462           } else { // if (gt_half_ulp)
1463             // res = +1 * 10^expmin
1464             res.w[1] = 0x0000000000000000ull;
1465             res.w[0] = 0x0000000000000001ull;
1466           }
1467           e4 = expmin;
1468         } else { // if (x0 > q4)
1469           // the second rounding is for 0.0...d(0)d(1)...d(q4-1) * 10^emin
1470           res.w[1] = 0;
1471           res.w[0] = 0;
1472           e4 = expmin;
1473           is_inexact_lt_midpoint = 1;
1474         }
1475         // avoid a double rounding error
1476         if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
1477             is_midpoint_lt_even) { // double rounding error upward
1478           // res = res - 1
1479           res.w[0]--;
1480           if (res.w[0] == 0xffffffffffffffffull)
1481             res.w[1]--;
1482           // Note: a double rounding error upward is not possible; for this
1483           // the result after the first rounding would have to be 99...95
1484           // (35 digits in all), possibly followed by a number of zeros; this
1485           // not possible for f * f + 0
1486           is_midpoint_lt_even = 0;
1487           is_inexact_lt_midpoint = 1;
1488         } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
1489             is_midpoint_gt_even) { // double rounding error downward
1490           // res = res + 1
1491           res.w[0]++;
1492           if (res.w[0] == 0)
1493             res.w[1]++;
1494           is_midpoint_gt_even = 0;
1495           is_inexact_gt_midpoint = 1;
1496         } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
1497         	   !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
1498           // if this second rounding was exact the result may still be
1499           // inexact because of the first rounding
1500           if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
1501             is_inexact_gt_midpoint = 1;
1502           }
1503           if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
1504             is_inexact_lt_midpoint = 1;
1505           }
1506         } else if (is_midpoint_gt_even &&
1507         	   (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
1508           // pulled up to a midpoint
1509           is_inexact_lt_midpoint = 1;
1510           is_inexact_gt_midpoint = 0;
1511           is_midpoint_lt_even = 0;
1512           is_midpoint_gt_even = 0;
1513         } else if (is_midpoint_lt_even &&
1514         	   (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
1515           // pulled down to a midpoint
1516           is_inexact_lt_midpoint = 0;
1517           is_inexact_gt_midpoint = 1;
1518           is_midpoint_lt_even = 0;
1519           is_midpoint_gt_even = 0;
1520         } else {
1521           ;
1522         }
1523       } else { // if e4 >= emin then q4 < P and the result is tiny and exact
1524         if (e3 < e4) {
1525           // if (e3 < e4) the preferred exponent is e3
1526           // return (C4 * 10^scale) * 10^(e4 - scale)
1527           // where scale = min (p34-q4, (e4 - e3))
1528           scale = p34 - q4;
1529           ind = e4 - e3;
1530           if (ind < scale)
1531             scale = ind;
1532           if (scale == 0) {
1533             ; // res and e4 are unchanged
1534           } else if (q4 <= 19) { // C4 fits in 64 bits
1535             if (scale <= 19) { // 10^scale fits in 64 bits
1536               // 64 x 64 res.w[0] * ten2k64[scale]
1537               __mul_64x64_to_128MACH (res, res.w[0], ten2k64[scale]);
1538             } else { // 10^scale fits in 128 bits
1539               // 64 x 128 res.w[0] * ten2k128[scale - 20]
1540               __mul_128x64_to_128 (res, res.w[0], ten2k128[scale - 20]);
1541             }
1542           } else { // res fits in 128 bits, but 10^scale must fit in 64 bits
1543             // 64 x 128 ten2k64[scale] * C3
1544             __mul_128x64_to_128 (res, ten2k64[scale], res);
1545           }
1546           // subtract scale from the exponent
1547           e4 = e4 - scale;
1548         }
1549       }
1550 
1551       // check for inexact result
1552       if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
1553           is_midpoint_lt_even || is_midpoint_gt_even) {
1554         // set the inexact flag and the underflow flag
1555         *pfpsf |= INEXACT_EXCEPTION;
1556         *pfpsf |= UNDERFLOW_EXCEPTION;
1557       }
1558       res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1];
1559       if (rnd_mode != ROUNDING_TO_NEAREST) {
1560         rounding_correction (rnd_mode,
1561         		     is_inexact_lt_midpoint,
1562         		     is_inexact_gt_midpoint,
1563         		     is_midpoint_lt_even, is_midpoint_gt_even,
1564         		     e4, &res, pfpsf);
1565       }
1566       *pfpsf |= save_fpsf;
1567       *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
1568       *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
1569       *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
1570       *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
1571       BID_SWAP128 (res);
1572       BID_RETURN (res)
1573     }
1574     // no overflow, and no underflow for rounding to nearest
1575     res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1];
1576 
1577     if (rnd_mode != ROUNDING_TO_NEAREST) {
1578       rounding_correction (rnd_mode,
1579         		   is_inexact_lt_midpoint,
1580         		   is_inexact_gt_midpoint,
1581         		   is_midpoint_lt_even, is_midpoint_gt_even,
1582         		   e4, &res, pfpsf);
1583       // if e4 = expmin && significand < 10^33 => result is tiny (for RD, RZ)
1584       if (e4 == expmin) {
1585         if ((res.w[1] & MASK_COEFF) < 0x0000314dc6448d93ull ||
1586             ((res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull &&
1587              res.w[0] < 0x38c15b0a00000000ull)) {
1588           is_tiny = 1;
1589         }
1590       }
1591     }
1592 
1593     if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
1594         is_midpoint_lt_even || is_midpoint_gt_even) {
1595       // set the inexact flag
1596       *pfpsf |= INEXACT_EXCEPTION;
1597       if (is_tiny)
1598         *pfpsf |= UNDERFLOW_EXCEPTION;
1599     }
1600 
1601     if ((*pfpsf & INEXACT_EXCEPTION) == 0) { // x * y is exact
1602       // need to ensure that the result has the preferred exponent
1603       p_exp = res.w[1] & MASK_EXP;
1604       if (z_exp < p_exp) { // the preferred exponent is z_exp
1605         // signficand of res in C3
1606         C3.w[1] = res.w[1] & MASK_COEFF;
1607         C3.w[0] = res.w[0];
1608         // the number of decimal digits of x * y is q4 <= 34
1609         // Note: the coefficient fits in 128 bits
1610 
1611         // return (C3 * 10^scale) * 10^(p_exp - scale)
1612         // where scale = min (p34-q4, (p_exp-z_exp) >> 49)
1613         scale = p34 - q4;
1614         ind = (p_exp - z_exp) >> 49;
1615         if (ind < scale)
1616           scale = ind;
1617         // subtract scale from the exponent
1618         p_exp = p_exp - ((UINT64) scale << 49);
1619         if (scale == 0) {
1620           ; // leave res unchanged
1621         } else if (q4 <= 19) { // x * y fits in 64 bits
1622           if (scale <= 19) { // 10^scale fits in 64 bits
1623             // 64 x 64 C3.w[0] * ten2k64[scale]
1624             __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
1625           } else { // 10^scale fits in 128 bits
1626             // 64 x 128 C3.w[0] * ten2k128[scale - 20]
1627             __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
1628           }
1629           res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1];
1630         } else { // x * y fits in 128 bits, but 10^scale must fit in 64 bits
1631           // 64 x 128 ten2k64[scale] * C3
1632           __mul_128x64_to_128 (res, ten2k64[scale], C3);
1633           res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1];
1634         }
1635       } // else leave the result as it is, because p_exp <= z_exp
1636     }
1637     *pfpsf |= save_fpsf;
1638     *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
1639     *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
1640     *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
1641     *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
1642     BID_SWAP128 (res);
1643     BID_RETURN (res)
1644   } // else we have f * f + f
1645 
1646   // continue with x = f, y = f, z = f
1647 
1648   delta = q3 + e3 - q4 - e4;
1649 delta_ge_zero:
1650   if (delta >= 0) {
1651 
1652     if (p34 <= delta - 1 ||	// Case (1')
1653         (p34 == delta && e3 + 6176 < p34 - q3)) { // Case (1''A)
1654       // check for overflow, which can occur only in Case (1')
1655       if ((q3 + e3) > (p34 + expmax) && p34 <= delta - 1) {
1656         // e3 > expmax implies p34 <= delta-1 and e3 > expmax is a necessary
1657         // condition for (q3 + e3) > (p34 + expmax)
1658         if (rnd_mode == ROUNDING_TO_NEAREST) {
1659           res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
1660           res.w[0] = 0x0000000000000000ull;
1661           *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
1662         } else {
1663           if (p_sign == z_sign) {
1664             is_inexact_lt_midpoint = 1;
1665           } else {
1666             is_inexact_gt_midpoint = 1;
1667           }
1668           // q3 <= p34; if (q3 < p34) scale C3 up by 10^(p34-q3)
1669           scale = p34 - q3;
1670           if (scale == 0) {
1671             res.w[1] = z_sign | C3.w[1];
1672             res.w[0] = C3.w[0];
1673           } else {
1674             if (q3 <= 19) { // C3 fits in 64 bits
1675               if (scale <= 19) { // 10^scale fits in 64 bits
1676         	// 64 x 64 C3.w[0] * ten2k64[scale]
1677         	__mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
1678               } else { // 10^scale fits in 128 bits
1679         	// 64 x 128 C3.w[0] * ten2k128[scale - 20]
1680         	__mul_128x64_to_128 (res, C3.w[0],
1681         			     ten2k128[scale - 20]);
1682               }
1683             } else { // C3 fits in 128 bits, but 10^scale must fit in 64 bits
1684               // 64 x 128 ten2k64[scale] * C3
1685               __mul_128x64_to_128 (res, ten2k64[scale], C3);
1686             }
1687             // the coefficient in res has q3 + scale = p34 digits
1688           }
1689           e3 = e3 - scale;
1690           res.w[1] = z_sign | res.w[1];
1691           rounding_correction (rnd_mode,
1692         		       is_inexact_lt_midpoint,
1693         		       is_inexact_gt_midpoint,
1694         		       is_midpoint_lt_even, is_midpoint_gt_even,
1695         		       e3, &res, pfpsf);
1696         }
1697         *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
1698         *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
1699         *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
1700         *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
1701         BID_SWAP128 (res);
1702         BID_RETURN (res)
1703       }
1704       // res = z
1705       if (q3 < p34) { // the preferred exponent is z_exp - (p34 - q3)
1706         // return (C3 * 10^scale) * 10^(z_exp - scale)
1707         // where scale = min (p34-q3, z_exp-EMIN)
1708         scale = p34 - q3;
1709         ind = e3 + 6176;
1710         if (ind < scale)
1711           scale = ind;
1712         if (scale == 0) {
1713           res.w[1] = C3.w[1];
1714           res.w[0] = C3.w[0];
1715         } else if (q3 <= 19) { // z fits in 64 bits
1716           if (scale <= 19) { // 10^scale fits in 64 bits
1717             // 64 x 64 C3.w[0] * ten2k64[scale]
1718             __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
1719           } else { // 10^scale fits in 128 bits
1720             // 64 x 128 C3.w[0] * ten2k128[scale - 20]
1721             __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
1722           }
1723         } else { // z fits in 128 bits, but 10^scale must fit in 64 bits
1724           // 64 x 128 ten2k64[scale] * C3
1725           __mul_128x64_to_128 (res, ten2k64[scale], C3);
1726         }
1727         // the coefficient in res has q3 + scale digits
1728         // subtract scale from the exponent
1729         z_exp = z_exp - ((UINT64) scale << 49);
1730         e3 = e3 - scale;
1731         res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
1732         if (scale + q3 < p34)
1733           *pfpsf |= UNDERFLOW_EXCEPTION;
1734       } else {
1735         scale = 0;
1736         res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | C3.w[1];
1737         res.w[0] = C3.w[0];
1738       }
1739 
1740       // use the following to avoid double rounding errors when operating on
1741       // mixed formats in rounding to nearest, and for correcting the result
1742       // if not rounding to nearest
1743       if ((p_sign != z_sign) && (delta == (q3 + scale + 1))) {
1744         // there is a gap of exactly one digit between the scaled C3 and C4
1745         // C3 * 10^ scale = 10^(q3+scale-1) <=> C3 = 10^(q3-1) is special case
1746         if ((q3 <= 19 && C3.w[0] != ten2k64[q3 - 1]) ||
1747             (q3 == 20 && (C3.w[1] != 0 || C3.w[0] != ten2k64[19])) ||
1748             (q3 >= 21 && (C3.w[1] != ten2k128[q3 - 21].w[1] ||
1749         		  C3.w[0] != ten2k128[q3 - 21].w[0]))) {
1750           // C3 * 10^ scale != 10^(q3-1)
1751           // if ((res.w[1] & MASK_COEFF) != 0x0000314dc6448d93ull ||
1752           // res.w[0] != 0x38c15b0a00000000ull) { // C3 * 10^scale != 10^33
1753           is_inexact_gt_midpoint = 1; // if (z_sign), set as if for abs. value
1754         } else { // if C3 * 10^scale = 10^(q3+scale-1)
1755           // ok from above e3 = (z_exp >> 49) - 6176;
1756           // the result is always inexact
1757           if (q4 == 1) {
1758             R64 = C4.w[0];
1759           } else {
1760             // if q4 > 1 then truncate C4 from q4 digits to 1 digit;
1761             // x = q4-1, 1 <= x <= 67 and check if this operation is exact
1762             if (q4 <= 18) { // 2 <= q4 <= 18
1763               round64_2_18 (q4, q4 - 1, C4.w[0], &R64, &incr_exp,
1764         		    &is_midpoint_lt_even, &is_midpoint_gt_even,
1765         		    &is_inexact_lt_midpoint,
1766         		    &is_inexact_gt_midpoint);
1767             } else if (q4 <= 38) {
1768               P128.w[1] = C4.w[1];
1769               P128.w[0] = C4.w[0];
1770               round128_19_38 (q4, q4 - 1, P128, &R128, &incr_exp,
1771         		      &is_midpoint_lt_even,
1772         		      &is_midpoint_gt_even,
1773         		      &is_inexact_lt_midpoint,
1774         		      &is_inexact_gt_midpoint);
1775               R64 = R128.w[0]; // one decimal digit
1776             } else if (q4 <= 57) {
1777               P192.w[2] = C4.w[2];
1778               P192.w[1] = C4.w[1];
1779               P192.w[0] = C4.w[0];
1780               round192_39_57 (q4, q4 - 1, P192, &R192, &incr_exp,
1781         		      &is_midpoint_lt_even,
1782         		      &is_midpoint_gt_even,
1783         		      &is_inexact_lt_midpoint,
1784         		      &is_inexact_gt_midpoint);
1785               R64 = R192.w[0]; // one decimal digit
1786             } else { // if (q4 <= 68)
1787               round256_58_76 (q4, q4 - 1, C4, &R256, &incr_exp,
1788         		      &is_midpoint_lt_even,
1789         		      &is_midpoint_gt_even,
1790         		      &is_inexact_lt_midpoint,
1791         		      &is_inexact_gt_midpoint);
1792               R64 = R256.w[0]; // one decimal digit
1793             }
1794             if (incr_exp) {
1795               R64 = 10;
1796             }
1797           }
1798           if (q4 == 1 && C4.w[0] == 5) {
1799             is_inexact_lt_midpoint = 0;
1800             is_inexact_gt_midpoint = 0;
1801             is_midpoint_lt_even = 1;
1802             is_midpoint_gt_even = 0;
1803           } else if ((e3 == expmin) ||
1804         	     R64 < 5 || (R64 == 5 && is_inexact_gt_midpoint)) {
1805             // result does not change
1806             is_inexact_lt_midpoint = 0;
1807             is_inexact_gt_midpoint = 1;
1808             is_midpoint_lt_even = 0;
1809             is_midpoint_gt_even = 0;
1810           } else {
1811             is_inexact_lt_midpoint = 1;
1812             is_inexact_gt_midpoint = 0;
1813             is_midpoint_lt_even = 0;
1814             is_midpoint_gt_even = 0;
1815             // result decremented is 10^(q3+scale) - 1
1816             if ((q3 + scale) <= 19) {
1817               res.w[1] = 0;
1818               res.w[0] = ten2k64[q3 + scale];
1819             } else { // if ((q3 + scale + 1) <= 35)
1820               res.w[1] = ten2k128[q3 + scale - 20].w[1];
1821               res.w[0] = ten2k128[q3 + scale - 20].w[0];
1822             }
1823             res.w[0] = res.w[0] - 1; // borrow never occurs
1824             z_exp = z_exp - EXP_P1;
1825             e3 = e3 - 1;
1826             res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1];
1827           }
1828           if (e3 == expmin) {
1829             if (R64 < 5 || (R64 == 5 && !is_inexact_lt_midpoint)) {
1830               ; // result not tiny (in round-to-nearest mode)
1831             } else {
1832               *pfpsf |= UNDERFLOW_EXCEPTION;
1833             }
1834           }
1835         } // end 10^(q3+scale-1)
1836         // set the inexact flag
1837         *pfpsf |= INEXACT_EXCEPTION;
1838       } else {
1839         if (p_sign == z_sign) {
1840           // if (z_sign), set as if for absolute value
1841           is_inexact_lt_midpoint = 1;
1842         } else { // if (p_sign != z_sign)
1843           // if (z_sign), set as if for absolute value
1844           is_inexact_gt_midpoint = 1;
1845         }
1846         *pfpsf |= INEXACT_EXCEPTION;
1847       }
1848       // the result is always inexact => set the inexact flag
1849       // Determine tininess:
1850       //    if (exp > expmin)
1851       //      the result is not tiny
1852       //    else // if exp = emin
1853       //      if (q3 + scale < p34)
1854       //        the result is tiny
1855       //      else // if (q3 + scale = p34)
1856       //        if (C3 * 10^scale > 10^33)
1857       //          the result is not tiny
1858       //        else // if C3 * 10^scale = 10^33
1859       //          if (xy * z > 0)
1860       //            the result is not tiny
1861       //          else // if (xy * z < 0)
1862       //            if (z > 0)
1863       //              if rnd_mode != RP
1864       //                the result is tiny
1865       //              else // if RP
1866       //                the result is not tiny
1867       //            else // if (z < 0)
1868       //              if rnd_mode != RM
1869       //                the result is tiny
1870       //              else // if RM
1871       //                the result is not tiny
1872       //              endif
1873       //            endif
1874       //          endif
1875       //        endif
1876       //      endif
1877       //    endif
1878       if ((e3 == expmin && (q3 + scale) < p34) ||
1879           (e3 == expmin && (q3 + scale) == p34 &&
1880           (res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull &&	// 10^33_high
1881           res.w[0] == 0x38c15b0a00000000ull &&	// 10^33_low
1882           z_sign != p_sign && ((!z_sign && rnd_mode != ROUNDING_UP) ||
1883           (z_sign && rnd_mode != ROUNDING_DOWN)))) {
1884         *pfpsf |= UNDERFLOW_EXCEPTION;
1885       }
1886       if (rnd_mode != ROUNDING_TO_NEAREST) {
1887         rounding_correction (rnd_mode,
1888         		     is_inexact_lt_midpoint,
1889         		     is_inexact_gt_midpoint,
1890         		     is_midpoint_lt_even, is_midpoint_gt_even,
1891         		     e3, &res, pfpsf);
1892       }
1893       *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
1894       *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
1895       *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
1896       *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
1897       BID_SWAP128 (res);
1898       BID_RETURN (res)
1899 
1900     } else if (p34 == delta) { // Case (1''B)
1901 
1902       // because Case (1''A) was treated above, e3 + 6176 >= p34 - q3
1903       // and C3 can be scaled up to p34 digits if needed
1904 
1905       // scale C3 to p34 digits if needed
1906       scale = p34 - q3; // 0 <= scale <= p34 - 1
1907       if (scale == 0) {
1908         res.w[1] = C3.w[1];
1909         res.w[0] = C3.w[0];
1910       } else if (q3 <= 19) { // z fits in 64 bits
1911         if (scale <= 19) { // 10^scale fits in 64 bits
1912           // 64 x 64 C3.w[0] * ten2k64[scale]
1913           __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
1914         } else { // 10^scale fits in 128 bits
1915           // 64 x 128 C3.w[0] * ten2k128[scale - 20]
1916           __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
1917         }
1918       } else { // z fits in 128 bits, but 10^scale must fit in 64 bits
1919         // 64 x 128 ten2k64[scale] * C3
1920         __mul_128x64_to_128 (res, ten2k64[scale], C3);
1921       }
1922       // subtract scale from the exponent
1923       z_exp = z_exp - ((UINT64) scale << 49);
1924       e3 = e3 - scale;
1925       // now z_sign, z_exp, and res correspond to a z scaled to p34 = 34 digits
1926 
1927       // determine whether x * y is less than, equal to, or greater than
1928       // 1/2 ulp (z)
1929       if (q4 <= 19) {
1930         if (C4.w[0] < midpoint64[q4 - 1]) { // < 1/2 ulp
1931           lt_half_ulp = 1;
1932         } else if (C4.w[0] == midpoint64[q4 - 1]) { // = 1/2 ulp
1933           eq_half_ulp = 1;
1934         } else { // > 1/2 ulp
1935           gt_half_ulp = 1;
1936         }
1937       } else if (q4 <= 38) {
1938         if (C4.w[2] == 0 && (C4.w[1] < midpoint128[q4 - 20].w[1] ||
1939             (C4.w[1] == midpoint128[q4 - 20].w[1] &&
1940             C4.w[0] < midpoint128[q4 - 20].w[0]))) { // < 1/2 ulp
1941           lt_half_ulp = 1;
1942         } else if (C4.w[2] == 0 && C4.w[1] == midpoint128[q4 - 20].w[1] &&
1943             C4.w[0] == midpoint128[q4 - 20].w[0]) { // = 1/2 ulp
1944           eq_half_ulp = 1;
1945         } else { // > 1/2 ulp
1946           gt_half_ulp = 1;
1947         }
1948       } else if (q4 <= 58) {
1949         if (C4.w[3] == 0 && (C4.w[2] < midpoint192[q4 - 39].w[2] ||
1950             (C4.w[2] == midpoint192[q4 - 39].w[2] &&
1951             C4.w[1] < midpoint192[q4 - 39].w[1]) ||
1952             (C4.w[2] == midpoint192[q4 - 39].w[2] &&
1953             C4.w[1] == midpoint192[q4 - 39].w[1] &&
1954             C4.w[0] < midpoint192[q4 - 39].w[0]))) { // < 1/2 ulp
1955           lt_half_ulp = 1;
1956         } else if (C4.w[3] == 0 && C4.w[2] == midpoint192[q4 - 39].w[2] &&
1957             C4.w[1] == midpoint192[q4 - 39].w[1] &&
1958             C4.w[0] == midpoint192[q4 - 39].w[0]) { // = 1/2 ulp
1959           eq_half_ulp = 1;
1960         } else { // > 1/2 ulp
1961           gt_half_ulp = 1;
1962         }
1963       } else {
1964         if (C4.w[3] < midpoint256[q4 - 59].w[3] ||
1965             (C4.w[3] == midpoint256[q4 - 59].w[3] &&
1966             C4.w[2] < midpoint256[q4 - 59].w[2]) ||
1967             (C4.w[3] == midpoint256[q4 - 59].w[3] &&
1968             C4.w[2] == midpoint256[q4 - 59].w[2] &&
1969             C4.w[1] < midpoint256[q4 - 59].w[1]) ||
1970             (C4.w[3] == midpoint256[q4 - 59].w[3] &&
1971             C4.w[2] == midpoint256[q4 - 59].w[2] &&
1972             C4.w[1] == midpoint256[q4 - 59].w[1] &&
1973             C4.w[0] < midpoint256[q4 - 59].w[0])) { // < 1/2 ulp
1974           lt_half_ulp = 1;
1975         } else if (C4.w[3] == midpoint256[q4 - 59].w[3] &&
1976             C4.w[2] == midpoint256[q4 - 59].w[2] &&
1977             C4.w[1] == midpoint256[q4 - 59].w[1] &&
1978             C4.w[0] == midpoint256[q4 - 59].w[0]) { // = 1/2 ulp
1979           eq_half_ulp = 1;
1980         } else { // > 1/2 ulp
1981           gt_half_ulp = 1;
1982         }
1983       }
1984 
1985       if (p_sign == z_sign) {
1986         if (lt_half_ulp) {
1987           res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
1988           // use the following to avoid double rounding errors when operating on
1989           // mixed formats in rounding to nearest
1990           is_inexact_lt_midpoint = 1; // if (z_sign), as if for absolute value
1991         } else if ((eq_half_ulp && (res.w[0] & 0x01)) || gt_half_ulp) {
1992           // add 1 ulp to the significand
1993           res.w[0]++;
1994           if (res.w[0] == 0x0ull)
1995             res.w[1]++;
1996           // check for rounding overflow, when coeff == 10^34
1997           if ((res.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull &&
1998               res.w[0] == 0x378d8e6400000000ull) { // coefficient = 10^34
1999             e3 = e3 + 1;
2000             // coeff = 10^33
2001             z_exp = ((UINT64) (e3 + 6176) << 49) & MASK_EXP;
2002             res.w[1] = 0x0000314dc6448d93ull;
2003             res.w[0] = 0x38c15b0a00000000ull;
2004           }
2005           // end add 1 ulp
2006           res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
2007           if (eq_half_ulp) {
2008             is_midpoint_lt_even = 1; // if (z_sign), as if for absolute value
2009           } else {
2010             is_inexact_gt_midpoint = 1; // if (z_sign), as if for absolute value
2011           }
2012         } else { // if (eq_half_ulp && !(res.w[0] & 0x01))
2013           // leave unchanged
2014           res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
2015           is_midpoint_gt_even = 1; // if (z_sign), as if for absolute value
2016         }
2017         // the result is always inexact, and never tiny
2018         // set the inexact flag
2019         *pfpsf |= INEXACT_EXCEPTION;
2020         // check for overflow
2021         if (e3 > expmax && rnd_mode == ROUNDING_TO_NEAREST) {
2022           res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
2023           res.w[0] = 0x0000000000000000ull;
2024           *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
2025           *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
2026           *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
2027           *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
2028           *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
2029           BID_SWAP128 (res);
2030           BID_RETURN (res)
2031         }
2032         if (rnd_mode != ROUNDING_TO_NEAREST) {
2033           rounding_correction (rnd_mode,
2034         		       is_inexact_lt_midpoint,
2035         		       is_inexact_gt_midpoint,
2036         		       is_midpoint_lt_even, is_midpoint_gt_even,
2037         		       e3, &res, pfpsf);
2038           z_exp = res.w[1] & MASK_EXP;
2039         }
2040       } else { // if (p_sign != z_sign)
2041         // consider two cases, because C3 * 10^scale = 10^33 is a special case
2042         if (res.w[1] != 0x0000314dc6448d93ull ||
2043             res.w[0] != 0x38c15b0a00000000ull) { // C3 * 10^scale != 10^33
2044           if (lt_half_ulp) {
2045             res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
2046             // use the following to avoid double rounding errors when operating
2047             // on mixed formats in rounding to nearest
2048             is_inexact_gt_midpoint = 1; // if (z_sign), as if for absolute value
2049           } else if ((eq_half_ulp && (res.w[0] & 0x01)) || gt_half_ulp) {
2050             // subtract 1 ulp from the significand
2051             res.w[0]--;
2052             if (res.w[0] == 0xffffffffffffffffull)
2053               res.w[1]--;
2054             res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
2055             if (eq_half_ulp) {
2056               is_midpoint_gt_even = 1; // if (z_sign), as if for absolute value
2057             } else {
2058               is_inexact_lt_midpoint = 1; //if(z_sign), as if for absolute value
2059             }
2060           } else { // if (eq_half_ulp && !(res.w[0] & 0x01))
2061             // leave unchanged
2062             res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
2063             is_midpoint_lt_even = 1; // if (z_sign), as if for absolute value
2064           }
2065           // the result is always inexact, and never tiny
2066           // check for overflow for RN
2067           if (e3 > expmax) {
2068             if (rnd_mode == ROUNDING_TO_NEAREST) {
2069               res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
2070               res.w[0] = 0x0000000000000000ull;
2071               *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
2072             } else {
2073               rounding_correction (rnd_mode,
2074         			   is_inexact_lt_midpoint,
2075         			   is_inexact_gt_midpoint,
2076         			   is_midpoint_lt_even,
2077         			   is_midpoint_gt_even, e3, &res,
2078         			   pfpsf);
2079             }
2080             *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
2081             *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
2082             *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
2083             *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
2084             BID_SWAP128 (res);
2085             BID_RETURN (res)
2086           }
2087           // set the inexact flag
2088           *pfpsf |= INEXACT_EXCEPTION;
2089           if (rnd_mode != ROUNDING_TO_NEAREST) {
2090             rounding_correction (rnd_mode,
2091         			 is_inexact_lt_midpoint,
2092         			 is_inexact_gt_midpoint,
2093         			 is_midpoint_lt_even,
2094         			 is_midpoint_gt_even, e3, &res, pfpsf);
2095           }
2096           z_exp = res.w[1] & MASK_EXP;
2097         } else { // if C3 * 10^scale = 10^33
2098           e3 = (z_exp >> 49) - 6176;
2099           if (e3 > expmin) {
2100             // the result is exact if exp > expmin and C4 = d*10^(q4-1),
2101             // where d = 1, 2, 3, ..., 9; it could be tiny too, but exact
2102             if (q4 == 1) {
2103               // if q4 = 1 the result is exact
2104               // result coefficient = 10^34 - C4
2105               res.w[1] = 0x0001ed09bead87c0ull;
2106               res.w[0] = 0x378d8e6400000000ull - C4.w[0];
2107               z_exp = z_exp - EXP_P1;
2108               e3 = e3 - 1;
2109               res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
2110             } else {
2111               // if q4 > 1 then truncate C4 from q4 digits to 1 digit;
2112               // x = q4-1, 1 <= x <= 67 and check if this operation is exact
2113               if (q4 <= 18) { // 2 <= q4 <= 18
2114         	round64_2_18 (q4, q4 - 1, C4.w[0], &R64, &incr_exp,
2115         		      &is_midpoint_lt_even,
2116         		      &is_midpoint_gt_even,
2117         		      &is_inexact_lt_midpoint,
2118         		      &is_inexact_gt_midpoint);
2119               } else if (q4 <= 38) {
2120         	P128.w[1] = C4.w[1];
2121         	P128.w[0] = C4.w[0];
2122         	round128_19_38 (q4, q4 - 1, P128, &R128, &incr_exp,
2123         			&is_midpoint_lt_even,
2124         			&is_midpoint_gt_even,
2125         			&is_inexact_lt_midpoint,
2126         			&is_inexact_gt_midpoint);
2127         	R64 = R128.w[0]; // one decimal digit
2128               } else if (q4 <= 57) {
2129         	P192.w[2] = C4.w[2];
2130         	P192.w[1] = C4.w[1];
2131         	P192.w[0] = C4.w[0];
2132         	round192_39_57 (q4, q4 - 1, P192, &R192, &incr_exp,
2133         			&is_midpoint_lt_even,
2134         			&is_midpoint_gt_even,
2135         			&is_inexact_lt_midpoint,
2136         			&is_inexact_gt_midpoint);
2137         	R64 = R192.w[0]; // one decimal digit
2138               } else { // if (q4 <= 68)
2139         	round256_58_76 (q4, q4 - 1, C4, &R256, &incr_exp,
2140         			&is_midpoint_lt_even,
2141         			&is_midpoint_gt_even,
2142         			&is_inexact_lt_midpoint,
2143         			&is_inexact_gt_midpoint);
2144         	R64 = R256.w[0]; // one decimal digit
2145               }
2146               if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
2147         	  !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
2148         	// the result is exact: 10^34 - R64
2149         	// incr_exp = 0 with certainty
2150         	z_exp = z_exp - EXP_P1;
2151         	e3 = e3 - 1;
2152         	res.w[1] =
2153         	  z_sign | (z_exp & MASK_EXP) | 0x0001ed09bead87c0ull;
2154         	res.w[0] = 0x378d8e6400000000ull - R64;
2155               } else {
2156         	// We want R64 to be the top digit of C4, but we actually
2157         	// obtained (C4 * 10^(-q4+1))RN; a correction may be needed,
2158         	// because the top digit is (C4 * 10^(-q4+1))RZ
2159         	// however, if incr_exp = 1 then R64 = 10 with certainty
2160         	if (incr_exp) {
2161         	  R64 = 10;
2162         	}
2163         	// the result is inexact as C4 has more than 1 significant digit
2164         	// and C3 * 10^scale = 10^33
2165         	// example of case that is treated here:
2166         	// 100...0 * 10^e3 - 0.41 * 10^e3 =
2167         	// 0999...9.59 * 10^e3 -> rounds to 99...96*10^(e3-1)
2168         	// note that (e3 > expmin}
2169         	// in order to round, subtract R64 from 10^34 and then compare
2170         	// C4 - R64 * 10^(q4-1) with 1/2 ulp
2171         	// calculate 10^34 - R64
2172         	res.w[1] = 0x0001ed09bead87c0ull;
2173         	res.w[0] = 0x378d8e6400000000ull - R64;
2174         	z_exp = z_exp - EXP_P1; // will be OR-ed with sign & significand
2175         	// calculate C4 - R64 * 10^(q4-1); this is a rare case and
2176         	// R64 is small, 1 <= R64 <= 9
2177         	e3 = e3 - 1;
2178         	if (is_inexact_lt_midpoint) {
2179         	  is_inexact_lt_midpoint = 0;
2180         	  is_inexact_gt_midpoint = 1;
2181         	} else if (is_inexact_gt_midpoint) {
2182         	  is_inexact_gt_midpoint = 0;
2183         	  is_inexact_lt_midpoint = 1;
2184         	} else if (is_midpoint_lt_even) {
2185         	  is_midpoint_lt_even = 0;
2186         	  is_midpoint_gt_even = 1;
2187         	} else if (is_midpoint_gt_even) {
2188         	  is_midpoint_gt_even = 0;
2189         	  is_midpoint_lt_even = 1;
2190         	} else {
2191         	  ;
2192         	}
2193         	// the result is always inexact, and never tiny
2194         	// check for overflow for RN
2195         	if (e3 > expmax) {
2196         	  if (rnd_mode == ROUNDING_TO_NEAREST) {
2197         	    res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
2198         	    res.w[0] = 0x0000000000000000ull;
2199         	    *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
2200         	  } else {
2201         	    rounding_correction (rnd_mode,
2202         				 is_inexact_lt_midpoint,
2203         				 is_inexact_gt_midpoint,
2204         				 is_midpoint_lt_even,
2205         				 is_midpoint_gt_even, e3, &res,
2206         				 pfpsf);
2207         	  }
2208         	  *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
2209         	  *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
2210         	  *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
2211         	  *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
2212         	  BID_SWAP128 (res);
2213         	  BID_RETURN (res)
2214         	}
2215         	// set the inexact flag
2216         	*pfpsf |= INEXACT_EXCEPTION;
2217         	res.w[1] =
2218         	  z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1];
2219         	if (rnd_mode != ROUNDING_TO_NEAREST) {
2220         	  rounding_correction (rnd_mode,
2221         			       is_inexact_lt_midpoint,
2222         			       is_inexact_gt_midpoint,
2223         			       is_midpoint_lt_even,
2224         			       is_midpoint_gt_even, e3, &res,
2225         			       pfpsf);
2226         	}
2227         	z_exp = res.w[1] & MASK_EXP;
2228               } // end result is inexact
2229             } // end q4 > 1
2230           } else { // if (e3 = emin)
2231             // if e3 = expmin the result is also tiny (the condition for
2232             // tininess is C4 > 050...0 [q4 digits] which is met because
2233             // the msd of C4 is not zero)
2234             // the result is tiny and inexact in all rounding modes;
2235             // it is either 100...0 or 0999...9 (use lt_half_ulp, eq_half_ulp,
2236             // gt_half_ulp to calculate)
2237             // if (lt_half_ulp || eq_half_ulp) res = 10^33 stays unchanged
2238 
2239             // p_sign != z_sign so swap gt_half_ulp and lt_half_ulp
2240             if (gt_half_ulp) { // res = 10^33 - 1
2241               res.w[1] = 0x0000314dc6448d93ull;
2242               res.w[0] = 0x38c15b09ffffffffull;
2243             } else {
2244               res.w[1] = 0x0000314dc6448d93ull;
2245               res.w[0] = 0x38c15b0a00000000ull;
2246             }
2247             res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
2248             *pfpsf |= UNDERFLOW_EXCEPTION; // inexact is set later
2249 
2250             if (eq_half_ulp) {
2251               is_midpoint_lt_even = 1; // if (z_sign), as if for absolute value
2252             } else if (lt_half_ulp) {
2253               is_inexact_gt_midpoint = 1; //if(z_sign), as if for absolute value
2254             } else { // if (gt_half_ulp)
2255               is_inexact_lt_midpoint = 1; //if(z_sign), as if for absolute value
2256             }
2257 
2258             if (rnd_mode != ROUNDING_TO_NEAREST) {
2259               rounding_correction (rnd_mode,
2260                   is_inexact_lt_midpoint,
2261                   is_inexact_gt_midpoint,
2262                   is_midpoint_lt_even,
2263                   is_midpoint_gt_even, e3, &res,
2264                   pfpsf);
2265               z_exp = res.w[1] & MASK_EXP;
2266             }
2267           } // end e3 = emin
2268           // set the inexact flag (if the result was not exact)
2269           if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
2270               is_midpoint_lt_even || is_midpoint_gt_even)
2271             *pfpsf |= INEXACT_EXCEPTION;
2272         } // end 10^33
2273       } // end if (p_sign != z_sign)
2274       res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
2275       *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
2276       *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
2277       *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
2278       *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
2279       BID_SWAP128 (res);
2280       BID_RETURN (res)
2281 
2282     } else if (((q3 <= delta && delta < p34 && p34 < delta + q4) || // Case (2)
2283         (q3 <= delta && delta + q4 <= p34) || // Case (3)
2284         (delta < q3 && p34 < delta + q4) || // Case (4)
2285         (delta < q3 && q3 <= delta + q4 && delta + q4 <= p34) || // Case (5)
2286         (delta + q4 < q3)) && // Case (6)
2287         !(delta <= 1 && p_sign != z_sign)) { // Case (2), (3), (4), (5) or (6)
2288 
2289       // the result has the sign of z
2290 
2291       if ((q3 <= delta && delta < p34 && p34 < delta + q4) || // Case (2)
2292           (delta < q3 && p34 < delta + q4)) { // Case (4)
2293         // round first the sum x * y + z with unbounded exponent
2294         // scale C3 up by scale = p34 - q3, 1 <= scale <= p34-1,
2295         // 1 <= scale <= 33
2296         // calculate res = C3 * 10^scale
2297         scale = p34 - q3;
2298         x0 = delta + q4 - p34;
2299       } else if (delta + q4 < q3) { // Case (6)
2300         // make Case (6) look like Case (3) or Case (5) with scale = 0
2301         // by scaling up C4 by 10^(q3 - delta - q4)
2302         scale = q3 - delta - q4; // 1 <= scale <= 33
2303         if (q4 <= 19) { // 1 <= scale <= 19; C4 fits in 64 bits
2304           if (scale <= 19) { // 10^scale fits in 64 bits
2305             // 64 x 64 C4.w[0] * ten2k64[scale]
2306             __mul_64x64_to_128MACH (P128, C4.w[0], ten2k64[scale]);
2307           } else { // 10^scale fits in 128 bits
2308             // 64 x 128 C4.w[0] * ten2k128[scale - 20]
2309             __mul_128x64_to_128 (P128, C4.w[0], ten2k128[scale - 20]);
2310           }
2311         } else { // C4 fits in 128 bits, but 10^scale must fit in 64 bits
2312           // 64 x 128 ten2k64[scale] * C4
2313           __mul_128x64_to_128 (P128, ten2k64[scale], C4);
2314         }
2315         C4.w[0] = P128.w[0];
2316         C4.w[1] = P128.w[1];
2317         // e4 does not need adjustment, as it is not used from this point on
2318         scale = 0;
2319         x0 = 0;
2320         // now Case (6) looks like Case (3) or Case (5) with scale = 0
2321       } else { // if Case (3) or Case (5)
2322         // Note: Case (3) is similar to Case (2), but scale differs and the
2323         // result is exact, unless it is tiny (so x0 = 0 when calculating the
2324         // result with unbounded exponent)
2325 
2326         // calculate first the sum x * y + z with unbounded exponent (exact)
2327         // scale C3 up by scale = delta + q4 - q3, 1 <= scale <= p34-1,
2328         // 1 <= scale <= 33
2329         // calculate res = C3 * 10^scale
2330         scale = delta + q4 - q3;
2331         x0 = 0;
2332         // Note: the comments which follow refer [mainly] to Case (2)]
2333       }
2334 
2335     case2_repeat:
2336       if (scale == 0) { // this could happen e.g. if we return to case2_repeat
2337         // or in Case (4)
2338         res.w[1] = C3.w[1];
2339         res.w[0] = C3.w[0];
2340       } else if (q3 <= 19) { // 1 <= scale <= 19; z fits in 64 bits
2341         if (scale <= 19) { // 10^scale fits in 64 bits
2342           // 64 x 64 C3.w[0] * ten2k64[scale]
2343           __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
2344         } else { // 10^scale fits in 128 bits
2345           // 64 x 128 C3.w[0] * ten2k128[scale - 20]
2346           __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
2347         }
2348       } else { // z fits in 128 bits, but 10^scale must fit in 64 bits
2349         // 64 x 128 ten2k64[scale] * C3
2350         __mul_128x64_to_128 (res, ten2k64[scale], C3);
2351       }
2352       // e3 is already calculated
2353       e3 = e3 - scale;
2354       // now res = C3 * 10^scale and e3 = e3 - scale
2355       // Note: C3 * 10^scale could be 10^34 if we returned to case2_repeat
2356       // because the result was too small
2357 
2358       // round C4 to nearest to q4 - x0 digits, where x0 = delta + q4 - p34,
2359       // 1 <= x0 <= min (q4 - 1, 2 * p34 - 1) <=> 1 <= x0 <= min (q4 - 1, 67)
2360       // Also: 1 <= q4 - x0 <= p34 -1 => 1 <= q4 - x0 <= 33 (so the result of
2361       // the rounding fits in 128 bits!)
2362       // x0 = delta + q4 - p34 (calculated before reaching case2_repeat)
2363       // because q3 + q4 - x0 <= P => x0 >= q3 + q4 - p34
2364       if (x0 == 0) { // this could happen only if we return to case2_repeat, or
2365         // for Case (3) or Case (6)
2366         R128.w[1] = C4.w[1];
2367         R128.w[0] = C4.w[0];
2368       } else if (q4 <= 18) {
2369         // 2 <= q4 <= 18, max(1, q3+q4-p34) <= x0 <= q4 - 1, 1 <= x0 <= 17
2370         round64_2_18 (q4, x0, C4.w[0], &R64, &incr_exp,
2371             &is_midpoint_lt_even, &is_midpoint_gt_even,
2372             &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
2373         if (incr_exp) {
2374           // R64 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 17
2375           R64 = ten2k64[q4 - x0];
2376         }
2377         R128.w[1] = 0;
2378         R128.w[0] = R64;
2379       } else if (q4 <= 38) {
2380         // 19 <= q4 <= 38, max(1, q3+q4-p34) <= x0 <= q4 - 1, 1 <= x0 <= 37
2381         P128.w[1] = C4.w[1];
2382         P128.w[0] = C4.w[0];
2383         round128_19_38 (q4, x0, P128, &R128, &incr_exp,
2384             &is_midpoint_lt_even, &is_midpoint_gt_even,
2385             &is_inexact_lt_midpoint,
2386             &is_inexact_gt_midpoint);
2387         if (incr_exp) {
2388           // R128 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 37
2389           if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
2390             R128.w[0] = ten2k64[q4 - x0];
2391             // R128.w[1] stays 0
2392           } else { // 20 <= q4 - x0 <= 37
2393             R128.w[0] = ten2k128[q4 - x0 - 20].w[0];
2394             R128.w[1] = ten2k128[q4 - x0 - 20].w[1];
2395           }
2396         }
2397       } else if (q4 <= 57) {
2398         // 38 <= q4 <= 57, max(1, q3+q4-p34) <= x0 <= q4 - 1, 5 <= x0 <= 56
2399         P192.w[2] = C4.w[2];
2400         P192.w[1] = C4.w[1];
2401         P192.w[0] = C4.w[0];
2402         round192_39_57 (q4, x0, P192, &R192, &incr_exp,
2403             &is_midpoint_lt_even, &is_midpoint_gt_even,
2404             &is_inexact_lt_midpoint,
2405             &is_inexact_gt_midpoint);
2406         // R192.w[2] is always 0
2407         if (incr_exp) {
2408           // R192 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 5, 1 <= q4 - x0 <= 52
2409           if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
2410             R192.w[0] = ten2k64[q4 - x0];
2411             // R192.w[1] stays 0
2412             // R192.w[2] stays 0
2413           } else { // 20 <= q4 - x0 <= 33
2414             R192.w[0] = ten2k128[q4 - x0 - 20].w[0];
2415             R192.w[1] = ten2k128[q4 - x0 - 20].w[1];
2416             // R192.w[2] stays 0
2417           }
2418         }
2419         R128.w[1] = R192.w[1];
2420         R128.w[0] = R192.w[0];
2421       } else {
2422         // 58 <= q4 <= 68, max(1, q3+q4-p34) <= x0 <= q4 - 1, 25 <= x0 <= 67
2423         round256_58_76 (q4, x0, C4, &R256, &incr_exp,
2424             &is_midpoint_lt_even, &is_midpoint_gt_even,
2425             &is_inexact_lt_midpoint,
2426             &is_inexact_gt_midpoint);
2427         // R256.w[3] and R256.w[2] are always 0
2428         if (incr_exp) {
2429           // R256 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 25, 1 <= q4 - x0 <= 43
2430           if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
2431             R256.w[0] = ten2k64[q4 - x0];
2432             // R256.w[1] stays 0
2433             // R256.w[2] stays 0
2434             // R256.w[3] stays 0
2435           } else { // 20 <= q4 - x0 <= 33
2436             R256.w[0] = ten2k128[q4 - x0 - 20].w[0];
2437             R256.w[1] = ten2k128[q4 - x0 - 20].w[1];
2438             // R256.w[2] stays 0
2439             // R256.w[3] stays 0
2440           }
2441         }
2442         R128.w[1] = R256.w[1];
2443         R128.w[0] = R256.w[0];
2444       }
2445       // now add C3 * 10^scale in res and the signed top (q4-x0) digits of C4,
2446       // rounded to nearest, which were copied into R128
2447       if (z_sign == p_sign) {
2448         lsb = res.w[0] & 0x01; // lsb of C3 * 10^scale
2449         // the sum can result in [up to] p34 or p34 + 1 digits
2450         res.w[0] = res.w[0] + R128.w[0];
2451         res.w[1] = res.w[1] + R128.w[1];
2452         if (res.w[0] < R128.w[0])
2453           res.w[1]++; // carry
2454         // if res > 10^34 - 1 need to increase x0 and decrease scale by 1
2455         if (res.w[1] > 0x0001ed09bead87c0ull ||
2456             (res.w[1] == 0x0001ed09bead87c0ull &&
2457              res.w[0] > 0x378d8e63ffffffffull)) {
2458           // avoid double rounding error
2459           is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
2460           is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
2461           is_midpoint_lt_even0 = is_midpoint_lt_even;
2462           is_midpoint_gt_even0 = is_midpoint_gt_even;
2463           is_inexact_lt_midpoint = 0;
2464           is_inexact_gt_midpoint = 0;
2465           is_midpoint_lt_even = 0;
2466           is_midpoint_gt_even = 0;
2467           P128.w[1] = res.w[1];
2468           P128.w[0] = res.w[0];
2469           round128_19_38 (35, 1, P128, &res, &incr_exp,
2470               &is_midpoint_lt_even, &is_midpoint_gt_even,
2471               &is_inexact_lt_midpoint,
2472               &is_inexact_gt_midpoint);
2473           // incr_exp is 0 with certainty in this case
2474           // avoid a double rounding error
2475           if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
2476               is_midpoint_lt_even) { // double rounding error upward
2477             // res = res - 1
2478             res.w[0]--;
2479             if (res.w[0] == 0xffffffffffffffffull)
2480               res.w[1]--;
2481             // Note: a double rounding error upward is not possible; for this
2482             // the result after the first rounding would have to be 99...95
2483             // (35 digits in all), possibly followed by a number of zeros; this
2484             // not possible in Cases (2)-(6) or (15)-(17) which may get here
2485             is_midpoint_lt_even = 0;
2486             is_inexact_lt_midpoint = 1;
2487           } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
2488               is_midpoint_gt_even) { // double rounding error downward
2489             // res = res + 1
2490             res.w[0]++;
2491             if (res.w[0] == 0)
2492               res.w[1]++;
2493             is_midpoint_gt_even = 0;
2494             is_inexact_gt_midpoint = 1;
2495           } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
2496         	     !is_inexact_lt_midpoint
2497         	     && !is_inexact_gt_midpoint) {
2498             // if this second rounding was exact the result may still be
2499             // inexact because of the first rounding
2500             if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
2501               is_inexact_gt_midpoint = 1;
2502             }
2503             if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
2504               is_inexact_lt_midpoint = 1;
2505             }
2506           } else if (is_midpoint_gt_even &&
2507         	     (is_inexact_gt_midpoint0
2508         	      || is_midpoint_lt_even0)) {
2509             // pulled up to a midpoint
2510             is_inexact_lt_midpoint = 1;
2511             is_inexact_gt_midpoint = 0;
2512             is_midpoint_lt_even = 0;
2513             is_midpoint_gt_even = 0;
2514           } else if (is_midpoint_lt_even &&
2515         	     (is_inexact_lt_midpoint0
2516         	      || is_midpoint_gt_even0)) {
2517             // pulled down to a midpoint
2518             is_inexact_lt_midpoint = 0;
2519             is_inexact_gt_midpoint = 1;
2520             is_midpoint_lt_even = 0;
2521             is_midpoint_gt_even = 0;
2522           } else {
2523             ;
2524           }
2525           // adjust exponent
2526           e3 = e3 + 1;
2527           if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
2528               !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
2529             if (is_midpoint_lt_even0 || is_midpoint_gt_even0 ||
2530         	is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0) {
2531               is_inexact_lt_midpoint = 1;
2532             }
2533           }
2534         } else {
2535           // this is the result rounded with unbounded exponent, unless a
2536           // correction is needed
2537           res.w[1] = res.w[1] & MASK_COEFF;
2538           if (lsb == 1) {
2539             if (is_midpoint_gt_even) {
2540               // res = res + 1
2541               is_midpoint_gt_even = 0;
2542               is_midpoint_lt_even = 1;
2543               res.w[0]++;
2544               if (res.w[0] == 0x0)
2545         	res.w[1]++;
2546               // check for rounding overflow
2547               if (res.w[1] == 0x0001ed09bead87c0ull &&
2548         	  res.w[0] == 0x378d8e6400000000ull) {
2549         	// res = 10^34 => rounding overflow
2550         	res.w[1] = 0x0000314dc6448d93ull;
2551         	res.w[0] = 0x38c15b0a00000000ull; // 10^33
2552         	e3++;
2553               }
2554             } else if (is_midpoint_lt_even) {
2555               // res = res - 1
2556               is_midpoint_lt_even = 0;
2557               is_midpoint_gt_even = 1;
2558               res.w[0]--;
2559               if (res.w[0] == 0xffffffffffffffffull)
2560         	res.w[1]--;
2561               // if the result is pure zero, the sign depends on the rounding
2562               // mode (x*y and z had opposite signs)
2563               if (res.w[1] == 0x0ull && res.w[0] == 0x0ull) {
2564         	if (rnd_mode != ROUNDING_DOWN)
2565         	  z_sign = 0x0000000000000000ull;
2566         	else
2567         	  z_sign = 0x8000000000000000ull;
2568         	// the exponent is max (e3, expmin)
2569         	res.w[1] = 0x0;
2570         	res.w[0] = 0x0;
2571         	*ptr_is_midpoint_lt_even = is_midpoint_lt_even;
2572         	*ptr_is_midpoint_gt_even = is_midpoint_gt_even;
2573         	*ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
2574         	*ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
2575         	BID_SWAP128 (res);
2576         	BID_RETURN (res)
2577               }
2578             } else {
2579               ;
2580             }
2581           }
2582         }
2583       } else { // if (z_sign != p_sign)
2584         lsb = res.w[0] & 0x01; // lsb of C3 * 10^scale; R128 contains rounded C4
2585         // used to swap rounding indicators if p_sign != z_sign
2586         // the sum can result in [up to] p34 or p34 - 1 digits
2587         tmp64 = res.w[0];
2588         res.w[0] = res.w[0] - R128.w[0];
2589         res.w[1] = res.w[1] - R128.w[1];
2590         if (res.w[0] > tmp64)
2591           res.w[1]--; // borrow
2592         // if res < 10^33 and exp > expmin need to decrease x0 and
2593         // increase scale by 1
2594         if (e3 > expmin && ((res.w[1] < 0x0000314dc6448d93ull ||
2595         		     (res.w[1] == 0x0000314dc6448d93ull &&
2596         		      res.w[0] < 0x38c15b0a00000000ull)) ||
2597         		    (is_inexact_lt_midpoint
2598         		     && res.w[1] == 0x0000314dc6448d93ull
2599         		     && res.w[0] == 0x38c15b0a00000000ull))
2600             && x0 >= 1) {
2601           x0 = x0 - 1;
2602           // first restore e3, otherwise it will be too small
2603           e3 = e3 + scale;
2604           scale = scale + 1;
2605           is_inexact_lt_midpoint = 0;
2606           is_inexact_gt_midpoint = 0;
2607           is_midpoint_lt_even = 0;
2608           is_midpoint_gt_even = 0;
2609           incr_exp = 0;
2610           goto case2_repeat;
2611         }
2612         // else this is the result rounded with unbounded exponent;
2613         // because the result has opposite sign to that of C4 which was
2614         // rounded, need to change the rounding indicators
2615         if (is_inexact_lt_midpoint) {
2616           is_inexact_lt_midpoint = 0;
2617           is_inexact_gt_midpoint = 1;
2618         } else if (is_inexact_gt_midpoint) {
2619           is_inexact_gt_midpoint = 0;
2620           is_inexact_lt_midpoint = 1;
2621         } else if (lsb == 0) {
2622           if (is_midpoint_lt_even) {
2623             is_midpoint_lt_even = 0;
2624             is_midpoint_gt_even = 1;
2625           } else if (is_midpoint_gt_even) {
2626             is_midpoint_gt_even = 0;
2627             is_midpoint_lt_even = 1;
2628           } else {
2629             ;
2630           }
2631         } else if (lsb == 1) {
2632           if (is_midpoint_lt_even) {
2633             // res = res + 1
2634             res.w[0]++;
2635             if (res.w[0] == 0x0)
2636               res.w[1]++;
2637             // check for rounding overflow
2638             if (res.w[1] == 0x0001ed09bead87c0ull &&
2639         	res.w[0] == 0x378d8e6400000000ull) {
2640               // res = 10^34 => rounding overflow
2641               res.w[1] = 0x0000314dc6448d93ull;
2642               res.w[0] = 0x38c15b0a00000000ull; // 10^33
2643               e3++;
2644             }
2645           } else if (is_midpoint_gt_even) {
2646             // res = res - 1
2647             res.w[0]--;
2648             if (res.w[0] == 0xffffffffffffffffull)
2649               res.w[1]--;
2650             // if the result is pure zero, the sign depends on the rounding
2651             // mode (x*y and z had opposite signs)
2652             if (res.w[1] == 0x0ull && res.w[0] == 0x0ull) {
2653               if (rnd_mode != ROUNDING_DOWN)
2654         	z_sign = 0x0000000000000000ull;
2655               else
2656         	z_sign = 0x8000000000000000ull;
2657               // the exponent is max (e3, expmin)
2658               res.w[1] = 0x0;
2659               res.w[0] = 0x0;
2660               *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
2661               *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
2662               *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
2663               *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
2664               BID_SWAP128 (res);
2665               BID_RETURN (res)
2666             }
2667           } else {
2668             ;
2669           }
2670         } else {
2671           ;
2672         }
2673       }
2674       // check for underflow
2675       if (e3 == expmin) { // and if significand < 10^33 => result is tiny
2676         if ((res.w[1] & MASK_COEFF) < 0x0000314dc6448d93ull ||
2677             ((res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull &&
2678              res.w[0] < 0x38c15b0a00000000ull)) {
2679           is_tiny = 1;
2680         }
2681       } else if (e3 < expmin) {
2682         // the result is tiny, so we must truncate more of res
2683         is_tiny = 1;
2684         x0 = expmin - e3;
2685         is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
2686         is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
2687         is_midpoint_lt_even0 = is_midpoint_lt_even;
2688         is_midpoint_gt_even0 = is_midpoint_gt_even;
2689         is_inexact_lt_midpoint = 0;
2690         is_inexact_gt_midpoint = 0;
2691         is_midpoint_lt_even = 0;
2692         is_midpoint_gt_even = 0;
2693         // determine the number of decimal digits in res
2694         if (res.w[1] == 0x0) {
2695           // between 1 and 19 digits
2696           for (ind = 1; ind <= 19; ind++) {
2697             if (res.w[0] < ten2k64[ind]) {
2698               break;
2699             }
2700           }
2701           // ind digits
2702         } else if (res.w[1] < ten2k128[0].w[1] ||
2703         	   (res.w[1] == ten2k128[0].w[1]
2704         	    && res.w[0] < ten2k128[0].w[0])) {
2705           // 20 digits
2706           ind = 20;
2707         } else { // between 21 and 38 digits
2708           for (ind = 1; ind <= 18; ind++) {
2709             if (res.w[1] < ten2k128[ind].w[1] ||
2710         	(res.w[1] == ten2k128[ind].w[1] &&
2711         	 res.w[0] < ten2k128[ind].w[0])) {
2712               break;
2713             }
2714           }
2715           // ind + 20 digits
2716           ind = ind + 20;
2717         }
2718 
2719         // at this point ind >= x0; because delta >= 2 on this path, the case
2720         // ind = x0 can occur only in Case (2) or case (3), when C3 has one
2721         // digit (q3 = 1) equal to 1 (C3 = 1), e3 is expmin (e3 = expmin),
2722         // the signs of x * y and z are opposite, and through cancellation
2723         // the most significant decimal digit in res has the weight
2724         // 10^(emin-1); however, it is clear that in this case the most
2725         // significant digit is 9, so the result before rounding is
2726         // 0.9... * 10^emin
2727         // Otherwise, ind > x0 because there are non-zero decimal digits in the
2728         // result with weight of at least 10^emin, and correction for underflow
2729         //  can be carried out using the round*_*_2_* () routines
2730         if (x0 == ind) { // the result before rounding is 0.9... * 10^emin
2731           res.w[1] = 0x0;
2732           res.w[0] = 0x1;
2733           is_inexact_gt_midpoint = 1;
2734         } else if (ind <= 18) { // check that 2 <= ind
2735           // 2 <= ind <= 18, 1 <= x0 <= 17
2736           round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
2737         		&is_midpoint_lt_even, &is_midpoint_gt_even,
2738         		&is_inexact_lt_midpoint,
2739         		&is_inexact_gt_midpoint);
2740           if (incr_exp) {
2741             // R64 = 10^(ind-x0), 1 <= ind - x0 <= ind - 1, 1 <= ind - x0 <= 17
2742             R64 = ten2k64[ind - x0];
2743           }
2744           res.w[1] = 0;
2745           res.w[0] = R64;
2746         } else if (ind <= 38) {
2747           // 19 <= ind <= 38
2748           P128.w[1] = res.w[1];
2749           P128.w[0] = res.w[0];
2750           round128_19_38 (ind, x0, P128, &res, &incr_exp,
2751         		  &is_midpoint_lt_even, &is_midpoint_gt_even,
2752         		  &is_inexact_lt_midpoint,
2753         		  &is_inexact_gt_midpoint);
2754           if (incr_exp) {
2755             // R128 = 10^(ind-x0), 1 <= ind - x0 <= ind - 1, 1 <= ind - x0 <= 37
2756             if (ind - x0 <= 19) { // 1 <= ind - x0 <= 19
2757               res.w[0] = ten2k64[ind - x0];
2758               // res.w[1] stays 0
2759             } else { // 20 <= ind - x0 <= 37
2760               res.w[0] = ten2k128[ind - x0 - 20].w[0];
2761               res.w[1] = ten2k128[ind - x0 - 20].w[1];
2762             }
2763           }
2764         }
2765         // avoid a double rounding error
2766         if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
2767             is_midpoint_lt_even) { // double rounding error upward
2768           // res = res - 1
2769           res.w[0]--;
2770           if (res.w[0] == 0xffffffffffffffffull)
2771             res.w[1]--;
2772           // Note: a double rounding error upward is not possible; for this
2773           // the result after the first rounding would have to be 99...95
2774           // (35 digits in all), possibly followed by a number of zeros; this
2775           // not possible in Cases (2)-(6) which may get here
2776           is_midpoint_lt_even = 0;
2777           is_inexact_lt_midpoint = 1;
2778         } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
2779             is_midpoint_gt_even) { // double rounding error downward
2780           // res = res + 1
2781           res.w[0]++;
2782           if (res.w[0] == 0)
2783             res.w[1]++;
2784           is_midpoint_gt_even = 0;
2785           is_inexact_gt_midpoint = 1;
2786         } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
2787         	   !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
2788           // if this second rounding was exact the result may still be
2789           // inexact because of the first rounding
2790           if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
2791             is_inexact_gt_midpoint = 1;
2792           }
2793           if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
2794             is_inexact_lt_midpoint = 1;
2795           }
2796         } else if (is_midpoint_gt_even &&
2797         	   (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
2798           // pulled up to a midpoint
2799           is_inexact_lt_midpoint = 1;
2800           is_inexact_gt_midpoint = 0;
2801           is_midpoint_lt_even = 0;
2802           is_midpoint_gt_even = 0;
2803         } else if (is_midpoint_lt_even &&
2804         	   (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
2805           // pulled down to a midpoint
2806           is_inexact_lt_midpoint = 0;
2807           is_inexact_gt_midpoint = 1;
2808           is_midpoint_lt_even = 0;
2809           is_midpoint_gt_even = 0;
2810         } else {
2811           ;
2812         }
2813         // adjust exponent
2814         e3 = e3 + x0;
2815         if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
2816             !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
2817           if (is_midpoint_lt_even0 || is_midpoint_gt_even0 ||
2818               is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0) {
2819             is_inexact_lt_midpoint = 1;
2820           }
2821         }
2822       } else {
2823         ; // not underflow
2824       }
2825       // check for inexact result
2826       if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
2827           is_midpoint_lt_even || is_midpoint_gt_even) {
2828         // set the inexact flag
2829         *pfpsf |= INEXACT_EXCEPTION;
2830         if (is_tiny)
2831           *pfpsf |= UNDERFLOW_EXCEPTION;
2832       }
2833       // now check for significand = 10^34 (may have resulted from going
2834       // back to case2_repeat)
2835       if (res.w[1] == 0x0001ed09bead87c0ull &&
2836           res.w[0] == 0x378d8e6400000000ull) { // if  res = 10^34
2837         res.w[1] = 0x0000314dc6448d93ull; // res = 10^33
2838         res.w[0] = 0x38c15b0a00000000ull;
2839         e3 = e3 + 1;
2840       }
2841       res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1];
2842       // check for overflow
2843       if (rnd_mode == ROUNDING_TO_NEAREST && e3 > expmax) {
2844         res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
2845         res.w[0] = 0x0000000000000000ull;
2846         *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
2847       }
2848       if (rnd_mode != ROUNDING_TO_NEAREST) {
2849         rounding_correction (rnd_mode,
2850         		     is_inexact_lt_midpoint,
2851         		     is_inexact_gt_midpoint,
2852         		     is_midpoint_lt_even, is_midpoint_gt_even,
2853         		     e3, &res, pfpsf);
2854       }
2855       *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
2856       *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
2857       *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
2858       *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
2859       BID_SWAP128 (res);
2860       BID_RETURN (res)
2861 
2862     } else {
2863 
2864       // we get here only if delta <= 1 in Cases (2), (3), (4), (5), or (6) and
2865       // the signs of x*y and z are opposite; in these cases massive
2866       // cancellation can occur, so it is better to scale either C3 or C4 and
2867       // to perform the subtraction before rounding; rounding is performed
2868       // next, depending on the number of decimal digits in the result and on
2869       // the exponent value
2870       // Note: overlow is not possible in this case
2871       // this is similar to Cases (15), (16), and (17)
2872 
2873       if (delta + q4 < q3) { // from Case (6)
2874         // Case (6) with 0<= delta <= 1 is similar to Cases (15), (16), and
2875         // (17) if we swap (C3, C4), (q3, q4), (e3, e4), (z_sign, p_sign)
2876         // and call add_and_round; delta stays positive
2877         // C4.w[3] = 0 and C4.w[2] = 0, so swap just the low part of C4 with C3
2878         P128.w[1] = C3.w[1];
2879         P128.w[0] = C3.w[0];
2880         C3.w[1] = C4.w[1];
2881         C3.w[0] = C4.w[0];
2882         C4.w[1] = P128.w[1];
2883         C4.w[0] = P128.w[0];
2884         ind = q3;
2885         q3 = q4;
2886         q4 = ind;
2887         ind = e3;
2888         e3 = e4;
2889         e4 = ind;
2890         tmp_sign = z_sign;
2891         z_sign = p_sign;
2892         p_sign = tmp_sign;
2893       } else { // from Cases (2), (3), (4), (5)
2894         // In Cases (2), (3), (4), (5) with 0 <= delta <= 1 C3 has to be
2895         // scaled up by q4 + delta - q3; this is the same as in Cases (15),
2896         // (16), and (17) if we just change the sign of delta
2897         delta = -delta;
2898       }
2899       add_and_round (q3, q4, e4, delta, p34, z_sign, p_sign, C3, C4,
2900         	     rnd_mode, &is_midpoint_lt_even,
2901         	     &is_midpoint_gt_even, &is_inexact_lt_midpoint,
2902         	     &is_inexact_gt_midpoint, pfpsf, &res);
2903       *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
2904       *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
2905       *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
2906       *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
2907       BID_SWAP128 (res);
2908       BID_RETURN (res)
2909 
2910     }
2911 
2912   } else { // if delta < 0
2913 
2914     delta = -delta;
2915 
2916     if (p34 < q4 && q4 <= delta) { // Case (7)
2917 
2918       // truncate C4 to p34 digits into res
2919       // x = q4-p34, 1 <= x <= 34 because 35 <= q4 <= 68
2920       x0 = q4 - p34;
2921       if (q4 <= 38) {
2922         P128.w[1] = C4.w[1];
2923         P128.w[0] = C4.w[0];
2924         round128_19_38 (q4, x0, P128, &res, &incr_exp,
2925         		&is_midpoint_lt_even, &is_midpoint_gt_even,
2926         		&is_inexact_lt_midpoint,
2927         		&is_inexact_gt_midpoint);
2928       } else if (q4 <= 57) { // 35 <= q4 <= 57
2929         P192.w[2] = C4.w[2];
2930         P192.w[1] = C4.w[1];
2931         P192.w[0] = C4.w[0];
2932         round192_39_57 (q4, x0, P192, &R192, &incr_exp,
2933         		&is_midpoint_lt_even, &is_midpoint_gt_even,
2934         		&is_inexact_lt_midpoint,
2935         		&is_inexact_gt_midpoint);
2936         res.w[0] = R192.w[0];
2937         res.w[1] = R192.w[1];
2938       } else { // if (q4 <= 68)
2939         round256_58_76 (q4, x0, C4, &R256, &incr_exp,
2940         		&is_midpoint_lt_even, &is_midpoint_gt_even,
2941         		&is_inexact_lt_midpoint,
2942         		&is_inexact_gt_midpoint);
2943         res.w[0] = R256.w[0];
2944         res.w[1] = R256.w[1];
2945       }
2946       e4 = e4 + x0;
2947       if (incr_exp) {
2948         e4 = e4 + 1;
2949       }
2950       if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
2951           !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
2952         // if C4 rounded to p34 digits is exact then the result is inexact,
2953         // in a way that depends on the signs of x * y and z
2954         if (p_sign == z_sign) {
2955           is_inexact_lt_midpoint = 1;
2956         } else { // if (p_sign != z_sign)
2957           if (res.w[1] != 0x0000314dc6448d93ull ||
2958               res.w[0] != 0x38c15b0a00000000ull) { // res != 10^33
2959             is_inexact_gt_midpoint = 1;
2960           } else { // res = 10^33 and exact is a special case
2961             // if C3 < 1/2 ulp then res = 10^33 and is_inexact_gt_midpoint = 1
2962             // if C3 = 1/2 ulp then res = 10^33 and is_midpoint_lt_even = 1
2963             // if C3 > 1/2 ulp then res = 10^34-1 and is_inexact_lt_midpoint = 1
2964             // Note: ulp is really ulp/10 (after borrow which propagates to msd)
2965             if (delta > p34 + 1) { // C3 < 1/2
2966               // res = 10^33, unchanged
2967               is_inexact_gt_midpoint = 1;
2968             } else { // if (delta == p34 + 1)
2969               if (q3 <= 19) {
2970         	if (C3.w[0] < midpoint64[q3 - 1]) { // C3 < 1/2 ulp
2971         	  // res = 10^33, unchanged
2972         	  is_inexact_gt_midpoint = 1;
2973         	} else if (C3.w[0] == midpoint64[q3 - 1]) { // C3 = 1/2 ulp
2974         	  // res = 10^33, unchanged
2975         	  is_midpoint_lt_even = 1;
2976         	} else { // if (C3.w[0] > midpoint64[q3-1]), C3 > 1/2 ulp
2977         	  res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
2978         	  res.w[0] = 0x378d8e63ffffffffull;
2979         	  e4 = e4 - 1;
2980         	  is_inexact_lt_midpoint = 1;
2981         	}
2982               } else { // if (20 <= q3 <=34)
2983         	if (C3.w[1] < midpoint128[q3 - 20].w[1] ||
2984                     (C3.w[1] == midpoint128[q3 - 20].w[1] &&
2985                     C3.w[0] < midpoint128[q3 - 20].w[0])) { // C3 < 1/2 ulp
2986         	  // res = 10^33, unchanged
2987         	  is_inexact_gt_midpoint = 1;
2988         	} else if (C3.w[1] == midpoint128[q3 - 20].w[1] &&
2989                     C3.w[0] == midpoint128[q3 - 20].w[0]) { // C3 = 1/2 ulp
2990         	  // res = 10^33, unchanged
2991         	  is_midpoint_lt_even = 1;
2992         	} else { // if (C3 > midpoint128[q3-20]), C3 > 1/2 ulp
2993         	  res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
2994         	  res.w[0] = 0x378d8e63ffffffffull;
2995         	  e4 = e4 - 1;
2996         	  is_inexact_lt_midpoint = 1;
2997         	}
2998               }
2999             }
3000           }
3001         }
3002       } else if (is_midpoint_lt_even) {
3003         if (z_sign != p_sign) {
3004           // needs correction: res = res - 1
3005           res.w[0] = res.w[0] - 1;
3006           if (res.w[0] == 0xffffffffffffffffull)
3007             res.w[1]--;
3008           // if it is (10^33-1)*10^e4 then the corect result is
3009           // (10^34-1)*10(e4-1)
3010           if (res.w[1] == 0x0000314dc6448d93ull &&
3011               res.w[0] == 0x38c15b09ffffffffull) {
3012             res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
3013             res.w[0] = 0x378d8e63ffffffffull;
3014             e4 = e4 - 1;
3015           }
3016           is_midpoint_lt_even = 0;
3017           is_inexact_lt_midpoint = 1;
3018         } else { // if (z_sign == p_sign)
3019           is_midpoint_lt_even = 0;
3020           is_inexact_gt_midpoint = 1;
3021         }
3022       } else if (is_midpoint_gt_even) {
3023         if (z_sign == p_sign) {
3024           // needs correction: res = res + 1 (cannot cross in the next binade)
3025           res.w[0] = res.w[0] + 1;
3026           if (res.w[0] == 0x0000000000000000ull)
3027             res.w[1]++;
3028           is_midpoint_gt_even = 0;
3029           is_inexact_gt_midpoint = 1;
3030         } else { // if (z_sign != p_sign)
3031           is_midpoint_gt_even = 0;
3032           is_inexact_lt_midpoint = 1;
3033         }
3034       } else {
3035         ; // the rounded result is already correct
3036       }
3037       // check for overflow
3038       if (rnd_mode == ROUNDING_TO_NEAREST && e4 > expmax) {
3039         res.w[1] = p_sign | 0x7800000000000000ull;
3040         res.w[0] = 0x0000000000000000ull;
3041         *pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
3042       } else { // no overflow or not RN
3043         p_exp = ((UINT64) (e4 + 6176) << 49);
3044         res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1];
3045       }
3046       if (rnd_mode != ROUNDING_TO_NEAREST) {
3047         rounding_correction (rnd_mode,
3048         		     is_inexact_lt_midpoint,
3049         		     is_inexact_gt_midpoint,
3050         		     is_midpoint_lt_even, is_midpoint_gt_even,
3051         		     e4, &res, pfpsf);
3052       }
3053       if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
3054           is_midpoint_lt_even || is_midpoint_gt_even) {
3055         // set the inexact flag
3056         *pfpsf |= INEXACT_EXCEPTION;
3057       }
3058       *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
3059       *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
3060       *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
3061       *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
3062       BID_SWAP128 (res);
3063       BID_RETURN (res)
3064 
3065     } else if ((q4 <= p34 && p34 <= delta) || // Case (8)
3066         (q4 <= delta && delta < p34 && p34 < delta + q3) || // Case (9)
3067         (q4 <= delta && delta + q3 <= p34) || // Case (10)
3068         (delta < q4 && q4 <= p34 && p34 < delta + q3) || // Case (13)
3069         (delta < q4 && q4 <= delta + q3 && delta + q3 <= p34) || // Case (14)
3070         (delta + q3 < q4 && q4 <= p34)) { // Case (18)
3071 
3072       // Case (8) is similar to Case (1), with C3 and C4 swapped
3073       // Case (9) is similar to Case (2), with C3 and C4 swapped
3074       // Case (10) is similar to Case (3), with C3 and C4 swapped
3075       // Case (13) is similar to Case (4), with C3 and C4 swapped
3076       // Case (14) is similar to Case (5), with C3 and C4 swapped
3077       // Case (18) is similar to Case (6), with C3 and C4 swapped
3078 
3079       // swap (C3, C4), (q3, q4), (e3, 34), (z_sign, p_sign), (z_exp, p_exp)
3080       // and go back to delta_ge_zero
3081       // C4.w[3] = 0 and C4.w[2] = 0, so swap just the low part of C4 with C3
3082       P128.w[1] = C3.w[1];
3083       P128.w[0] = C3.w[0];
3084       C3.w[1] = C4.w[1];
3085       C3.w[0] = C4.w[0];
3086       C4.w[1] = P128.w[1];
3087       C4.w[0] = P128.w[0];
3088       ind = q3;
3089       q3 = q4;
3090       q4 = ind;
3091       ind = e3;
3092       e3 = e4;
3093       e4 = ind;
3094       tmp_sign = z_sign;
3095       z_sign = p_sign;
3096       p_sign = tmp_sign;
3097       tmp.ui64 = z_exp;
3098       z_exp = p_exp;
3099       p_exp = tmp.ui64;
3100       goto delta_ge_zero;
3101 
3102     } else if ((p34 <= delta && delta < q4 && q4 < delta + q3) || // Case (11)
3103                (delta < p34 && p34 < q4 && q4 < delta + q3)) { // Case (12)
3104 
3105       // round C3 to nearest to q3 - x0 digits, where x0 = e4 - e3,
3106       // 1 <= x0 <= q3 - 1 <= p34 - 1
3107       x0 = e4 - e3; // or x0 = delta + q3 - q4
3108       if (q3 <= 18) { // 2 <= q3 <= 18
3109         round64_2_18 (q3, x0, C3.w[0], &R64, &incr_exp,
3110         	      &is_midpoint_lt_even, &is_midpoint_gt_even,
3111         	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
3112         // C3.w[1] = 0;
3113         C3.w[0] = R64;
3114       } else if (q3 <= 38) {
3115         round128_19_38 (q3, x0, C3, &R128, &incr_exp,
3116         		&is_midpoint_lt_even, &is_midpoint_gt_even,
3117         		&is_inexact_lt_midpoint,
3118         		&is_inexact_gt_midpoint);
3119         C3.w[1] = R128.w[1];
3120         C3.w[0] = R128.w[0];
3121       }
3122       // the rounded result has q3 - x0 digits
3123       // we want the exponent to be e4, so if incr_exp = 1 then
3124       // multiply the rounded result by 10 - it will still fit in 113 bits
3125       if (incr_exp) {
3126         // 64 x 128 -> 128
3127         P128.w[1] = C3.w[1];
3128         P128.w[0] = C3.w[0];
3129         __mul_64x128_to_128 (C3, ten2k64[1], P128);
3130       }
3131       e3 = e3 + x0; // this is e4
3132       // now add/subtract the 256-bit C4 and the new (and shorter) 128-bit C3;
3133       // the result will have the sign of x * y; the exponent is e4
3134       R256.w[3] = 0;
3135       R256.w[2] = 0;
3136       R256.w[1] = C3.w[1];
3137       R256.w[0] = C3.w[0];
3138       if (p_sign == z_sign) { // R256 = C4 + R256
3139         add256 (C4, R256, &R256);
3140       } else { // if (p_sign != z_sign) { // R256 = C4 - R256
3141         sub256 (C4, R256, &R256); // the result cannot be pure zero
3142         // because the result has opposite sign to that of R256 which was
3143         // rounded, need to change the rounding indicators
3144         lsb = C4.w[0] & 0x01;
3145         if (is_inexact_lt_midpoint) {
3146           is_inexact_lt_midpoint = 0;
3147           is_inexact_gt_midpoint = 1;
3148         } else if (is_inexact_gt_midpoint) {
3149           is_inexact_gt_midpoint = 0;
3150           is_inexact_lt_midpoint = 1;
3151         } else if (lsb == 0) {
3152           if (is_midpoint_lt_even) {
3153             is_midpoint_lt_even = 0;
3154             is_midpoint_gt_even = 1;
3155           } else if (is_midpoint_gt_even) {
3156             is_midpoint_gt_even = 0;
3157             is_midpoint_lt_even = 1;
3158           } else {
3159             ;
3160           }
3161         } else if (lsb == 1) {
3162           if (is_midpoint_lt_even) {
3163             // res = res + 1
3164             R256.w[0]++;
3165             if (R256.w[0] == 0x0) {
3166               R256.w[1]++;
3167               if (R256.w[1] == 0x0) {
3168         	R256.w[2]++;
3169         	if (R256.w[2] == 0x0) {
3170         	  R256.w[3]++;
3171         	}
3172               }
3173             }
3174             // no check for rounding overflow - R256 was a difference
3175           } else if (is_midpoint_gt_even) {
3176             // res = res - 1
3177             R256.w[0]--;
3178             if (R256.w[0] == 0xffffffffffffffffull) {
3179               R256.w[1]--;
3180               if (R256.w[1] == 0xffffffffffffffffull) {
3181         	R256.w[2]--;
3182         	if (R256.w[2] == 0xffffffffffffffffull) {
3183         	  R256.w[3]--;
3184         	}
3185               }
3186             }
3187           } else {
3188             ;
3189           }
3190         } else {
3191           ;
3192         }
3193       }
3194       // determine the number of decimal digits in R256
3195       ind = nr_digits256 (R256); // ind >= p34
3196       // if R256 is sum, then ind > p34; if R256 is a difference, then
3197       // ind >= p34; this means that we can calculate the result rounded to
3198       // the destination precision, with unbounded exponent, starting from R256
3199       // and using the indicators from the rounding of C3 to avoid a double
3200       // rounding error
3201 
3202       if (ind < p34) {
3203         ;
3204       } else if (ind == p34) {
3205         // the result rounded to the destination precision with
3206         // unbounded exponent
3207         // is (-1)^p_sign * R256 * 10^e4
3208         res.w[1] = R256.w[1];
3209         res.w[0] = R256.w[0];
3210       } else { // if (ind > p34)
3211         // if more than P digits, round to nearest to P digits
3212         // round R256 to p34 digits
3213         x0 = ind - p34; // 1 <= x0 <= 34 as 35 <= ind <= 68
3214         // save C3 rounding indicators to help avoid double rounding error
3215         is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
3216         is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
3217         is_midpoint_lt_even0 = is_midpoint_lt_even;
3218         is_midpoint_gt_even0 = is_midpoint_gt_even;
3219         // initialize rounding indicators
3220         is_inexact_lt_midpoint = 0;
3221         is_inexact_gt_midpoint = 0;
3222         is_midpoint_lt_even = 0;
3223         is_midpoint_gt_even = 0;
3224         // round to p34 digits; the result fits in 113 bits
3225         if (ind <= 38) {
3226           P128.w[1] = R256.w[1];
3227           P128.w[0] = R256.w[0];
3228           round128_19_38 (ind, x0, P128, &R128, &incr_exp,
3229         		  &is_midpoint_lt_even, &is_midpoint_gt_even,
3230         		  &is_inexact_lt_midpoint,
3231         		  &is_inexact_gt_midpoint);
3232         } else if (ind <= 57) {
3233           P192.w[2] = R256.w[2];
3234           P192.w[1] = R256.w[1];
3235           P192.w[0] = R256.w[0];
3236           round192_39_57 (ind, x0, P192, &R192, &incr_exp,
3237         		  &is_midpoint_lt_even, &is_midpoint_gt_even,
3238         		  &is_inexact_lt_midpoint,
3239         		  &is_inexact_gt_midpoint);
3240           R128.w[1] = R192.w[1];
3241           R128.w[0] = R192.w[0];
3242         } else { // if (ind <= 68)
3243           round256_58_76 (ind, x0, R256, &R256, &incr_exp,
3244         		  &is_midpoint_lt_even, &is_midpoint_gt_even,
3245         		  &is_inexact_lt_midpoint,
3246         		  &is_inexact_gt_midpoint);
3247           R128.w[1] = R256.w[1];
3248           R128.w[0] = R256.w[0];
3249         }
3250         // the rounded result has p34 = 34 digits
3251         e4 = e4 + x0 + incr_exp;
3252 
3253         res.w[1] = R128.w[1];
3254         res.w[0] = R128.w[0];
3255 
3256         // avoid a double rounding error
3257         if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
3258             is_midpoint_lt_even) { // double rounding error upward
3259           // res = res - 1
3260           res.w[0]--;
3261           if (res.w[0] == 0xffffffffffffffffull)
3262             res.w[1]--;
3263           is_midpoint_lt_even = 0;
3264           is_inexact_lt_midpoint = 1;
3265           // Note: a double rounding error upward is not possible; for this
3266           // the result after the first rounding would have to be 99...95
3267           // (35 digits in all), possibly followed by a number of zeros; this
3268           // not possible in Cases (2)-(6) or (15)-(17) which may get here
3269           // if this is 10^33 - 1 make it 10^34 - 1 and decrement exponent
3270           if (res.w[1] == 0x0000314dc6448d93ull &&
3271             res.w[0] == 0x38c15b09ffffffffull) { // 10^33 - 1
3272             res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
3273             res.w[0] = 0x378d8e63ffffffffull;
3274             e4--;
3275           }
3276         } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
3277             is_midpoint_gt_even) { // double rounding error downward
3278           // res = res + 1
3279           res.w[0]++;
3280           if (res.w[0] == 0)
3281             res.w[1]++;
3282           is_midpoint_gt_even = 0;
3283           is_inexact_gt_midpoint = 1;
3284         } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
3285         	   !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
3286           // if this second rounding was exact the result may still be
3287           // inexact because of the first rounding
3288           if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
3289             is_inexact_gt_midpoint = 1;
3290           }
3291           if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
3292             is_inexact_lt_midpoint = 1;
3293           }
3294         } else if (is_midpoint_gt_even &&
3295         	   (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
3296           // pulled up to a midpoint
3297           is_inexact_lt_midpoint = 1;
3298           is_inexact_gt_midpoint = 0;
3299           is_midpoint_lt_even = 0;
3300           is_midpoint_gt_even = 0;
3301         } else if (is_midpoint_lt_even &&
3302         	   (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
3303           // pulled down to a midpoint
3304           is_inexact_lt_midpoint = 0;
3305           is_inexact_gt_midpoint = 1;
3306           is_midpoint_lt_even = 0;
3307           is_midpoint_gt_even = 0;
3308         } else {
3309           ;
3310         }
3311       }
3312 
3313       // determine tininess
3314       if (rnd_mode == ROUNDING_TO_NEAREST) {
3315         if (e4 < expmin) {
3316           is_tiny = 1; // for other rounding modes apply correction
3317         }
3318       } else {
3319         // for RM, RP, RZ, RA apply correction in order to determine tininess
3320         // but do not save the result; apply the correction to
3321         // (-1)^p_sign * res * 10^0
3322         P128.w[1] = p_sign | 0x3040000000000000ull | res.w[1];
3323         P128.w[0] = res.w[0];
3324         rounding_correction (rnd_mode,
3325         		     is_inexact_lt_midpoint,
3326         		     is_inexact_gt_midpoint,
3327         		     is_midpoint_lt_even, is_midpoint_gt_even,
3328         		     0, &P128, pfpsf);
3329         scale = ((P128.w[1] & MASK_EXP) >> 49) - 6176; // -1, 0, or +1
3330         // the number of digits in the significand is p34 = 34
3331         if (e4 + scale < expmin) {
3332           is_tiny = 1;
3333         }
3334       }
3335 
3336       // the result rounded to the destination precision with unbounded exponent
3337       // is (-1)^p_sign * res * 10^e4
3338       res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1]; // RN
3339       // res.w[0] unchanged;
3340       // Note: res is correct only if expmin <= e4 <= expmax
3341       ind = p34; // the number of decimal digits in the signifcand of res
3342 
3343       // at this point we have the result rounded with unbounded exponent in
3344       // res and we know its tininess:
3345       // res = (-1)^p_sign * significand * 10^e4,
3346       // where q (significand) = ind = p34
3347       // Note: res is correct only if expmin <= e4 <= expmax
3348 
3349       // check for overflow if RN
3350       if (rnd_mode == ROUNDING_TO_NEAREST
3351           && (ind + e4) > (p34 + expmax)) {
3352         res.w[1] = p_sign | 0x7800000000000000ull;
3353         res.w[0] = 0x0000000000000000ull;
3354         *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
3355         *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
3356         *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
3357         *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
3358         *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
3359         BID_SWAP128 (res);
3360         BID_RETURN (res)
3361       } // else not overflow or not RN, so continue
3362 
3363       // from this point on this is similar to the last part of the computation
3364       // for Cases (15), (16), (17)
3365 
3366       // if (e4 >= expmin) we have the result rounded with bounded exponent
3367       if (e4 < expmin) {
3368         x0 = expmin - e4; // x0 >= 1; the number of digits to chop off of res
3369         // where the result rounded [at most] once is
3370         //   (-1)^p_sign * significand_res * 10^e4
3371 
3372         // avoid double rounding error
3373         is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
3374         is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
3375         is_midpoint_lt_even0 = is_midpoint_lt_even;
3376         is_midpoint_gt_even0 = is_midpoint_gt_even;
3377         is_inexact_lt_midpoint = 0;
3378         is_inexact_gt_midpoint = 0;
3379         is_midpoint_lt_even = 0;
3380         is_midpoint_gt_even = 0;
3381 
3382         if (x0 > ind) {
3383           // nothing is left of res when moving the decimal point left x0 digits
3384           is_inexact_lt_midpoint = 1;
3385           res.w[1] = p_sign | 0x0000000000000000ull;
3386           res.w[0] = 0x0000000000000000ull;
3387           e4 = expmin;
3388         } else if (x0 == ind) { // 1 <= x0 = ind <= p34 = 34
3389           // this is <, =, or > 1/2 ulp
3390           // compare the ind-digit value in the significand of res with
3391           // 1/2 ulp = 5*10^(ind-1), i.e. determine whether it is
3392           // less than, equal to, or greater than 1/2 ulp (significand of res)
3393           R128.w[1] = res.w[1] & MASK_COEFF;
3394           R128.w[0] = res.w[0];
3395           if (ind <= 19) {
3396             if (R128.w[0] < midpoint64[ind - 1]) { // < 1/2 ulp
3397               lt_half_ulp = 1;
3398               is_inexact_lt_midpoint = 1;
3399             } else if (R128.w[0] == midpoint64[ind - 1]) { // = 1/2 ulp
3400               eq_half_ulp = 1;
3401               is_midpoint_gt_even = 1;
3402             } else { // > 1/2 ulp
3403               gt_half_ulp = 1;
3404               is_inexact_gt_midpoint = 1;
3405             }
3406           } else { // if (ind <= 38)
3407             if (R128.w[1] < midpoint128[ind - 20].w[1] ||
3408                 (R128.w[1] == midpoint128[ind - 20].w[1] &&
3409                 R128.w[0] < midpoint128[ind - 20].w[0])) { // < 1/2 ulp
3410               lt_half_ulp = 1;
3411               is_inexact_lt_midpoint = 1;
3412             } else if (R128.w[1] == midpoint128[ind - 20].w[1] &&
3413                 R128.w[0] == midpoint128[ind - 20].w[0]) { // = 1/2 ulp
3414               eq_half_ulp = 1;
3415               is_midpoint_gt_even = 1;
3416             } else { // > 1/2 ulp
3417               gt_half_ulp = 1;
3418               is_inexact_gt_midpoint = 1;
3419             }
3420           }
3421           if (lt_half_ulp || eq_half_ulp) {
3422             // res = +0.0 * 10^expmin
3423             res.w[1] = 0x0000000000000000ull;
3424             res.w[0] = 0x0000000000000000ull;
3425           } else { // if (gt_half_ulp)
3426             // res = +1 * 10^expmin
3427             res.w[1] = 0x0000000000000000ull;
3428             res.w[0] = 0x0000000000000001ull;
3429           }
3430           res.w[1] = p_sign | res.w[1];
3431           e4 = expmin;
3432         } else { // if (1 <= x0 <= ind - 1 <= 33)
3433           // round the ind-digit result to ind - x0 digits
3434 
3435           if (ind <= 18) { // 2 <= ind <= 18
3436             round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
3437         		  &is_midpoint_lt_even, &is_midpoint_gt_even,
3438         		  &is_inexact_lt_midpoint,
3439         		  &is_inexact_gt_midpoint);
3440             res.w[1] = 0x0;
3441             res.w[0] = R64;
3442           } else if (ind <= 38) {
3443             P128.w[1] = res.w[1] & MASK_COEFF;
3444             P128.w[0] = res.w[0];
3445             round128_19_38 (ind, x0, P128, &res, &incr_exp,
3446         		    &is_midpoint_lt_even, &is_midpoint_gt_even,
3447         		    &is_inexact_lt_midpoint,
3448         		    &is_inexact_gt_midpoint);
3449           }
3450           e4 = e4 + x0; // expmin
3451           // we want the exponent to be expmin, so if incr_exp = 1 then
3452           // multiply the rounded result by 10 - it will still fit in 113 bits
3453           if (incr_exp) {
3454             // 64 x 128 -> 128
3455             P128.w[1] = res.w[1] & MASK_COEFF;
3456             P128.w[0] = res.w[0];
3457             __mul_64x128_to_128 (res, ten2k64[1], P128);
3458           }
3459           res.w[1] =
3460             p_sign | ((UINT64) (e4 + 6176) << 49) | (res.
3461         					     w[1] & MASK_COEFF);
3462           // avoid a double rounding error
3463           if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
3464                 is_midpoint_lt_even) { // double rounding error upward
3465             // res = res - 1
3466             res.w[0]--;
3467             if (res.w[0] == 0xffffffffffffffffull)
3468               res.w[1]--;
3469             // Note: a double rounding error upward is not possible; for this
3470             // the result after the first rounding would have to be 99...95
3471             // (35 digits in all), possibly followed by a number of zeros; this
3472             // not possible in this underflow case
3473             is_midpoint_lt_even = 0;
3474             is_inexact_lt_midpoint = 1;
3475           } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
3476                 is_midpoint_gt_even) { // double rounding error downward
3477             // res = res + 1
3478             res.w[0]++;
3479             if (res.w[0] == 0)
3480               res.w[1]++;
3481             is_midpoint_gt_even = 0;
3482             is_inexact_gt_midpoint = 1;
3483           } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
3484         	     !is_inexact_lt_midpoint
3485         	     && !is_inexact_gt_midpoint) {
3486             // if this second rounding was exact the result may still be
3487             // inexact because of the first rounding
3488             if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
3489               is_inexact_gt_midpoint = 1;
3490             }
3491             if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
3492               is_inexact_lt_midpoint = 1;
3493             }
3494           } else if (is_midpoint_gt_even &&
3495         	     (is_inexact_gt_midpoint0
3496         	      || is_midpoint_lt_even0)) {
3497             // pulled up to a midpoint
3498             is_inexact_lt_midpoint = 1;
3499             is_inexact_gt_midpoint = 0;
3500             is_midpoint_lt_even = 0;
3501             is_midpoint_gt_even = 0;
3502           } else if (is_midpoint_lt_even &&
3503         	     (is_inexact_lt_midpoint0
3504         	      || is_midpoint_gt_even0)) {
3505             // pulled down to a midpoint
3506             is_inexact_lt_midpoint = 0;
3507             is_inexact_gt_midpoint = 1;
3508             is_midpoint_lt_even = 0;
3509             is_midpoint_gt_even = 0;
3510           } else {
3511             ;
3512           }
3513         }
3514       }
3515       // res contains the correct result
3516       // apply correction if not rounding to nearest
3517       if (rnd_mode != ROUNDING_TO_NEAREST) {
3518         rounding_correction (rnd_mode,
3519         		     is_inexact_lt_midpoint,
3520         		     is_inexact_gt_midpoint,
3521         		     is_midpoint_lt_even, is_midpoint_gt_even,
3522         		     e4, &res, pfpsf);
3523       }
3524       if (is_midpoint_lt_even || is_midpoint_gt_even ||
3525           is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
3526         // set the inexact flag
3527         *pfpsf |= INEXACT_EXCEPTION;
3528         if (is_tiny)
3529           *pfpsf |= UNDERFLOW_EXCEPTION;
3530       }
3531       *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
3532       *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
3533       *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
3534       *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
3535       BID_SWAP128 (res);
3536       BID_RETURN (res)
3537 
3538     } else if ((p34 <= delta && delta + q3 <= q4) || // Case (15)
3539         (delta < p34 && p34 < delta + q3 && delta + q3 <= q4) || //Case (16)
3540         (delta + q3 <= p34 && p34 < q4)) { // Case (17)
3541 
3542       // calculate first the result rounded to the destination precision, with
3543       // unbounded exponent
3544 
3545       add_and_round (q3, q4, e4, delta, p34, z_sign, p_sign, C3, C4,
3546               rnd_mode, &is_midpoint_lt_even,
3547               &is_midpoint_gt_even, &is_inexact_lt_midpoint,
3548               &is_inexact_gt_midpoint, pfpsf, &res);
3549       *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
3550       *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
3551       *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
3552       *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
3553       BID_SWAP128 (res);
3554       BID_RETURN (res)
3555 
3556     } else {
3557       ;
3558     }
3559 
3560   } // end if delta < 0
3561 
3562   *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
3563   *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
3564   *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
3565   *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
3566   BID_SWAP128 (res);
3567   BID_RETURN (res)
3568 
3569 }
3570 
3571 
3572 #if DECIMAL_CALL_BY_REFERENCE
3573 void
3574 bid128_fma (UINT128 * pres, UINT128 * px, UINT128 * py, UINT128 * pz
3575             _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3576             _EXC_INFO_PARAM) {
3577   UINT128 x = *px, y = *py, z = *pz;
3578 #if !DECIMAL_GLOBAL_ROUNDING
3579   unsigned int rnd_mode = *prnd_mode;
3580 #endif
3581 #else
3582 UINT128
3583 bid128_fma (UINT128 x, UINT128 y, UINT128 z
3584             _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3585             _EXC_INFO_PARAM) {
3586 #endif
3587   int is_midpoint_lt_even, is_midpoint_gt_even,
3588     is_inexact_lt_midpoint, is_inexact_gt_midpoint;
3589   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
3590 
3591 #if DECIMAL_CALL_BY_REFERENCE
3592   bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3593         	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
3594         	  &res, &x, &y, &z
3595         	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3596         	  _EXC_INFO_ARG);
3597 #else
3598   res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3599         		&is_inexact_lt_midpoint,
3600         		&is_inexact_gt_midpoint, x, y,
3601         		z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3602         		_EXC_INFO_ARG);
3603 #endif
3604   BID_RETURN (res);
3605 }
3606 
3607 
3608 #if DECIMAL_CALL_BY_REFERENCE
3609 void
3610 bid128ddd_fma (UINT128 * pres, UINT64 * px, UINT64 * py, UINT64 * pz
3611                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3612                _EXC_INFO_PARAM) {
3613   UINT64 x = *px, y = *py, z = *pz;
3614 #if !DECIMAL_GLOBAL_ROUNDING
3615   unsigned int rnd_mode = *prnd_mode;
3616 #endif
3617 #else
3618 UINT128
3619 bid128ddd_fma (UINT64 x, UINT64 y, UINT64 z
3620                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3621                _EXC_INFO_PARAM) {
3622 #endif
3623   int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
3624     is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
3625   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
3626   UINT128 x1, y1, z1;
3627 
3628 #if DECIMAL_CALL_BY_REFERENCE
3629   bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3630   bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3631   bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3632   bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3633         	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
3634         	  &res, &x1, &y1, &z1
3635         	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3636         	  _EXC_INFO_ARG);
3637 #else
3638   x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3639   y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3640   z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3641   res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3642         		&is_inexact_lt_midpoint,
3643         		&is_inexact_gt_midpoint, x1, y1,
3644         		z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3645         		_EXC_INFO_ARG);
3646 #endif
3647   BID_RETURN (res);
3648 }
3649 
3650 
3651 #if DECIMAL_CALL_BY_REFERENCE
3652 void
3653 bid128ddq_fma (UINT128 * pres, UINT64 * px, UINT64 * py, UINT128 * pz
3654                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3655                _EXC_INFO_PARAM) {
3656   UINT64 x = *px, y = *py;
3657   UINT128 z = *pz;
3658 #if !DECIMAL_GLOBAL_ROUNDING
3659   unsigned int rnd_mode = *prnd_mode;
3660 #endif
3661 #else
3662 UINT128
3663 bid128ddq_fma (UINT64 x, UINT64 y, UINT128 z
3664                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3665                _EXC_INFO_PARAM) {
3666 #endif
3667   int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
3668     is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
3669   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
3670   UINT128 x1, y1;
3671 
3672 #if DECIMAL_CALL_BY_REFERENCE
3673   bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3674   bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3675   bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3676         	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
3677         	  &res, &x1, &y1, &z
3678         	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3679         	  _EXC_INFO_ARG);
3680 #else
3681   x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3682   y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3683   res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3684         		&is_inexact_lt_midpoint,
3685         		&is_inexact_gt_midpoint, x1, y1,
3686         		z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3687         		_EXC_INFO_ARG);
3688 #endif
3689   BID_RETURN (res);
3690 }
3691 
3692 
3693 #if DECIMAL_CALL_BY_REFERENCE
3694 void
3695 bid128dqd_fma (UINT128 * pres, UINT64 * px, UINT128 * py, UINT64 * pz
3696                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3697                _EXC_INFO_PARAM) {
3698   UINT64 x = *px, z = *pz;
3699 #if !DECIMAL_GLOBAL_ROUNDING
3700   unsigned int rnd_mode = *prnd_mode;
3701 #endif
3702 #else
3703 UINT128
3704 bid128dqd_fma (UINT64 x, UINT128 y, UINT64 z
3705                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3706                _EXC_INFO_PARAM) {
3707 #endif
3708   int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
3709     is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
3710   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
3711   UINT128 x1, z1;
3712 
3713 #if DECIMAL_CALL_BY_REFERENCE
3714   bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3715   bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3716   bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3717         	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
3718         	  &res, &x1, py, &z1
3719         	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3720         	  _EXC_INFO_ARG);
3721 #else
3722   x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3723   z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3724   res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3725         		&is_inexact_lt_midpoint,
3726         		&is_inexact_gt_midpoint, x1, y,
3727         		z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3728         		_EXC_INFO_ARG);
3729 #endif
3730   BID_RETURN (res);
3731 }
3732 
3733 
3734 #if DECIMAL_CALL_BY_REFERENCE
3735 void
3736 bid128dqq_fma (UINT128 * pres, UINT64 * px, UINT128 * py, UINT128 * pz
3737                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3738                _EXC_INFO_PARAM) {
3739   UINT64 x = *px;
3740 #if !DECIMAL_GLOBAL_ROUNDING
3741   unsigned int rnd_mode = *prnd_mode;
3742 #endif
3743 #else
3744 UINT128
3745 bid128dqq_fma (UINT64 x, UINT128 y, UINT128 z
3746                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3747                _EXC_INFO_PARAM) {
3748 #endif
3749   int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
3750     is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
3751   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
3752   UINT128 x1;
3753 
3754 #if DECIMAL_CALL_BY_REFERENCE
3755   bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3756   bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3757         	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
3758         	  &res, &x1, py, pz
3759         	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3760         	  _EXC_INFO_ARG);
3761 #else
3762   x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3763   res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3764         		&is_inexact_lt_midpoint,
3765         		&is_inexact_gt_midpoint, x1, y,
3766         		z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3767         		_EXC_INFO_ARG);
3768 #endif
3769   BID_RETURN (res);
3770 }
3771 
3772 
3773 #if DECIMAL_CALL_BY_REFERENCE
3774 void
3775 bid128qdd_fma (UINT128 * pres, UINT128 * px, UINT64 * py, UINT64 * pz
3776                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3777                _EXC_INFO_PARAM) {
3778   UINT64 y = *py, z = *pz;
3779 #if !DECIMAL_GLOBAL_ROUNDING
3780   unsigned int rnd_mode = *prnd_mode;
3781 #endif
3782 #else
3783 UINT128
3784 bid128qdd_fma (UINT128 x, UINT64 y, UINT64 z
3785                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3786                _EXC_INFO_PARAM) {
3787 #endif
3788   int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
3789     is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
3790   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
3791   UINT128 y1, z1;
3792 
3793 #if DECIMAL_CALL_BY_REFERENCE
3794   bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3795   bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3796   bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3797         	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
3798         	  &res, px, &y1, &z1
3799         	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3800         	  _EXC_INFO_ARG);
3801 #else
3802   y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3803   z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3804   res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3805         		&is_inexact_lt_midpoint,
3806         		&is_inexact_gt_midpoint, x, y1,
3807         		z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3808         		_EXC_INFO_ARG);
3809 #endif
3810   BID_RETURN (res);
3811 }
3812 
3813 
3814 #if DECIMAL_CALL_BY_REFERENCE
3815 void
3816 bid128qdq_fma (UINT128 * pres, UINT128 * px, UINT64 * py, UINT128 * pz
3817                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3818                _EXC_INFO_PARAM) {
3819   UINT64 y = *py;
3820 #if !DECIMAL_GLOBAL_ROUNDING
3821   unsigned int rnd_mode = *prnd_mode;
3822 #endif
3823 #else
3824 UINT128
3825 bid128qdq_fma (UINT128 x, UINT64 y, UINT128 z
3826                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3827                _EXC_INFO_PARAM) {
3828 #endif
3829   int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
3830     is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
3831   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
3832   UINT128 y1;
3833 
3834 #if DECIMAL_CALL_BY_REFERENCE
3835   bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3836   bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3837         	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
3838         	  &res, px, &y1, pz
3839         	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3840         	  _EXC_INFO_ARG);
3841 #else
3842   y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3843   res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3844         		&is_inexact_lt_midpoint,
3845         		&is_inexact_gt_midpoint, x, y1,
3846         		z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3847         		_EXC_INFO_ARG);
3848 #endif
3849   BID_RETURN (res);
3850 }
3851 
3852 
3853 #if DECIMAL_CALL_BY_REFERENCE
3854 void
3855 bid128qqd_fma (UINT128 * pres, UINT128 * px, UINT128 * py, UINT64 * pz
3856                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3857                _EXC_INFO_PARAM) {
3858   UINT64 z = *pz;
3859 #if !DECIMAL_GLOBAL_ROUNDING
3860   unsigned int rnd_mode = *prnd_mode;
3861 #endif
3862 #else
3863 UINT128
3864 bid128qqd_fma (UINT128 x, UINT128 y, UINT64 z
3865                _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3866                _EXC_INFO_PARAM) {
3867 #endif
3868   int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
3869     is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
3870   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
3871   UINT128 z1;
3872 
3873 #if DECIMAL_CALL_BY_REFERENCE
3874   bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3875   bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3876         	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
3877         	  &res, px, py, &z1
3878         	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3879         	  _EXC_INFO_ARG);
3880 #else
3881   z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3882   res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
3883         		&is_inexact_lt_midpoint,
3884         		&is_inexact_gt_midpoint, x, y,
3885         		z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3886         		_EXC_INFO_ARG);
3887 #endif
3888   BID_RETURN (res);
3889 }
3890 
3891 // Note: bid128qqq_fma is represented by bid128_fma
3892 
3893 // Note: bid64ddd_fma is represented by bid64_fma
3894 
3895 #if DECIMAL_CALL_BY_REFERENCE
3896 void
3897 bid64ddq_fma (UINT64 * pres, UINT64 * px, UINT64 * py, UINT128 * pz
3898               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3899               _EXC_INFO_PARAM) {
3900   UINT64 x = *px, y = *py;
3901 #if !DECIMAL_GLOBAL_ROUNDING
3902   unsigned int rnd_mode = *prnd_mode;
3903 #endif
3904 #else
3905 UINT64
3906 bid64ddq_fma (UINT64 x, UINT64 y, UINT128 z
3907               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3908               _EXC_INFO_PARAM) {
3909 #endif
3910   UINT64 res1 = 0xbaddbaddbaddbaddull;
3911   UINT128 x1, y1;
3912 
3913 #if DECIMAL_CALL_BY_REFERENCE
3914   bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3915   bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3916   bid64qqq_fma (&res1, &x1, &y1, pz
3917         	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3918         	_EXC_INFO_ARG);
3919 #else
3920   x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3921   y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3922   res1 = bid64qqq_fma (x1, y1, z
3923         	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3924         	       _EXC_INFO_ARG);
3925 #endif
3926   BID_RETURN (res1);
3927 }
3928 
3929 
3930 #if DECIMAL_CALL_BY_REFERENCE
3931 void
3932 bid64dqd_fma (UINT64 * pres, UINT64 * px, UINT128 * py, UINT64 * pz
3933               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3934               _EXC_INFO_PARAM) {
3935   UINT64 x = *px, z = *pz;
3936 #if !DECIMAL_GLOBAL_ROUNDING
3937   unsigned int rnd_mode = *prnd_mode;
3938 #endif
3939 #else
3940 UINT64
3941 bid64dqd_fma (UINT64 x, UINT128 y, UINT64 z
3942               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3943               _EXC_INFO_PARAM) {
3944 #endif
3945   UINT64 res1 = 0xbaddbaddbaddbaddull;
3946   UINT128 x1, z1;
3947 
3948 #if DECIMAL_CALL_BY_REFERENCE
3949   bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3950   bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3951   bid64qqq_fma (&res1, &x1, py, &z1
3952         	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3953         	_EXC_INFO_ARG);
3954 #else
3955   x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3956   z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3957   res1 = bid64qqq_fma (x1, y, z1
3958         	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3959         	       _EXC_INFO_ARG);
3960 #endif
3961   BID_RETURN (res1);
3962 }
3963 
3964 
3965 #if DECIMAL_CALL_BY_REFERENCE
3966 void
3967 bid64dqq_fma (UINT64 * pres, UINT64 * px, UINT128 * py, UINT128 * pz
3968               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3969               _EXC_INFO_PARAM) {
3970   UINT64 x = *px;
3971 #if !DECIMAL_GLOBAL_ROUNDING
3972   unsigned int rnd_mode = *prnd_mode;
3973 #endif
3974 #else
3975 UINT64
3976 bid64dqq_fma (UINT64 x, UINT128 y, UINT128 z
3977               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
3978               _EXC_INFO_PARAM) {
3979 #endif
3980   UINT64 res1 = 0xbaddbaddbaddbaddull;
3981   UINT128 x1;
3982 
3983 #if DECIMAL_CALL_BY_REFERENCE
3984   bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3985   bid64qqq_fma (&res1, &x1, py, pz
3986         	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3987         	_EXC_INFO_ARG);
3988 #else
3989   x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
3990   res1 = bid64qqq_fma (x1, y, z
3991         	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
3992         	       _EXC_INFO_ARG);
3993 #endif
3994   BID_RETURN (res1);
3995 }
3996 
3997 
3998 #if DECIMAL_CALL_BY_REFERENCE
3999 void
4000 bid64qdd_fma (UINT64 * pres, UINT128 * px, UINT64 * py, UINT64 * pz
4001               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
4002               _EXC_INFO_PARAM) {
4003   UINT64 y = *py, z = *pz;
4004 #if !DECIMAL_GLOBAL_ROUNDING
4005   unsigned int rnd_mode = *prnd_mode;
4006 #endif
4007 #else
4008 UINT64
4009 bid64qdd_fma (UINT128 x, UINT64 y, UINT64 z
4010               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
4011               _EXC_INFO_PARAM) {
4012 #endif
4013   UINT64 res1 = 0xbaddbaddbaddbaddull;
4014   UINT128 y1, z1;
4015 
4016 #if DECIMAL_CALL_BY_REFERENCE
4017   bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
4018   bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
4019   bid64qqq_fma (&res1, px, &y1, &z1
4020         	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
4021         	_EXC_INFO_ARG);
4022 #else
4023   y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
4024   z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
4025   res1 = bid64qqq_fma (x, y1, z1
4026         	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
4027         	       _EXC_INFO_ARG);
4028 #endif
4029   BID_RETURN (res1);
4030 }
4031 
4032 
4033 #if DECIMAL_CALL_BY_REFERENCE
4034 void
4035 bid64qdq_fma (UINT64 * pres, UINT128 * px, UINT64 * py, UINT128 * pz
4036               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
4037               _EXC_INFO_PARAM) {
4038   UINT64 y = *py;
4039 #if !DECIMAL_GLOBAL_ROUNDING
4040   unsigned int rnd_mode = *prnd_mode;
4041 #endif
4042 #else
4043 UINT64
4044 bid64qdq_fma (UINT128 x, UINT64 y, UINT128 z
4045               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
4046               _EXC_INFO_PARAM) {
4047 #endif
4048   UINT64 res1 = 0xbaddbaddbaddbaddull;
4049   UINT128 y1;
4050 
4051 #if DECIMAL_CALL_BY_REFERENCE
4052   bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
4053   bid64qqq_fma (&res1, px, &y1, pz
4054         	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
4055         	_EXC_INFO_ARG);
4056 #else
4057   y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
4058   res1 = bid64qqq_fma (x, y1, z
4059         	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
4060         	       _EXC_INFO_ARG);
4061 #endif
4062   BID_RETURN (res1);
4063 }
4064 
4065 
4066 #if DECIMAL_CALL_BY_REFERENCE
4067 void
4068 bid64qqd_fma (UINT64 * pres, UINT128 * px, UINT128 * py, UINT64 * pz
4069               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
4070               _EXC_INFO_PARAM) {
4071   UINT64 z = *pz;
4072 #if !DECIMAL_GLOBAL_ROUNDING
4073   unsigned int rnd_mode = *prnd_mode;
4074 #endif
4075 #else
4076 UINT64
4077 bid64qqd_fma (UINT128 x, UINT128 y, UINT64 z
4078               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
4079               _EXC_INFO_PARAM) {
4080 #endif
4081   UINT64 res1 = 0xbaddbaddbaddbaddull;
4082   UINT128 z1;
4083 
4084 #if DECIMAL_CALL_BY_REFERENCE
4085   bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
4086   bid64qqq_fma (&res1, px, py, &z1
4087         	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
4088         	_EXC_INFO_ARG);
4089 #else
4090   z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
4091   res1 = bid64qqq_fma (x, y, z1
4092         	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
4093         	       _EXC_INFO_ARG);
4094 #endif
4095   BID_RETURN (res1);
4096 }
4097 
4098 
4099 #if DECIMAL_CALL_BY_REFERENCE
4100 void
4101 bid64qqq_fma (UINT64 * pres, UINT128 * px, UINT128 * py, UINT128 * pz
4102               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
4103               _EXC_INFO_PARAM) {
4104   UINT128 x = *px, y = *py, z = *pz;
4105 #if !DECIMAL_GLOBAL_ROUNDING
4106   unsigned int rnd_mode = *prnd_mode;
4107 #endif
4108 #else
4109 UINT64
4110 bid64qqq_fma (UINT128 x, UINT128 y, UINT128 z
4111               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
4112               _EXC_INFO_PARAM) {
4113 #endif
4114   int is_midpoint_lt_even0 = 0, is_midpoint_gt_even0 = 0,
4115     is_inexact_lt_midpoint0 = 0, is_inexact_gt_midpoint0 = 0;
4116   int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
4117     is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
4118   int incr_exp;
4119   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
4120   UINT128 res128 = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
4121   UINT64 res1 = 0xbaddbaddbaddbaddull;
4122   unsigned int save_fpsf; // needed because of the call to bid128_ext_fma
4123   UINT64 sign;
4124   UINT64 exp;
4125   int unbexp;
4126   UINT128 C;
4127   BID_UI64DOUBLE tmp;
4128   int nr_bits;
4129   int q, x0;
4130   int scale;
4131   int lt_half_ulp = 0, eq_half_ulp = 0;
4132 
4133   // Note: for rounding modes other than RN or RA, the result can be obtained
4134   // by rounding first to BID128 and then to BID64
4135 
4136   save_fpsf = *pfpsf; // sticky bits - caller value must be preserved
4137   *pfpsf = 0;
4138 
4139 #if DECIMAL_CALL_BY_REFERENCE
4140   bid128_ext_fma (&is_midpoint_lt_even0, &is_midpoint_gt_even0,
4141         	  &is_inexact_lt_midpoint0, &is_inexact_gt_midpoint0,
4142         	  &res, &x, &y, &z
4143         	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
4144         	  _EXC_INFO_ARG);
4145 #else
4146   res = bid128_ext_fma (&is_midpoint_lt_even0, &is_midpoint_gt_even0,
4147         		&is_inexact_lt_midpoint0,
4148         		&is_inexact_gt_midpoint0, x, y,
4149         		z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
4150         		_EXC_INFO_ARG);
4151 #endif
4152 
4153   if ((rnd_mode == ROUNDING_DOWN) || (rnd_mode == ROUNDING_UP) ||
4154       (rnd_mode == ROUNDING_TO_ZERO) || // no double rounding error is possible
4155       ((res.w[HIGH_128W] & MASK_NAN) == MASK_NAN) || //res=QNaN (cannot be SNaN)
4156       ((res.w[HIGH_128W] & MASK_ANY_INF) == MASK_INF)) { // result is infinity
4157 #if DECIMAL_CALL_BY_REFERENCE
4158     bid128_to_bid64 (&res1, &res _RND_MODE_ARG _EXC_FLAGS_ARG);
4159 #else
4160     res1 = bid128_to_bid64 (res _RND_MODE_ARG _EXC_FLAGS_ARG);
4161 #endif
4162     // determine the unbiased exponent of the result
4163     unbexp = ((res1 >> 53) & 0x3ff) - 398;
4164 
4165     // if subnormal, res1  must have exp = -398
4166     // if tiny and inexact set underflow and inexact status flags
4167     if (!((res1 & MASK_NAN) == MASK_NAN) &&	// res1 not NaN
4168         (unbexp == -398)
4169         && ((res1 & MASK_BINARY_SIG1) < 1000000000000000ull)
4170         && (is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0
4171             || is_midpoint_lt_even0 || is_midpoint_gt_even0)) {
4172       // set the inexact flag and the underflow flag
4173       *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION);
4174     } else if (is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0 ||
4175                is_midpoint_lt_even0 || is_midpoint_gt_even0) {
4176       // set the inexact flag and the underflow flag
4177       *pfpsf |= INEXACT_EXCEPTION;
4178     }
4179 
4180     *pfpsf |= save_fpsf;
4181     BID_RETURN (res1);
4182   } // else continue, and use rounding to nearest to round to 16 digits
4183 
4184   // at this point the result is rounded to nearest (even or away) to 34 digits
4185   // (or less if exact), and it is zero or finite non-zero canonical [sub]normal
4186   sign = res.w[HIGH_128W] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
4187   exp = res.w[HIGH_128W] & MASK_EXP; // biased and shifted left 49 bits
4188   unbexp = (exp >> 49) - 6176;
4189   C.w[1] = res.w[HIGH_128W] & MASK_COEFF;
4190   C.w[0] = res.w[LOW_128W];
4191 
4192   if ((C.w[1] == 0x0 && C.w[0] == 0x0) ||	// result is zero
4193       (unbexp <= (-398 - 35)) || (unbexp >= (369 + 16))) {
4194       // clear under/overflow
4195 #if DECIMAL_CALL_BY_REFERENCE
4196     bid128_to_bid64 (&res1, &res _RND_MODE_ARG _EXC_FLAGS_ARG);
4197 #else
4198     res1 = bid128_to_bid64 (res _RND_MODE_ARG _EXC_FLAGS_ARG);
4199 #endif
4200     *pfpsf |= save_fpsf;
4201     BID_RETURN (res1);
4202   } // else continue
4203 
4204   // -398 - 34 <= unbexp <= 369 + 15
4205   if (rnd_mode == ROUNDING_TIES_AWAY) {
4206     // apply correction, if needed, to make the result rounded to nearest-even
4207     if (is_midpoint_gt_even) {
4208       // res = res - 1
4209       res1--; // res1 is now even
4210     } // else the result is already correctly rounded to nearest-even
4211   }
4212   // at this point the result is finite, non-zero canonical normal or subnormal,
4213   // and in most cases overflow or underflow will not occur
4214 
4215   // determine the number of digits q in the result
4216   // q = nr. of decimal digits in x
4217   // determine first the nr. of bits in x
4218   if (C.w[1] == 0) {
4219     if (C.w[0] >= 0x0020000000000000ull) { // x >= 2^53
4220       // split the 64-bit value in two 32-bit halves to avoid rounding errors
4221       if (C.w[0] >= 0x0000000100000000ull) { // x >= 2^32
4222         tmp.d = (double) (C.w[0] >> 32); // exact conversion
4223         nr_bits =
4224           33 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
4225       } else { // x < 2^32
4226         tmp.d = (double) (C.w[0]); // exact conversion
4227         nr_bits =
4228           1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
4229       }
4230     } else { // if x < 2^53
4231       tmp.d = (double) C.w[0]; // exact conversion
4232       nr_bits =
4233         1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
4234     }
4235   } else { // C.w[1] != 0 => nr. bits = 64 + nr_bits (C.w[1])
4236     tmp.d = (double) C.w[1]; // exact conversion
4237     nr_bits =
4238       65 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
4239   }
4240   q = nr_digits[nr_bits - 1].digits;
4241   if (q == 0) {
4242     q = nr_digits[nr_bits - 1].digits1;
4243     if (C.w[1] > nr_digits[nr_bits - 1].threshold_hi ||
4244         (C.w[1] == nr_digits[nr_bits - 1].threshold_hi &&
4245          C.w[0] >= nr_digits[nr_bits - 1].threshold_lo))
4246       q++;
4247   }
4248   // if q > 16, round to nearest even to 16 digits (but for underflow it may
4249   // have to be truncated even more)
4250   if (q > 16) {
4251     x0 = q - 16;
4252     if (q <= 18) {
4253       round64_2_18 (q, x0, C.w[0], &res1, &incr_exp,
4254         	    &is_midpoint_lt_even, &is_midpoint_gt_even,
4255         	    &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
4256     } else { // 19 <= q <= 34
4257       round128_19_38 (q, x0, C, &res128, &incr_exp,
4258         	      &is_midpoint_lt_even, &is_midpoint_gt_even,
4259         	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
4260       res1 = res128.w[0]; // the result fits in 64 bits
4261     }
4262     unbexp = unbexp + x0;
4263     if (incr_exp)
4264       unbexp++;
4265     q = 16; // need to set in case denormalization is necessary
4266   } else {
4267     // the result does not require a second rounding (and it must have
4268     // been exact in the first rounding, since q <= 16)
4269     res1 = C.w[0];
4270   }
4271 
4272   // avoid a double rounding error
4273   if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
4274       is_midpoint_lt_even) { // double rounding error upward
4275     // res = res - 1
4276     res1--; // res1 becomes odd
4277     is_midpoint_lt_even = 0;
4278     is_inexact_lt_midpoint = 1;
4279     if (res1 == 0x00038d7ea4c67fffull) { // 10^15 - 1
4280       res1 = 0x002386f26fc0ffffull; // 10^16 - 1
4281       unbexp--;
4282     }
4283   } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
4284       is_midpoint_gt_even) { // double rounding error downward
4285     // res = res + 1
4286     res1++; // res1 becomes odd (so it cannot be 10^16)
4287     is_midpoint_gt_even = 0;
4288     is_inexact_gt_midpoint = 1;
4289   } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
4290              !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
4291     // if this second rounding was exact the result may still be
4292     // inexact because of the first rounding
4293     if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
4294       is_inexact_gt_midpoint = 1;
4295     }
4296     if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
4297       is_inexact_lt_midpoint = 1;
4298     }
4299   } else if (is_midpoint_gt_even &&
4300              (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
4301     // pulled up to a midpoint
4302     is_inexact_lt_midpoint = 1;
4303     is_inexact_gt_midpoint = 0;
4304     is_midpoint_lt_even = 0;
4305     is_midpoint_gt_even = 0;
4306   } else if (is_midpoint_lt_even &&
4307              (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
4308     // pulled down to a midpoint
4309     is_inexact_lt_midpoint = 0;
4310     is_inexact_gt_midpoint = 1;
4311     is_midpoint_lt_even = 0;
4312     is_midpoint_gt_even = 0;
4313   } else {
4314     ;
4315   }
4316   // this is the result rounded correctly to nearest even, with unbounded exp.
4317 
4318   // check for overflow
4319   if (q + unbexp > P16 + expmax16) {
4320     res1 = sign | 0x7800000000000000ull;
4321     *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
4322     *pfpsf |= save_fpsf;
4323     BID_RETURN (res1)
4324   } else if (unbexp > expmax16) { // q + unbexp <= P16 + expmax16
4325     // not overflow; the result must be exact, and we can multiply res1 by
4326     // 10^(unbexp - expmax16) and the product will fit in 16 decimal digits
4327     scale = unbexp - expmax16;
4328     res1 = res1 * ten2k64[scale]; // res1 * 10^scale
4329     unbexp = expmax16; // unbexp - scale
4330   } else {
4331     ; // continue
4332   }
4333 
4334   // check for underflow
4335   if (q + unbexp < P16 + expmin16) {
4336     if (unbexp < expmin16) {
4337       // we must truncate more of res
4338       x0 = expmin16 - unbexp; // x0 >= 1
4339       is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
4340       is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
4341       is_midpoint_lt_even0 = is_midpoint_lt_even;
4342       is_midpoint_gt_even0 = is_midpoint_gt_even;
4343       is_inexact_lt_midpoint = 0;
4344       is_inexact_gt_midpoint = 0;
4345       is_midpoint_lt_even = 0;
4346       is_midpoint_gt_even = 0;
4347       // the number of decimal digits in res1 is q
4348       if (x0 < q) { // 1 <= x0 <= q-1 => round res to q - x0 digits
4349         // 2 <= q <= 16, 1 <= x0 <= 15
4350         round64_2_18 (q, x0, res1, &res1, &incr_exp,
4351         	      &is_midpoint_lt_even, &is_midpoint_gt_even,
4352         	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
4353         if (incr_exp) {
4354           // res1 = 10^(q-x0), 1 <= q - x0 <= q - 1, 1 <= q - x0 <= 15
4355           res1 = ten2k64[q - x0];
4356         }
4357         unbexp = unbexp + x0; // expmin16
4358       } else if (x0 == q) {
4359         // the second rounding is for 0.d(0)d(1)...d(q-1) * 10^emin
4360         // determine relationship with 1/2 ulp
4361         // q <= 16
4362         if (res1 < midpoint64[q - 1]) { // < 1/2 ulp
4363           lt_half_ulp = 1;
4364           is_inexact_lt_midpoint = 1;
4365         } else if (res1 == midpoint64[q - 1]) { // = 1/2 ulp
4366           eq_half_ulp = 1;
4367           is_midpoint_gt_even = 1;
4368         } else { // > 1/2 ulp
4369           // gt_half_ulp = 1;
4370           is_inexact_gt_midpoint = 1;
4371         }
4372         if (lt_half_ulp || eq_half_ulp) {
4373           // res = +0.0 * 10^expmin16
4374           res1 = 0x0000000000000000ull;
4375         } else { // if (gt_half_ulp)
4376           // res = +1 * 10^expmin16
4377           res1 = 0x0000000000000001ull;
4378         }
4379         unbexp = expmin16;
4380       } else { // if (x0 > q)
4381         // the second rounding is for 0.0...d(0)d(1)...d(q-1) * 10^emin
4382         res1 = 0x0000000000000000ull;
4383         unbexp = expmin16;
4384         is_inexact_lt_midpoint = 1;
4385       }
4386       // avoid a double rounding error
4387       if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) &&
4388           is_midpoint_lt_even) { // double rounding error upward
4389         // res = res - 1
4390         res1--; // res1 becomes odd
4391         is_midpoint_lt_even = 0;
4392         is_inexact_lt_midpoint = 1;
4393       } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) &&
4394           is_midpoint_gt_even) { // double rounding error downward
4395         // res = res + 1
4396         res1++; // res1 becomes odd
4397         is_midpoint_gt_even = 0;
4398         is_inexact_gt_midpoint = 1;
4399       } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
4400         	 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
4401         // if this rounding was exact the result may still be
4402         // inexact because of the previous roundings
4403         if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
4404           is_inexact_gt_midpoint = 1;
4405         }
4406         if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
4407           is_inexact_lt_midpoint = 1;
4408         }
4409       } else if (is_midpoint_gt_even &&
4410         	 (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
4411         // pulled up to a midpoint
4412         is_inexact_lt_midpoint = 1;
4413         is_inexact_gt_midpoint = 0;
4414         is_midpoint_lt_even = 0;
4415         is_midpoint_gt_even = 0;
4416       } else if (is_midpoint_lt_even &&
4417         	 (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
4418         // pulled down to a midpoint
4419         is_inexact_lt_midpoint = 0;
4420         is_inexact_gt_midpoint = 1;
4421         is_midpoint_lt_even = 0;
4422         is_midpoint_gt_even = 0;
4423       } else {
4424         ;
4425       }
4426     }
4427     // else if unbexp >= emin then q < P (because q + unbexp < P16 + expmin16)
4428     // and the result is tiny and exact
4429 
4430     // check for inexact result
4431     if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
4432         is_midpoint_lt_even || is_midpoint_gt_even ||
4433         is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0 ||
4434         is_midpoint_lt_even0 || is_midpoint_gt_even0) {
4435       // set the inexact flag and the underflow flag
4436       *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION);
4437     }
4438   } else if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
4439              is_midpoint_lt_even || is_midpoint_gt_even) {
4440     *pfpsf |= INEXACT_EXCEPTION;
4441   }
4442   // this is the result rounded correctly to nearest, with bounded exponent
4443 
4444   if (rnd_mode == ROUNDING_TIES_AWAY && is_midpoint_gt_even) { // correction
4445     // res = res + 1
4446     res1++; // res1 is now odd
4447   } // else the result is already correct
4448 
4449   // assemble the result
4450   if (res1 < 0x0020000000000000ull) { // res < 2^53
4451     res1 = sign | ((UINT64) (unbexp + 398) << 53) | res1;
4452   } else { // res1 >= 2^53
4453     res1 = sign | MASK_STEERING_BITS |
4454       ((UINT64) (unbexp + 398) << 51) | (res1 & MASK_BINARY_SIG2);
4455   }
4456   *pfpsf |= save_fpsf;
4457   BID_RETURN (res1);
4458 }
4459