1 /* Copyright (C) 2007-2024 Free Software Foundation, Inc.
3 This file is part of GCC.
5 GCC is free software; you can redistribute it and/or modify it under
6 the terms of the GNU General Public License as published by the Free
7 Software Foundation; either version 3, or (at your option) any later
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 Under Section 7 of GPL version 3, you are granted additional
16 permissions described in the GCC Runtime Library Exception, version
17 3.1, as published by the Free Software Foundation.
19 You should have received a copy of the GNU General Public License and
20 a copy of the GCC Runtime Library Exception along with this program;
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22 <http://www.gnu.org/licenses/>. */
24 /*****************************************************************************
26 *****************************************************************************
28 * Algorithm description:
30 * if(exponent_a < exponent_b)
32 * diff_expon = exponent_a - exponent_b
35 * if(coefficient_a*10^diff_expon guaranteed below 2^62)
36 * S = sign_a*coefficient_a*10^diff_expon + sign_b*coefficient_b
38 * return get_BID64(sign(S),exponent_b,|S|)
40 * determine number of extra digits in S (1, 2, or 3)
41 * return rounded result
42 * else // large exponent difference
43 * if(number_digits(coefficient_a*10^diff_expon) +/- 10^16)
44 * guaranteed the same as
45 * number_digits(coefficient_a*10^diff_expon) )
46 * S = normalize(coefficient_a + (sign_a^sign_b)*10^(16-diff_expon))
47 * corr = 10^16 + (sign_a^sign_b)*coefficient_b
48 * corr*10^exponent_b is rounded so it aligns with S*10^exponent_S
49 * return get_BID64(sign_a,exponent(S),S+rounded(corr))
51 * add sign_a*coefficient_a*10^diff_expon, sign_b*coefficient_b
52 * in 128-bit integer arithmetic, then round to 16 decimal digits
55 ****************************************************************************/
57 #include "bid_internal.h"
59 #if DECIMAL_CALL_BY_REFERENCE
60 void bid64_add (UINT64
* pres
, UINT64
* px
,
62 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
65 UINT64
bid64_add (UINT64 x
,
66 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
67 _EXC_MASKS_PARAM _EXC_INFO_PARAM
);
70 #if DECIMAL_CALL_BY_REFERENCE
73 bid64_sub (UINT64
* pres
, UINT64
* px
,
75 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
78 #if !DECIMAL_GLOBAL_ROUNDING
79 _IDEC_round rnd_mode
= *prnd_mode
;
81 // check if y is not NaN
82 if (((y
& NAN_MASK64
) != NAN_MASK64
))
83 y
^= 0x8000000000000000ull
;
85 &y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
92 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
93 _EXC_MASKS_PARAM _EXC_INFO_PARAM
) {
94 // check if y is not NaN
95 if (((y
& NAN_MASK64
) != NAN_MASK64
))
96 y
^= 0x8000000000000000ull
;
99 y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
106 #if DECIMAL_CALL_BY_REFERENCE
109 bid64_add (UINT64
* pres
, UINT64
* px
,
111 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
118 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
119 _EXC_MASKS_PARAM _EXC_INFO_PARAM
) {
122 UINT128 CA
, CT
, CT_new
;
123 UINT64 sign_x
, sign_y
, coefficient_x
, coefficient_y
, C64_new
;
124 UINT64 valid_x
, valid_y
;
126 UINT64 sign_a
, sign_b
, coefficient_a
, coefficient_b
, sign_s
, sign_ab
,
128 UINT64 saved_ca
, saved_cb
, C0_64
, C64
, remainder_h
, T1
, carry
, tmp
;
130 int exponent_x
, exponent_y
, exponent_a
, exponent_b
, diff_dec_expon
;
131 int bin_expon_ca
, extra_digits
, amount
, scale_k
, scale_ca
;
132 unsigned rmode
, status
;
134 #if DECIMAL_CALL_BY_REFERENCE
135 #if !DECIMAL_GLOBAL_ROUNDING
136 _IDEC_round rnd_mode
= *prnd_mode
;
142 valid_x
= unpack_BID64 (&sign_x
, &exponent_x
, &coefficient_x
, x
);
143 valid_y
= unpack_BID64 (&sign_y
, &exponent_y
, &coefficient_y
, y
);
145 // unpack arguments, check for NaN or Infinity
150 if ((x
& NAN_MASK64
) == NAN_MASK64
) {
151 #ifdef SET_STATUS_FLAGS
152 if (((x
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
153 || ((y
& SNAN_MASK64
) == SNAN_MASK64
))
154 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
156 res
= coefficient_x
& QUIET_MASK64
;
160 if ((x
& INFINITY_MASK64
) == INFINITY_MASK64
) {
162 if (((y
& NAN_MASK64
) == INFINITY_MASK64
)) {
163 if (sign_x
== (y
& 0x8000000000000000ull
)) {
169 #ifdef SET_STATUS_FLAGS
170 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
177 if (((y
& NAN_MASK64
) == NAN_MASK64
)) {
178 res
= coefficient_y
& QUIET_MASK64
;
179 #ifdef SET_STATUS_FLAGS
180 if (((y
& SNAN_MASK64
) == SNAN_MASK64
))
181 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
185 // otherwise return +/-Inf
193 if (((y
& INFINITY_MASK64
) != INFINITY_MASK64
) && coefficient_y
) {
194 if (exponent_y
<= exponent_x
) {
204 if (((y
& INFINITY_MASK64
) == INFINITY_MASK64
)) {
205 #ifdef SET_STATUS_FLAGS
206 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
207 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
209 res
= coefficient_y
& QUIET_MASK64
;
213 if (!coefficient_x
) { // x==0
214 if (exponent_x
<= exponent_y
)
215 res
= ((UINT64
) exponent_x
) << 53;
217 res
= ((UINT64
) exponent_y
) << 53;
218 if (sign_x
== sign_y
)
220 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
221 #ifndef IEEE_ROUND_NEAREST
222 if (rnd_mode
== ROUNDING_DOWN
&& sign_x
!= sign_y
)
223 res
|= 0x8000000000000000ull
;
227 } else if (exponent_y
>= exponent_x
) {
232 // sort arguments by exponent
233 if (exponent_x
< exponent_y
) {
235 exponent_a
= exponent_y
;
236 coefficient_a
= coefficient_y
;
238 exponent_b
= exponent_x
;
239 coefficient_b
= coefficient_x
;
242 exponent_a
= exponent_x
;
243 coefficient_a
= coefficient_x
;
245 exponent_b
= exponent_y
;
246 coefficient_b
= coefficient_y
;
249 // exponent difference
250 diff_dec_expon
= exponent_a
- exponent_b
;
252 /* get binary coefficients of x and y */
254 //--- get number of bits in the coefficients of x and y ---
256 // version 2 (original)
257 tempx
.d
= (double) coefficient_a
;
258 bin_expon_ca
= ((tempx
.i
& MASK_BINARY_EXPONENT
) >> 52) - 0x3ff;
260 if (diff_dec_expon
> MAX_FORMAT_DIGITS
) {
261 // normalize a to a 16-digit coefficient
263 scale_ca
= estimate_decimal_digits
[bin_expon_ca
];
264 if (coefficient_a
>= power10_table_128
[scale_ca
].w
[0])
267 scale_k
= 16 - scale_ca
;
269 coefficient_a
*= power10_table_128
[scale_k
].w
[0];
271 diff_dec_expon
-= scale_k
;
272 exponent_a
-= scale_k
;
274 /* get binary coefficients of x and y */
276 //--- get number of bits in the coefficients of x and y ---
277 tempx
.d
= (double) coefficient_a
;
278 bin_expon_ca
= ((tempx
.i
& MASK_BINARY_EXPONENT
) >> 52) - 0x3ff;
280 if (diff_dec_expon
> MAX_FORMAT_DIGITS
) {
281 #ifdef SET_STATUS_FLAGS
283 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
287 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
288 #ifndef IEEE_ROUND_NEAREST
289 if (((rnd_mode
) & 3) && coefficient_b
) // not ROUNDING_TO_NEAREST
294 coefficient_a
-= ((((SINT64
) sign_a
) >> 63) | 1);
295 if (coefficient_a
< 1000000000000000ull) {
297 coefficient_a
= 9999999999999999ull;
298 } else if (coefficient_a
>= 10000000000000000ull) {
300 coefficient_a
= 1000000000000000ull;
306 coefficient_a
+= ((((SINT64
) sign_a
) >> 63) | 1);
307 if (coefficient_a
< 1000000000000000ull) {
309 coefficient_a
= 9999999999999999ull;
310 } else if (coefficient_a
>= 10000000000000000ull) {
312 coefficient_a
= 1000000000000000ull;
317 if (sign_a
!= sign_b
) {
319 if (coefficient_a
< 1000000000000000ull) {
321 coefficient_a
= 9999999999999999ull;
329 // check special case here
330 if ((coefficient_a
== 1000000000000000ull)
331 && (diff_dec_expon
== MAX_FORMAT_DIGITS
+ 1)
333 && (coefficient_b
> 5000000000000000ull)) {
334 coefficient_a
= 9999999999999999ull;
339 fast_get_BID64_check_OF (sign_a
, exponent_a
, coefficient_a
,
344 // test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62
345 if (bin_expon_ca
+ estimate_bin_expon
[diff_dec_expon
] < 60) {
346 // coefficient_a*10^(exponent_a-exponent_b)<2^63
348 // multiply by 10^(exponent_a-exponent_b)
349 coefficient_a
*= power10_table_128
[diff_dec_expon
].w
[0];
352 sign_b
= ((SINT64
) sign_b
) >> 63;
353 // apply sign to coeff. of b
354 coefficient_b
= (coefficient_b
+ sign_b
) ^ sign_b
;
356 // apply sign to coefficient a
357 sign_a
= ((SINT64
) sign_a
) >> 63;
358 coefficient_a
= (coefficient_a
+ sign_a
) ^ sign_a
;
360 coefficient_a
+= coefficient_b
;
362 sign_s
= ((SINT64
) coefficient_a
) >> 63;
363 coefficient_a
= (coefficient_a
+ sign_s
) ^ sign_s
;
364 sign_s
&= 0x8000000000000000ull
;
366 // coefficient_a < 10^16 ?
367 if (coefficient_a
< power10_table_128
[MAX_FORMAT_DIGITS
].w
[0]) {
368 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
369 #ifndef IEEE_ROUND_NEAREST
370 if (rnd_mode
== ROUNDING_DOWN
&& (!coefficient_a
)
372 sign_s
= 0x8000000000000000ull
;
375 res
= very_fast_get_BID64 (sign_s
, exponent_b
, coefficient_a
);
378 // otherwise rounding is necessary
380 // already know coefficient_a<10^19
381 // coefficient_a < 10^17 ?
382 if (coefficient_a
< power10_table_128
[17].w
[0])
384 else if (coefficient_a
< power10_table_128
[18].w
[0])
389 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
390 #ifndef IEEE_ROUND_NEAREST
392 if (sign_s
&& (unsigned) (rmode
- 1) < 2)
400 coefficient_a
+= round_const_table
[rmode
][extra_digits
];
402 // get P*(2^M[extra_digits])/10^extra_digits
403 __mul_64x64_to_128 (CT
, coefficient_a
,
404 reciprocals10_64
[extra_digits
]);
406 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
407 amount
= short_recip_scale
[extra_digits
];
408 C64
= CT
.w
[1] >> amount
;
411 // coefficient_a*10^(exponent_a-exponent_b) is large
414 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
415 #ifndef IEEE_ROUND_NEAREST
417 if (sign_s
&& (unsigned) (rmode
- 1) < 2)
426 // check whether we can take faster path
427 scale_ca
= estimate_decimal_digits
[bin_expon_ca
];
429 sign_ab
= sign_a
^ sign_b
;
430 sign_ab
= ((SINT64
) sign_ab
) >> 63;
432 // T1 = 10^(16-diff_dec_expon)
433 T1
= power10_table_128
[16 - diff_dec_expon
].w
[0];
435 // get number of digits in coefficient_a
436 if (coefficient_a
>= power10_table_128
[scale_ca
].w
[0]) {
440 scale_k
= 16 - scale_ca
;
443 saved_ca
= coefficient_a
- T1
;
445 (SINT64
) saved_ca
*(SINT64
) power10_table_128
[scale_k
].w
[0];
446 extra_digits
= diff_dec_expon
- scale_k
;
449 saved_cb
= (coefficient_b
+ sign_ab
) ^ sign_ab
;
450 // add 10^16 and rounding constant
452 saved_cb
+ 10000000000000000ull +
453 round_const_table
[rmode
][extra_digits
];
455 // get P*(2^M[extra_digits])/10^extra_digits
456 __mul_64x64_to_128 (CT
, coefficient_b
,
457 reciprocals10_64
[extra_digits
]);
459 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
460 amount
= short_recip_scale
[extra_digits
];
461 C0_64
= CT
.w
[1] >> amount
;
463 // result coefficient
464 C64
= C0_64
+ coefficient_a
;
465 // filter out difficult (corner) cases
466 // this test ensures the number of digits in coefficient_a does not change
467 // after adding (the appropriately scaled and rounded) coefficient_b
468 if ((UINT64
) (C64
- 1000000000000000ull - 1) >
469 9000000000000000ull - 2) {
470 if (C64
>= 10000000000000000ull) {
471 // result has more than 16 digits
473 // must divide coeff_a by 10
474 saved_ca
= saved_ca
+ T1
;
475 __mul_64x64_to_128 (CA
, saved_ca
, 0x3333333333333334ull
);
476 //reciprocals10_64[1]);
477 coefficient_a
= CA
.w
[1] >> 1;
479 saved_ca
- (coefficient_a
<< 3) - (coefficient_a
<< 1);
480 coefficient_a
= coefficient_a
- T1
;
482 saved_cb
+= rem_a
* power10_table_128
[diff_dec_expon
].w
[0];
485 (SINT64
) (saved_ca
- T1
-
486 (T1
<< 3)) * (SINT64
) power10_table_128
[scale_k
-
491 saved_cb
+ 100000000000000000ull +
492 round_const_table
[rmode
][extra_digits
];
494 // get P*(2^M[extra_digits])/10^extra_digits
495 __mul_64x64_to_128 (CT
, coefficient_b
,
496 reciprocals10_64
[extra_digits
]);
498 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
499 amount
= short_recip_scale
[extra_digits
];
500 C0_64
= CT
.w
[1] >> amount
;
502 // result coefficient
503 C64
= C0_64
+ coefficient_a
;
504 } else if (C64
<= 1000000000000000ull) {
505 // less than 16 digits in result
507 (SINT64
) saved_ca
*(SINT64
) power10_table_128
[scale_k
+
512 (saved_cb
<< 3) + (saved_cb
<< 1) + 100000000000000000ull +
513 round_const_table
[rmode
][extra_digits
];
515 // get P*(2^M[extra_digits])/10^extra_digits
516 __mul_64x64_to_128 (CT_new
, coefficient_b
,
517 reciprocals10_64
[extra_digits
]);
519 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
520 amount
= short_recip_scale
[extra_digits
];
521 C0_64
= CT_new
.w
[1] >> amount
;
523 // result coefficient
524 C64_new
= C0_64
+ coefficient_a
;
525 if (C64_new
< 10000000000000000ull) {
527 #ifdef SET_STATUS_FLAGS
538 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
539 #ifndef IEEE_ROUND_NEAREST
540 if (rmode
== 0) //ROUNDING_TO_NEAREST
543 // check whether fractional part of initial_P/10^extra_digits is
545 // this is the same as fractional part of
546 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
549 remainder_h
= CT
.w
[1] << (64 - amount
);
551 // test whether fractional part is 0
552 if (!remainder_h
&& (CT
.w
[0] < reciprocals10_64
[extra_digits
])) {
558 #ifdef SET_STATUS_FLAGS
559 status
= INEXACT_EXCEPTION
;
562 remainder_h
= CT
.w
[1] << (64 - amount
);
565 case ROUNDING_TO_NEAREST
:
566 case ROUNDING_TIES_AWAY
:
567 // test whether fractional part is 0
568 if ((remainder_h
== 0x8000000000000000ull
)
569 && (CT
.w
[0] < reciprocals10_64
[extra_digits
]))
570 status
= EXACT_STATUS
;
573 case ROUNDING_TO_ZERO
:
574 if (!remainder_h
&& (CT
.w
[0] < reciprocals10_64
[extra_digits
]))
575 status
= EXACT_STATUS
;
576 //if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y;
580 __add_carry_out (tmp
, carry
, CT
.w
[0],
581 reciprocals10_64
[extra_digits
]);
582 if ((remainder_h
>> (64 - amount
)) + carry
>=
583 (((UINT64
) 1) << amount
))
584 status
= EXACT_STATUS
;
587 __set_status_flags (pfpsf
, status
);
592 fast_get_BID64_check_OF (sign_s
, exponent_b
+ extra_digits
, C64
,