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