libgcc/config/libbid/bid64_add.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 /*****************************************************************************
  30  *    BID64 add
  31  *****************************************************************************
  32  *
  33  *  Algorithm description:
  34  *
  35  *   if(exponent_a < exponent_b)
  36  *       switch a, b
  37  *   diff_expon = exponent_a - exponent_b
  38  *   if(diff_expon > 16)
  39  *      return normalize(a)
  40  *   if(coefficient_a*10^diff_expon guaranteed below 2^62)
  41  *       S = sign_a*coefficient_a*10^diff_expon + sign_b*coefficient_b
  42  *       if(|S|<10^16)
  43  *           return get_BID64(sign(S),exponent_b,|S|)
  44  *       else
  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))
  55  *       else
  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
  58  *
  59  *
  60  ****************************************************************************/
  61
  62 #include "bid_internal.h"
  63
  64 #if DECIMAL_CALL_BY_REFERENCE
  65 void bid64_add (UINT64 * pres, UINT64 * px,
  66                 UINT64 *
  67                 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  68                 _EXC_INFO_PARAM);
  69 #else
  70 UINT64 bid64_add (UINT64 x,
  71                   UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
  72                   _EXC_MASKS_PARAM _EXC_INFO_PARAM);
  73 #endif
  74
  75 #if DECIMAL_CALL_BY_REFERENCE
  76
  77 void
  78 bid64_sub (UINT64 * pres, UINT64 * px,
  79            UINT64 *
  80            py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  81            _EXC_INFO_PARAM) {
  82   UINT64 y = *py;
  83 #if !DECIMAL_GLOBAL_ROUNDING
  84   _IDEC_round rnd_mode = *prnd_mode;
  85 #endif
  86   // check if y is not NaN
  87   if (((y & NAN_MASK64) != NAN_MASK64))
  88     y ^= 0x8000000000000000ull;
  89   bid64_add (pres, px,
  90              &y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
  91              _EXC_INFO_ARG);
  92 }
  93 #else
  94
  95 UINT64
  96 bid64_sub (UINT64 x,
  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;
 102
 103   return bid64_add (x,
 104                     y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 105                     _EXC_INFO_ARG);
 106 }
 107 #endif
 108
 109
 110
 111 #if DECIMAL_CALL_BY_REFERENCE
 112
 113 void
 114 bid64_add (UINT64 * pres, UINT64 * px,
 115            UINT64 *
 116            py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 117            _EXC_INFO_PARAM) {
 118   UINT64 x, y;
 119 #else
 120
 121 UINT64
 122 bid64_add (UINT64 x,
 123            UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
 124            _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
 125 #endif
 126
 127   UINT128 CA, CT, CT_new;
 128   UINT64 sign_x, sign_y, coefficient_x, coefficient_y, C64_new;
 129   UINT64 valid_x, valid_y;
 130   UINT64 res;
 131   UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab,
 132     rem_a;
 133   UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp;
 134   int_double tempx;
 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;
 138
 139 #if DECIMAL_CALL_BY_REFERENCE
 140 #if !DECIMAL_GLOBAL_ROUNDING
 141   _IDEC_round rnd_mode = *prnd_mode;
 142 #endif
 143   x = *px;
 144   y = *py;
 145 #endif
 146
 147   valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
 148   valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
 149
 150   // unpack arguments, check for NaN or Infinity
 151   if (!valid_x) {
 152     // x is Inf. or NaN
 153
 154     // test if x is NaN
 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);
 160 #endif
 161       res = coefficient_x & QUIET_MASK64;
 162       BID_RETURN (res);
 163     }
 164     // x is Infinity?
 165     if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
 166       // check if y is Inf
 167       if (((y & NAN_MASK64) == INFINITY_MASK64)) {
 168         if (sign_x == (y & 0x8000000000000000ull)) {
 169           res = coefficient_x;
 170           BID_RETURN (res);
 171         }
 172         // return NaN
 173         {
 174 #ifdef SET_STATUS_FLAGS
 175           __set_status_flags (pfpsf, INVALID_EXCEPTION);
 176 #endif
 177           res = NAN_MASK64;
 178           BID_RETURN (res);
 179         }
 180       }
 181       // check if y is NaN
 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);
 187 #endif
 188         BID_RETURN (res);
 189       }
 190       // otherwise return +/-Inf
 191       {
 192         res = coefficient_x;
 193         BID_RETURN (res);
 194       }
 195     }
 196     // x is 0
 197     {
 198       if (((y & INFINITY_MASK64) != INFINITY_MASK64) && coefficient_y) {
 199         if (exponent_y <= exponent_x) {
 200           res = y;
 201           BID_RETURN (res);
 202         }
 203       }
 204     }
 205
 206   }
 207   if (!valid_y) {
 208     // y is Inf. or NaN?
 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);
 213 #endif
 214       res = coefficient_y & QUIET_MASK64;
 215       BID_RETURN (res);
 216     }
 217     // y is 0
 218     if (!coefficient_x) {       // x==0
 219       if (exponent_x <= exponent_y)
 220         res = ((UINT64) exponent_x) << 53;
 221       else
 222         res = ((UINT64) exponent_y) << 53;
 223       if (sign_x == sign_y)
 224         res |= sign_x;
 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;
 229 #endif
 230 #endif
 231       BID_RETURN (res);
 232     } else if (exponent_y >= exponent_x) {
 233       res = x;
 234       BID_RETURN (res);
 235     }
 236   }
 237   // sort arguments by exponent
 238   if (exponent_x < exponent_y) {
 239     sign_a = sign_y;
 240     exponent_a = exponent_y;
 241     coefficient_a = coefficient_y;
 242     sign_b = sign_x;
 243     exponent_b = exponent_x;
 244     coefficient_b = coefficient_x;
 245   } else {
 246     sign_a = sign_x;
 247     exponent_a = exponent_x;
 248     coefficient_a = coefficient_x;
 249     sign_b = sign_y;
 250     exponent_b = exponent_y;
 251     coefficient_b = coefficient_y;
 252   }
 253
 254   // exponent difference
 255   diff_dec_expon = exponent_a - exponent_b;
 256
 257   /* get binary coefficients of x and y */
 258
 259   //--- get number of bits in the coefficients of x and y ---
 260
 261   // version 2 (original)
 262   tempx.d = (double) coefficient_a;
 263   bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
 264
 265   if (diff_dec_expon > MAX_FORMAT_DIGITS) {
 266     // normalize a to a 16-digit coefficient
 267
 268     scale_ca = estimate_decimal_digits[bin_expon_ca];
 269     if (coefficient_a >= power10_table_128[scale_ca].w[0])
 270       scale_ca++;
 271
 272     scale_k = 16 - scale_ca;
 273
 274     coefficient_a *= power10_table_128[scale_k].w[0];
 275
 276     diff_dec_expon -= scale_k;
 277     exponent_a -= scale_k;
 278
 279     /* get binary coefficients of x and y */
 280
 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;
 284
 285     if (diff_dec_expon > MAX_FORMAT_DIGITS) {
 286 #ifdef SET_STATUS_FLAGS
 287       if (coefficient_b) {
 288         __set_status_flags (pfpsf, INEXACT_EXCEPTION);
 289       }
 290 #endif
 291
 292 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 293 #ifndef IEEE_ROUND_NEAREST
 294       if (((rnd_mode) & 3) && coefficient_b)    // not ROUNDING_TO_NEAREST
 295       {
 296         switch (rnd_mode) {
 297         case ROUNDING_DOWN:
 298           if (sign_b) {
 299             coefficient_a -= ((((SINT64) sign_a) >> 63) | 1);
 300             if (coefficient_a < 1000000000000000ull) {
 301               exponent_a--;
 302               coefficient_a = 9999999999999999ull;
 303             } else if (coefficient_a >= 10000000000000000ull) {
 304               exponent_a++;
 305               coefficient_a = 1000000000000000ull;
 306             }
 307           }
 308           break;
 309         case ROUNDING_UP:
 310           if (!sign_b) {
 311             coefficient_a += ((((SINT64) sign_a) >> 63) | 1);
 312             if (coefficient_a < 1000000000000000ull) {
 313               exponent_a--;
 314               coefficient_a = 9999999999999999ull;
 315             } else if (coefficient_a >= 10000000000000000ull) {
 316               exponent_a++;
 317               coefficient_a = 1000000000000000ull;
 318             }
 319           }
 320           break;
 321         default:        // RZ
 322           if (sign_a != sign_b) {
 323             coefficient_a--;
 324             if (coefficient_a < 1000000000000000ull) {
 325               exponent_a--;
 326               coefficient_a = 9999999999999999ull;
 327             }
 328           }
 329           break;
 330         }
 331       } else
 332 #endif
 333 #endif
 334         // check special case here
 335         if ((coefficient_a == 1000000000000000ull)
 336             && (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
 337             && (sign_a ^ sign_b)
 338             && (coefficient_b > 5000000000000000ull)) {
 339         coefficient_a = 9999999999999999ull;
 340         exponent_a--;
 341       }
 342
 343       res =
 344         fast_get_BID64_check_OF (sign_a, exponent_a, coefficient_a,
 345                                  rnd_mode, pfpsf);
 346       BID_RETURN (res);
 347     }
 348   }
 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
 352
 353     // multiply by 10^(exponent_a-exponent_b)
 354     coefficient_a *= power10_table_128[diff_dec_expon].w[0];
 355
 356     // sign mask
 357     sign_b = ((SINT64) sign_b) >> 63;
 358     // apply sign to coeff. of b
 359     coefficient_b = (coefficient_b + sign_b) ^ sign_b;
 360
 361     // apply sign to coefficient a
 362     sign_a = ((SINT64) sign_a) >> 63;
 363     coefficient_a = (coefficient_a + sign_a) ^ sign_a;
 364
 365     coefficient_a += coefficient_b;
 366     // get sign
 367     sign_s = ((SINT64) coefficient_a) >> 63;
 368     coefficient_a = (coefficient_a + sign_s) ^ sign_s;
 369     sign_s &= 0x8000000000000000ull;
 370
 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)
 376           && sign_a != sign_b)
 377         sign_s = 0x8000000000000000ull;
 378 #endif
 379 #endif
 380       res = very_fast_get_BID64 (sign_s, exponent_b, coefficient_a);
 381       BID_RETURN (res);
 382     }
 383     // otherwise rounding is necessary
 384
 385     // already know coefficient_a<10^19
 386     // coefficient_a < 10^17 ?
 387     if (coefficient_a < power10_table_128[17].w[0])
 388       extra_digits = 1;
 389     else if (coefficient_a < power10_table_128[18].w[0])
 390       extra_digits = 2;
 391     else
 392       extra_digits = 3;
 393
 394 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 395 #ifndef IEEE_ROUND_NEAREST
 396     rmode = rnd_mode;
 397     if (sign_s && (unsigned) (rmode - 1) < 2)
 398       rmode = 3 - rmode;
 399 #else
 400     rmode = 0;
 401 #endif
 402 #else
 403     rmode = 0;
 404 #endif
 405     coefficient_a += round_const_table[rmode][extra_digits];
 406
 407     // get P*(2^M[extra_digits])/10^extra_digits
 408     __mul_64x64_to_128 (CT, coefficient_a,
 409                         reciprocals10_64[extra_digits]);
 410
 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;
 414
 415   } else {
 416     // coefficient_a*10^(exponent_a-exponent_b) is large
 417     sign_s = sign_a;
 418
 419 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 420 #ifndef IEEE_ROUND_NEAREST
 421     rmode = rnd_mode;
 422     if (sign_s && (unsigned) (rmode - 1) < 2)
 423       rmode = 3 - rmode;
 424 #else
 425     rmode = 0;
 426 #endif
 427 #else
 428     rmode = 0;
 429 #endif
 430
 431     // check whether we can take faster path
 432     scale_ca = estimate_decimal_digits[bin_expon_ca];
 433
 434     sign_ab = sign_a ^ sign_b;
 435     sign_ab = ((SINT64) sign_ab) >> 63;
 436
 437     // T1 = 10^(16-diff_dec_expon)
 438     T1 = power10_table_128[16 - diff_dec_expon].w[0];
 439
 440     // get number of digits in coefficient_a
 441     if (coefficient_a >= power10_table_128[scale_ca].w[0]) {
 442       scale_ca++;
 443     }
 444
 445     scale_k = 16 - scale_ca;
 446
 447     // addition
 448     saved_ca = coefficient_a - T1;
 449     coefficient_a =
 450       (SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0];
 451     extra_digits = diff_dec_expon - scale_k;
 452
 453     // apply sign
 454     saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
 455     // add 10^16 and rounding constant
 456     coefficient_b =
 457       saved_cb + 10000000000000000ull +
 458       round_const_table[rmode][extra_digits];
 459
 460     // get P*(2^M[extra_digits])/10^extra_digits
 461     __mul_64x64_to_128 (CT, coefficient_b,
 462                         reciprocals10_64[extra_digits]);
 463
 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;
 467
 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
 477         if (!scale_k) {
 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;
 483           rem_a =
 484             saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
 485           coefficient_a = coefficient_a - T1;
 486
 487           saved_cb += rem_a * power10_table_128[diff_dec_expon].w[0];
 488         } else
 489           coefficient_a =
 490             (SINT64) (saved_ca - T1 -
 491                       (T1 << 3)) * (SINT64) power10_table_128[scale_k -
 492                                                               1].w[0];
 493
 494         extra_digits++;
 495         coefficient_b =
 496           saved_cb + 100000000000000000ull +
 497           round_const_table[rmode][extra_digits];
 498
 499         // get P*(2^M[extra_digits])/10^extra_digits
 500         __mul_64x64_to_128 (CT, coefficient_b,
 501                             reciprocals10_64[extra_digits]);
 502
 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;
 506
 507         // result coefficient
 508         C64 = C0_64 + coefficient_a;
 509       } else if (C64 <= 1000000000000000ull) {
 510         // less than 16 digits in result
 511         coefficient_a =
 512           (SINT64) saved_ca *(SINT64) power10_table_128[scale_k +
 513                                                         1].w[0];
 514         //extra_digits --;
 515         exponent_b--;
 516         coefficient_b =
 517           (saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull +
 518           round_const_table[rmode][extra_digits];
 519
 520         // get P*(2^M[extra_digits])/10^extra_digits
 521         __mul_64x64_to_128 (CT_new, coefficient_b,
 522                             reciprocals10_64[extra_digits]);
 523
 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;
 527
 528         // result coefficient
 529         C64_new = C0_64 + coefficient_a;
 530         if (C64_new < 10000000000000000ull) {
 531           C64 = C64_new;
 532 #ifdef SET_STATUS_FLAGS
 533           CT = CT_new;
 534 #endif
 535         } else
 536           exponent_b++;
 537       }
 538
 539     }
 540
 541   }
 542
 543 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 544 #ifndef IEEE_ROUND_NEAREST
 545   if (rmode == 0)       //ROUNDING_TO_NEAREST
 546 #endif
 547     if (C64 & 1) {
 548       // check whether fractional part of initial_P/10^extra_digits is
 549       // exactly .5
 550       // this is the same as fractional part of
 551       //      (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
 552
 553       // get remainder
 554       remainder_h = CT.w[1] << (64 - amount);
 555
 556       // test whether fractional part is 0
 557       if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
 558         C64--;
 559       }
 560     }
 561 #endif
 562
 563 #ifdef SET_STATUS_FLAGS
 564   status = INEXACT_EXCEPTION;
 565
 566   // get remainder
 567   remainder_h = CT.w[1] << (64 - amount);
 568
 569   switch (rmode) {
 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;
 576     break;
 577   case ROUNDING_DOWN:
 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;
 582     break;
 583   default:
 584     // round up
 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;
 590     break;
 591   }
 592   __set_status_flags (pfpsf, status);
 593
 594 #endif
 595
 596   res =
 597     fast_get_BID64_check_OF (sign_s, exponent_b + extra_digits, C64,
 598                              rnd_mode, pfpsf);
 599   BID_RETURN (res);
 600 }