1 /* Copyright (C) 2007 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 2, or (at your option) any later
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
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
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
29 /*****************************************************************************
31 *****************************************************************************
33 * Algorithm description:
35 * if(exponent_a < exponent_b)
37 * diff_expon = exponent_a - exponent_b
40 * if(coefficient_a*10^diff_expon guaranteed below 2^62)
41 * S = sign_a*coefficient_a*10^diff_expon + sign_b*coefficient_b
43 * return get_BID64(sign(S),exponent_b,|S|)
45 * determine number of extra digits in S (1, 2, or 3)
46 * return rounded result
47 * else // large exponent difference
48 * if(number_digits(coefficient_a*10^diff_expon) +/- 10^16)
49 * guaranteed the same as
50 * number_digits(coefficient_a*10^diff_expon) )
51 * S = normalize(coefficient_a + (sign_a^sign_b)*10^(16-diff_expon))
52 * corr = 10^16 + (sign_a^sign_b)*coefficient_b
53 * corr*10^exponent_b is rounded so it aligns with S*10^exponent_S
54 * return get_BID64(sign_a,exponent(S),S+rounded(corr))
56 * add sign_a*coefficient_a*10^diff_expon, sign_b*coefficient_b
57 * in 128-bit integer arithmetic, then round to 16 decimal digits
60 ****************************************************************************/
62 #include "bid_internal.h"
64 #if DECIMAL_CALL_BY_REFERENCE
65 void bid64_add (UINT64
* pres
, UINT64
* px
,
67 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
70 UINT64
bid64_add (UINT64 x
,
71 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
72 _EXC_MASKS_PARAM _EXC_INFO_PARAM
);
75 #if DECIMAL_CALL_BY_REFERENCE
78 bid64_sub (UINT64
* pres
, UINT64
* px
,
80 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
83 #if !DECIMAL_GLOBAL_ROUNDING
84 _IDEC_round rnd_mode
= *prnd_mode
;
86 // check if y is not NaN
87 if (((y
& NAN_MASK64
) != NAN_MASK64
))
88 y
^= 0x8000000000000000ull
;
90 &y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
97 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
98 _EXC_MASKS_PARAM _EXC_INFO_PARAM
) {
99 // check if y is not NaN
100 if (((y
& NAN_MASK64
) != NAN_MASK64
))
101 y
^= 0x8000000000000000ull
;
104 y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
111 #if DECIMAL_CALL_BY_REFERENCE
114 bid64_add (UINT64
* pres
, UINT64
* px
,
116 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
123 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
124 _EXC_MASKS_PARAM _EXC_INFO_PARAM
) {
127 UINT128 CA
, CT
, CT_new
;
128 UINT64 sign_x
, sign_y
, coefficient_x
, coefficient_y
, C64_new
;
129 UINT64 valid_x
, valid_y
;
131 UINT64 sign_a
, sign_b
, coefficient_a
, coefficient_b
, sign_s
, sign_ab
,
133 UINT64 saved_ca
, saved_cb
, C0_64
, C64
, remainder_h
, T1
, carry
, tmp
;
135 int exponent_x
, exponent_y
, exponent_a
, exponent_b
, diff_dec_expon
;
136 int bin_expon_ca
, extra_digits
, amount
, scale_k
, scale_ca
;
137 unsigned rmode
, status
;
139 #if DECIMAL_CALL_BY_REFERENCE
140 #if !DECIMAL_GLOBAL_ROUNDING
141 _IDEC_round rnd_mode
= *prnd_mode
;
147 valid_x
= unpack_BID64 (&sign_x
, &exponent_x
, &coefficient_x
, x
);
148 valid_y
= unpack_BID64 (&sign_y
, &exponent_y
, &coefficient_y
, y
);
150 // unpack arguments, check for NaN or Infinity
155 if ((x
& NAN_MASK64
) == NAN_MASK64
) {
156 #ifdef SET_STATUS_FLAGS
157 if (((x
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
158 || ((y
& SNAN_MASK64
) == SNAN_MASK64
))
159 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
161 res
= coefficient_x
& QUIET_MASK64
;
165 if ((x
& INFINITY_MASK64
) == INFINITY_MASK64
) {
167 if (((y
& NAN_MASK64
) == INFINITY_MASK64
)) {
168 if (sign_x
== (y
& 0x8000000000000000ull
)) {
174 #ifdef SET_STATUS_FLAGS
175 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
182 if (((y
& NAN_MASK64
) == NAN_MASK64
)) {
183 res
= coefficient_y
& QUIET_MASK64
;
184 #ifdef SET_STATUS_FLAGS
185 if (((y
& SNAN_MASK64
) == SNAN_MASK64
))
186 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
190 // otherwise return +/-Inf
198 if (((y
& INFINITY_MASK64
) != INFINITY_MASK64
) && coefficient_y
) {
199 if (exponent_y
<= exponent_x
) {
209 if (((y
& INFINITY_MASK64
) == INFINITY_MASK64
)) {
210 #ifdef SET_STATUS_FLAGS
211 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
212 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
214 res
= coefficient_y
& QUIET_MASK64
;
218 if (!coefficient_x
) { // x==0
219 if (exponent_x
<= exponent_y
)
220 res
= ((UINT64
) exponent_x
) << 53;
222 res
= ((UINT64
) exponent_y
) << 53;
223 if (sign_x
== sign_y
)
225 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
226 #ifndef IEEE_ROUND_NEAREST
227 if (rnd_mode
== ROUNDING_DOWN
&& sign_x
!= sign_y
)
228 res
|= 0x8000000000000000ull
;
232 } else if (exponent_y
>= exponent_x
) {
237 // sort arguments by exponent
238 if (exponent_x
< exponent_y
) {
240 exponent_a
= exponent_y
;
241 coefficient_a
= coefficient_y
;
243 exponent_b
= exponent_x
;
244 coefficient_b
= coefficient_x
;
247 exponent_a
= exponent_x
;
248 coefficient_a
= coefficient_x
;
250 exponent_b
= exponent_y
;
251 coefficient_b
= coefficient_y
;
254 // exponent difference
255 diff_dec_expon
= exponent_a
- exponent_b
;
257 /* get binary coefficients of x and y */
259 //--- get number of bits in the coefficients of x and y ---
261 // version 2 (original)
262 tempx
.d
= (double) coefficient_a
;
263 bin_expon_ca
= ((tempx
.i
& MASK_BINARY_EXPONENT
) >> 52) - 0x3ff;
265 if (diff_dec_expon
> MAX_FORMAT_DIGITS
) {
266 // normalize a to a 16-digit coefficient
268 scale_ca
= estimate_decimal_digits
[bin_expon_ca
];
269 if (coefficient_a
>= power10_table_128
[scale_ca
].w
[0])
272 scale_k
= 16 - scale_ca
;
274 coefficient_a
*= power10_table_128
[scale_k
].w
[0];
276 diff_dec_expon
-= scale_k
;
277 exponent_a
-= scale_k
;
279 /* get binary coefficients of x and y */
281 //--- get number of bits in the coefficients of x and y ---
282 tempx
.d
= (double) coefficient_a
;
283 bin_expon_ca
= ((tempx
.i
& MASK_BINARY_EXPONENT
) >> 52) - 0x3ff;
285 if (diff_dec_expon
> MAX_FORMAT_DIGITS
) {
286 #ifdef SET_STATUS_FLAGS
288 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
292 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
293 #ifndef IEEE_ROUND_NEAREST
294 if (((rnd_mode
) & 3) && coefficient_b
) // not ROUNDING_TO_NEAREST
299 coefficient_a
-= ((((SINT64
) sign_a
) >> 63) | 1);
300 if (coefficient_a
< 1000000000000000ull) {
302 coefficient_a
= 9999999999999999ull;
303 } else if (coefficient_a
>= 10000000000000000ull) {
305 coefficient_a
= 1000000000000000ull;
311 coefficient_a
+= ((((SINT64
) sign_a
) >> 63) | 1);
312 if (coefficient_a
< 1000000000000000ull) {
314 coefficient_a
= 9999999999999999ull;
315 } else if (coefficient_a
>= 10000000000000000ull) {
317 coefficient_a
= 1000000000000000ull;
322 if (sign_a
!= sign_b
) {
324 if (coefficient_a
< 1000000000000000ull) {
326 coefficient_a
= 9999999999999999ull;
334 // check special case here
335 if ((coefficient_a
== 1000000000000000ull)
336 && (diff_dec_expon
== MAX_FORMAT_DIGITS
+ 1)
338 && (coefficient_b
> 5000000000000000ull)) {
339 coefficient_a
= 9999999999999999ull;
344 fast_get_BID64_check_OF (sign_a
, exponent_a
, coefficient_a
,
349 // test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62
350 if (bin_expon_ca
+ estimate_bin_expon
[diff_dec_expon
] < 60) {
351 // coefficient_a*10^(exponent_a-exponent_b)<2^63
353 // multiply by 10^(exponent_a-exponent_b)
354 coefficient_a
*= power10_table_128
[diff_dec_expon
].w
[0];
357 sign_b
= ((SINT64
) sign_b
) >> 63;
358 // apply sign to coeff. of b
359 coefficient_b
= (coefficient_b
+ sign_b
) ^ sign_b
;
361 // apply sign to coefficient a
362 sign_a
= ((SINT64
) sign_a
) >> 63;
363 coefficient_a
= (coefficient_a
+ sign_a
) ^ sign_a
;
365 coefficient_a
+= coefficient_b
;
367 sign_s
= ((SINT64
) coefficient_a
) >> 63;
368 coefficient_a
= (coefficient_a
+ sign_s
) ^ sign_s
;
369 sign_s
&= 0x8000000000000000ull
;
371 // coefficient_a < 10^16 ?
372 if (coefficient_a
< power10_table_128
[MAX_FORMAT_DIGITS
].w
[0]) {
373 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
374 #ifndef IEEE_ROUND_NEAREST
375 if (rnd_mode
== ROUNDING_DOWN
&& (!coefficient_a
)
377 sign_s
= 0x8000000000000000ull
;
380 res
= very_fast_get_BID64 (sign_s
, exponent_b
, coefficient_a
);
383 // otherwise rounding is necessary
385 // already know coefficient_a<10^19
386 // coefficient_a < 10^17 ?
387 if (coefficient_a
< power10_table_128
[17].w
[0])
389 else if (coefficient_a
< power10_table_128
[18].w
[0])
394 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
395 #ifndef IEEE_ROUND_NEAREST
397 if (sign_s
&& (unsigned) (rmode
- 1) < 2)
405 coefficient_a
+= round_const_table
[rmode
][extra_digits
];
407 // get P*(2^M[extra_digits])/10^extra_digits
408 __mul_64x64_to_128 (CT
, coefficient_a
,
409 reciprocals10_64
[extra_digits
]);
411 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
412 amount
= short_recip_scale
[extra_digits
];
413 C64
= CT
.w
[1] >> amount
;
416 // coefficient_a*10^(exponent_a-exponent_b) is large
419 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
420 #ifndef IEEE_ROUND_NEAREST
422 if (sign_s
&& (unsigned) (rmode
- 1) < 2)
431 // check whether we can take faster path
432 scale_ca
= estimate_decimal_digits
[bin_expon_ca
];
434 sign_ab
= sign_a
^ sign_b
;
435 sign_ab
= ((SINT64
) sign_ab
) >> 63;
437 // T1 = 10^(16-diff_dec_expon)
438 T1
= power10_table_128
[16 - diff_dec_expon
].w
[0];
440 // get number of digits in coefficient_a
441 if (coefficient_a
>= power10_table_128
[scale_ca
].w
[0]) {
445 scale_k
= 16 - scale_ca
;
448 saved_ca
= coefficient_a
- T1
;
450 (SINT64
) saved_ca
*(SINT64
) power10_table_128
[scale_k
].w
[0];
451 extra_digits
= diff_dec_expon
- scale_k
;
454 saved_cb
= (coefficient_b
+ sign_ab
) ^ sign_ab
;
455 // add 10^16 and rounding constant
457 saved_cb
+ 10000000000000000ull +
458 round_const_table
[rmode
][extra_digits
];
460 // get P*(2^M[extra_digits])/10^extra_digits
461 __mul_64x64_to_128 (CT
, coefficient_b
,
462 reciprocals10_64
[extra_digits
]);
464 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
465 amount
= short_recip_scale
[extra_digits
];
466 C0_64
= CT
.w
[1] >> amount
;
468 // result coefficient
469 C64
= C0_64
+ coefficient_a
;
470 // filter out difficult (corner) cases
471 // this test ensures the number of digits in coefficient_a does not change
472 // after adding (the appropriately scaled and rounded) coefficient_b
473 if ((UINT64
) (C64
- 1000000000000000ull - 1) >
474 9000000000000000ull - 2) {
475 if (C64
>= 10000000000000000ull) {
476 // result has more than 16 digits
478 // must divide coeff_a by 10
479 saved_ca
= saved_ca
+ T1
;
480 __mul_64x64_to_128 (CA
, saved_ca
, 0x3333333333333334ull
);
481 //reciprocals10_64[1]);
482 coefficient_a
= CA
.w
[1] >> 1;
484 saved_ca
- (coefficient_a
<< 3) - (coefficient_a
<< 1);
485 coefficient_a
= coefficient_a
- T1
;
487 saved_cb
+= rem_a
* power10_table_128
[diff_dec_expon
].w
[0];
490 (SINT64
) (saved_ca
- T1
-
491 (T1
<< 3)) * (SINT64
) power10_table_128
[scale_k
-
496 saved_cb
+ 100000000000000000ull +
497 round_const_table
[rmode
][extra_digits
];
499 // get P*(2^M[extra_digits])/10^extra_digits
500 __mul_64x64_to_128 (CT
, coefficient_b
,
501 reciprocals10_64
[extra_digits
]);
503 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
504 amount
= short_recip_scale
[extra_digits
];
505 C0_64
= CT
.w
[1] >> amount
;
507 // result coefficient
508 C64
= C0_64
+ coefficient_a
;
509 } else if (C64
<= 1000000000000000ull) {
510 // less than 16 digits in result
512 (SINT64
) saved_ca
*(SINT64
) power10_table_128
[scale_k
+
517 (saved_cb
<< 3) + (saved_cb
<< 1) + 100000000000000000ull +
518 round_const_table
[rmode
][extra_digits
];
520 // get P*(2^M[extra_digits])/10^extra_digits
521 __mul_64x64_to_128 (CT_new
, coefficient_b
,
522 reciprocals10_64
[extra_digits
]);
524 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
525 amount
= short_recip_scale
[extra_digits
];
526 C0_64
= CT_new
.w
[1] >> amount
;
528 // result coefficient
529 C64_new
= C0_64
+ coefficient_a
;
530 if (C64_new
< 10000000000000000ull) {
532 #ifdef SET_STATUS_FLAGS
543 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
544 #ifndef IEEE_ROUND_NEAREST
545 if (rmode
== 0) //ROUNDING_TO_NEAREST
548 // check whether fractional part of initial_P/10^extra_digits is
550 // this is the same as fractional part of
551 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
554 remainder_h
= CT
.w
[1] << (64 - amount
);
556 // test whether fractional part is 0
557 if (!remainder_h
&& (CT
.w
[0] < reciprocals10_64
[extra_digits
])) {
563 #ifdef SET_STATUS_FLAGS
564 status
= INEXACT_EXCEPTION
;
567 remainder_h
= CT
.w
[1] << (64 - amount
);
570 case ROUNDING_TO_NEAREST
:
571 case ROUNDING_TIES_AWAY
:
572 // test whether fractional part is 0
573 if ((remainder_h
== 0x8000000000000000ull
)
574 && (CT
.w
[0] < reciprocals10_64
[extra_digits
]))
575 status
= EXACT_STATUS
;
578 case ROUNDING_TO_ZERO
:
579 if (!remainder_h
&& (CT
.w
[0] < reciprocals10_64
[extra_digits
]))
580 status
= EXACT_STATUS
;
581 //if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y;
585 __add_carry_out (tmp
, carry
, CT
.w
[0],
586 reciprocals10_64
[extra_digits
]);
587 if ((remainder_h
>> (64 - amount
)) + carry
>=
588 (((UINT64
) 1) << amount
))
589 status
= EXACT_STATUS
;
592 __set_status_flags (pfpsf
, status
);
597 fast_get_BID64_check_OF (sign_s
, exponent_b
+ extra_digits
, C64
,