]> git.ipfire.org Git - thirdparty/gcc.git/blob - libgcc/config/libbid/bid128_mul.c
projects / thirdparty / gcc.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Merged with libbbid branch at revision 126349.
[thirdparty/gcc.git] / libgcc / config / libbid / bid128_mul.c
1 /* Copyright (C) 2007 Free Software Foundation, Inc.
2
3 This file is part of GCC.
4
5 GCC is free software; you can redistribute it and/or modify it under
6 the terms of the GNU General Public License as published by the Free
7 Software Foundation; either version 2, or (at your option) any later
8 version.
9
10 In addition to the permissions in the GNU General Public License, the
11 Free Software Foundation gives you unlimited permission to link the
12 compiled version of this file into combinations with other programs,
13 and to distribute those combinations without any restriction coming
14 from the use of this file. (The General Public License restrictions
15 do apply in other respects; for example, they cover modification of
16 the file, and distribution when not linked into a combine
17 executable.)
18
19 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
20 WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with GCC; see the file COPYING. If not, write to the Free
26 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
27 02110-1301, USA. */
28
29 #include "bid_internal.h"
30
31 #if DECIMAL_CALL_BY_REFERENCE
32 void
33 __bid128_mul (UINT128 * pres, UINT128 * px,
34 UINT128 *
35 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
36 _EXC_INFO_PARAM) {
37 UINT128 x = *px, y = *py;
38
39 #if !DECIMAL_GLOBAL_ROUNDING
40 unsigned int rnd_mode = *prnd_mode;
41
42 #endif
43 #else
44 UINT128
45 __bid128_mul (UINT128 x,
46 UINT128 y _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
47 _EXC_INFO_PARAM) {
48
49 #endif
50 UINT128 res;
51 UINT64 x_sign, y_sign, sign;
52 UINT64 x_exp, y_exp;
53
54 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
55 // Note: C2.w[1], C2.w[0] represent y_signif_hi, y_signif_lo (all are UINT64)
56 UINT64 tmp64, tmp64A;
57 BID_UI64DOUBLE tmp1, tmp2;
58 int x_nr_bits, y_nr_bits;
59 int q1, q2, ind, shift;
60 UINT128 C1, C2;
61 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
62 UINT384 fstar;
63 int q;
64 UINT128 P128, R128; // for underflow path
65 UINT192 P192, R192; // for underflow path
66 UINT256 C, P256, R256;
67 UINT384 P384;
68 UINT512 P512;
69 int incr_exp = 0; // for underflow path
70 int incr_exp1 = 0; // for underflow path
71 int tmp_fpa = 0; // if possible underflow and q>=34, use to undo the rounding
72 UINT64 C1_hi, C2_hi;
73 UINT64 C1_lo, C2_lo;
74 int is_inexact = 0;
75 int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0;
76 int is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
77 int is_midpoint_lt_even1 = 0, is_midpoint_gt_even1 = 0;
78 int is_inexact_lt_midpoint1 = 0, is_inexact_gt_midpoint1 = 0;
79 int is_overflow = 0;
80 int no_underflow = 0;
81
82 // unpack the arguments
83 // unpack x
84 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
85 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
86 C1_hi = x.w[1] & MASK_COEFF;
87 C1_lo = x.w[0];
88
89 // unpack y
90 y_sign = y.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
91 y_exp = y.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
92 C2_hi = y.w[1] & MASK_COEFF;
93 C2_lo = y.w[0];
94 sign = x_sign ^ y_sign;
95
96 // check for NaN or Infinity
97 if (((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL)
98 || ((y.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
99
100 // x is special or y is special
101 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
102 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
103 // set invalid flag
104 *pfpsf |= INVALID_EXCEPTION;
105
106 // return quiet (x)
107 res.w[1] = x.w[1] & 0xfdffffffffffffffull;
108 res.w[0] = x.w[0];
109 } else { // x is QNaN
110 if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
111 // set invalid flag
112 *pfpsf |= INVALID_EXCEPTION;
113 }
114 // return x
115 res.w[1] = x.w[1];
116 res.w[0] = x.w[0];
117 }
118 BID_RETURN (res);
119 } else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
120 if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
121 // set invalid flag
122 *pfpsf |= INVALID_EXCEPTION;
123
124 // return quiet (y)
125 res.w[1] = y.w[1] & 0xfdffffffffffffffull;
126 res.w[0] = y.w[0];
127 } else { // y is QNaN
128 // return y
129 res.w[1] = y.w[1];
130 res.w[0] = y.w[0];
131 }
132 BID_RETURN (res);
133 } else { // neither x nor y is NaN; at least one is infinity
134 if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x is infinity
135 if (((y.w[1] & MASK_ANY_INF) == MASK_INF) || (C2_hi != 0x0ull)
136 || (C2_lo != 0x0ull)) {
137
138 // y is infinity OR y is finite
139 // if same sign, return +inf otherwise return -inf
140 if (!sign) {
141 res.w[1] = 0x7800000000000000ull; // +inf
142 res.w[0] = 0x0000000000000000ull;
143 } else { // x and y are infinities of opposite signs
144 res.w[1] = 0xf800000000000000ull; // -inf
145 res.w[0] = 0x0000000000000000ull;
146 }
147 } else { // if y is 0
148 // set invalid flag
149 *pfpsf |= INVALID_EXCEPTION;
150
151 // return QNaN Indefinite
152 res.w[1] = 0x7c00000000000000ull;
153 res.w[0] = 0x0000000000000000ull;
154 }
155 } else { // x is not NaN or infinity, so y must be infinity
156 if ((C1_hi != 0x0ull) || (C1_lo != 0x0ull)) {
157
158 // x is finite
159 // if same sign, return +inf otherwise return -inf
160 if (!sign) {
161 res.w[1] = 0x7800000000000000ull; // +inf
162 res.w[0] = 0x0000000000000000ull;
163 } else { // y and x are of opposite signs
164 res.w[1] = 0xf800000000000000ull; // -inf
165 res.w[0] = 0x0000000000000000ull;
166 }
167 } else { // if x is 0
168 // set invalid flag
169 *pfpsf |= INVALID_EXCEPTION;
170
171 // return QNaN Indefinite
172 res.w[1] = 0x7c00000000000000ull;
173 res.w[0] = 0x0000000000000000ull;
174 }
175 }
176 BID_RETURN (res);
177 }
178 }
179 // test for non-canonical values:
180 // - values whose encoding begins with x00, x01, or x10 and whose
181 // coefficient is larger than 10^34 -1, or
182 // - values whose encoding begins with x1100, x1101, x1110 (if NaNs
183 // and infinitis were eliminated already this test is reduced to
184 // checking for x10x)
185
186 // test for non-canonical values of the argument x
187 if ((((C1_hi > 0x0001ed09bead87c0ull) ||
188 ((C1_hi == 0x0001ed09bead87c0ull) && (C1_lo > 0x378d8e63ffffffffull))) &&
189 ((x.w[1] & 0x6000000000000000ull) != 0x6000000000000000ull)) ||
190 ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
191 // check for the case where the exponent is shifted right by 2 bits!
192 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {
193 x_exp = (x.w[1] << 2) & MASK_EXP; // same position as for G[0]G[1] != 11
194 }
195 x.w[1] = x.w[1] & 0x8000000000000000ull; // preserve the sign bit
196 x.w[0] = 0;
197 C1_hi = 0;
198 C1_lo = 0;
199 }
200 // test for non-canonical values of the argument y
201 if ((((C2_hi > 0x0001ed09bead87c0ull)
202 || ((C2_hi == 0x0001ed09bead87c0ull)
203 && (C2_lo > 0x378d8e63ffffffffull)))
204 && ((y.w[1] & 0x6000000000000000ull) != 0x6000000000000000ull))
205 || ((y.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
206
207 // check for the case where the exponent is shifted right by 2 bits!
208 if ((y.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {
209 y_exp = (y.w[1] << 2) & MASK_EXP; // same position as for G[0]G[1] != 11
210 }
211 y.w[1] = y.w[1] & 0x8000000000000000ull; // preserve the sign bit
212 y.w[0] = 0;
213 C2_hi = 0;
214 C2_lo = 0;
215 }
216 if (((C1_hi == 0x0ull) && (C1_lo == 0x0ull)) || ((C2_hi == 0x0ull)
217 && (C2_lo == 0x0ull))) {
218
219 // x is 0 and y is not special OR y is 0 and x is not special
220 // if same sign, return +0 otherwise return -0
221 ind = (x_exp >> 49) + (y_exp >> 49) - 6176;
222 if (ind < 0)
223 ind = 0;
224 if (ind > 0x2fff)
225 ind = 0x2fff; // 6111 + 6176
226 if ((x.w[1] & MASK_SIGN) == (y.w[1] & MASK_SIGN)) {
227 res.w[1] = 0x0000000000000000ull | ((UINT64) ind << 49); // +0.0
228 res.w[0] = 0x0000000000000000ull;
229 } else { // opposite signs
230 res.w[1] = 0x8000000000000000ull | ((UINT64) ind << 49); // -0.0
231 res.w[0] = 0x0000000000000000ull;
232 }
233 BID_RETURN (res);
234 } else { // x and y are not special and are not zero
235 // unpack x
236 C1.w[1] = C1_hi;
237 C1.w[0] = C1_lo;
238
239 // q1 = nr. of decimal digits in x
240 // determine first the nr. of bits in x
241 if (C1.w[1] == 0) {
242 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
243 // split the 64-bit value in two 32-bit halves to avoid rounding errors
244 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
245 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
246 x_nr_bits =
247 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
248 } else { // x < 2^32
249 tmp1.d = (double) (C1.w[0]); // exact conversion
250 x_nr_bits =
251 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
252 }} else { // if x < 2^53
253 tmp1.d = (double) C1.w[0]; // exact conversion
254 x_nr_bits =
255 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
256 }} else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
257 tmp1.d = (double) C1.w[1]; // exact conversion
258 x_nr_bits =
259 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
260 } q1 = __bid_nr_digits[x_nr_bits - 1].digits;
261 if (q1 == 0) {
262 q1 = __bid_nr_digits[x_nr_bits - 1].digits1;
263 if (C1.w[1] > __bid_nr_digits[x_nr_bits - 1].threshold_hi
264 || (C1.w[1] == __bid_nr_digits[x_nr_bits - 1].threshold_hi
265 && C1.w[0] >= __bid_nr_digits[x_nr_bits - 1].threshold_lo))
266 q1++;
267 }
268 C2.w[1] = C2_hi;
269 C2.w[0] = C2_lo;
270 if (C2.w[1] == 0) {
271 if (C2.w[0] >= 0x0020000000000000ull) { // y >= 2^53
272 // split the 64-bit value in two 32-bit halves to avoid rounding errors
273 if (C2.w[0] >= 0x0000000100000000ull) { // y >= 2^32
274 tmp2.d = (double) (C2.w[0] >> 32); // exact conversion
275 y_nr_bits =
276 32 + ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff);
277 } else { // y < 2^32
278 tmp2.d = (double) C2.w[0]; // exact conversion
279 y_nr_bits =
280 ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff);
281 }} else { // if y < 2^53
282 tmp2.d = (double) C2.w[0]; // exact conversion
283 y_nr_bits =
284 ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff);
285 }} else { // C2.w[1] != 0 => nr. bits = 64 + nr_bits (C2.w[1])
286 tmp2.d = (double) C2.w[1]; // exact conversion
287 y_nr_bits =
288 64 + ((((unsigned int) (tmp2.ui64 >> 52)) & 0x7ff) - 0x3ff);
289 } q2 = __bid_nr_digits[y_nr_bits].digits;
290 if (q2 == 0) {
291 q2 = __bid_nr_digits[y_nr_bits].digits1;
292 if (C2.w[1] > __bid_nr_digits[y_nr_bits].threshold_hi
293 || (C2.w[1] == __bid_nr_digits[y_nr_bits].threshold_hi
294 && C2.w[0] >= __bid_nr_digits[y_nr_bits].threshold_lo))
295 q2++;
296 }
297 // the exact product has either q1 + q2 - 1 or q1 + q2 decimal digits
298 // where 2 <= q1 + q2 <= 68
299 // calculate C' = C1 * C2 and determine q
300 C.w[3] = C.w[2] = C.w[1] = C.w[0] = 0;
301 if (q1 + q2 <= 19) { // if 2 <= q1 + q2 <= 19, C' = C1 * C2 fits in 64 bits
302 C.w[0] = C1.w[0] * C2.w[0];
303
304 // if C' < 10^(q1+q2-1) then q = q1 + q2 - 1 else q = q1 + q2
305 if (C.w[0] < __bid_ten2k64[q1 + q2 - 1])
306 q = q1 + q2 - 1; // q in [1, 18]
307 else
308 q = q1 + q2; // q in [2, 19]
309 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
310 } else if (q1 + q2 == 20) { // C' = C1 * C2 fits in 64 or 128 bits
311 // q1 <= 19 and q2 <= 19 so both C1 and C2 fit in 64 bits
312 __mul_64x64_to_128MACH (C, C1.w[0], C2.w[0]);
313
314 // if C' < 10^(q1+q2-1) = 10^19 then q = q1+q2-1 = 19 else q = q1+q2 = 20
315 if (C.w[1] == 0 && C.w[0] < __bid_ten2k64[19]) { // 19 = q1+q2-1
316 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
317 q = 19; // 19 = q1 + q2 - 1
318 } else {
319
320 // if (C.w[1] == 0)
321 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
322 // else
323 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
324 q = 20; // 20 = q1 + q2
325 }
326 } else if (q1 + q2 <= 38) { // 21 <= q1 + q2 <= 38
327 // C' = C1 * C2 fits in 64 or 128 bits
328 // (64 bits possibly, but only when q1 + q2 = 21 and C' has 20 digits)
329 // at least one of C1, C2 has at most 19 decimal digits & fits in 64 bits
330 if (q1 <= 19) {
331 __mul_128x64_to_128 (C, C1.w[0], C2);
332 } else { // q2 <= 19
333 __mul_128x64_to_128 (C, C2.w[0], C1);
334 }
335
336 // if C' < 10^(q1+q2-1) then q = q1 + q2 - 1 else q = q1 + q2
337 if (C.w[1] < __bid_ten2k128[q1 + q2 - 21].w[1]
338 || (C.w[1] == __bid_ten2k128[q1 + q2 - 21].w[1]
339 && C.w[0] < __bid_ten2k128[q1 + q2 - 21].w[0])) {
340
341 // if (C.w[1] == 0) // q = 20, necessarily
342 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
343 // else
344 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
345 q = q1 + q2 - 1; // q in [20, 37]
346 } else {
347
348 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
349 q = q1 + q2; // q in [21, 38]
350 }
351 } else if (q1 + q2 == 39) { // C' = C1 * C2 fits in 128 or 192 bits
352 // both C1 and C2 fit in 128 bits (actually in 113 bits)
353 // may replace this by 128x128_to192
354 __mul_128x128_to_256 (C, C1, C2); // C.w[3] is 0
355 // if C' < 10^(q1+q2-1) = 10^38 then q = q1+q2-1 = 38 else q = q1+q2 = 39
356 if (C.w[2] == 0 && (C.w[1] < __bid_ten2k128[18].w[1] ||
357 (C.w[1] == __bid_ten2k128[18].w[1] && C.w[0] < __bid_ten2k128[18].w[0]))) {
358 // 18 = 38 - 20 = q1+q2-1 - 20
359 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
360 q = 38; // 38 = q1 + q2 - 1
361 } else {
362
363 // if (C.w[2] == 0)
364 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
365 // else
366 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
367 q = 39; // 39 = q1 + q2
368 }
369 } else if (q1 + q2 <= 57) { // 40 <= q1 + q2 <= 57
370 // C' = C1 * C2 fits in 128 or 192 bits
371 // (128 bits possibly, but only when q1 + q2 = 40 and C' has 39 digits)
372 // both C1 and C2 fit in 128 bits (actually in 113 bits); at most one
373 // may fit in 64 bits
374 if (C1.w[1] == 0) { // C1 fits in 64 bits
375 // __mul_64x128_full (REShi64, RESlo128, A64, B128)
376 __mul_64x128_full (C.w[2], C, C1.w[0], C2);
377 } else if (C2.w[1] == 0) { // C2 fits in 64 bits
378 // __mul_64x128_full (REShi64, RESlo128, A64, B128)
379 __mul_64x128_full (C.w[2], C, C2.w[0], C1);
380 } else { // both C1 and C2 require 128 bits
381 // may use __mul_128x128_to_192 (C.w[2], C.w[0], C2.w[0], C1);
382 __mul_128x128_to_256 (C, C1, C2); // C.w[3] = 0
383 }
384
385 // if C' < 10^(q1+q2-1) then q = q1 + q2 - 1 else q = q1 + q2
386 if (C.w[2] < __bid_ten2k256[q1 + q2 - 40].w[2]
387 || (C.w[2] == __bid_ten2k256[q1 + q2 - 40].w[2]
388 && (C.w[1] < __bid_ten2k256[q1 + q2 - 40].w[1]
389 || (C.w[1] == __bid_ten2k256[q1 + q2 - 40].w[1]
390 && C.w[0] < __bid_ten2k256[q1 + q2 - 40].w[0])))) {
391
392 // if (C.w[2] == 0) // q = 39, necessarily
393 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
394 // else
395 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
396 q = q1 + q2 - 1; // q in [39, 56]
397 } else {
398
399 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
400 q = q1 + q2; // q in [40, 57]
401 }
402 } else if (q1 + q2 == 58) { // C' = C1 * C2 fits in 192 or 256 bits
403 // both C1 and C2 fit in 128 bits (actually in 113 bits); at most one
404 // may fit in 64 bits
405 if (C1.w[1] == 0) { // C1 * C2 will fit in 192 bits
406 __mul_64x128_full (C.w[2], C, C1.w[0], C2); // may use 64x128_to_192
407 } else if (C2.w[1] == 0) { // C1 * C2 will fit in 192 bits
408 __mul_64x128_full (C.w[2], C, C2.w[0], C1); // may use 64x128_to_192
409 } else { // C1 * C2 will fit in 192 bits or in 256 bits
410 __mul_128x128_to_256 (C, C1, C2);
411 }
412
413 // if C' < 10^(q1+q2-1) = 10^57 then q = q1+q2-1 = 57 else q = q1+q2 = 58
414 if (C.w[3] == 0 && (C.w[2] < __bid_ten2k256[18].w[2] ||
415 (C.w[2] == __bid_ten2k256[18].w[2] && (C.w[1] < __bid_ten2k256[18].w[1] ||
416 (C.w[1] == __bid_ten2k256[18].w[1] && C.w[0] < __bid_ten2k256[18].w[0]))))) {
417 // 18 = 57 - 39 = q1+q2-1 - 39
418 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
419 q = 57; // 57 = q1 + q2 - 1
420 } else {
421
422 // if (C.w[3] == 0)
423 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
424 // else
425 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 256;
426 q = 58; // 58 = q1 + q2
427 }
428 } else { // if 59 <= q1 + q2 <= 68
429 // C' = C1 * C2 fits in 192 or 256 bits
430 // (192 bits possibly, but only when q1 + q2 = 59 and C' has 58 digits)
431 // both C1 and C2 fit in 128 bits (actually in 113 bits); none fits in
432 // 64 bits
433 // may use __mul_128x128_to_192 (C.w[2], C.w[0], C2.w[0], C1);
434 __mul_128x128_to_256 (C, C1, C2); // C.w[3] = 0
435 // if C' < 10^(q1+q2-1) then q = q1 + q2 - 1 else q = q1 + q2
436 if (C.w[3] < __bid_ten2k256[q1 + q2 - 40].w[3]
437 || (C.w[3] == __bid_ten2k256[q1 + q2 - 40].w[3]
438 && (C.w[2] < __bid_ten2k256[q1 + q2 - 40].w[2]
439 || (C.w[2] == __bid_ten2k256[q1 + q2 - 40].w[2]
440 && (C.w[1] < __bid_ten2k256[q1 + q2 - 40].w[1]
441 || (C.w[1] == __bid_ten2k256[q1 + q2 - 40].w[1]
442 && C.w[0] < __bid_ten2k256[q1 + q2 - 40].w[0])))))) {
443
444 // if (C.w[3] == 0) // q = 58, necessarily
445 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
446 // else
447 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 256;
448 q = q1 + q2 - 1; // q in [58, 67]
449 } else {
450
451 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 256;
452 q = q1 + q2; // q in [59, 68]
453 }
454 }
455 if (((UINT64) q << 49) + x_exp + y_exp <
456 ((UINT64) P34 << 49) + EXP_MIN + BIN_EXP_BIAS) {
457
458 // possible underflow
459 // q + ex + ey < P34 + EMIN <=> q - P34 < EMIN - ex - ey <=> q - P34 < ind
460 goto _underflow_path;
461 }
462 if (q <= 34) { // 2 <= q <= 34 the result is exact, and fits in 113 bits
463 tmp64 = x_exp + y_exp;
464 if (tmp64 > EXP_MAX + BIN_EXP_BIAS) { // possible overflow
465 ind = (tmp64 - EXP_MAX - BIN_EXP_BIAS) >> 49;
466 if (ind > 34 - q) { // overflow in all rounding modes
467 // |res| >= 10^p * 10^emax = 10^(p-1) * 10^(emax+1)
468 // assemble the result
469 if (rnd_mode == ROUNDING_TO_NEAREST
470 || rnd_mode == ROUNDING_TIES_AWAY) {
471 res.w[1] = sign | 0x7800000000000000ull;
472 res.w[0] = 0x0ull;
473 } else if (rnd_mode == ROUNDING_DOWN) {
474 if (sign) { // res = -inf
475 res.w[1] = 0xf800000000000000ull;
476 res.w[0] = 0x0ull;
477 } else { // res = +MAXFP
478 res.w[1] = 0x5fffed09bead87c0ull;
479 res.w[0] = 0x378d8e63ffffffffull;
480 }
481 } else if (rnd_mode == ROUNDING_UP) {
482 if (sign) { // res = -MAXFP
483 res.w[1] = 0xdfffed09bead87c0ull;
484 res.w[0] = 0x378d8e63ffffffffull;
485 } else { // res = +inf
486 res.w[1] = 0x7800000000000000ull;
487 res.w[0] = 0x0ull;
488 }
489 } else { // if (rnd_mode == ROUNDING_TO_ZERO)
490 // |res| = (10^34 - 1) * 10^6111 = +MAXFP
491 res.w[1] = sign | 0x5fffed09bead87c0ull;
492 res.w[0] = 0x378d8e63ffffffffull;
493 }
494
495 // set the inexact flag
496 *pfpsf |= INEXACT_EXCEPTION;
497
498 // set the overflow flag
499 *pfpsf |= OVERFLOW_EXCEPTION;
500
501 // is_overflow = 1;
502 BID_RETURN (res);
503 } else { // tmp64 > EXP_MAX + BIN_EXP_BIAS but
504 // ind = ((tmp64-EXP_MAX-BIN_EXP_BIAS)>>49) <= 34 - q
505 // the exponent will be the maximum exponent
506 // multiply C by 10^ind; the result fits in 34 digits
507 if (ind <= 19) { // multiply by __bid_ten2k64[ind]
508 if (q <= 19) { // 64x64 -> 128
509 __mul_64x64_to_128MACH (C, C.w[0], __bid_ten2k64[ind]);
510 } else { // 128 x 64 -> 128
511 // may optimize to multiply 64 x 128 -> 128
512 __mul_64x128_full (tmp64, C, __bid_ten2k64[ind], C);
513 }
514 } else { // if 20 <= ind <= 32 multiply by __bid_ten2k128[ind - 20]
515 // it must be that C.w[1] = 0, as C < 10^14
516 // may optimize to multiply 64 x 128 -> 128
517 __mul_64x128_full (tmp64, C, C.w[0], __bid_ten2k128[ind - 20]);
518 }
519 res.w[0] = C.w[0];
520 res.w[1] = C.w[1];
521 res.w[1] |= EXP_MAX; // EXP MAX
522 }
523 } else {
524 res.w[0] = C.w[0];
525 res.w[1] = C.w[1];
526 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
527 }
528 res.w[1] |= sign;
529 } else if (q <= 38) { // 35 <= q <= 38; exact coefficient fits in 128 bits
530 // C = C + 1/2 * 10^x where the result C fits in 127 bits
531 ind = q - 35;
532 tmp64 = C.w[0];
533 C.w[0] = C.w[0] + __bid_midpoint64[ind];
534 if (C.w[0] < tmp64)
535 C.w[1]++;
536
537 // x = q - p = q - 34, 1 <= x <= 4
538 // kx = 10^(-x) = __bid_ten2mk128M[ind]
539 // C* = (C + 1/2 * 10^x) * 10^(-x)
540 // the approximation of 10^(-x) was rounded up to 128 bits
541 __mul_128x128_to_256 (P256, C, __bid_ten2mk128M[ind]);
542 Cstar.w[1] = P256.w[3];
543 Cstar.w[0] = P256.w[2];
544 fstar.w[2] = Cstar.w[0] & __bid_maskhigh128M[ind]; // fstar.w[3|4|5]=0
545 fstar.w[1] = P256.w[1];
546 fstar.w[0] = P256.w[0];
547
548 // calculate C* and f*
549 // C* is actually floor(C*) in this case
550 // C* and f* need shifting and masking, as shown by
551 // __bid_shiftright128M[] and __bid_maskhigh128M[]
552 // the top Ex bits of 10^(-x) are T* = __bid_ten2mk128truncM[ind], e.g.
553 // if x=1, T*=__bid_ten2mk128truncM[0]=0xcccccccccccccccccccccccccccccccc
554 // if (0 < f* < 10^(-x)) then the result is a midpoint
555 // if floor(C*) is even then C* = floor(C*) - logical right
556 // shift; C* has p decimal digits, correct by Prop. 1)
557 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
558 // shift; C* has p decimal digits, correct by Pr. 1)
559 // else
560 // C* = floor(C*) (logical right shift; C has p decimal digits,
561 // correct by Property 1)
562 // n = C* * 10^(e+x)
563
564 // shift right C* by Ex-128 = __bid_shiftright128M[ind]
565 shift = __bid_shiftright128M[ind]; // 3 <= shift <= 13
566 Cstar.w[0] = (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
567 Cstar.w[1] = (Cstar.w[1] >> shift);
568
569 // determine inexactness of the rounding of C*
570 // if (0 < f* - 1/2 < 10^(-x)) then
571 // the result is exact
572 // else // if (f* - 1/2 > T*) then
573 // the result is inexact
574 if (fstar.w[2] > __bid_one_half128M[ind]
575 || (fstar.w[2] == __bid_one_half128M[ind]
576 && (fstar.w[1] || fstar.w[0]))) {
577
578 // f* > 1/2 and the result may be exact
579 // Calculate f* - 1/2
580 tmp64 = fstar.w[2] - __bid_one_half128M[ind];
581 if (tmp64 || fstar.w[1] > __bid_ten2mk128truncM[ind].w[1] ||
582 (fstar.w[1] == __bid_ten2mk128truncM[ind].w[1] &&
583 fstar.w[0] > __bid_ten2mk128truncM[ind].w[0])) { // f* - 1/2 > 10^(-x)
584 // set the inexact flag
585 *pfpsf |= INEXACT_EXCEPTION;
586 is_inexact_lt_midpoint = 1;
587 } // else the result is exact
588 } else { // the result is inexact; f2* <= 1/2
589 // set the inexact flag
590 *pfpsf |= INEXACT_EXCEPTION;
591 tmp_fpa = 1;
592 is_inexact_gt_midpoint = 1;
593 }
594
595 // check for midpoints (could do this before determining inexactness)
596 if ((fstar.w[2] == 0) && (fstar.w[1] || fstar.w[0])
597 && (fstar.w[1] < __bid_ten2mk128truncM[ind].w[1]
598 || (fstar.w[1] == __bid_ten2mk128truncM[ind].w[1]
599 && fstar.w[0] <= __bid_ten2mk128truncM[ind].w[0]))) {
600
601 // the result is a midpoint
602 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
603 // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
604 Cstar.w[0]--; // Cstar.w[0] is now even
605 if (tmp_fpa == 1)
606 tmp_fpa = 0;
607 is_midpoint_gt_even = 1;
608 is_inexact_lt_midpoint = 0;
609 is_inexact_gt_midpoint = 0;
610 } else { // else MP in [ODD, EVEN]
611 is_midpoint_lt_even = 1;
612 is_inexact_lt_midpoint = 0;
613 is_inexact_gt_midpoint = 0;
614 }
615 }
616 // check for rounding overflow
617 if (Cstar.w[1] == 0x0001ed09bead87c0ull &&
618 Cstar.w[0] == 0x378d8e6400000000ull) { // if Cstar = 10^34
619 tmp64 = x_exp + y_exp + ((UINT64) (ind + 2) << 49);
620 Cstar.w[1] = 0x0000314dc6448d93ull; // Cstar = 10^33
621 Cstar.w[0] = 0x38c15b0a00000000ull;
622
623 // if rounding overflow made the exponent equal to emin, set underflow
624 if (tmp64 == EXP_MIN + BIN_EXP_BIAS)
625 *pfpsf |= UNDERFLOW_EXCEPTION;
626 } else { // 10^33 <= Cstar <= 10^34 - 1
627 tmp64 = x_exp + y_exp + ((UINT64) (ind + 1) << 49); // ind+1 = q-34
628 }
629 if (tmp64 >= EXP_MAX + BIN_EXP_BIAS) { // possibble overflow
630 // exp >= emax for the result rounded to nearest even
631 if (rnd_mode == ROUNDING_TO_NEAREST
632 || rnd_mode == ROUNDING_TIES_AWAY) {
633 if (tmp64 > EXP_MAX + BIN_EXP_BIAS) {
634
635 // |res| >= 10^(p-1) * 10^(emax+1) <=> exp >= emax+1
636 res.w[1] = sign | 0x7800000000000000ull; // +/-inf
637 res.w[0] = 0x0ull;
638 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
639 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
640 is_overflow = 1;
641 } else { // not overflow
642 res.w[0] = Cstar.w[0];
643 res.w[1] = Cstar.w[1];
644 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
645 }
646 } else if (rnd_mode == ROUNDING_DOWN) {
647 if (!sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
648 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
649 Cstar.w[1] == 0x0000314dc6448d93ull &&
650 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
651 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
652
653 // res = +MAXFP
654 res.w[1] = 0x5fffed09bead87c0ull;
655 res.w[0] = 0x378d8e63ffffffffull; // (10^34-1) * 10^emax
656 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
657 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
658 is_overflow = 1;
659 } else if (sign && ((tmp64 > EXP_MAX + BIN_EXP_BIAS) ||
660 ((tmp64 == EXP_MAX + BIN_EXP_BIAS) &&
661 Cstar.w[1] == 0x0001ed09bead87c0ull &&
662 Cstar.w[0] == 0x378d8e63ffffffffull && // (10^34-1) * 10^emax
663 is_inexact_lt_midpoint))) {
664 res.w[1] = 0xf800000000000000ull; // -inf
665 res.w[0] = 0x0ull;
666 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
667 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
668 is_overflow = 1;
669 } else { // not overflow
670 res.w[0] = Cstar.w[0];
671 res.w[1] = Cstar.w[1];
672 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
673 }
674 } else if (rnd_mode == ROUNDING_UP) {
675 if (sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
676 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
677 Cstar.w[1] == 0x0000314dc6448d93ull &&
678 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
679 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
680 // res = -MAXFP
681 res.w[1] = 0xdfffed09bead87c0ull;
682 res.w[0] = 0x378d8e63ffffffffull; // -(10^34-1) * 10^emax
683 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
684 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
685 is_overflow = 1;
686 } else if (!sign && ((tmp64 > EXP_MAX + BIN_EXP_BIAS) ||
687 ((tmp64 == EXP_MAX + BIN_EXP_BIAS) &&
688 Cstar.w[1] == 0x0001ed09bead87c0ull &&
689 Cstar.w[0] == 0x378d8e63ffffffffull && // (10^34-1) * 10^emax
690 is_inexact_lt_midpoint))) {
691 res.w[1] = 0x7800000000000000ull; // inf
692 res.w[0] = 0x0ull;
693 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
694 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
695 is_overflow = 1;
696 } else { // not overflow
697 res.w[0] = Cstar.w[0];
698 res.w[1] = Cstar.w[1];
699 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
700 }
701 } else { // if (rnd_mode == ROUNDING_TO_ZERO)
702 if (!sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
703 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
704 Cstar.w[1] == 0x0000314dc6448d93ull &&
705 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
706 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
707 // res = +MAXFP
708 res.w[1] = 0x5fffed09bead87c0ull;
709 res.w[0] = 0x378d8e63ffffffffull; // (10^34-1) * 10^emax
710 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
711 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
712 is_overflow = 1;
713 } else if (sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
714 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
715 Cstar.w[1] == 0x0000314dc6448d93ull &&
716 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
717 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
718 // res = -MAXFP
719 res.w[1] = 0xdfffed09bead87c0ull;
720 res.w[0] = 0x378d8e63ffffffffull; // -(10^34-1) * 10^emax
721 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
722 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
723 is_overflow = 1;
724 } else { // not overflow
725 res.w[0] = Cstar.w[0];
726 res.w[1] = Cstar.w[1];
727 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
728 }
729 }
730 if (is_overflow) { // return for overflow
731 // set the inexact flag
732 *pfpsf |= INEXACT_EXCEPTION;
733
734 // set the overflow flag
735 *pfpsf |= OVERFLOW_EXCEPTION;
736
737 // is_overflow = 1;
738 BID_RETURN (res);
739 }
740 } else {
741 res.w[0] = Cstar.w[0];
742 res.w[1] = Cstar.w[1];
743 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
744 }
745 res.w[1] |= sign;
746 } else if (q <= 57) { // 39 <= q <= 57; exact coefficient takes 128-192 bits
747 // C = C + 1/2 * 10^x where the result C fits in 190 bits
748 // (10^57 - 1 + 1/2 * 10^23 can be represented with 190 bits)
749 // x = q - p = q - 34, 5 <= x <= 23
750 // kx = 10^(-x) = __bid_ten2mk192M[ind]
751 // C* = (C + 1/2 * 10^x) * 10^(-x)
752 // the approximation of 10^(-x) was rounded up to 192 bits
753 ind = q - 39; // 0 <= ind <= 18
754 tmp64 = C.w[0];
755 tmp64A = C.w[1];
756
757 // Note:
758 // if 5 <= x <= 19 <=> 0 <= ind <= 14 then
759 // f* has 256 bits
760 // else // if 20 <= x <= 23 <=> 15 <= ind <= 18 then
761 // f* has 320 bits
762 if (ind <= 14) { // x - 1 = q - 35 = ind + 4 <= 18
763 // add one 64-bit word
764 C.w[0] = C.w[0] + __bid_midpoint64[ind + 4];
765 if (C.w[0] < tmp64)
766 C.w[1]++;
767 if (C.w[1] < tmp64A)
768 C.w[2]++;
769 __mul_192x192_to_384 (P384, C, __bid_ten2mk192M[ind])
770 // calculate C* and f*; C* is actually floor(C*) in this case
771 // C* and f* need shifting and masking, as shown by
772 // __bid_shiftright192M[] and __bid_maskhigh192M[]
773 // C* has 128 bits; P384.w[5], P384.w[4], P384.w[3] need to be
774 // shifted right by Ex-192 = __bid_shiftright192M[ind]
775 shift = __bid_shiftright192M[ind]; // 16 <= shift <= 63
776 Cstar.w[0] = (P384.w[3] >> shift) | (P384.w[4] << (64 - shift));
777 Cstar.w[1] = (P384.w[4] >> shift) | (P384.w[5] << (64 - shift));
778
779 // f* has 256 bits
780 fstar.w[3] = P384.w[3] & __bid_maskhigh192M[ind];
781 fstar.w[2] = P384.w[2];
782 fstar.w[1] = P384.w[1];
783 fstar.w[0] = P384.w[0];
784
785 // the top Ex bits of 10^(-x) are T* = __bid_ten2mk192truncM[ind], e.g.
786 // if x=5, T* = __bid_ten2mk192truncM[0] =
787 // 0xa7c5ac471b4784230fcf80dc33721d53cddd6e04c0592103
788 // if (0 < f* < 10^(-x)) then the result is a midpoint
789 // if floor(C*) is even then C* = floor(C*) - logical right
790 // shift; C* has p decimal digits, correct by Prop. 1)
791 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
792 // shift; C* has p decimal digits, correct by Pr. 1)
793 // else
794 // C* = floor(C*) (logical right shift; C has p decimal digits,
795 // correct by Property 1)
796 // n = C* * 10^(e+x)
797
798 // determine inexactness of the rounding of C*
799 // if (0 < f* - 1/2 < T* ~= 10^(-x)) then
800 // the result is exact
801 // else // if (f* - 1/2 >= T*) then
802 // the result is inexact
803 if (fstar.w[3] > __bid_one_half192M[ind]
804 || (fstar.w[3] == __bid_one_half192M[ind]
805 && (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
806
807 // f* > 1/2 and the result may be exact
808 // Calculate f* - 1/2
809 tmp64 = fstar.w[3] - __bid_one_half192M[ind];
810 if (tmp64 || fstar.w[2] > __bid_ten2mk192truncM[ind].w[2] ||
811 (fstar.w[2] == __bid_ten2mk192truncM[ind].w[2] &&
812 fstar.w[1] > __bid_ten2mk192truncM[ind].w[1]) ||
813 (fstar.w[2] == __bid_ten2mk192truncM[ind].w[2] &&
814 fstar.w[1] == __bid_ten2mk192truncM[ind].w[1] &&
815 fstar.w[0] > __bid_ten2mk192truncM[ind].w[0])) { // f* - 1/2 > 10^(-x)
816 // set the inexact flag
817 *pfpsf |= INEXACT_EXCEPTION;
818 is_inexact_lt_midpoint = 1;
819 } // else the result is exact
820 } else { // the result is inexact; f2* <= 1/2
821 // set the inexact flag
822 *pfpsf |= INEXACT_EXCEPTION;
823 tmp_fpa = 1;
824 is_inexact_gt_midpoint = 1;
825 }
826
827 // check for midpoints (could do this before determining inexactness)
828 if ((fstar.w[3] == 0)
829 && (fstar.w[2] || fstar.w[1] || fstar.w[0])
830 && (fstar.w[2] < __bid_ten2mk192truncM[ind].w[2]
831 || (fstar.w[2] == __bid_ten2mk192truncM[ind].w[2]
832 && fstar.w[1] < __bid_ten2mk192truncM[ind].w[1])
833 || (fstar.w[2] == __bid_ten2mk192truncM[ind].w[2]
834 && fstar.w[1] == __bid_ten2mk192truncM[ind].w[1]
835 && fstar.w[0] <= __bid_ten2mk192truncM[ind].w[0]))) {
836
837 // the result is a midpoint
838 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
839 // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
840 Cstar.w[0]--; // Cstar.w[0] is now even
841 if (tmp_fpa == 1)
842 tmp_fpa = 0;
843 is_midpoint_gt_even = 1;
844 is_inexact_lt_midpoint = 0;
845 is_inexact_gt_midpoint = 0;
846 } else { // else MP in [ODD, EVEN]
847 is_midpoint_lt_even = 1;
848 is_inexact_lt_midpoint = 0;
849 is_inexact_gt_midpoint = 0;
850 }
851 }
852 } else { // if ind >= 15 <=> x - 1 = q - 35 = ind + 4 >= 19
853 // add two 64-bit words
854 C.w[0] = C.w[0] + __bid_midpoint128[ind - 15].w[0];
855 C.w[1] = C.w[1] + __bid_midpoint128[ind - 15].w[1];
856 if (C.w[0] < tmp64)
857 C.w[1]++;
858 if (C.w[1] < tmp64A)
859 C.w[2]++;
860 __mul_192x192_to_384 (P384, C, __bid_ten2mk192M[ind])
861 // calculate C* and f*; C* is actually floor(C*) in this case
862 // C* and f* need shifting and masking, as shown by
863 // __bid_shiftright192M[] and __bid_maskhigh192M[]
864 // C* has 128 bits; P384.w[5], P384.w[4], need to be
865 // shifted right by Ex-256 = __bid_shiftright192M[ind]
866 shift = __bid_shiftright192M[ind]; // 2 <= shift <= 12
867 Cstar.w[0] = (P384.w[4] >> shift) | (P384.w[5] << (64 - shift));
868 Cstar.w[1] = (P384.w[5] >> shift);
869
870 // f* has 320 bits
871 fstar.w[4] = P384.w[4] & __bid_maskhigh192M[ind];
872 fstar.w[3] = P384.w[3];
873 fstar.w[2] = P384.w[2];
874 fstar.w[1] = P384.w[1];
875 fstar.w[0] = P384.w[0];
876
877 // the top Ex bits of 10^(-x) are T* = __bid_ten2mk192truncM[ind], e.g.
878 // if x=23, T* = __bid_ten2mk192truncM[18] =
879 // 0xc16d9a0095928a2775b7053c0f1782938d6f439b43088650
880 // if (0 < f* < 10^(-x)) then the result is a midpoint
881 // if floor(C*) is even then C* = floor(C*) - logical right
882 // shift; C* has p decimal digits, correct by Prop. 1)
883 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
884 // shift; C* has p decimal digits, correct by Pr. 1)
885 // else
886 // C* = floor(C*) (logical right shift; C has p decimal digits,
887 // correct by Property 1)
888 // n = C* * 10^(e+x)
889
890 // determine inexactness of the rounding of C*
891 // if (0 < f* - 1/2 < T* ~= 10^(-x)) then
892 // the result is exact
893 // else // if (f* - 1/2 >= T*) then
894 // the result is inexact
895 if (fstar.w[4] > __bid_one_half192M[ind]
896 || (fstar.w[4] == __bid_one_half192M[ind]
897 && (fstar.w[3] || fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
898
899 // f* > 1/2 and the result may be exact
900 // Calculate f* - 1/2
901 tmp64 = fstar.w[4] - __bid_one_half192M[ind];
902 if (tmp64 || fstar.w[3] || fstar.w[2] > __bid_ten2mk192truncM[ind].w[2] ||
903 (fstar.w[2] == __bid_ten2mk192truncM[ind].w[2] &&
904 fstar.w[1] > __bid_ten2mk192truncM[ind].w[1]) ||
905 (fstar.w[2] == __bid_ten2mk192truncM[ind].w[2] &&
906 fstar.w[1] == __bid_ten2mk192truncM[ind].w[1] &&
907 fstar.w[0] > __bid_ten2mk192truncM[ind].w[0])) { // f* - 1/2 > 10^(-x)
908 // set the inexact flag
909 *pfpsf |= INEXACT_EXCEPTION;
910 is_inexact_lt_midpoint = 1;
911 } // else the result is exact
912 } else { // the result is inexact; f2* <= 1/2
913 // set the inexact flag
914 *pfpsf |= INEXACT_EXCEPTION;
915 tmp_fpa = 1;
916 is_inexact_gt_midpoint = 1;
917 }
918
919 // check for midpoints (could do this before determining inexactness)
920 if ((fstar.w[4] == 0) && (fstar.w[3] == 0)
921 && (fstar.w[2] || fstar.w[1] || fstar.w[0])
922 && (fstar.w[2] < __bid_ten2mk192truncM[ind].w[2]
923 || (fstar.w[2] == __bid_ten2mk192truncM[ind].w[2]
924 && fstar.w[1] < __bid_ten2mk192truncM[ind].w[1])
925 || (fstar.w[2] == __bid_ten2mk192truncM[ind].w[2]
926 && fstar.w[1] == __bid_ten2mk192truncM[ind].w[1]
927 && fstar.w[0] <= __bid_ten2mk192truncM[ind].w[0]))) {
928
929 // the result is a midpoint
930 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
931 // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
932 Cstar.w[0]--; // Cstar.w[0] is now even
933 if (tmp_fpa == 1)
934 tmp_fpa = 0;
935 is_midpoint_gt_even = 1;
936 is_inexact_lt_midpoint = 0;
937 is_inexact_gt_midpoint = 0;
938 } else { // else MP in [ODD, EVEN]
939 is_midpoint_lt_even = 1;
940 is_inexact_lt_midpoint = 0;
941 is_inexact_gt_midpoint = 0;
942 }
943 }
944 }
945
946 // check for rounding overflow
947 if (Cstar.w[1] == 0x0001ed09bead87c0ull &&
948 Cstar.w[0] == 0x378d8e6400000000ull) { // if Cstar = 10^34
949 tmp64 = x_exp + y_exp + ((UINT64) (ind + 6) << 49);
950 Cstar.w[1] = 0x0000314dc6448d93ull; // Cstar = 10^33
951 Cstar.w[0] = 0x38c15b0a00000000ull;
952 } else { // 10^33 <= Cstar <= 10^34 - 1
953 tmp64 = x_exp + y_exp + ((UINT64) (ind + 5) << 49); // ind+5 = q-34
954 }
955 if (tmp64 >= EXP_MAX + BIN_EXP_BIAS) { // possibble overflow
956 // exp >= emax for the result rounded to nearest even
957 if (rnd_mode == ROUNDING_TO_NEAREST
958 || rnd_mode == ROUNDING_TIES_AWAY) {
959 if (tmp64 > EXP_MAX + BIN_EXP_BIAS) {
960
961 // |res| >= 10^(p-1) * 10^(emax+1) <=> exp >= emax+1
962 res.w[1] = sign | 0x7800000000000000ull; // +/-inf
963 res.w[0] = 0x0ull;
964 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
965 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
966 is_overflow = 1;
967 } else { // not overflow
968 res.w[0] = Cstar.w[0];
969 res.w[1] = Cstar.w[1];
970 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
971 }
972 } else if (rnd_mode == ROUNDING_DOWN) {
973 if (!sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
974 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
975 Cstar.w[1] == 0x0000314dc6448d93ull &&
976 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
977 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
978 // res = +MAXFP
979 res.w[1] = 0x5fffed09bead87c0ull;
980 res.w[0] = 0x378d8e63ffffffffull; // (10^34-1) * 10^emax
981 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
982 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
983 is_overflow = 1;
984 } else if (sign && ((tmp64 > EXP_MAX + BIN_EXP_BIAS) ||
985 ((tmp64 == EXP_MAX + BIN_EXP_BIAS) &&
986 Cstar.w[1] == 0x0001ed09bead87c0ull &&
987 Cstar.w[0] == 0x378d8e63ffffffffull && // (10^34-1) * 10^emax
988 is_inexact_lt_midpoint))) {
989 res.w[1] = 0xf800000000000000ull; // -inf
990 res.w[0] = 0x0ull;
991 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
992 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
993 is_overflow = 1;
994 } else { // not overflow
995 res.w[0] = Cstar.w[0];
996 res.w[1] = Cstar.w[1];
997 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
998 }
999 } else if (rnd_mode == ROUNDING_UP) {
1000 if (sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
1001 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
1002 Cstar.w[1] == 0x0000314dc6448d93ull &&
1003 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
1004 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
1005 // res = -MAXFP
1006 res.w[1] = 0xdfffed09bead87c0ull;
1007 res.w[0] = 0x378d8e63ffffffffull; // -(10^34-1) * 10^emax
1008 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1009 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1010 is_overflow = 1;
1011 } else if (!sign && ((tmp64 > EXP_MAX + BIN_EXP_BIAS) ||
1012 ((tmp64 == EXP_MAX + BIN_EXP_BIAS) &&
1013 Cstar.w[1] == 0x0001ed09bead87c0ull &&
1014 Cstar.w[0] == 0x378d8e63ffffffffull && // (10^34-1) * 10^emax
1015 is_inexact_lt_midpoint))) {
1016 res.w[1] = 0x7800000000000000ull; // inf
1017 res.w[0] = 0x0ull;
1018 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1019 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1020 is_overflow = 1;
1021 } else { // not overflow
1022 res.w[0] = Cstar.w[0];
1023 res.w[1] = Cstar.w[1];
1024 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
1025 }
1026 } else { // if (rnd_mode == ROUNDING_TO_ZERO)
1027 if (!sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
1028 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
1029 Cstar.w[1] == 0x0000314dc6448d93ull &&
1030 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
1031 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
1032 // res = +MAXFP
1033 res.w[1] = 0x5fffed09bead87c0ull;
1034 res.w[0] = 0x378d8e63ffffffffull; // (10^34-1) * 10^emax
1035 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1036 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1037 is_overflow = 1;
1038 } else if (sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
1039 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
1040 Cstar.w[1] == 0x0000314dc6448d93ull &&
1041 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
1042 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
1043 // res = -MAXFP
1044 res.w[1] = 0xdfffed09bead87c0ull;
1045 res.w[0] = 0x378d8e63ffffffffull; // -(10^34-1) * 10^emax
1046 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1047 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1048 is_overflow = 1;
1049 } else { // not overflow
1050 res.w[0] = Cstar.w[0];
1051 res.w[1] = Cstar.w[1];
1052 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
1053 }
1054 }
1055 if (is_overflow) { // return for overflow
1056 // set the inexact flag
1057 *pfpsf |= INEXACT_EXCEPTION;
1058
1059 // set the overflow flag
1060 *pfpsf |= OVERFLOW_EXCEPTION;
1061
1062 // is_overflow = 1;
1063 BID_RETURN (res)}
1064 } else {
1065 res.w[0] = Cstar.w[0];
1066 res.w[1] = Cstar.w[1];
1067 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
1068 }
1069 res.w[1] |= sign;
1070 } else { // if (58 <= q <= 68) exact coefficient takes 192-226 bits
1071 // C = C + 1/2 * 10^x where the result C fits in 226 bits
1072 // (10^68 - 1 + 1/2 * 10^34 can be represented with 226 bits)
1073 // x = q - p = q - 34, 24 <= x <= 34
1074 // kx = 10^(-x) = __bid_ten2mk256M[ind]
1075 // C* = (C + 1/2 * 10^x) * 10^(-x)
1076 // the approximation of 10^(-x) was rounded up to 256 bits
1077 ind = q - 58; // 0 <= ind <= 10
1078 tmp64 = C.w[0];
1079 tmp64A = C.w[1];
1080
1081 // Note:
1082 // f* has 384 bits (more than 320 bits)
1083 // x - 1 = q - 35 = ind + 23
1084 // add two 64-bit words; e.g. for ind=0 <=> q=58, add 1/2*10^24
1085 C.w[0] = C.w[0] + __bid_midpoint128[ind + 4].w[0];
1086 C.w[1] = C.w[1] + __bid_midpoint128[ind + 4].w[1];
1087 if (C.w[0] < tmp64)
1088 C.w[1]++;
1089 if (C.w[1] < tmp64A)
1090 C.w[2]++;
1091 if (C.w[2] == 0)
1092 C.w[3]++;
1093 __mul_256x256_to_512 (P512, C, __bid_ten2mk256M[ind])
1094 // calculate C* and f*; C* is actually floor(C*) in this case
1095 // C* and f* need shifting and masking, as shown by
1096 // __bid_shiftright256M[] and __bid_maskhigh256M[]
1097 // C* has 128 bits; P512.w[7], P512.w[6], P512.w[5] need to be
1098 // shifted right by Ex-320 = __bid_shiftright256M[ind]
1099 shift = __bid_shiftright256M[ind]; // 15 <= shift <= 48
1100 if (shift == 32) {
1101 Cstar.w[0] =
1102 ((P512.w[5] >> 31) >> 1) | ((P512.w[6] << 31) << 1);
1103 Cstar.w[1] =
1104 ((P512.w[6] >> 31) >> 1) | ((P512.w[7] << 31) << 1);
1105 } else {
1106 Cstar.w[0] = (P512.w[5] >> shift) | (P512.w[6] << (64 - shift));
1107 Cstar.w[1] = (P512.w[6] >> shift) | (P512.w[7] << (64 - shift));
1108 }
1109 // f* has 384 bits
1110 fstar.w[5] = P512.w[5] & __bid_maskhigh256M[ind];
1111 fstar.w[4] = P512.w[4];
1112 fstar.w[3] = P512.w[3];
1113 fstar.w[2] = P512.w[2];
1114 fstar.w[1] = P512.w[1];
1115 fstar.w[0] = P512.w[0];
1116
1117 // the top Ex bits of 10^(-x) are T* = __bid_ten2mk256truncM[ind], e.g.
1118 // if x=24, T* = __bid_ten2mk256truncM[0] =
1119 // 0x9abe14cd44753b52c4926a9672793542d78c3615cf3a050cf23472530ce6e3ec =~
1120 // 10^(-24) * 2^335
1121 // if (0 < f* < 10^(-x)) then the result is a midpoint
1122 // if floor(C*) is even then C* = floor(C*) - logical right
1123 // shift; C* has p decimal digits, correct by Prop. 1)
1124 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
1125 // shift; C* has p decimal digits, correct by Pr. 1)
1126 // else
1127 // C* = floor(C*) (logical right shift; C has p decimal digits,
1128 // correct by Property 1)
1129 // n = C* * 10^(e+x)
1130
1131 // determine inexactness of the rounding of C*
1132 // if (0 < f* - 1/2 < T* ~= 10^(-x)) then
1133 // the result is exact
1134 // else // if (f* - 1/2 >= T*) then
1135 // the result is inexact
1136 if (fstar.w[5] > __bid_one_half256M[ind]
1137 || (fstar.w[5] == __bid_one_half256M[ind]
1138 && (fstar.w[4] || fstar.w[3] || fstar.w[2] || fstar.w[1]
1139 || fstar.w[0]))) {
1140
1141 // f* > 1/2 and the result may be exact
1142 // Calculate f* - 1/2
1143 tmp64 = fstar.w[5] - __bid_one_half256M[ind]; // tmp64 >= 0
1144 if (tmp64 || fstar.w[4] || fstar.w[3] > __bid_ten2mk256truncM[ind].w[3] ||
1145 (fstar.w[3] == __bid_ten2mk256truncM[ind].w[3] &&
1146 fstar.w[2] > __bid_ten2mk256truncM[ind].w[2]) ||
1147 (fstar.w[3] == __bid_ten2mk256truncM[ind].w[3] &&
1148 fstar.w[2] == __bid_ten2mk256truncM[ind].w[2] &&
1149 fstar.w[1] > __bid_ten2mk256truncM[ind].w[1]) ||
1150 (fstar.w[3] == __bid_ten2mk256truncM[ind].w[3] &&
1151 fstar.w[2] == __bid_ten2mk256truncM[ind].w[2] &&
1152 fstar.w[1] == __bid_ten2mk256truncM[ind].w[1] &&
1153 fstar.w[0] > __bid_ten2mk256truncM[ind].w[0])) { // f* - 1/2 > 10^(-x)
1154 // set the inexact flag
1155 *pfpsf |= INEXACT_EXCEPTION;
1156 is_inexact_lt_midpoint = 1;
1157 } // else the result is exact
1158 } else { // the result is inexact; f2* <= 1/2
1159 // set the inexact flag
1160 *pfpsf |= INEXACT_EXCEPTION;
1161 tmp_fpa = 1;
1162 is_inexact_gt_midpoint = 1;
1163 }
1164
1165 // check for midpoints (could do this before determining inexactness)
1166 if ((fstar.w[5] == 0) && (fstar.w[4] == 0)
1167 && (fstar.w[3] || fstar.w[2] || fstar.w[1] || fstar.w[0])
1168 && (fstar.w[3] < __bid_ten2mk256truncM[ind].w[3]
1169 || (fstar.w[3] == __bid_ten2mk256truncM[ind].w[3]
1170 && fstar.w[2] < __bid_ten2mk256truncM[ind].w[2])
1171 || (fstar.w[3] == __bid_ten2mk256truncM[ind].w[3]
1172 && fstar.w[2] == __bid_ten2mk256truncM[ind].w[2]
1173 && fstar.w[1] < __bid_ten2mk256truncM[ind].w[1])
1174 || (fstar.w[3] == __bid_ten2mk256truncM[ind].w[3]
1175 && fstar.w[2] == __bid_ten2mk256truncM[ind].w[2]
1176 && fstar.w[1] == __bid_ten2mk256truncM[ind].w[1]
1177 && fstar.w[0] <= __bid_ten2mk256truncM[ind].w[1]))) {
1178
1179 // the result is a midpoint
1180 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1181 // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
1182 Cstar.w[0]--; // Cstar.w[0] is now even
1183 if (tmp_fpa == 1)
1184 tmp_fpa = 0;
1185 is_midpoint_gt_even = 1;
1186 is_inexact_lt_midpoint = 0;
1187 is_inexact_gt_midpoint = 0;
1188 } else { // else MP in [ODD, EVEN]
1189 is_midpoint_lt_even = 1;
1190 is_inexact_lt_midpoint = 0;
1191 is_inexact_gt_midpoint = 0;
1192 }
1193 }
1194 // check for rounding overflow
1195 if (Cstar.w[1] == 0x0001ed09bead87c0ull &&
1196 Cstar.w[0] == 0x378d8e6400000000ull) { // if Cstar = 10^34
1197 tmp64 = x_exp + y_exp + ((UINT64) (ind + 25) << 49);
1198 Cstar.w[1] = 0x0000314dc6448d93ull; // Cstar = 10^33
1199 Cstar.w[0] = 0x38c15b0a00000000ull;
1200 } else { // 10^33 <= Cstar <= 10^34 - 1
1201 tmp64 = x_exp + y_exp + ((UINT64) (ind + 24) << 49); // ind+24 = q-34
1202 }
1203 if (tmp64 >= EXP_MAX + BIN_EXP_BIAS) { // possibble overflow
1204 // exp >= emax for the result rounded to nearest even
1205 if (rnd_mode == ROUNDING_TO_NEAREST
1206 || rnd_mode == ROUNDING_TIES_AWAY) {
1207 if (tmp64 > EXP_MAX + BIN_EXP_BIAS) {
1208
1209 // |res| >= 10^(p-1) * 10^(emax+1) <=> exp >= emax+1
1210 res.w[1] = sign | 0x7800000000000000ull; // +/-inf
1211 res.w[0] = 0x0ull;
1212 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1213 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1214 is_overflow = 1;
1215 } else { // not overflow
1216 res.w[0] = Cstar.w[0];
1217 res.w[1] = Cstar.w[1];
1218 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
1219 }
1220 } else if (rnd_mode == ROUNDING_DOWN) {
1221 if (!sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
1222 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
1223 Cstar.w[1] == 0x0000314dc6448d93ull &&
1224 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
1225 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
1226 // res = +MAXFP
1227 res.w[1] = 0x5fffed09bead87c0ull;
1228 res.w[0] = 0x378d8e63ffffffffull; // (10^34-1) * 10^emax
1229 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1230 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1231 is_overflow = 1;
1232 } else if (sign && ((tmp64 > EXP_MAX + BIN_EXP_BIAS) ||
1233 ((tmp64 == EXP_MAX + BIN_EXP_BIAS) &&
1234 Cstar.w[1] == 0x0001ed09bead87c0ull &&
1235 Cstar.w[0] == 0x378d8e63ffffffffull && // (10^34-1) * 10^emax
1236 is_inexact_lt_midpoint))) {
1237 res.w[1] = 0xf800000000000000ull; // -inf
1238 res.w[0] = 0x0ull;
1239 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1240 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1241 is_overflow = 1;
1242 } else { // not overflow
1243 res.w[0] = Cstar.w[0];
1244 res.w[1] = Cstar.w[1];
1245 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
1246 }
1247 } else if (rnd_mode == ROUNDING_UP) {
1248 if (sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
1249 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
1250 Cstar.w[1] == 0x0000314dc6448d93ull &&
1251 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
1252 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
1253 // res = -MAXFP
1254 res.w[1] = 0xdfffed09bead87c0ull;
1255 res.w[0] = 0x378d8e63ffffffffull; // -(10^34-1) * 10^emax
1256 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1257 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1258 is_overflow = 1;
1259 } else if (!sign && ((tmp64 > EXP_MAX + BIN_EXP_BIAS) ||
1260 ((tmp64 == EXP_MAX + BIN_EXP_BIAS) &&
1261 Cstar.w[1] == 0x0001ed09bead87c0ull &&
1262 Cstar.w[0] == 0x378d8e63ffffffffull && // (10^34-1) * 10^emax
1263 is_inexact_lt_midpoint))) {
1264 res.w[1] = 0x7800000000000000ull; // inf
1265 res.w[0] = 0x0ull;
1266 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1267 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1268 is_overflow = 1;
1269 } else { // not overflow
1270 res.w[0] = Cstar.w[0];
1271 res.w[1] = Cstar.w[1];
1272 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
1273 }
1274 } else { // if (rnd_mode == ROUNDING_TO_ZERO)
1275 if (!sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
1276 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
1277 Cstar.w[1] == 0x0000314dc6448d93ull &&
1278 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
1279 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
1280 // res = +MAXFP
1281 res.w[1] = 0x5fffed09bead87c0ull;
1282 res.w[0] = 0x378d8e63ffffffffull; // (10^34-1) * 10^emax
1283 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1284 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1285 is_overflow = 1;
1286 } else if (sign && (tmp64 > EXP_MAX + BIN_EXP_BIAS) &&
1287 !(tmp64 == EXP_MAX + BIN_EXP_BIAS + EXP_P1 &&
1288 Cstar.w[1] == 0x0000314dc6448d93ull &&
1289 Cstar.w[0] == 0x38c15b0a00000000ull && // 10^33 * 10^(emax+1)
1290 (is_midpoint_lt_even || is_inexact_gt_midpoint))) {
1291 // res = -MAXFP
1292 res.w[1] = 0xdfffed09bead87c0ull;
1293 res.w[0] = 0x378d8e63ffffffffull; // -(10^34-1) * 10^emax
1294 *pfpsf |= INEXACT_EXCEPTION; // set the inexact flag
1295 *pfpsf |= OVERFLOW_EXCEPTION; // set the overflow flag
1296 is_overflow = 1;
1297 } else { // not overflow
1298 res.w[0] = Cstar.w[0];
1299 res.w[1] = Cstar.w[1];
1300 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
1301 }
1302 }
1303 if (is_overflow) { // return for overflow
1304 // set the inexact flag
1305 *pfpsf |= INEXACT_EXCEPTION;
1306
1307 // set the overflow flag
1308 *pfpsf |= OVERFLOW_EXCEPTION;
1309
1310 // is_overflow = 1;
1311 BID_RETURN (res);
1312 }
1313 } else {
1314 res.w[0] = Cstar.w[0];
1315 res.w[1] = Cstar.w[1];
1316 res.w[1] |= (tmp64 - BIN_EXP_BIAS);
1317 }
1318 res.w[1] |= sign;
1319 }
1320
1321 // general correction from RN to RA, RM, RP, RZ
1322 if (rnd_mode != ROUNDING_TO_NEAREST && !is_overflow) { // overflow is solved
1323 x_exp = res.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1324 C1_hi = res.w[1] & MASK_COEFF;
1325 C1_lo = res.w[0];
1326 if ((!sign && ((rnd_mode == ROUNDING_UP && is_inexact_lt_midpoint) ||
1327 ((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_UP) &&
1328 is_midpoint_gt_even))) ||
1329 (sign && ((rnd_mode == ROUNDING_DOWN && is_inexact_lt_midpoint) ||
1330 ((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_DOWN) &&
1331 is_midpoint_gt_even)))) {
1332
1333 // C1 = C1 + 1
1334 C1_lo = C1_lo + 1;
1335 if (C1_lo == 0) { // rounding overflow in the low 64 bits
1336 C1_hi = C1_hi + 1;
1337 if (C1_hi == 0x0001ed09bead87c0ull
1338 && C1_lo == 0x378d8e6400000000ull) {
1339
1340 // C1 = 10^34 => rounding overflow
1341 C1_hi = 0x0000314dc6448d93ull;
1342 C1_lo = 0x38c15b0a00000000ull; // 10^33
1343 x_exp = x_exp + EXP_P1;
1344 }
1345 }
1346 } else if ((is_midpoint_lt_even || is_inexact_gt_midpoint)
1347 && ((sign && (rnd_mode == ROUNDING_UP ||
1348 rnd_mode == ROUNDING_TO_ZERO)) ||
1349 (!sign && (rnd_mode == ROUNDING_DOWN ||
1350 rnd_mode == ROUNDING_TO_ZERO)))) {
1351
1352 // C1 = C1 - 1
1353 C1_lo = C1_lo - 1;
1354 if (C1_lo == 0xffffffffffffffffull)
1355 C1_hi--;
1356
1357 // check if we crossed into the lower decade
1358 if (C1_hi == 0x0000314dc6448d93ull && C1_lo == 0x38c15b09ffffffffull) {
1359 // 10^33 - 1
1360 C1_hi = 0x0001ed09bead87c0ull; // 10^34 - 1
1361 C1_lo = 0x378d8e63ffffffffull;
1362 x_exp = x_exp - EXP_P1; // no underflow (TO CHECK)
1363 }
1364 } else {
1365 ; // exact, the result is already correct
1366 }
1367
1368 // assemble the result
1369 res.w[1] = x_exp | C1_hi;
1370 res.w[0] = C1_lo;
1371 }
1372 res.w[1] |= sign;
1373 BID_RETURN (res);
1374 }
1375 _underflow_path:
1376 // got here because q - P34 < ind where ind = EMIN - ex - ey
1377 // q is the number of digits in C; ind is the [positive] exponent of the
1378 // negative power of 10 that must multiply C in order to make the result's
1379 // exponent equal to e_min - P34 + 1 = -6176
1380 ind =
1381 (int) (((SINT64) EXP_MIN + (SINT64) BIN_EXP_BIAS - (SINT64) x_exp -
1382 (SINT64) y_exp) >> 49);
1383
1384 // q - P34 < ind => -P34 + 1 < ind => -P34 + 2 <= ind
1385 // ind = EMIN - ex - ey < -6176 + 6176 + 6176 = 6176
1386 if (q < ind) { // q - ind < 0; result rounds to 0 when rounding to nearest
1387 // set the inexact and underflow flags
1388 *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION);
1389 res.w[1] = EXP_MIN; // EXP_MIN = 0x0
1390 res.w[0] = 0x0;
1391 if (rnd_mode != ROUNDING_TO_NEAREST) {
1392 if ((rnd_mode == ROUNDING_DOWN && sign) ||
1393 (rnd_mode == ROUNDING_UP && !sign))
1394 res.w[0] = 0x0000000000000001ull;
1395 }
1396 } else if (q == ind) { // q - ind = 0; result rounds to 0 or +/-1*10^EMIN
1397 // set the inexact and underflow flags
1398 *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION);
1399
1400 // if C <= 5*10^(q-1) then C = 0 else C = 1
1401 if (q <= 19) {
1402 if (C.w[0] == __bid_midpoint64[q - 1]) { // C = 0.5 * 10^emin
1403 if (rnd_mode == ROUNDING_TO_NEAREST || (rnd_mode == ROUNDING_DOWN
1404 && !sign) || (rnd_mode == ROUNDING_UP && sign)
1405 || rnd_mode == ROUNDING_TO_ZERO) {
1406 res.w[1] = EXP_MIN;
1407 res.w[0] = 0;
1408 } else {
1409 res.w[1] = EXP_MIN;
1410 res.w[0] = 1;
1411 }
1412 } else if (C.w[0] < __bid_midpoint64[q - 1]) { // C < 0.5 * 10^emin
1413 if (rnd_mode == ROUNDING_TO_NEAREST
1414 || rnd_mode == ROUNDING_TIES_AWAY
1415 || (rnd_mode == ROUNDING_DOWN && !sign)
1416 || (rnd_mode == ROUNDING_UP && sign)
1417 || rnd_mode == ROUNDING_TO_ZERO) {
1418 res.w[1] = EXP_MIN;
1419 res.w[0] = 0;
1420 } else {
1421 res.w[1] = EXP_MIN;
1422 res.w[0] = 1;
1423 }
1424 } else { // C > 0.5 * 10^emin
1425 if (rnd_mode == ROUNDING_TO_NEAREST
1426 || rnd_mode == ROUNDING_TIES_AWAY
1427 || (rnd_mode == ROUNDING_DOWN && sign)
1428 || (rnd_mode == ROUNDING_UP && !sign)) {
1429 res.w[1] = EXP_MIN;
1430 res.w[0] = 1;
1431 } else {
1432 res.w[1] = EXP_MIN;
1433 res.w[0] = 0;
1434 }
1435 }
1436 } else if (q <= 38) { // 20 <= q <= 38
1437 // if q <= P34 = 34 the exact result rounded to P34 digits with unbounded
1438 // exponent will have an exponent smaller than e_min; otherwise if
1439 // 35 <= q <= 38, it depends
1440 if (C.w[1] == __bid_midpoint128[q - 20].w[1] &&
1441 C.w[0] == __bid_midpoint128[q - 20].w[0]) { // C = 0.5 * 10^emin
1442 if (rnd_mode == ROUNDING_TO_NEAREST || (rnd_mode == ROUNDING_DOWN
1443 && !sign) || (rnd_mode == ROUNDING_UP && sign)
1444 || rnd_mode == ROUNDING_TO_ZERO) {
1445 res.w[1] = EXP_MIN;
1446 res.w[0] = 0;
1447 } else {
1448 res.w[1] = EXP_MIN;
1449 res.w[0] = 1;
1450 }
1451 } else if (C.w[1] < __bid_midpoint128[q - 20].w[1] ||
1452 (C.w[1] == __bid_midpoint128[q - 20].w[1] &&
1453 C.w[0] < __bid_midpoint128[q - 20].w[0])) { // C < 0.5 * 10^emin
1454 if (rnd_mode == ROUNDING_TO_NEAREST
1455 || rnd_mode == ROUNDING_TIES_AWAY
1456 || (rnd_mode == ROUNDING_DOWN && !sign)
1457 || (rnd_mode == ROUNDING_UP && sign)
1458 || rnd_mode == ROUNDING_TO_ZERO) {
1459 res.w[1] = EXP_MIN;
1460 res.w[0] = 0;
1461 } else {
1462 res.w[1] = EXP_MIN;
1463 res.w[0] = 1;
1464 }
1465 } else { // C > 0.5 * 10^emin
1466 if (rnd_mode == ROUNDING_TO_NEAREST
1467 || rnd_mode == ROUNDING_TIES_AWAY
1468 || (rnd_mode == ROUNDING_DOWN && sign)
1469 || (rnd_mode == ROUNDING_UP && !sign)) {
1470 res.w[1] = EXP_MIN;
1471 res.w[0] = 1;
1472 } else {
1473 res.w[1] = EXP_MIN;
1474 res.w[0] = 0;
1475 }
1476 }
1477 } else if (q <= 58) { // 39 <= q <= 58
1478 // Note: for q = 58 C may take 256 bits, so need to test C.w[3]
1479 if (C.w[3] == 0x0 && C.w[2] == __bid_midpoint192[q - 39].w[2] &&
1480 C.w[1] == __bid_midpoint192[q - 39].w[1] &&
1481 C.w[0] == __bid_midpoint192[q - 39].w[0]) { // C = 0.5 * 10^emin
1482 if (rnd_mode == ROUNDING_TO_NEAREST || (rnd_mode == ROUNDING_DOWN
1483 && !sign) || (rnd_mode == ROUNDING_UP && sign)
1484 || rnd_mode == ROUNDING_TO_ZERO) {
1485 res.w[1] = EXP_MIN;
1486 res.w[0] = 0;
1487 } else {
1488 res.w[1] = EXP_MIN;
1489 res.w[0] = 1;
1490 }
1491 } else if ((C.w[3] == 0x0 && C.w[2] < __bid_midpoint192[q - 39].w[2]) ||
1492 (C.w[3] == 0x0 && C.w[2] == __bid_midpoint192[q - 39].w[2] &&
1493 C.w[1] < __bid_midpoint192[q - 39].w[1]) || (C.w[3] == 0x0 &&
1494 C.w[2] == __bid_midpoint192[q - 39].w[2] &&
1495 C.w[1] == __bid_midpoint192[q - 39].w[1] &&
1496 C.w[0] < __bid_midpoint192[q - 39].w[0])) { // C < 0.5 * 10^emin
1497 if (rnd_mode == ROUNDING_TO_NEAREST
1498 || rnd_mode == ROUNDING_TIES_AWAY
1499 || (rnd_mode == ROUNDING_DOWN && !sign)
1500 || (rnd_mode == ROUNDING_UP && sign)
1501 || rnd_mode == ROUNDING_TO_ZERO) {
1502 res.w[1] = EXP_MIN;
1503 res.w[0] = 0;
1504 } else {
1505 res.w[1] = EXP_MIN;
1506 res.w[0] = 1;
1507 }
1508 } else { // C > 0.5 * 10^emin
1509 if (rnd_mode == ROUNDING_TO_NEAREST
1510 || rnd_mode == ROUNDING_TIES_AWAY
1511 || (rnd_mode == ROUNDING_DOWN && sign)
1512 || (rnd_mode == ROUNDING_UP && !sign)) {
1513 res.w[1] = EXP_MIN;
1514 res.w[0] = 1;
1515 } else {
1516 res.w[1] = EXP_MIN;
1517 res.w[0] = 0;
1518 }
1519 }
1520 } else { // if (q <= 68), i.e. 59 <= q <= 68
1521 if (C.w[3] == __bid_midpoint256[q - 59].w[3] &&
1522 C.w[2] == __bid_midpoint256[q - 59].w[2] &&
1523 C.w[1] == __bid_midpoint256[q - 59].w[1] &&
1524 C.w[0] == __bid_midpoint256[q - 59].w[0]) { // C = 0.5 * 10^emin
1525 if (rnd_mode == ROUNDING_TO_NEAREST || (rnd_mode == ROUNDING_DOWN
1526 && !sign) || (rnd_mode == ROUNDING_UP && sign)
1527 || rnd_mode == ROUNDING_TO_ZERO) {
1528 res.w[1] = EXP_MIN;
1529 res.w[0] = 0;
1530 } else {
1531 res.w[1] = EXP_MIN;
1532 res.w[0] = 1;
1533 }
1534 } else if (C.w[3] < __bid_midpoint256[q - 59].w[3] ||
1535 (C.w[3] == __bid_midpoint256[q - 59].w[3] &&
1536 C.w[2] < __bid_midpoint256[q - 59].w[2]) ||
1537 (C.w[3] == __bid_midpoint256[q - 59].w[3] &&
1538 C.w[2] == __bid_midpoint256[q - 59].w[2] &&
1539 C.w[1] < __bid_midpoint256[q - 59].w[1]) ||
1540 (C.w[3] == __bid_midpoint256[q - 59].w[3] &&
1541 C.w[2] == __bid_midpoint256[q - 59].w[2] &&
1542 C.w[1] == __bid_midpoint256[q - 59].w[1] &&
1543 C.w[0] < __bid_midpoint256[q - 59].w[0])) { // C < 0.5 * 10^emin
1544 if (rnd_mode == ROUNDING_TO_NEAREST
1545 || rnd_mode == ROUNDING_TIES_AWAY
1546 || (rnd_mode == ROUNDING_DOWN && !sign)
1547 || (rnd_mode == ROUNDING_UP && sign)
1548 || rnd_mode == ROUNDING_TO_ZERO) {
1549 res.w[1] = EXP_MIN;
1550 res.w[0] = 0;
1551 } else {
1552 res.w[1] = EXP_MIN;
1553 res.w[0] = 1;
1554 }
1555 } else { // C > 0.5 * 10^emin
1556 if (rnd_mode == ROUNDING_TO_NEAREST
1557 || rnd_mode == ROUNDING_TIES_AWAY
1558 || (rnd_mode == ROUNDING_DOWN && sign)
1559 || (rnd_mode == ROUNDING_UP && !sign)) {
1560 res.w[1] = EXP_MIN;
1561 res.w[0] = 1;
1562 } else {
1563 res.w[1] = EXP_MIN;
1564 res.w[0] = 0;
1565 }
1566 }
1567 }
1568 } else { // if 0 < q - ind < P34 <=> 1 <= q - ind <= P34 - 1 = 33
1569 // In general -P34 + 2 <= ind <= 6176 => -P34 + 2 <= ind < q =>
1570 // -P34 + 2 <= ind <= q - 1
1571 if (rnd_mode != ROUNDING_TO_NEAREST) {
1572 is_inexact_lt_midpoint = 0;
1573 is_inexact_gt_midpoint = 0;
1574 is_midpoint_lt_even = 0;
1575 is_midpoint_gt_even = 0;
1576 }
1577 if (ind <= 0) { // 0 <= -ind
1578 // the result is exact
1579 res.w[1] = (x_exp + y_exp - BIN_EXP_BIAS) | C.w[1];
1580 res.w[0] = C.w[0];
1581
1582 // because the result is exact the U and I status flags are not set
1583 } else {
1584
1585 // if ind > 0 <=> 1 <= ind <= q - 1; must remove ind digits
1586 // from C, which may have up to 68 digits; note that q >= ind + 1 >= 2
1587 // Note: there is no underflow in some cases when the coefficient of
1588 // the result is 10^33 or 10^33 - 1
1589 if (q <= 18) { // 2 <= q <= 18
1590 __bid_round64_2_18 (q, ind, C.w[0], &res.w[0], &incr_exp,
1591 &is_midpoint_lt_even, &is_midpoint_gt_even,
1592 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
1593 if (incr_exp) {
1594
1595 // multiply by 10; this cannot be 10^33
1596 __mul_64x64_to_128MACH (res, res.w[0], __bid_ten2k64[1]);
1597 res.w[1] |= (UINT64) EXP_MIN;
1598 } else { // underflow
1599 res.w[1] = (UINT64) EXP_MIN;
1600 }
1601 if (is_midpoint_lt_even || is_midpoint_gt_even
1602 || is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
1603
1604 // set the inexact and underflow flags
1605 *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION);
1606 }
1607 } else if (q <= 38) { // 19 <= q <= 38
1608 P128.w[1] = C.w[1];
1609 P128.w[0] = C.w[0];
1610 __bid_round128_19_38 (q, ind, P128, &res, &incr_exp,
1611 &is_midpoint_lt_even, &is_midpoint_gt_even,
1612 &is_inexact_lt_midpoint,
1613 &is_inexact_gt_midpoint);
1614 if (incr_exp) {
1615
1616 // multiply by 10 and check is this is 10^33, because in that case
1617 // it is possible that this is not underflow
1618 if (q - ind <= 19) {
1619 __mul_64x64_to_128MACH (res, res.w[0], __bid_ten2k64[1]);
1620 } else { // if 20 <= q - ind
1621 __mul_128x64_to_128 (res, __bid_ten2k64[1], res);
1622 }
1623 if ((q - ind + 1) == P34) { // the result is 10^(P34-1)
1624 // if the result rounded directly to P34 digits is the same, then
1625 // there is no underflow
1626 __bid_round128_19_38 (q, ind - 1, P128, &R128, &incr_exp1,
1627 &is_midpoint_lt_even1,
1628 &is_midpoint_gt_even1,
1629 &is_inexact_lt_midpoint1,
1630 &is_inexact_gt_midpoint1);
1631 if (res.w[1] == R128.w[1] && res.w[0] == R128.w[0]) {
1632 no_underflow = 1;
1633 }
1634 }
1635 // res.w[1] |= (UINT64)EXP_MIN; // redundant
1636 } else { // underflow
1637 // res.w[1] = (UINT64)EXP_MIN | res.w[1]; // redundant
1638 }
1639 if (is_midpoint_lt_even || is_midpoint_gt_even
1640 || is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
1641
1642 // set the inexact and underflow flags
1643 *pfpsf |= INEXACT_EXCEPTION;
1644 is_inexact = 1;
1645 if (!no_underflow)
1646 *pfpsf |= UNDERFLOW_EXCEPTION;
1647 }
1648 } else if (q <= 57) { // 39 <= q <= 57
1649 P192.w[2] = C.w[2];
1650 P192.w[1] = C.w[1];
1651 P192.w[0] = C.w[0];
1652 __bid_round192_39_57 (q, ind, P192, &R192, &incr_exp,
1653 &is_midpoint_lt_even, &is_midpoint_gt_even,
1654 &is_inexact_lt_midpoint,
1655 &is_inexact_gt_midpoint);
1656 if (incr_exp) {
1657
1658 // multiply by 10 and check is this is 10^33, because in that case
1659 // it is possible that this is not underflow
1660 res.w[1] = R192.w[1]; // res has q - ind digits
1661 res.w[0] = R192.w[0];
1662 if (q - ind <= 19) {
1663 __mul_64x64_to_128MACH (res, res.w[0], __bid_ten2k64[1]);
1664 } else { // if 20 <= q - ind
1665 __mul_128x64_to_128 (res, __bid_ten2k64[1], res);
1666 }
1667 if ((q - ind + 1) == P34) { // the result is 10^(P34-1)
1668 // if the result rounded directly to P34 digits is the same, then
1669 // there is no underflow
1670 __bid_round192_39_57 (q, ind - 1, P192, &R192, &incr_exp1,
1671 &is_midpoint_lt_even1,
1672 &is_midpoint_gt_even1,
1673 &is_inexact_lt_midpoint1,
1674 &is_inexact_gt_midpoint1);
1675 if (res.w[1] == R192.w[1] && res.w[0] == R192.w[0]) {
1676 no_underflow = 1;
1677 }
1678 }
1679 // res.w[1] |= (UINT64)EXP_MIN; // redundant
1680 } else { // underflow
1681 res.w[1] = (UINT64) EXP_MIN | R192.w[1];
1682 res.w[0] = R192.w[0];
1683 }
1684 if (is_midpoint_lt_even || is_midpoint_gt_even
1685 || is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
1686
1687 // set the inexact and underflow flags
1688 *pfpsf |= INEXACT_EXCEPTION;
1689 is_inexact = 1;
1690 if (!no_underflow)
1691 *pfpsf |= UNDERFLOW_EXCEPTION;
1692 }
1693 } else if (q <= 76) { // 58 <= q <= 76 (actually 58 <= q <= 68)
1694 P256.w[3] = C.w[3];
1695 P256.w[2] = C.w[2];
1696 P256.w[1] = C.w[1];
1697 P256.w[0] = C.w[0];
1698 __bid_round256_58_76 (q, ind, P256, &R256, &incr_exp,
1699 &is_midpoint_lt_even, &is_midpoint_gt_even,
1700 &is_inexact_lt_midpoint,
1701 &is_inexact_gt_midpoint);
1702 if (incr_exp) {
1703
1704 // multiply by 10 and check is this is 10^33, because in that case
1705 // it is possible that this is not underflow
1706 res.w[1] = R256.w[1]; // res has q - ind digits
1707 res.w[0] = R256.w[0];
1708 if (q - ind <= 19) {
1709 __mul_64x64_to_128MACH (res, res.w[0], __bid_ten2k64[1]);
1710 } else { // if 20 <= q - ind
1711 __mul_128x64_to_128 (res, __bid_ten2k64[1], res);
1712 }
1713 if ((q - ind + 1) == P34) { // the result is 10^(P34-1)
1714 // if the result rounded directly to P34 digits is the same, then
1715 // there is no underflow
1716 __bid_round256_58_76 (q, ind - 1, P256, &R256, &incr_exp1,
1717 &is_midpoint_lt_even1,
1718 &is_midpoint_gt_even1,
1719 &is_inexact_lt_midpoint1,
1720 &is_inexact_gt_midpoint1);
1721 if (res.w[1] == R256.w[1] && res.w[0] == R256.w[0]) {
1722 no_underflow = 1;
1723 }
1724 }
1725 // res.w[1] |= (UINT64)EXP_MIN; // redundant
1726 } else { // underflow
1727 res.w[1] = (UINT64) EXP_MIN | R256.w[1];
1728 res.w[0] = R256.w[0];
1729 }
1730 if (is_midpoint_lt_even || is_midpoint_gt_even
1731 || is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
1732
1733 // set the inexact and underflow flags
1734 *pfpsf |= INEXACT_EXCEPTION;
1735 is_inexact = 1;
1736 if (!no_underflow)
1737 *pfpsf |= UNDERFLOW_EXCEPTION;
1738 }
1739 }
1740 }
1741
1742 // general correction from RN to RA, RM, RP, RZ
1743 if (rnd_mode != ROUNDING_TO_NEAREST) {
1744 x_exp = res.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1745 // this must be e_min
1746 C1_hi = res.w[1] & MASK_COEFF;
1747 C1_lo = res.w[0];
1748 if ((!sign && ((rnd_mode == ROUNDING_UP && is_inexact_lt_midpoint) ||
1749 ((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_UP) &&
1750 is_midpoint_gt_even))) ||
1751 (sign && ((rnd_mode == ROUNDING_DOWN && is_inexact_lt_midpoint) ||
1752 ((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_DOWN) &&
1753 is_midpoint_gt_even)))) {
1754
1755 // C1 = C1 + 1
1756 C1_lo = C1_lo + 1;
1757 if (C1_lo == 0) { // rounding overflow in the low 64 bits
1758 C1_hi = C1_hi + 1;
1759 if (C1_hi == 0x0001ed09bead87c0ull
1760 && C1_lo == 0x378d8e6400000000ull) {
1761
1762 // C1 = 10^34 => rounding overflow (not possible) TO CHECK
1763 C1_hi = 0x0000314dc6448d93ull;
1764 C1_lo = 0x38c15b0a00000000ull; // 10^33
1765 x_exp = x_exp + EXP_P1; // this must be e_min
1766 }
1767 }
1768 } else if ((is_midpoint_lt_even || is_inexact_gt_midpoint) &&
1769 ((sign &&
1770 (rnd_mode == ROUNDING_UP || rnd_mode == ROUNDING_TO_ZERO)) ||
1771 (!sign &&
1772 (rnd_mode == ROUNDING_DOWN || rnd_mode == ROUNDING_TO_ZERO)))) {
1773
1774 // C1 = C1 - 1 (the exponent is emin already)
1775 C1_lo = C1_lo - 1;
1776 if (C1_lo == 0xffffffffffffffffull)
1777 C1_hi--;
1778
1779 // cannot cross into the lower decade anymore, but the result can be 0
1780 } else {
1781 ; // exact, the result is already correct
1782 }
1783
1784 // no overflow is possible
1785 // assemble the result
1786 res.w[1] = x_exp | C1_hi;
1787 res.w[0] = C1_lo;
1788
1789 // Now fix the case where the general rounding routine returned a non-tiny
1790 // result, but after the correction for rounding modes other than to
1791 // nearest, the result is less in magnitude than 100...0[34] * 10^(-6176)
1792 // (this is due to the fact that the general rounding routine works only
1793 // with rounding to nearest)
1794 if (is_inexact && (x_exp == EXP_MIN)
1795 && (C1_hi < 0x0000314dc6448d93ull
1796 || (C1_hi == 0x0000314dc6448d93ull
1797 && C1_lo < 0x38c15b0a00000000ull))) {
1798 *pfpsf |= UNDERFLOW_EXCEPTION;
1799 }
1800 }
1801 }
1802 res.w[1] |= sign;
1803 BID_RETURN (res);
1804 }
Mirror of https://gcc.gnu.org/git/gcc.git