libgcc/config/libbid/bid64_mul.c

   1 /* Copyright (C) 2007-2024 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 3, or (at your option) any later
   8 version.
   9
  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
  13 for more details.
  14
  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.
  18
  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/>.  */
  23
  24 /*****************************************************************************
  25  *    BID64 multiply
  26  *****************************************************************************
  27  *
  28  *  Algorithm description:
  29  *
  30  *  if(number_digits(coefficient_x)+number_digits(coefficient_y) guaranteed
  31  *       below 16)
  32  *      return get_BID64(sign_x^sign_y, exponent_x + exponent_y - dec_bias,
  33  *                     coefficient_x*coefficient_y)
  34  *  else
  35  *      get long product: coefficient_x*coefficient_y
  36  *      determine number of digits to round off (extra_digits)
  37  *      rounding is performed as a 128x128-bit multiplication by
  38  *         2^M[extra_digits]/10^extra_digits, followed by a shift
  39  *         M[extra_digits] is sufficiently large for required accuracy
  40  *
  41  ****************************************************************************/
  42
  43 #include "bid_internal.h"
  44
  45 #if DECIMAL_CALL_BY_REFERENCE
  46
  47 void
  48 bid64_mul (UINT64 * pres, UINT64 * px,
  49            UINT64 *
  50            py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  51            _EXC_INFO_PARAM) {
  52   UINT64 x, y;
  53 #else
  54
  55 UINT64
  56 bid64_mul (UINT64 x,
  57            UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
  58            _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  59 #endif
  60   UINT128 P, PU, C128, Q_high, Q_low, Stemp;
  61   UINT64 sign_x, sign_y, coefficient_x, coefficient_y;
  62   UINT64 C64, remainder_h, carry, CY, res;
  63   UINT64 valid_x, valid_y;
  64   int_double tempx, tempy;
  65   int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy,
  66     bin_expon_product;
  67   int rmode, digits_p, bp, amount, amount2, final_exponent, round_up;
  68   unsigned status, uf_status;
  69
  70 #if DECIMAL_CALL_BY_REFERENCE
  71 #if !DECIMAL_GLOBAL_ROUNDING
  72   _IDEC_round rnd_mode = *prnd_mode;
  73 #endif
  74   x = *px;
  75   y = *py;
  76 #endif
  77
  78   valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
  79   valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
  80
  81   // unpack arguments, check for NaN or Infinity
  82   if (!valid_x) {
  83
  84 #ifdef SET_STATUS_FLAGS
  85     if ((y & SNAN_MASK64) == SNAN_MASK64)       // y is sNaN
  86       __set_status_flags (pfpsf, INVALID_EXCEPTION);
  87 #endif
  88     // x is Inf. or NaN
  89
  90     // test if x is NaN
  91     if ((x & NAN_MASK64) == NAN_MASK64) {
  92 #ifdef SET_STATUS_FLAGS
  93       if ((x & SNAN_MASK64) == SNAN_MASK64)     // sNaN
  94         __set_status_flags (pfpsf, INVALID_EXCEPTION);
  95 #endif
  96       BID_RETURN (coefficient_x & QUIET_MASK64);
  97     }
  98     // x is Infinity?
  99     if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
 100       // check if y is 0
 101       if (((y & INFINITY_MASK64) != INFINITY_MASK64)
 102           && !coefficient_y) {
 103 #ifdef SET_STATUS_FLAGS
 104         __set_status_flags (pfpsf, INVALID_EXCEPTION);
 105 #endif
 106         // y==0 , return NaN
 107         BID_RETURN (NAN_MASK64);
 108       }
 109       // check if y is NaN
 110       if ((y & NAN_MASK64) == NAN_MASK64)
 111         // y==NaN , return NaN
 112         BID_RETURN (coefficient_y & QUIET_MASK64);
 113       // otherwise return +/-Inf
 114       BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64);
 115     }
 116     // x is 0
 117     if (((y & INFINITY_MASK64) != INFINITY_MASK64)) {
 118       if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64)
 119         exponent_y = ((UINT32) (y >> 51)) & 0x3ff;
 120       else
 121         exponent_y = ((UINT32) (y >> 53)) & 0x3ff;
 122       sign_y = y & 0x8000000000000000ull;
 123
 124       exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS;
 125       if (exponent_x > DECIMAL_MAX_EXPON_64)
 126         exponent_x = DECIMAL_MAX_EXPON_64;
 127       else if (exponent_x < 0)
 128         exponent_x = 0;
 129       BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53));
 130     }
 131   }
 132   if (!valid_y) {
 133     // y is Inf. or NaN
 134
 135     // test if y is NaN
 136     if ((y & NAN_MASK64) == NAN_MASK64) {
 137 #ifdef SET_STATUS_FLAGS
 138       if ((y & SNAN_MASK64) == SNAN_MASK64)     // sNaN
 139         __set_status_flags (pfpsf, INVALID_EXCEPTION);
 140 #endif
 141       BID_RETURN (coefficient_y & QUIET_MASK64);
 142     }
 143     // y is Infinity?
 144     if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
 145       // check if x is 0
 146       if (!coefficient_x) {
 147         __set_status_flags (pfpsf, INVALID_EXCEPTION);
 148         // x==0, return NaN
 149         BID_RETURN (NAN_MASK64);
 150       }
 151       // otherwise return +/-Inf
 152       BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64);
 153     }
 154     // y is 0
 155     exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS;
 156     if (exponent_x > DECIMAL_MAX_EXPON_64)
 157       exponent_x = DECIMAL_MAX_EXPON_64;
 158     else if (exponent_x < 0)
 159       exponent_x = 0;
 160     BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53));
 161   }
 162   //--- get number of bits in the coefficients of x and y ---
 163   // version 2 (original)
 164   tempx.d = (double) coefficient_x;
 165   bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52);
 166   tempy.d = (double) coefficient_y;
 167   bin_expon_cy = ((tempy.i & MASK_BINARY_EXPONENT) >> 52);
 168
 169   // magnitude estimate for coefficient_x*coefficient_y is
 170   //        2^(unbiased_bin_expon_cx + unbiased_bin_expon_cx)
 171   bin_expon_product = bin_expon_cx + bin_expon_cy;
 172
 173   // check if coefficient_x*coefficient_y<2^(10*k+3)
 174   // equivalent to unbiased_bin_expon_cx + unbiased_bin_expon_cx < 10*k+1
 175   if (bin_expon_product < UPPER_EXPON_LIMIT + 2 * BINARY_EXPONENT_BIAS) {
 176     //  easy multiply
 177     C64 = coefficient_x * coefficient_y;
 178
 179     res =
 180       get_BID64_small_mantissa (sign_x ^ sign_y,
 181                                 exponent_x + exponent_y -
 182                                 DECIMAL_EXPONENT_BIAS, C64, rnd_mode,
 183                                 pfpsf);
 184     BID_RETURN (res);
 185   } else {
 186     uf_status = 0;
 187     // get 128-bit product: coefficient_x*coefficient_y
 188     __mul_64x64_to_128 (P, coefficient_x, coefficient_y);
 189
 190     // tighten binary range of P:  leading bit is 2^bp
 191     // unbiased_bin_expon_product <= bp <= unbiased_bin_expon_product+1
 192     bin_expon_product -= 2 * BINARY_EXPONENT_BIAS;
 193
 194     __tight_bin_range_128 (bp, P, bin_expon_product);
 195
 196     // get number of decimal digits in the product
 197     digits_p = estimate_decimal_digits[bp];
 198     if (!(__unsigned_compare_gt_128 (power10_table_128[digits_p], P)))
 199       digits_p++;       // if power10_table_128[digits_p] <= P
 200
 201     // determine number of decimal digits to be rounded out
 202     extra_digits = digits_p - MAX_FORMAT_DIGITS;
 203     final_exponent =
 204       exponent_x + exponent_y + extra_digits - DECIMAL_EXPONENT_BIAS;
 205
 206 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 207 #ifndef IEEE_ROUND_NEAREST
 208     rmode = rnd_mode;
 209     if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
 210       rmode = 3 - rmode;
 211 #else
 212     rmode = 0;
 213 #endif
 214 #else
 215     rmode = 0;
 216 #endif
 217
 218     round_up = 0;
 219     if (((unsigned) final_exponent) >= 3 * 256) {
 220       if (final_exponent < 0) {
 221         // underflow
 222         if (final_exponent + 16 < 0) {
 223           res = sign_x ^ sign_y;
 224           __set_status_flags (pfpsf,
 225                               UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION);
 226           if (rmode == ROUNDING_UP)
 227             res |= 1;
 228           BID_RETURN (res);
 229         }
 230
 231         uf_status = UNDERFLOW_EXCEPTION;
 232         if (final_exponent == -1) {
 233           __add_128_64 (PU, P, round_const_table[rmode][extra_digits]);
 234           if (__unsigned_compare_ge_128
 235               (PU, power10_table_128[extra_digits + 16]))
 236             uf_status = 0;
 237         }
 238         extra_digits -= final_exponent;
 239         final_exponent = 0;
 240
 241         if (extra_digits > 17) {
 242           __mul_128x128_full (Q_high, Q_low, P, reciprocals10_128[16]);
 243
 244           amount = recip_scale[16];
 245           __shr_128 (P, Q_high, amount);
 246
 247           // get sticky bits
 248           amount2 = 64 - amount;
 249           remainder_h = 0;
 250           remainder_h--;
 251           remainder_h >>= amount2;
 252           remainder_h = remainder_h & Q_high.w[0];
 253
 254           extra_digits -= 16;
 255           if (remainder_h || (Q_low.w[1] > reciprocals10_128[16].w[1]
 256                               || (Q_low.w[1] ==
 257                                   reciprocals10_128[16].w[1]
 258                                   && Q_low.w[0] >=
 259                                   reciprocals10_128[16].w[0]))) {
 260             round_up = 1;
 261             __set_status_flags (pfpsf,
 262                                 UNDERFLOW_EXCEPTION |
 263                                 INEXACT_EXCEPTION);
 264             P.w[0] = (P.w[0] << 3) + (P.w[0] << 1);
 265             P.w[0] |= 1;
 266             extra_digits++;
 267           }
 268         }
 269       } else {
 270         res =
 271           fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent,
 272                                    1000000000000000ull, rnd_mode,
 273                                    pfpsf);
 274         BID_RETURN (res);
 275       }
 276     }
 277
 278
 279     if (extra_digits > 0) {
 280       // will divide by 10^(digits_p - 16)
 281
 282       // add a constant to P, depending on rounding mode
 283       // 0.5*10^(digits_p - 16) for round-to-nearest
 284       __add_128_64 (P, P, round_const_table[rmode][extra_digits]);
 285
 286       // get P*(2^M[extra_digits])/10^extra_digits
 287       __mul_128x128_full (Q_high, Q_low, P,
 288                           reciprocals10_128[extra_digits]);
 289
 290       // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
 291       amount = recip_scale[extra_digits];
 292       __shr_128 (C128, Q_high, amount);
 293
 294       C64 = __low_64 (C128);
 295
 296 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 297 #ifndef IEEE_ROUND_NEAREST
 298       if (rmode == 0)   //ROUNDING_TO_NEAREST
 299 #endif
 300         if ((C64 & 1) && !round_up) {
 301           // check whether fractional part of initial_P/10^extra_digits
 302           // is exactly .5
 303           // this is the same as fractional part of
 304           // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
 305
 306           // get remainder
 307           remainder_h = Q_high.w[0] << (64 - amount);
 308
 309           // test whether fractional part is 0
 310           if (!remainder_h
 311               && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
 312                   || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
 313                       && Q_low.w[0] <
 314                       reciprocals10_128[extra_digits].w[0]))) {
 315             C64--;
 316           }
 317         }
 318 #endif
 319
 320 #ifdef SET_STATUS_FLAGS
 321       status = INEXACT_EXCEPTION | uf_status;
 322
 323       // get remainder
 324       remainder_h = Q_high.w[0] << (64 - amount);
 325
 326       switch (rmode) {
 327       case ROUNDING_TO_NEAREST:
 328       case ROUNDING_TIES_AWAY:
 329         // test whether fractional part is 0
 330         if (remainder_h == 0x8000000000000000ull
 331             && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
 332                 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
 333                     && Q_low.w[0] <
 334                     reciprocals10_128[extra_digits].w[0])))
 335           status = EXACT_STATUS;
 336         break;
 337       case ROUNDING_DOWN:
 338       case ROUNDING_TO_ZERO:
 339         if (!remainder_h
 340             && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
 341                 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
 342                     && Q_low.w[0] <
 343                     reciprocals10_128[extra_digits].w[0])))
 344           status = EXACT_STATUS;
 345         break;
 346       default:
 347         // round up
 348         __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
 349                          reciprocals10_128[extra_digits].w[0]);
 350         __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
 351                             reciprocals10_128[extra_digits].w[1], CY);
 352         if ((remainder_h >> (64 - amount)) + carry >=
 353             (((UINT64) 1) << amount))
 354           status = EXACT_STATUS;
 355       }
 356
 357       __set_status_flags (pfpsf, status);
 358 #endif
 359
 360       // convert to BID and return
 361       res =
 362         fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent, C64,
 363                                  rmode, pfpsf);
 364       BID_RETURN (res);
 365     }
 366     // go to convert_format and exit
 367     C64 = __low_64 (P);
 368     res =
 369       get_BID64 (sign_x ^ sign_y,
 370                  exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64,
 371                  rmode, pfpsf);
 372     BID_RETURN (res);
 373   }
 374 }