]> git.ipfire.org Git - thirdparty/gcc.git/blame - libgcc/config/libbid/bid64_add.c
Update copyright years.
[thirdparty/gcc.git] / libgcc / config / libbid / bid64_add.c
CommitLineData
8d9254fc 1/* Copyright (C) 2007-2020 Free Software Foundation, Inc.
200359e8
L
2
3This file is part of GCC.
4
5GCC is free software; you can redistribute it and/or modify it under
6the terms of the GNU General Public License as published by the Free
748086b7 7Software Foundation; either version 3, or (at your option) any later
200359e8
L
8version.
9
200359e8
L
10GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11WARRANTY; without even the implied warranty of MERCHANTABILITY or
12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13for more details.
14
748086b7
JJ
15Under Section 7 of GPL version 3, you are granted additional
16permissions described in the GCC Runtime Library Exception, version
173.1, as published by the Free Software Foundation.
18
19You should have received a copy of the GNU General Public License and
20a copy of the GCC Runtime Library Exception along with this program;
21see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22<http://www.gnu.org/licenses/>. */
200359e8
L
23
24/*****************************************************************************
25 * BID64 add
26 *****************************************************************************
27 *
28 * Algorithm description:
29 *
30 * if(exponent_a < exponent_b)
31 * switch a, b
32 * diff_expon = exponent_a - exponent_b
33 * if(diff_expon > 16)
34 * return normalize(a)
35 * if(coefficient_a*10^diff_expon guaranteed below 2^62)
36 * S = sign_a*coefficient_a*10^diff_expon + sign_b*coefficient_b
37 * if(|S|<10^16)
38 * return get_BID64(sign(S),exponent_b,|S|)
39 * else
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))
50 * else
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
53 *
54 *
55 ****************************************************************************/
56
57#include "bid_internal.h"
58
59#if DECIMAL_CALL_BY_REFERENCE
b2a00c89 60void bid64_add (UINT64 * pres, UINT64 * px,
200359e8
L
61 UINT64 *
62 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
63 _EXC_INFO_PARAM);
64#else
b2a00c89
L
65UINT64 bid64_add (UINT64 x,
66 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
67 _EXC_MASKS_PARAM _EXC_INFO_PARAM);
200359e8
L
68#endif
69
70#if DECIMAL_CALL_BY_REFERENCE
71
72void
b2a00c89 73bid64_sub (UINT64 * pres, UINT64 * px,
200359e8
L
74 UINT64 *
75 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
76 _EXC_INFO_PARAM) {
77 UINT64 y = *py;
78#if !DECIMAL_GLOBAL_ROUNDING
79 _IDEC_round rnd_mode = *prnd_mode;
80#endif
81 // check if y is not NaN
82 if (((y & NAN_MASK64) != NAN_MASK64))
83 y ^= 0x8000000000000000ull;
b2a00c89 84 bid64_add (pres, px,
200359e8
L
85 &y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
86 _EXC_INFO_ARG);
87}
88#else
89
90UINT64
b2a00c89 91bid64_sub (UINT64 x,
200359e8
L
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;
97
b2a00c89 98 return bid64_add (x,
200359e8
L
99 y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
100 _EXC_INFO_ARG);
101}
102#endif
103
104
105
106#if DECIMAL_CALL_BY_REFERENCE
107
108void
b2a00c89 109bid64_add (UINT64 * pres, UINT64 * px,
200359e8
L
110 UINT64 *
111 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
112 _EXC_INFO_PARAM) {
113 UINT64 x, y;
114#else
115
116UINT64
b2a00c89 117bid64_add (UINT64 x,
200359e8
L
118 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
119 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
120#endif
121
122 UINT128 CA, CT, CT_new;
123 UINT64 sign_x, sign_y, coefficient_x, coefficient_y, C64_new;
124 UINT64 valid_x, valid_y;
125 UINT64 res;
126 UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab,
127 rem_a;
128 UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp;
129 int_double tempx;
b2a00c89 130 int exponent_x, exponent_y, exponent_a, exponent_b, diff_dec_expon;
200359e8
L
131 int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
132 unsigned rmode, status;
133
134#if DECIMAL_CALL_BY_REFERENCE
135#if !DECIMAL_GLOBAL_ROUNDING
136 _IDEC_round rnd_mode = *prnd_mode;
137#endif
138 x = *px;
139 y = *py;
140#endif
141
142 valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
143 valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
144
145 // unpack arguments, check for NaN or Infinity
146 if (!valid_x) {
147 // x is Inf. or NaN
148
149 // test if x is NaN
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);
155#endif
b2a00c89 156 res = coefficient_x & QUIET_MASK64;
200359e8
L
157 BID_RETURN (res);
158 }
159 // x is Infinity?
160 if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
161 // check if y is Inf
162 if (((y & NAN_MASK64) == INFINITY_MASK64)) {
163 if (sign_x == (y & 0x8000000000000000ull)) {
b2a00c89 164 res = coefficient_x;
200359e8
L
165 BID_RETURN (res);
166 }
167 // return NaN
168 {
169#ifdef SET_STATUS_FLAGS
170 __set_status_flags (pfpsf, INVALID_EXCEPTION);
171#endif
172 res = NAN_MASK64;
173 BID_RETURN (res);
174 }
175 }
176 // check if y is NaN
177 if (((y & NAN_MASK64) == NAN_MASK64)) {
b2a00c89 178 res = coefficient_y & QUIET_MASK64;
200359e8
L
179#ifdef SET_STATUS_FLAGS
180 if (((y & SNAN_MASK64) == SNAN_MASK64))
181 __set_status_flags (pfpsf, INVALID_EXCEPTION);
182#endif
183 BID_RETURN (res);
184 }
185 // otherwise return +/-Inf
186 {
b2a00c89 187 res = coefficient_x;
200359e8
L
188 BID_RETURN (res);
189 }
190 }
191 // x is 0
192 {
b2a00c89 193 if (((y & INFINITY_MASK64) != INFINITY_MASK64) && coefficient_y) {
200359e8
L
194 if (exponent_y <= exponent_x) {
195 res = y;
196 BID_RETURN (res);
197 }
198 }
199 }
200
201 }
202 if (!valid_y) {
203 // y is Inf. or NaN?
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);
208#endif
b2a00c89 209 res = coefficient_y & QUIET_MASK64;
200359e8
L
210 BID_RETURN (res);
211 }
212 // y is 0
b2a00c89 213 if (!coefficient_x) { // x==0
200359e8 214 if (exponent_x <= exponent_y)
b2a00c89 215 res = ((UINT64) exponent_x) << 53;
200359e8 216 else
b2a00c89 217 res = ((UINT64) exponent_y) << 53;
200359e8
L
218 if (sign_x == sign_y)
219 res |= sign_x;
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;
224#endif
225#endif
226 BID_RETURN (res);
227 } else if (exponent_y >= exponent_x) {
228 res = x;
229 BID_RETURN (res);
230 }
231 }
232 // sort arguments by exponent
233 if (exponent_x < exponent_y) {
234 sign_a = sign_y;
235 exponent_a = exponent_y;
236 coefficient_a = coefficient_y;
237 sign_b = sign_x;
238 exponent_b = exponent_x;
239 coefficient_b = coefficient_x;
240 } else {
241 sign_a = sign_x;
242 exponent_a = exponent_x;
243 coefficient_a = coefficient_x;
244 sign_b = sign_y;
245 exponent_b = exponent_y;
246 coefficient_b = coefficient_y;
247 }
248
249 // exponent difference
250 diff_dec_expon = exponent_a - exponent_b;
251
252 /* get binary coefficients of x and y */
253
254 //--- get number of bits in the coefficients of x and y ---
255
256 // version 2 (original)
257 tempx.d = (double) coefficient_a;
258 bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
259
260 if (diff_dec_expon > MAX_FORMAT_DIGITS) {
261 // normalize a to a 16-digit coefficient
262
b2a00c89
L
263 scale_ca = estimate_decimal_digits[bin_expon_ca];
264 if (coefficient_a >= power10_table_128[scale_ca].w[0])
200359e8
L
265 scale_ca++;
266
267 scale_k = 16 - scale_ca;
268
b2a00c89 269 coefficient_a *= power10_table_128[scale_k].w[0];
200359e8
L
270
271 diff_dec_expon -= scale_k;
272 exponent_a -= scale_k;
273
274 /* get binary coefficients of x and y */
275
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;
279
280 if (diff_dec_expon > MAX_FORMAT_DIGITS) {
281#ifdef SET_STATUS_FLAGS
282 if (coefficient_b) {
283 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
284 }
285#endif
286
287#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
288#ifndef IEEE_ROUND_NEAREST
289 if (((rnd_mode) & 3) && coefficient_b) // not ROUNDING_TO_NEAREST
290 {
291 switch (rnd_mode) {
292 case ROUNDING_DOWN:
293 if (sign_b) {
294 coefficient_a -= ((((SINT64) sign_a) >> 63) | 1);
295 if (coefficient_a < 1000000000000000ull) {
296 exponent_a--;
297 coefficient_a = 9999999999999999ull;
298 } else if (coefficient_a >= 10000000000000000ull) {
299 exponent_a++;
300 coefficient_a = 1000000000000000ull;
301 }
302 }
303 break;
304 case ROUNDING_UP:
305 if (!sign_b) {
306 coefficient_a += ((((SINT64) sign_a) >> 63) | 1);
307 if (coefficient_a < 1000000000000000ull) {
308 exponent_a--;
309 coefficient_a = 9999999999999999ull;
310 } else if (coefficient_a >= 10000000000000000ull) {
311 exponent_a++;
312 coefficient_a = 1000000000000000ull;
313 }
314 }
315 break;
316 default: // RZ
317 if (sign_a != sign_b) {
318 coefficient_a--;
319 if (coefficient_a < 1000000000000000ull) {
320 exponent_a--;
321 coefficient_a = 9999999999999999ull;
322 }
323 }
324 break;
325 }
326 } else
327#endif
328#endif
329 // check special case here
330 if ((coefficient_a == 1000000000000000ull)
331 && (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
b2a00c89
L
332 && (sign_a ^ sign_b)
333 && (coefficient_b > 5000000000000000ull)) {
200359e8
L
334 coefficient_a = 9999999999999999ull;
335 exponent_a--;
336 }
337
338 res =
339 fast_get_BID64_check_OF (sign_a, exponent_a, coefficient_a,
340 rnd_mode, pfpsf);
341 BID_RETURN (res);
342 }
343 }
344 // test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62
b2a00c89 345 if (bin_expon_ca + estimate_bin_expon[diff_dec_expon] < 60) {
200359e8
L
346 // coefficient_a*10^(exponent_a-exponent_b)<2^63
347
348 // multiply by 10^(exponent_a-exponent_b)
b2a00c89 349 coefficient_a *= power10_table_128[diff_dec_expon].w[0];
200359e8
L
350
351 // sign mask
352 sign_b = ((SINT64) sign_b) >> 63;
353 // apply sign to coeff. of b
354 coefficient_b = (coefficient_b + sign_b) ^ sign_b;
355
356 // apply sign to coefficient a
357 sign_a = ((SINT64) sign_a) >> 63;
358 coefficient_a = (coefficient_a + sign_a) ^ sign_a;
359
360 coefficient_a += coefficient_b;
361 // get sign
362 sign_s = ((SINT64) coefficient_a) >> 63;
363 coefficient_a = (coefficient_a + sign_s) ^ sign_s;
364 sign_s &= 0x8000000000000000ull;
365
366 // coefficient_a < 10^16 ?
b2a00c89 367 if (coefficient_a < power10_table_128[MAX_FORMAT_DIGITS].w[0]) {
200359e8
L
368#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
369#ifndef IEEE_ROUND_NEAREST
370 if (rnd_mode == ROUNDING_DOWN && (!coefficient_a)
371 && sign_a != sign_b)
372 sign_s = 0x8000000000000000ull;
373#endif
374#endif
375 res = very_fast_get_BID64 (sign_s, exponent_b, coefficient_a);
376 BID_RETURN (res);
377 }
378 // otherwise rounding is necessary
379
380 // already know coefficient_a<10^19
381 // coefficient_a < 10^17 ?
b2a00c89 382 if (coefficient_a < power10_table_128[17].w[0])
200359e8 383 extra_digits = 1;
b2a00c89 384 else if (coefficient_a < power10_table_128[18].w[0])
200359e8
L
385 extra_digits = 2;
386 else
387 extra_digits = 3;
388
389#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
390#ifndef IEEE_ROUND_NEAREST
391 rmode = rnd_mode;
392 if (sign_s && (unsigned) (rmode - 1) < 2)
393 rmode = 3 - rmode;
394#else
395 rmode = 0;
396#endif
397#else
398 rmode = 0;
399#endif
b2a00c89 400 coefficient_a += round_const_table[rmode][extra_digits];
200359e8
L
401
402 // get P*(2^M[extra_digits])/10^extra_digits
403 __mul_64x64_to_128 (CT, coefficient_a,
b2a00c89 404 reciprocals10_64[extra_digits]);
200359e8
L
405
406 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
b2a00c89 407 amount = short_recip_scale[extra_digits];
200359e8
L
408 C64 = CT.w[1] >> amount;
409
410 } else {
411 // coefficient_a*10^(exponent_a-exponent_b) is large
412 sign_s = sign_a;
413
414#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
415#ifndef IEEE_ROUND_NEAREST
416 rmode = rnd_mode;
417 if (sign_s && (unsigned) (rmode - 1) < 2)
418 rmode = 3 - rmode;
419#else
420 rmode = 0;
421#endif
422#else
423 rmode = 0;
424#endif
425
426 // check whether we can take faster path
b2a00c89 427 scale_ca = estimate_decimal_digits[bin_expon_ca];
200359e8
L
428
429 sign_ab = sign_a ^ sign_b;
430 sign_ab = ((SINT64) sign_ab) >> 63;
431
432 // T1 = 10^(16-diff_dec_expon)
b2a00c89 433 T1 = power10_table_128[16 - diff_dec_expon].w[0];
200359e8
L
434
435 // get number of digits in coefficient_a
b2a00c89 436 if (coefficient_a >= power10_table_128[scale_ca].w[0]) {
200359e8
L
437 scale_ca++;
438 }
439
440 scale_k = 16 - scale_ca;
441
442 // addition
443 saved_ca = coefficient_a - T1;
444 coefficient_a =
b2a00c89 445 (SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0];
200359e8
L
446 extra_digits = diff_dec_expon - scale_k;
447
448 // apply sign
449 saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
450 // add 10^16 and rounding constant
451 coefficient_b =
452 saved_cb + 10000000000000000ull +
b2a00c89 453 round_const_table[rmode][extra_digits];
200359e8
L
454
455 // get P*(2^M[extra_digits])/10^extra_digits
456 __mul_64x64_to_128 (CT, coefficient_b,
b2a00c89 457 reciprocals10_64[extra_digits]);
200359e8
L
458
459 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
b2a00c89 460 amount = short_recip_scale[extra_digits];
200359e8
L
461 C0_64 = CT.w[1] >> amount;
462
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
b2a00c89
L
468 if ((UINT64) (C64 - 1000000000000000ull - 1) >
469 9000000000000000ull - 2) {
200359e8
L
470 if (C64 >= 10000000000000000ull) {
471 // result has more than 16 digits
472 if (!scale_k) {
473 // must divide coeff_a by 10
474 saved_ca = saved_ca + T1;
b2a00c89
L
475 __mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull);
476 //reciprocals10_64[1]);
200359e8
L
477 coefficient_a = CA.w[1] >> 1;
478 rem_a =
479 saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
480 coefficient_a = coefficient_a - T1;
481
b2a00c89 482 saved_cb += rem_a * power10_table_128[diff_dec_expon].w[0];
200359e8
L
483 } else
484 coefficient_a =
485 (SINT64) (saved_ca - T1 -
b2a00c89 486 (T1 << 3)) * (SINT64) power10_table_128[scale_k -
200359e8
L
487 1].w[0];
488
489 extra_digits++;
490 coefficient_b =
491 saved_cb + 100000000000000000ull +
b2a00c89 492 round_const_table[rmode][extra_digits];
200359e8
L
493
494 // get P*(2^M[extra_digits])/10^extra_digits
495 __mul_64x64_to_128 (CT, coefficient_b,
b2a00c89 496 reciprocals10_64[extra_digits]);
200359e8
L
497
498 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
b2a00c89 499 amount = short_recip_scale[extra_digits];
200359e8
L
500 C0_64 = CT.w[1] >> amount;
501
502 // result coefficient
503 C64 = C0_64 + coefficient_a;
504 } else if (C64 <= 1000000000000000ull) {
505 // less than 16 digits in result
506 coefficient_a =
b2a00c89 507 (SINT64) saved_ca *(SINT64) power10_table_128[scale_k +
200359e8
L
508 1].w[0];
509 //extra_digits --;
510 exponent_b--;
511 coefficient_b =
512 (saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull +
b2a00c89 513 round_const_table[rmode][extra_digits];
200359e8
L
514
515 // get P*(2^M[extra_digits])/10^extra_digits
516 __mul_64x64_to_128 (CT_new, coefficient_b,
b2a00c89 517 reciprocals10_64[extra_digits]);
200359e8
L
518
519 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
b2a00c89 520 amount = short_recip_scale[extra_digits];
200359e8
L
521 C0_64 = CT_new.w[1] >> amount;
522
523 // result coefficient
524 C64_new = C0_64 + coefficient_a;
525 if (C64_new < 10000000000000000ull) {
526 C64 = C64_new;
527#ifdef SET_STATUS_FLAGS
528 CT = CT_new;
529#endif
530 } else
531 exponent_b++;
532 }
533
534 }
535
536 }
537
538#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
539#ifndef IEEE_ROUND_NEAREST
540 if (rmode == 0) //ROUNDING_TO_NEAREST
541#endif
542 if (C64 & 1) {
543 // check whether fractional part of initial_P/10^extra_digits is
544 // exactly .5
545 // this is the same as fractional part of
546 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
547
548 // get remainder
549 remainder_h = CT.w[1] << (64 - amount);
550
551 // test whether fractional part is 0
b2a00c89 552 if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
200359e8
L
553 C64--;
554 }
555 }
556#endif
557
558#ifdef SET_STATUS_FLAGS
559 status = INEXACT_EXCEPTION;
560
561 // get remainder
562 remainder_h = CT.w[1] << (64 - amount);
563
564 switch (rmode) {
565 case ROUNDING_TO_NEAREST:
566 case ROUNDING_TIES_AWAY:
567 // test whether fractional part is 0
568 if ((remainder_h == 0x8000000000000000ull)
b2a00c89 569 && (CT.w[0] < reciprocals10_64[extra_digits]))
200359e8
L
570 status = EXACT_STATUS;
571 break;
572 case ROUNDING_DOWN:
573 case ROUNDING_TO_ZERO:
b2a00c89 574 if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits]))
200359e8
L
575 status = EXACT_STATUS;
576 //if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y;
577 break;
578 default:
579 // round up
580 __add_carry_out (tmp, carry, CT.w[0],
b2a00c89 581 reciprocals10_64[extra_digits]);
200359e8
L
582 if ((remainder_h >> (64 - amount)) + carry >=
583 (((UINT64) 1) << amount))
584 status = EXACT_STATUS;
585 break;
586 }
587 __set_status_flags (pfpsf, status);
588
589#endif
590
591 res =
592 fast_get_BID64_check_OF (sign_s, exponent_b + extra_digits, C64,
593 rnd_mode, pfpsf);
594 BID_RETURN (res);
595}