1 /* Copyright (C) 2007-2018 Free Software Foundation, Inc.
2
3 This file is part of GCC.
4
5 GCC is free software; you can redistribute it and/or modify it under
6 the terms of the GNU General Public License as published by the Free
7 Software Foundation; either version 3, or (at your option) any later
8 version.
9
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 for more details.
14
15 Under Section 7 of GPL version 3, you are granted additional
16 permissions described in the GCC Runtime Library Exception, version
17 3.1, as published by the Free Software Foundation.
18
19 You should have received a copy of the GNU General Public License and
20 a copy of the GCC Runtime Library Exception along with this program;
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22 <http://www.gnu.org/licenses/>. */
23
24 /*****************************************************************************
25 *
26 * 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