1 /******************************************************************************
2 Copyright (c) 2007-2011, Intel Corp.
3 All rights reserved.
4
5 Redistribution and use in source and binary forms, with or without
6 modification, are permitted provided that the following conditions are met:
7
8 * Redistributions of source code must retain the above copyright notice,
9 this list of conditions and the following disclaimer.
10 * Redistributions in binary form must reproduce the above copyright
11 notice, this list of conditions and the following disclaimer in the
12 documentation and/or other materials provided with the distribution.
13 * Neither the name of Intel Corporation nor the names of its contributors
14 may be used to endorse or promote products derived from this software
15 without specific prior written permission.
16
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
27 THE POSSIBILITY OF SUCH DAMAGE.
28 ******************************************************************************/
29
30 /*****************************************************************************
31 *
32 * Helper add functions (for fma)
33 *
34 * __BID_INLINE__ BID_UINT64 bid_get_add64(
35 * BID_UINT64 sign_x, int exponent_x, BID_UINT64 coefficient_x,
36 * BID_UINT64 sign_y, int exponent_y, BID_UINT64 coefficient_y,
37 * int rounding_mode)
38 *
39 * __BID_INLINE__ BID_UINT64 bid_get_add128(
40 * BID_UINT64 sign_x, int exponent_x, BID_UINT64 coefficient_x,
41 * BID_UINT64 sign_y, int final_exponent_y, BID_UINT128 CY,
42 * int extra_digits, int rounding_mode)
43 *
44 *****************************************************************************
45 *
46 * Algorithm description:
47 *
48 * bid_get_add64: same as BID64 add, but arguments are unpacked and there
49 * are no special case checks
50 *
51 * bid_get_add128: add 64-bit coefficient to 128-bit product (which contains
52 * 16+extra_digits decimal digits),
53 * return BID64 result
54 * - the exponents are compared and the two coefficients are
55 * properly aligned for addition/subtraction
56 * - multiple paths are needed
57 * - final result exponent is calculated and the lower term is
58 * rounded first if necessary, to avoid manipulating
59 * coefficients longer than 128 bits
60 *
61 ****************************************************************************/
62
63 #ifndef _INLINE_BID_ADD_H_
64 #define _INLINE_BID_ADD_H_
65
66 #include "bid_internal.h"
67
68 #define MAX_FORMAT_DIGITS 16
69 #define DECIMAL_EXPONENT_BIAS 398
70 #define MASK_BINARY_EXPONENT 0x7ff0000000000000ull
71 #define BINARY_EXPONENT_BIAS 0x3ff
72 #define UPPER_EXPON_LIMIT 51
73
74 ///////////////////////////////////////////////////////////////////////
75 //
76 // bid_get_add64() is essentially the same as bid_add(), except that
77 // the arguments are unpacked
78 //
79 //////////////////////////////////////////////////////////////////////
80 __BID_INLINE__ BID_UINT64
bid_get_add64(BID_UINT64 sign_x,int exponent_x,BID_UINT64 coefficient_x,BID_UINT64 sign_y,int exponent_y,BID_UINT64 coefficient_y,int rounding_mode,unsigned * fpsc)81 bid_get_add64 (BID_UINT64 sign_x, int exponent_x, BID_UINT64 coefficient_x,
82 BID_UINT64 sign_y, int exponent_y, BID_UINT64 coefficient_y,
83 int rounding_mode, unsigned *fpsc) {
84 BID_UINT128 CA, CT, CT_new;
85 BID_UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab,
86 rem_a;
87 BID_UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp,
88 C64_new;
89 int_double tempx;
90 int exponent_a, exponent_b, diff_dec_expon;
91 int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
92 unsigned rmode, status;
93
94 // sort arguments by exponent
95 if (exponent_x <= exponent_y) {
96 sign_a = sign_y;
97 exponent_a = exponent_y;
98 coefficient_a = coefficient_y;
99 sign_b = sign_x;
100 exponent_b = exponent_x;
101 coefficient_b = coefficient_x;
102 } else {
103 sign_a = sign_x;
104 exponent_a = exponent_x;
105 coefficient_a = coefficient_x;
106 sign_b = sign_y;
107 exponent_b = exponent_y;
108 coefficient_b = coefficient_y;
109 }
110
111 // exponent difference
112 diff_dec_expon = exponent_a - exponent_b;
113
114 /* get binary coefficients of x and y */
115
116 //--- get number of bits in the coefficients of x and y ---
117
118 tempx.d = (double) coefficient_a;
119 bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
120
121 if (!coefficient_a) {
122 return get_BID64 (sign_b, exponent_b, coefficient_b, rounding_mode,
123 fpsc);
124 }
125 if (diff_dec_expon > MAX_FORMAT_DIGITS) {
126 // normalize a to a 16-digit coefficient
127
128 scale_ca = bid_estimate_decimal_digits[bin_expon_ca];
129 if (coefficient_a >= bid_power10_table_128[scale_ca].w[0])
130 scale_ca++;
131
132 scale_k = 16 - scale_ca;
133
134 coefficient_a *= bid_power10_table_128[scale_k].w[0];
135
136 diff_dec_expon -= scale_k;
137 exponent_a -= scale_k;
138
139 /* get binary coefficients of x and y */
140
141 //--- get number of bits in the coefficients of x and y ---
142 tempx.d = (double) coefficient_a;
143 bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
144
145 if (diff_dec_expon > MAX_FORMAT_DIGITS) {
146 #ifdef BID_SET_STATUS_FLAGS
147 if (coefficient_b) {
148 __set_status_flags (fpsc, BID_INEXACT_EXCEPTION);
149 }
150 #endif
151
152 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
153 #ifndef IEEE_ROUND_NEAREST
154 if (((rounding_mode) & 3) && coefficient_b) // not BID_ROUNDING_TO_NEAREST
155 {
156 switch (rounding_mode) {
157 case BID_ROUNDING_DOWN:
158 if (sign_b) {
159 coefficient_a -= ((((BID_SINT64) sign_a) >> 63) | 1);
160 if (coefficient_a < 1000000000000000ull) {
161 exponent_a--;
162 coefficient_a = 9999999999999999ull;
163 } else if (coefficient_a >= 10000000000000000ull) {
164 exponent_a++;
165 coefficient_a = 1000000000000000ull;
166 }
167 }
168 break;
169 case BID_ROUNDING_UP:
170 if (!sign_b) {
171 coefficient_a += ((((BID_SINT64) sign_a) >> 63) | 1);
172 if (coefficient_a < 1000000000000000ull) {
173 exponent_a--;
174 coefficient_a = 9999999999999999ull;
175 } else if (coefficient_a >= 10000000000000000ull) {
176 exponent_a++;
177 coefficient_a = 1000000000000000ull;
178 }
179 }
180 break;
181 default: // RZ
182 if (sign_a != sign_b) {
183 coefficient_a--;
184 if (coefficient_a < 1000000000000000ull) {
185 exponent_a--;
186 coefficient_a = 9999999999999999ull;
187 }
188 }
189 break;
190 }
191 } else
192 #endif
193 #endif
194 // check special case here
195 if ((coefficient_a == 1000000000000000ull)
196 && (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
197 && (sign_a ^ sign_b)
198 && (coefficient_b > 5000000000000000ull)) {
199 coefficient_a = 9999999999999999ull;
200 exponent_a--;
201 }
202
203 return get_BID64 (sign_a, exponent_a, coefficient_a,
204 rounding_mode, fpsc);
205 }
206 }
207 // test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62
208 if (bin_expon_ca + bid_estimate_bin_expon[diff_dec_expon] < 60) {
209 // coefficient_a*10^(exponent_a-exponent_b)<2^63
210
211 // multiply by 10^(exponent_a-exponent_b)
212 coefficient_a *= bid_power10_table_128[diff_dec_expon].w[0];
213
214 // sign mask
215 sign_b = ((BID_SINT64) sign_b) >> 63;
216 // apply sign to coeff. of b
217 coefficient_b = (coefficient_b + sign_b) ^ sign_b;
218
219 // apply sign to coefficient a
220 sign_a = ((BID_SINT64) sign_a) >> 63;
221 coefficient_a = (coefficient_a + sign_a) ^ sign_a;
222
223 coefficient_a += coefficient_b;
224 // get sign
225 sign_s = ((BID_SINT64) coefficient_a) >> 63;
226 coefficient_a = (coefficient_a + sign_s) ^ sign_s;
227 sign_s &= 0x8000000000000000ull;
228
229 // coefficient_a < 10^16 ?
230 if (coefficient_a < bid_power10_table_128[MAX_FORMAT_DIGITS].w[0]) {
231 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
232 #ifndef IEEE_ROUND_NEAREST
233 if (rounding_mode == BID_ROUNDING_DOWN && (!coefficient_a)
234 && sign_a != sign_b)
235 sign_s = 0x8000000000000000ull;
236 #endif
237 #endif
238 return get_BID64 (sign_s, exponent_b, coefficient_a,
239 rounding_mode, fpsc);
240 }
241 // otherwise rounding is necessary
242
243 // already know coefficient_a<10^19
244 // coefficient_a < 10^17 ?
245 if (coefficient_a < bid_power10_table_128[17].w[0])
246 extra_digits = 1;
247 else if (coefficient_a < bid_power10_table_128[18].w[0])
248 extra_digits = 2;
249 else
250 extra_digits = 3;
251
252 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
253 #ifndef IEEE_ROUND_NEAREST
254 rmode = rounding_mode;
255 if (sign_s && (unsigned) (rmode - 1) < 2)
256 rmode = 3 - rmode;
257 #else
258 rmode = 0;
259 #endif
260 #else
261 rmode = 0;
262 #endif
263 coefficient_a += bid_round_const_table[rmode][extra_digits];
264
265 // get P*(2^M[extra_digits])/10^extra_digits
266 __mul_64x64_to_128 (CT, coefficient_a,
267 bid_reciprocals10_64[extra_digits]);
268
269 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
270 amount = bid_short_recip_scale[extra_digits];
271 C64 = CT.w[1] >> amount;
272
273 } else {
274 // coefficient_a*10^(exponent_a-exponent_b) is large
275 sign_s = sign_a;
276
277 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
278 #ifndef IEEE_ROUND_NEAREST
279 rmode = rounding_mode;
280 if (sign_s && (unsigned) (rmode - 1) < 2)
281 rmode = 3 - rmode;
282 #else
283 rmode = 0;
284 #endif
285 #else
286 rmode = 0;
287 #endif
288
289 // check whether we can take faster path
290 scale_ca = bid_estimate_decimal_digits[bin_expon_ca];
291
292 sign_ab = sign_a ^ sign_b;
293 sign_ab = ((BID_SINT64) sign_ab) >> 63;
294
295 // T1 = 10^(16-diff_dec_expon)
296 T1 = bid_power10_table_128[16 - diff_dec_expon].w[0];
297
298 // get number of digits in coefficient_a
299 //P_ca = bid_power10_table_128[scale_ca].w[0];
300 //P_ca_m1 = bid_power10_table_128[scale_ca-1].w[0];
301 if (coefficient_a >= bid_power10_table_128[scale_ca].w[0]) {
302 scale_ca++;
303 //P_ca_m1 = P_ca;
304 //P_ca = bid_power10_table_128[scale_ca].w[0];
305 }
306
307 scale_k = 16 - scale_ca;
308
309 // apply sign
310 //Ts = (T1 + sign_ab) ^ sign_ab;
311
312 // test range of ca
313 //X = coefficient_a + Ts - P_ca_m1;
314
315 // addition
316 saved_ca = coefficient_a - T1;
317 coefficient_a =
318 (BID_SINT64) saved_ca *(BID_SINT64) bid_power10_table_128[scale_k].w[0];
319 extra_digits = diff_dec_expon - scale_k;
320
321 // apply sign
322 saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
323 // add 10^16 and rounding constant
324 coefficient_b =
325 saved_cb + 10000000000000000ull +
326 bid_round_const_table[rmode][extra_digits];
327
328 // get P*(2^M[extra_digits])/10^extra_digits
329 __mul_64x64_to_128 (CT, coefficient_b,
330 bid_reciprocals10_64[extra_digits]);
331
332 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
333 amount = bid_short_recip_scale[extra_digits];
334 C0_64 = CT.w[1] >> amount;
335
336 // result coefficient
337 C64 = C0_64 + coefficient_a;
338 // filter out difficult (corner) cases
339 // the following test is equivalent to
340 // ( (initial_coefficient_a + Ts) < P_ca &&
341 // (initial_coefficient_a + Ts) > P_ca_m1 ),
342 // which ensures the number of digits in coefficient_a does not change
343 // after adding (the appropriately scaled and rounded) coefficient_b
344 if ((BID_UINT64) (C64 - 1000000000000000ull - 1) >
345 9000000000000000ull - 2) {
346 if (C64 >= 10000000000000000ull) {
347 // result has more than 16 digits
348 if (!scale_k) {
349 // must divide coeff_a by 10
350 saved_ca = saved_ca + T1;
351 __mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull);
352 //reciprocals10_64[1]);
353 coefficient_a = CA.w[1] >> 1;
354 rem_a =
355 saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
356 coefficient_a = coefficient_a - T1;
357
358 saved_cb +=
359 /*90000000000000000 */ +rem_a *
360 bid_power10_table_128[diff_dec_expon].w[0];
361 } else
362 coefficient_a =
363 (BID_SINT64) (saved_ca - T1 -
364 (T1 << 3)) * (BID_SINT64) bid_power10_table_128[scale_k -
365 1].w[0];
366
367 extra_digits++;
368 coefficient_b =
369 saved_cb + 100000000000000000ull +
370 bid_round_const_table[rmode][extra_digits];
371
372 // get P*(2^M[extra_digits])/10^extra_digits
373 __mul_64x64_to_128 (CT, coefficient_b,
374 bid_reciprocals10_64[extra_digits]);
375
376 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
377 amount = bid_short_recip_scale[extra_digits];
378 C0_64 = CT.w[1] >> amount;
379
380 // result coefficient
381 C64 = C0_64 + coefficient_a;
382 } else if (C64 <= 1000000000000000ull) {
383 // less than 16 digits in result
384 coefficient_a =
385 (BID_SINT64) saved_ca *(BID_SINT64) bid_power10_table_128[scale_k +
386 1].w[0];
387 //extra_digits --;
388 exponent_b--;
389 coefficient_b =
390 (saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull +
391 bid_round_const_table[rmode][extra_digits];
392
393 // get P*(2^M[extra_digits])/10^extra_digits
394 __mul_64x64_to_128 (CT_new, coefficient_b,
395 bid_reciprocals10_64[extra_digits]);
396
397 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
398 amount = bid_short_recip_scale[extra_digits];
399 C0_64 = CT_new.w[1] >> amount;
400
401 // result coefficient
402 C64_new = C0_64 + coefficient_a;
403 if (C64_new < 10000000000000000ull) {
404 C64 = C64_new;
405 #ifdef BID_SET_STATUS_FLAGS
406 CT = CT_new;
407 #endif
408 } else
409 exponent_b++;
410 }
411
412 }
413
414 }
415
416 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
417 #ifndef IEEE_ROUND_NEAREST
418 if (rmode == 0) //BID_ROUNDING_TO_NEAREST
419 #endif
420 if (C64 & 1) {
421 // check whether fractional part of initial_P/10^extra_digits
422 // is exactly .5
423 // this is the same as fractional part of
424 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
425
426 // get remainder
427 remainder_h = CT.w[1] << (64 - amount);
428
429 // test whether fractional part is 0
430 if (!remainder_h && (CT.w[0] < bid_reciprocals10_64[extra_digits])) {
431 C64--;
432 }
433 }
434 #endif
435
436 #ifdef BID_SET_STATUS_FLAGS
437 status = BID_INEXACT_EXCEPTION;
438
439 // get remainder
440 remainder_h = CT.w[1] << (64 - amount);
441
442 switch (rmode) {
443 case BID_ROUNDING_TO_NEAREST:
444 case BID_ROUNDING_TIES_AWAY:
445 // test whether fractional part is 0
446 if ((remainder_h == 0x8000000000000000ull)
447 && (CT.w[0] < bid_reciprocals10_64[extra_digits]))
448 status = BID_EXACT_STATUS;
449 break;
450 case BID_ROUNDING_DOWN:
451 case BID_ROUNDING_TO_ZERO:
452 if (!remainder_h && (CT.w[0] < bid_reciprocals10_64[extra_digits]))
453 status = BID_EXACT_STATUS;
454 break;
455 default:
456 // round up
457 __add_carry_out (tmp, carry, CT.w[0],
458 bid_reciprocals10_64[extra_digits]);
459 if ((remainder_h >> (64 - amount)) + carry >=
460 (((BID_UINT64) 1) << amount))
461 status = BID_EXACT_STATUS;
462 break;
463 }
464 __set_status_flags (fpsc, status);
465
466 #endif
467
468 return get_BID64 (sign_s, exponent_b + extra_digits, C64,
469 rounding_mode, fpsc);
470 }
471
472
473 ///////////////////////////////////////////////////////////////////
474 // round 128-bit coefficient and return result in BID64 format
475 // do not worry about midpoint cases
476 //////////////////////////////////////////////////////////////////
477 static BID_UINT64
__bid_simple_round64_sticky(BID_UINT64 sign,int exponent,BID_UINT128 P,int extra_digits,int rounding_mode,unsigned * fpsc)478 __bid_simple_round64_sticky (BID_UINT64 sign, int exponent, BID_UINT128 P,
479 int extra_digits, int rounding_mode,
480 unsigned *fpsc) {
481 BID_UINT128 Q_high, Q_low, C128;
482 BID_UINT64 C64;
483 int amount, rmode;
484
485 rmode = rounding_mode;
486 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
487 #ifndef IEEE_ROUND_NEAREST
488 if (sign && (unsigned) (rmode - 1) < 2)
489 rmode = 3 - rmode;
490 #endif
491 #endif
492 __add_128_64 (P, P, bid_round_const_table[rmode][extra_digits]);
493
494 // get P*(2^M[extra_digits])/10^extra_digits
495 __mul_128x128_full (Q_high, Q_low, P,
496 bid_reciprocals10_128[extra_digits]);
497
498 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
499 amount = bid_recip_scale[extra_digits];
500 __shr_128 (C128, Q_high, amount);
501
502 C64 = __low_64 (C128);
503
504 #ifdef BID_SET_STATUS_FLAGS
505
506 __set_status_flags (fpsc, BID_INEXACT_EXCEPTION);
507
508 #endif
509
510 return get_BID64 (sign, exponent, C64, rounding_mode, fpsc);
511 }
512
513 ///////////////////////////////////////////////////////////////////
514 // round 128-bit coefficient and return result in BID64 format
515 ///////////////////////////////////////////////////////////////////
516 static BID_UINT64
__bid_full_round64(BID_UINT64 sign,int exponent,BID_UINT128 P,int extra_digits,int rounding_mode,unsigned * fpsc)517 __bid_full_round64 (BID_UINT64 sign, int exponent, BID_UINT128 P,
518 int extra_digits, int rounding_mode,
519 unsigned *fpsc) {
520 BID_UINT128 Q_high, Q_low, C128, Stemp;
521 #ifdef BID_SET_STATUS_FLAGS
522 #if DECIMAL_TINY_DETECTION_AFTER_ROUNDING
523 BID_UINT128 PU;
524 #endif
525 #endif
526 BID_UINT64 remainder_h, C64, carry, CY;
527 int amount, amount2, rmode, status = 0;
528
529 if (exponent < 0) {
530 if (exponent >= -16 && (extra_digits + exponent < 0)) {
531 extra_digits = -exponent;
532 #ifdef BID_SET_STATUS_FLAGS
533 #if DECIMAL_TINY_DETECTION_AFTER_ROUNDING
534 if (extra_digits > 0) {
535 rmode = rounding_mode;
536 if (sign && (unsigned) (rmode - 1) < 2)
537 rmode = 3 - rmode;
538 __add_128_128 (PU, P,
539 bid_round_const_table_128[rmode][extra_digits]);
540 if (__unsigned_compare_gt_128
541 (bid_power10_table_128[extra_digits + 15], PU))
542 status = BID_UNDERFLOW_EXCEPTION;
543 }
544 #else
545 status = BID_UNDERFLOW_EXCEPTION;
546 #endif
547 #endif
548 }
549 }
550
551 if (extra_digits > 0) {
552 exponent += extra_digits;
553 rmode = rounding_mode;
554 if (sign && (unsigned) (rmode - 1) < 2)
555 rmode = 3 - rmode;
556 __add_128_128 (P, P, bid_round_const_table_128[rmode][extra_digits]);
557
558 // get P*(2^M[extra_digits])/10^extra_digits
559 __mul_128x128_full (Q_high, Q_low, P,
560 bid_reciprocals10_128[extra_digits]);
561
562 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
563 amount = bid_recip_scale[extra_digits];
564 __shr_128_long (C128, Q_high, amount);
565
566 C64 = __low_64 (C128);
567
568 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
569 #ifndef IEEE_ROUND_NEAREST
570 if (rmode == 0) //BID_ROUNDING_TO_NEAREST
571 #endif
572 if (C64 & 1) {
573 // check whether fractional part of initial_P/10^extra_digits
574 // is exactly .5
575
576 // get remainder
577 amount2 = 64 - amount;
578 remainder_h = 0;
579 remainder_h--;
580 remainder_h >>= amount2;
581 remainder_h = remainder_h & Q_high.w[0];
582
583 if (!remainder_h
584 && (Q_low.w[1] < bid_reciprocals10_128[extra_digits].w[1]
585 || (Q_low.w[1] == bid_reciprocals10_128[extra_digits].w[1]
586 && Q_low.w[0] <
587 bid_reciprocals10_128[extra_digits].w[0]))) {
588 C64--;
589 }
590 }
591 #endif
592
593 #ifdef BID_SET_STATUS_FLAGS
594 status |= BID_INEXACT_EXCEPTION;
595
596 // get remainder
597 remainder_h = Q_high.w[0] << (64 - amount);
598
599 switch (rmode) {
600 case BID_ROUNDING_TO_NEAREST:
601 case BID_ROUNDING_TIES_AWAY:
602 // test whether fractional part is 0
603 if (remainder_h == 0x8000000000000000ull
604 && (Q_low.w[1] < bid_reciprocals10_128[extra_digits].w[1]
605 || (Q_low.w[1] == bid_reciprocals10_128[extra_digits].w[1]
606 && Q_low.w[0] <
607 bid_reciprocals10_128[extra_digits].w[0])))
608 status = BID_EXACT_STATUS;
609 break;
610 case BID_ROUNDING_DOWN:
611 case BID_ROUNDING_TO_ZERO:
612 if (!remainder_h
613 && (Q_low.w[1] < bid_reciprocals10_128[extra_digits].w[1]
614 || (Q_low.w[1] == bid_reciprocals10_128[extra_digits].w[1]
615 && Q_low.w[0] <
616 bid_reciprocals10_128[extra_digits].w[0])))
617 status = BID_EXACT_STATUS;
618 break;
619 default:
620 // round up
621 __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
622 bid_reciprocals10_128[extra_digits].w[0]);
623 __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
624 bid_reciprocals10_128[extra_digits].w[1], CY);
625 if ((remainder_h >> (64 - amount)) + carry >=
626 (((BID_UINT64) 1) << amount))
627 status = BID_EXACT_STATUS;
628 }
629
630 __set_status_flags (fpsc, status);
631
632 #endif
633 } else {
634 C64 = P.w[0];
635 if (!C64) {
636 sign = 0;
637 if (rounding_mode == BID_ROUNDING_DOWN)
638 sign = 0x8000000000000000ull;
639 }
640 }
641 return get_BID64 (sign, exponent, C64, rounding_mode, fpsc);
642 }
643
644 /////////////////////////////////////////////////////////////////////////////////
645 // round 192-bit coefficient (P, remainder_P) and return result in BID64 format
646 // the lowest 64 bits (remainder_P) are used for midpoint checking only
647 ////////////////////////////////////////////////////////////////////////////////
648 static BID_UINT64
__bid_full_round64_remainder(BID_UINT64 sign,int exponent,BID_UINT128 P,int extra_digits,BID_UINT64 remainder_P,int rounding_mode,unsigned * fpsc,unsigned uf_status)649 __bid_full_round64_remainder (BID_UINT64 sign, int exponent, BID_UINT128 P,
650 int extra_digits, BID_UINT64 remainder_P,
651 int rounding_mode, unsigned *fpsc,
652 unsigned uf_status) {
653 BID_UINT128 Q_high, Q_low, C128, Stemp;
654 BID_UINT64 remainder_h, C64, carry, CY;
655 int amount, amount2, rmode, status = uf_status;
656
657 rmode = rounding_mode;
658 if (sign && (unsigned) (rmode - 1) < 2)
659 rmode = 3 - rmode;
660 if (rmode == BID_ROUNDING_UP && remainder_P) {
661 P.w[0]++;
662 if (!P.w[0])
663 P.w[1]++;
664 }
665
666 if (extra_digits) {
667 __add_128_64 (P, P, bid_round_const_table[rmode][extra_digits]);
668
669 // get P*(2^M[extra_digits])/10^extra_digits
670 __mul_128x128_full (Q_high, Q_low, P,
671 bid_reciprocals10_128[extra_digits]);
672
673 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
674 amount = bid_recip_scale[extra_digits];
675 __shr_128 (C128, Q_high, amount);
676
677 C64 = __low_64 (C128);
678
679 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
680 #ifndef IEEE_ROUND_NEAREST
681 if (rmode == 0) //BID_ROUNDING_TO_NEAREST
682 #endif
683 if (!remainder_P && (C64 & 1)) {
684 // check whether fractional part of initial_P/10^extra_digits
685 // is exactly .5
686
687 // get remainder
688 amount2 = 64 - amount;
689 remainder_h = 0;
690 remainder_h--;
691 remainder_h >>= amount2;
692 remainder_h = remainder_h & Q_high.w[0];
693
694 if (!remainder_h
695 && (Q_low.w[1] < bid_reciprocals10_128[extra_digits].w[1]
696 || (Q_low.w[1] == bid_reciprocals10_128[extra_digits].w[1]
697 && Q_low.w[0] <
698 bid_reciprocals10_128[extra_digits].w[0]))) {
699 C64--;
700 }
701 }
702 #endif
703
704 #ifdef BID_SET_STATUS_FLAGS
705 status |= BID_INEXACT_EXCEPTION;
706
707 if (!remainder_P) {
708 // get remainder
709 remainder_h = Q_high.w[0] << (64 - amount);
710
711 switch (rmode) {
712 case BID_ROUNDING_TO_NEAREST:
713 case BID_ROUNDING_TIES_AWAY:
714 // test whether fractional part is 0
715 if (remainder_h == 0x8000000000000000ull
716 && (Q_low.w[1] < bid_reciprocals10_128[extra_digits].w[1]
717 || (Q_low.w[1] == bid_reciprocals10_128[extra_digits].w[1]
718 && Q_low.w[0] <
719 bid_reciprocals10_128[extra_digits].w[0])))
720 status = BID_EXACT_STATUS;
721 break;
722 case BID_ROUNDING_DOWN:
723 case BID_ROUNDING_TO_ZERO:
724 if (!remainder_h
725 && (Q_low.w[1] < bid_reciprocals10_128[extra_digits].w[1]
726 || (Q_low.w[1] == bid_reciprocals10_128[extra_digits].w[1]
727 && Q_low.w[0] <
728 bid_reciprocals10_128[extra_digits].w[0])))
729 status = BID_EXACT_STATUS;
730 break;
731 default:
732 // round up
733 __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
734 bid_reciprocals10_128[extra_digits].w[0]);
735 __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
736 bid_reciprocals10_128[extra_digits].w[1], CY);
737 if ((remainder_h >> (64 - amount)) + carry >=
738 (((BID_UINT64) 1) << amount))
739 status = BID_EXACT_STATUS;
740 }
741 }
742 __set_status_flags (fpsc, status);
743
744 #endif
745 } else {
746 C64 = P.w[0];
747 #ifdef BID_SET_STATUS_FLAGS
748 if (remainder_P) {
749 __set_status_flags (fpsc, uf_status | BID_INEXACT_EXCEPTION);
750 }
751 #endif
752 }
753
754 return get_BID64 (sign, exponent + extra_digits, C64, rounding_mode,
755 fpsc);
756 }
757
758
759 ///////////////////////////////////////////////////////////////////
760 // get P/10^extra_digits
761 // result fits in 64 bits
762 ///////////////////////////////////////////////////////////////////
763 __BID_INLINE__ BID_UINT64
__truncate(BID_UINT128 P,int extra_digits)764 __truncate (BID_UINT128 P, int extra_digits)
765 // extra_digits <= 16
766 {
767 BID_UINT128 Q_high, Q_low, C128;
768 BID_UINT64 C64;
769 int amount;
770
771 // get P*(2^M[extra_digits])/10^extra_digits
772 __mul_128x128_full (Q_high, Q_low, P,
773 bid_reciprocals10_128[extra_digits]);
774
775 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
776 amount = bid_recip_scale[extra_digits];
777 __shr_128 (C128, Q_high, amount);
778
779 C64 = __low_64 (C128);
780
781 return C64;
782 }
783
784
785 ///////////////////////////////////////////////////////////////////
786 // return number of decimal digits in 128-bit value X
787 ///////////////////////////////////////////////////////////////////
788 __BID_INLINE__ int
__get_dec_digits64(BID_UINT128 X)789 __get_dec_digits64 (BID_UINT128 X) {
790 int_double tempx;
791 int digits_x, bin_expon_cx;
792
793 if (!X.w[1]) {
794 if(!X.w[0]) return 0;
795 //--- get number of bits in the coefficients of x and y ---
796 tempx.d = (double) X.w[0];
797 bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
798 // get number of decimal digits in the coeff_x
799 digits_x = bid_estimate_decimal_digits[bin_expon_cx];
800 if (X.w[0] >= bid_power10_table_128[digits_x].w[0])
801 digits_x++;
802 return digits_x;
803 }
804 tempx.d = (double) X.w[1];
805 bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
806 // get number of decimal digits in the coeff_x
807 digits_x = bid_estimate_decimal_digits[bin_expon_cx + 64];
808 if (__unsigned_compare_ge_128 (X, bid_power10_table_128[digits_x]))
809 digits_x++;
810
811 return digits_x;
812 }
813
814
815 ////////////////////////////////////////////////////////////////////////////////
816 //
817 // add 64-bit coefficient to 128-bit coefficient, return result in BID64 format
818 //
819 ////////////////////////////////////////////////////////////////////////////////
820 __BID_INLINE__ BID_UINT64
bid_get_add128(BID_UINT64 sign_x,int exponent_x,BID_UINT64 coefficient_x,BID_UINT64 sign_y,int final_exponent_y,BID_UINT128 CY,int extra_digits,int rounding_mode,unsigned * fpsc)821 bid_get_add128 (BID_UINT64 sign_x, int exponent_x, BID_UINT64 coefficient_x,
822 BID_UINT64 sign_y, int final_exponent_y, BID_UINT128 CY,
823 int extra_digits, int rounding_mode, unsigned *fpsc) {
824 BID_UINT128 CY_L, CX, FS, F, CT, ST, T2;
825 BID_UINT64 CYh, CY0L, T, S, coefficient_y, remainder_y;
826 BID_SINT64 D = 0;
827 int_double tempx;
828 int diff_dec_expon, extra_digits2, exponent_y, status;
829 int extra_dx, diff_dec2, bin_expon_cx, digits_x, rmode;
830
831 // CY has more than 16 decimal digits
832
833 exponent_y = final_exponent_y - extra_digits;
834
835 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
836 rounding_mode = 0;
837 #endif
838 #ifdef IEEE_ROUND_NEAREST
839 rounding_mode = 0;
840 #endif
841
842 if (exponent_x > exponent_y) {
843 // normalize x
844 //--- get number of bits in the coefficients of x and y ---
845 tempx.d = (double) coefficient_x;
846 bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
847 // get number of decimal digits in the coeff_x
848 digits_x = bid_estimate_decimal_digits[bin_expon_cx];
849 if (coefficient_x >= bid_power10_table_128[digits_x].w[0])
850 digits_x++;
851
852 extra_dx = 16 - digits_x;
853 coefficient_x *= bid_power10_table_128[extra_dx].w[0];
854 if ((sign_x ^ sign_y) && (coefficient_x == 1000000000000000ull)) {
855 extra_dx++;
856 coefficient_x = 10000000000000000ull;
857 }
858 exponent_x -= extra_dx;
859
860 if (exponent_x > exponent_y) {
861
862 // exponent_x > exponent_y
863 diff_dec_expon = exponent_x - exponent_y;
864
865 if (exponent_x <= final_exponent_y + 1) {
866 __mul_64x64_to_128 (CX, coefficient_x,
867 bid_power10_table_128[diff_dec_expon].w[0]);
868
869 if (sign_x == sign_y) {
870 __add_128_128 (CT, CY, CX);
871 if ((exponent_x >
872 final_exponent_y) /*&& (final_exponent_y>0) */ )
873 extra_digits++;
874 if (__unsigned_compare_ge_128
875 (CT, bid_power10_table_128[16 + extra_digits]))
876 extra_digits++;
877 } else {
878 __sub_128_128 (CT, CY, CX);
879 if (((BID_SINT64) CT.w[1]) < 0) {
880 CT.w[0] = 0 - CT.w[0];
881 CT.w[1] = 0 - CT.w[1];
882 if (CT.w[0])
883 CT.w[1]--;
884 sign_y = sign_x;
885 } else if (!(CT.w[1] | CT.w[0])) {
886 sign_y =
887 (rounding_mode !=
888 BID_ROUNDING_DOWN) ? 0 : 0x8000000000000000ull;
889 }
890 if ((exponent_x + 1 >=
891 final_exponent_y) /*&& (final_exponent_y>=0) */ ) {
892 extra_digits = __get_dec_digits64 (CT) - 16;
893 if (extra_digits <= 0) {
894 if (!CT.w[0] && rounding_mode == BID_ROUNDING_DOWN)
895 sign_y = 0x8000000000000000ull;
896 return get_BID64 (sign_y, exponent_y, CT.w[0],
897 rounding_mode, fpsc);
898 }
899 } else
900 if (__unsigned_compare_gt_128
901 (bid_power10_table_128[15 + extra_digits], CT))
902 extra_digits--;
903 }
904
905 return __bid_full_round64 (sign_y, exponent_y, CT, extra_digits,
906 rounding_mode, fpsc);
907 }
908 // diff_dec2+extra_digits is the number of digits to eliminate from
909 // argument CY
910 diff_dec2 = exponent_x - final_exponent_y;
911
912 if (diff_dec2 >= 17) {
913 #ifndef IEEE_ROUND_NEAREST
914 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
915 if ((rounding_mode) & 3) {
916 switch (rounding_mode) {
917 case BID_ROUNDING_UP:
918 if (!sign_y) {
919 D = ((BID_SINT64) (sign_x ^ sign_y)) >> 63;
920 D = D + D + 1;
921 coefficient_x += D;
922 }
923 break;
924 case BID_ROUNDING_DOWN:
925 if (sign_y) {
926 D = ((BID_SINT64) (sign_x ^ sign_y)) >> 63;
927 D = D + D + 1;
928 coefficient_x += D;
929 }
930 break;
931 case BID_ROUNDING_TO_ZERO:
932 if (sign_y != sign_x) {
933 D = 0 - 1;
934 coefficient_x += D;
935 }
936 break;
937 }
938 if (coefficient_x < 1000000000000000ull) {
939 coefficient_x -= D;
940 coefficient_x =
941 D + (coefficient_x << 1) + (coefficient_x << 3);
942 exponent_x--;
943 }
944 }
945 #endif
946 #endif
947 #ifdef BID_SET_STATUS_FLAGS
948 if (CY.w[1] | CY.w[0])
949 __set_status_flags (fpsc, BID_INEXACT_EXCEPTION);
950 #endif
951 return get_BID64 (sign_x, exponent_x, coefficient_x,
952 rounding_mode, fpsc);
953 }
954 // here exponent_x <= 16+final_exponent_y
955
956 // truncate CY to 16 dec. digits
957 CYh = __truncate (CY, extra_digits);
958
959 // get remainder
960 T = bid_power10_table_128[extra_digits].w[0];
961 __mul_64x64_to_64 (CY0L, CYh, T);
962
963 remainder_y = CY.w[0] - CY0L;
964
965 // align coeff_x, CYh
966 __mul_64x64_to_128 (CX, coefficient_x,
967 bid_power10_table_128[diff_dec2].w[0]);
968
969 if (sign_x == sign_y) {
970 __add_128_64 (CT, CX, CYh);
971 if (__unsigned_compare_ge_128
972 (CT, bid_power10_table_128[16 + diff_dec2]))
973 diff_dec2++;
974 } else {
975 if (remainder_y)
976 CYh++;
977 __sub_128_64 (CT, CX, CYh);
978 if (__unsigned_compare_gt_128
979 (bid_power10_table_128[15 + diff_dec2], CT))
980 diff_dec2--;
981 }
982
983 return __bid_full_round64_remainder (sign_x, final_exponent_y, CT,
984 diff_dec2, remainder_y,
985 rounding_mode, fpsc, 0);
986 }
987 }
988 // Here (exponent_x <= exponent_y)
989 {
990 diff_dec_expon = exponent_y - exponent_x;
991
992 if (diff_dec_expon > MAX_FORMAT_DIGITS) {
993 rmode = rounding_mode;
994
995 if ((sign_x ^ sign_y)) {
996 if (!CY.w[0])
997 CY.w[1]--;
998 CY.w[0]--;
999 if (__unsigned_compare_gt_128
1000 (bid_power10_table_128[15 + extra_digits], CY)) {
1001 if (rmode & 3) {
1002 extra_digits--;
1003 final_exponent_y--;
1004 } else {
1005 CY.w[0] = 1000000000000000ull;
1006 CY.w[1] = 0;
1007 extra_digits = 0;
1008 }
1009 }
1010 }
1011 __scale128_10 (CY, CY);
1012 extra_digits++;
1013 CY.w[0] |= 1;
1014
1015 return __bid_simple_round64_sticky (sign_y, final_exponent_y, CY,
1016 extra_digits, rmode, fpsc);
1017 }
1018 // apply sign to coeff_x
1019 sign_x ^= sign_y;
1020 sign_x = ((BID_SINT64) sign_x) >> 63;
1021 CX.w[0] = (coefficient_x + sign_x) ^ sign_x;
1022 CX.w[1] = sign_x;
1023
1024 // check whether CY (rounded to 16 digits) and CX have
1025 // any digits in the same position
1026 diff_dec2 = final_exponent_y - exponent_x;
1027
1028 if (diff_dec2 <= 17) {
1029 // align CY to 10^ex
1030 S = bid_power10_table_128[diff_dec_expon].w[0];
1031 __mul_64x128_short (CY_L, S, CY);
1032
1033 __add_128_128 (ST, CY_L, CX);
1034 extra_digits2 = __get_dec_digits64 (ST) - 16;
1035 return __bid_full_round64 (sign_y, exponent_x, ST, extra_digits2,
1036 rounding_mode, fpsc);
1037 }
1038 // truncate CY to 16 dec. digits
1039 CYh = __truncate (CY, extra_digits);
1040
1041 // get remainder
1042 T = bid_power10_table_128[extra_digits].w[0];
1043 __mul_64x64_to_64 (CY0L, CYh, T);
1044
1045 coefficient_y = CY.w[0] - CY0L;
1046 // add rounding constant
1047 rmode = rounding_mode;
1048 if (sign_y && (unsigned) (rmode - 1) < 2)
1049 rmode = 3 - rmode;
1050 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1051 #ifndef IEEE_ROUND_NEAREST
1052 if (!(rmode & 3)) //BID_ROUNDING_TO_NEAREST
1053 #endif
1054 #endif
1055 {
1056 coefficient_y += bid_round_const_table[rmode][extra_digits];
1057 }
1058 // align coefficient_y, coefficient_x
1059 S = bid_power10_table_128[diff_dec_expon].w[0];
1060 __mul_64x64_to_128 (F, coefficient_y, S);
1061
1062 // fraction
1063 __add_128_128 (FS, F, CX);
1064
1065 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1066 #ifndef IEEE_ROUND_NEAREST
1067 if (rmode == 0) //BID_ROUNDING_TO_NEAREST
1068 #endif
1069 {
1070 // rounding code, here RN_EVEN
1071 // 10^(extra_digits+diff_dec_expon)
1072 T2 = bid_power10_table_128[diff_dec_expon + extra_digits];
1073 if (__unsigned_compare_gt_128 (FS, T2)
1074 || ((CYh & 1) && __test_equal_128 (FS, T2))) {
1075 CYh++;
1076 __sub_128_128 (FS, FS, T2);
1077 }
1078 }
1079 #endif
1080 #ifndef IEEE_ROUND_NEAREST
1081 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1082 if (rmode == 4) //BID_ROUNDING_TO_NEAREST
1083 #endif
1084 {
1085 // rounding code, here RN_AWAY
1086 // 10^(extra_digits+diff_dec_expon)
1087 T2 = bid_power10_table_128[diff_dec_expon + extra_digits];
1088 if (__unsigned_compare_ge_128 (FS, T2)) {
1089 CYh++;
1090 __sub_128_128 (FS, FS, T2);
1091 }
1092 }
1093 #endif
1094 #ifndef IEEE_ROUND_NEAREST
1095 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1096 switch (rmode) {
1097 case BID_ROUNDING_DOWN:
1098 case BID_ROUNDING_TO_ZERO:
1099 if ((BID_SINT64) FS.w[1] < 0) {
1100 CYh--;
1101 if (CYh < 1000000000000000ull) {
1102 CYh = 9999999999999999ull;
1103 final_exponent_y--;
1104 }
1105 } else {
1106 T2 = bid_power10_table_128[diff_dec_expon + extra_digits];
1107 if (__unsigned_compare_ge_128 (FS, T2)) {
1108 CYh++;
1109 __sub_128_128 (FS, FS, T2);
1110 }
1111 }
1112 break;
1113 case BID_ROUNDING_UP:
1114 if ((BID_SINT64) FS.w[1] < 0)
1115 break;
1116 T2 = bid_power10_table_128[diff_dec_expon + extra_digits];
1117 if (__unsigned_compare_gt_128 (FS, T2)) {
1118 CYh += 2;
1119 __sub_128_128 (FS, FS, T2);
1120 } else if ((FS.w[1] == T2.w[1]) && (FS.w[0] == T2.w[0])) {
1121 CYh++;
1122 FS.w[1] = FS.w[0] = 0;
1123 } else if (FS.w[1] | FS.w[0])
1124 CYh++;
1125 break;
1126 }
1127 #endif
1128 #endif
1129
1130 #ifdef BID_SET_STATUS_FLAGS
1131 status = BID_INEXACT_EXCEPTION;
1132 #ifndef IEEE_ROUND_NEAREST
1133 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1134 if (!(rmode & 3))
1135 #endif
1136 #endif
1137 {
1138 // RN modes
1139 if ((FS.w[1] ==
1140 bid_round_const_table_128[0][diff_dec_expon + extra_digits].w[1])
1141 && (FS.w[0] ==
1142 bid_round_const_table_128[0][diff_dec_expon +
1143 extra_digits].w[0]))
1144 status = BID_EXACT_STATUS;
1145 }
1146 #ifndef IEEE_ROUND_NEAREST
1147 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1148 else if (!FS.w[1] && !FS.w[0])
1149 status = BID_EXACT_STATUS;
1150 #endif
1151 #endif
1152
1153 __set_status_flags (fpsc, status);
1154 #endif
1155
1156 return get_BID64 (sign_y, final_exponent_y, CYh, rounding_mode,
1157 fpsc);
1158 }
1159
1160 }
1161
1162 //////////////////////////////////////////////////////////////////////////
1163 //
1164 // If coefficient_z is less than 16 digits long, normalize to 16 digits
1165 //
1166 /////////////////////////////////////////////////////////////////////////
1167 static BID_UINT64
BID_normalize(BID_UINT64 sign_z,int exponent_z,BID_UINT64 coefficient_z,BID_UINT64 round_dir,int round_flag,int rounding_mode,unsigned * fpsc)1168 BID_normalize (BID_UINT64 sign_z, int exponent_z,
1169 BID_UINT64 coefficient_z, BID_UINT64 round_dir, int round_flag,
1170 int rounding_mode, unsigned *fpsc) {
1171 BID_SINT64 D;
1172 int_double tempx;
1173 int digits_z, bin_expon, scale, rmode;
1174
1175 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1176 #ifndef IEEE_ROUND_NEAREST
1177 rmode = rounding_mode;
1178 if (sign_z && (unsigned) (rmode - 1) < 2)
1179 rmode = 3 - rmode;
1180 #else
1181 if (coefficient_z >= bid_power10_table_128[15].w[0])
1182 return z;
1183 #endif
1184 #endif
1185
1186 //--- get number of bits in the coefficients of x and y ---
1187 tempx.d = (double) coefficient_z;
1188 bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
1189 // get number of decimal digits in the coeff_x
1190 digits_z = bid_estimate_decimal_digits[bin_expon];
1191 if (coefficient_z >= bid_power10_table_128[digits_z].w[0])
1192 digits_z++;
1193
1194 scale = 16 - digits_z;
1195 exponent_z -= scale;
1196 if (exponent_z < 0) {
1197 scale += exponent_z;
1198 exponent_z = 0;
1199 }
1200 coefficient_z *= bid_power10_table_128[scale].w[0];
1201
1202 #ifdef BID_SET_STATUS_FLAGS
1203 if (round_flag) {
1204 __set_status_flags (fpsc, BID_INEXACT_EXCEPTION);
1205 if (coefficient_z < 1000000000000000ull)
1206 __set_status_flags (fpsc, BID_UNDERFLOW_EXCEPTION);
1207 else if ((coefficient_z == 1000000000000000ull) && !exponent_z
1208 && ((BID_SINT64) (round_dir ^ sign_z) < 0) && round_flag
1209 #if DECIMAL_TINY_DETECTION_AFTER_ROUNDING
1210 && (rmode == BID_ROUNDING_DOWN || rmode == BID_ROUNDING_TO_ZERO)
1211 #endif
1212 )
1213 __set_status_flags (fpsc, BID_UNDERFLOW_EXCEPTION);
1214 }
1215 #endif
1216
1217 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1218 #ifndef IEEE_ROUND_NEAREST
1219 if (round_flag && (rmode & 3)) {
1220 D = round_dir ^ sign_z;
1221
1222 if (rmode == BID_ROUNDING_UP) {
1223 if (D >= 0)
1224 coefficient_z++;
1225 } else {
1226 if (D < 0)
1227 coefficient_z--;
1228 if (coefficient_z < 1000000000000000ull && exponent_z) {
1229 coefficient_z = 9999999999999999ull;
1230 exponent_z--;
1231 }
1232 }
1233 }
1234 #endif
1235 #endif
1236
1237 return get_BID64 (sign_z, exponent_z, coefficient_z, rounding_mode,
1238 fpsc);
1239 }
1240
1241
1242 //////////////////////////////////////////////////////////////////////////
1243 //
1244 // 0*10^ey + cz*10^ez, ey<ez
1245 //
1246 //////////////////////////////////////////////////////////////////////////
1247
1248 __BID_INLINE__ BID_UINT64
add_zero64(int exponent_y,BID_UINT64 sign_z,int exponent_z,BID_UINT64 coefficient_z,unsigned * prounding_mode,unsigned * fpsc)1249 add_zero64 (int exponent_y, BID_UINT64 sign_z, int exponent_z,
1250 BID_UINT64 coefficient_z, unsigned *prounding_mode,
1251 unsigned *fpsc) {
1252 int_double tempx;
1253 int bin_expon, scale_k, scale_cz;
1254 int diff_expon;
1255
1256 diff_expon = exponent_z - exponent_y;
1257
1258 tempx.d = (double) coefficient_z;
1259 bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
1260 scale_cz = bid_estimate_decimal_digits[bin_expon];
1261 if (coefficient_z >= bid_power10_table_128[scale_cz].w[0])
1262 scale_cz++;
1263
1264 scale_k = 16 - scale_cz;
1265 if (diff_expon < scale_k)
1266 scale_k = diff_expon;
1267 coefficient_z *= bid_power10_table_128[scale_k].w[0];
1268
1269 return get_BID64 (sign_z, exponent_z - scale_k, coefficient_z,
1270 *prounding_mode, fpsc);
1271 }
1272 #endif
1273