| /* Copyright (C) 2007-2017 Free Software Foundation, Inc. |
| |
| This file is part of GCC. |
| |
| GCC is free software; you can redistribute it and/or modify it under |
| the terms of the GNU General Public License as published by the Free |
| Software Foundation; either version 3, or (at your option) any later |
| version. |
| |
| GCC is distributed in the hope that it will be useful, but WITHOUT ANY |
| WARRANTY; without even the implied warranty of MERCHANTABILITY or |
| FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
| for more details. |
| |
| Under Section 7 of GPL version 3, you are granted additional |
| permissions described in the GCC Runtime Library Exception, version |
| 3.1, as published by the Free Software Foundation. |
| |
| You should have received a copy of the GNU General Public License and |
| a copy of the GCC Runtime Library Exception along with this program; |
| see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |
| <http://www.gnu.org/licenses/>. */ |
| |
| /***************************************************************************** |
| * BID64 multiply |
| ***************************************************************************** |
| * |
| * Algorithm description: |
| * |
| * if(number_digits(coefficient_x)+number_digits(coefficient_y) guaranteed |
| * below 16) |
| * return get_BID64(sign_x^sign_y, exponent_x + exponent_y - dec_bias, |
| * coefficient_x*coefficient_y) |
| * else |
| * get long product: coefficient_x*coefficient_y |
| * determine number of digits to round off (extra_digits) |
| * rounding is performed as a 128x128-bit multiplication by |
| * 2^M[extra_digits]/10^extra_digits, followed by a shift |
| * M[extra_digits] is sufficiently large for required accuracy |
| * |
| ****************************************************************************/ |
| |
| #include "bid_internal.h" |
| |
| #if DECIMAL_CALL_BY_REFERENCE |
| |
| void |
| bid64_mul (UINT64 * pres, UINT64 * px, |
| UINT64 * |
| py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
| _EXC_INFO_PARAM) { |
| UINT64 x, y; |
| #else |
| |
| UINT64 |
| bid64_mul (UINT64 x, |
| UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM |
| _EXC_MASKS_PARAM _EXC_INFO_PARAM) { |
| #endif |
| UINT128 P, PU, C128, Q_high, Q_low, Stemp; |
| UINT64 sign_x, sign_y, coefficient_x, coefficient_y; |
| UINT64 C64, remainder_h, carry, CY, res; |
| UINT64 valid_x, valid_y; |
| int_double tempx, tempy; |
| int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy, |
| bin_expon_product; |
| int rmode, digits_p, bp, amount, amount2, final_exponent, round_up; |
| unsigned status, uf_status; |
| |
| #if DECIMAL_CALL_BY_REFERENCE |
| #if !DECIMAL_GLOBAL_ROUNDING |
| _IDEC_round rnd_mode = *prnd_mode; |
| #endif |
| x = *px; |
| y = *py; |
| #endif |
| |
| valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x); |
| valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y); |
| |
| // unpack arguments, check for NaN or Infinity |
| if (!valid_x) { |
| |
| #ifdef SET_STATUS_FLAGS |
| if ((y & SNAN_MASK64) == SNAN_MASK64) // y is sNaN |
| __set_status_flags (pfpsf, INVALID_EXCEPTION); |
| #endif |
| // x is Inf. or NaN |
| |
| // test if x is NaN |
| if ((x & NAN_MASK64) == NAN_MASK64) { |
| #ifdef SET_STATUS_FLAGS |
| if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN |
| __set_status_flags (pfpsf, INVALID_EXCEPTION); |
| #endif |
| BID_RETURN (coefficient_x & QUIET_MASK64); |
| } |
| // x is Infinity? |
| if ((x & INFINITY_MASK64) == INFINITY_MASK64) { |
| // check if y is 0 |
| if (((y & INFINITY_MASK64) != INFINITY_MASK64) |
| && !coefficient_y) { |
| #ifdef SET_STATUS_FLAGS |
| __set_status_flags (pfpsf, INVALID_EXCEPTION); |
| #endif |
| // y==0 , return NaN |
| BID_RETURN (NAN_MASK64); |
| } |
| // check if y is NaN |
| if ((y & NAN_MASK64) == NAN_MASK64) |
| // y==NaN , return NaN |
| BID_RETURN (coefficient_y & QUIET_MASK64); |
| // otherwise return +/-Inf |
| BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64); |
| } |
| // x is 0 |
| if (((y & INFINITY_MASK64) != INFINITY_MASK64)) { |
| if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64) |
| exponent_y = ((UINT32) (y >> 51)) & 0x3ff; |
| else |
| exponent_y = ((UINT32) (y >> 53)) & 0x3ff; |
| sign_y = y & 0x8000000000000000ull; |
| |
| exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS; |
| if (exponent_x > DECIMAL_MAX_EXPON_64) |
| exponent_x = DECIMAL_MAX_EXPON_64; |
| else if (exponent_x < 0) |
| exponent_x = 0; |
| BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53)); |
| } |
| } |
| if (!valid_y) { |
| // y is Inf. or NaN |
| |
| // test if y is NaN |
| if ((y & NAN_MASK64) == NAN_MASK64) { |
| #ifdef SET_STATUS_FLAGS |
| if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN |
| __set_status_flags (pfpsf, INVALID_EXCEPTION); |
| #endif |
| BID_RETURN (coefficient_y & QUIET_MASK64); |
| } |
| // y is Infinity? |
| if ((y & INFINITY_MASK64) == INFINITY_MASK64) { |
| // check if x is 0 |
| if (!coefficient_x) { |
| __set_status_flags (pfpsf, INVALID_EXCEPTION); |
| // x==0, return NaN |
| BID_RETURN (NAN_MASK64); |
| } |
| // otherwise return +/-Inf |
| BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64); |
| } |
| // y is 0 |
| exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS; |
| if (exponent_x > DECIMAL_MAX_EXPON_64) |
| exponent_x = DECIMAL_MAX_EXPON_64; |
| else if (exponent_x < 0) |
| exponent_x = 0; |
| BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53)); |
| } |
| //--- get number of bits in the coefficients of x and y --- |
| // version 2 (original) |
| tempx.d = (double) coefficient_x; |
| bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52); |
| tempy.d = (double) coefficient_y; |
| bin_expon_cy = ((tempy.i & MASK_BINARY_EXPONENT) >> 52); |
| |
| // magnitude estimate for coefficient_x*coefficient_y is |
| // 2^(unbiased_bin_expon_cx + unbiased_bin_expon_cx) |
| bin_expon_product = bin_expon_cx + bin_expon_cy; |
| |
| // check if coefficient_x*coefficient_y<2^(10*k+3) |
| // equivalent to unbiased_bin_expon_cx + unbiased_bin_expon_cx < 10*k+1 |
| if (bin_expon_product < UPPER_EXPON_LIMIT + 2 * BINARY_EXPONENT_BIAS) { |
| // easy multiply |
| C64 = coefficient_x * coefficient_y; |
| |
| res = |
| get_BID64_small_mantissa (sign_x ^ sign_y, |
| exponent_x + exponent_y - |
| DECIMAL_EXPONENT_BIAS, C64, rnd_mode, |
| pfpsf); |
| BID_RETURN (res); |
| } else { |
| uf_status = 0; |
| // get 128-bit product: coefficient_x*coefficient_y |
| __mul_64x64_to_128 (P, coefficient_x, coefficient_y); |
| |
| // tighten binary range of P: leading bit is 2^bp |
| // unbiased_bin_expon_product <= bp <= unbiased_bin_expon_product+1 |
| bin_expon_product -= 2 * BINARY_EXPONENT_BIAS; |
| |
| __tight_bin_range_128 (bp, P, bin_expon_product); |
| |
| // get number of decimal digits in the product |
| digits_p = estimate_decimal_digits[bp]; |
| if (!(__unsigned_compare_gt_128 (power10_table_128[digits_p], P))) |
| digits_p++; // if power10_table_128[digits_p] <= P |
| |
| // determine number of decimal digits to be rounded out |
| extra_digits = digits_p - MAX_FORMAT_DIGITS; |
| final_exponent = |
| exponent_x + exponent_y + extra_digits - DECIMAL_EXPONENT_BIAS; |
| |
| #ifndef IEEE_ROUND_NEAREST_TIES_AWAY |
| #ifndef IEEE_ROUND_NEAREST |
| rmode = rnd_mode; |
| if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2) |
| rmode = 3 - rmode; |
| #else |
| rmode = 0; |
| #endif |
| #else |
| rmode = 0; |
| #endif |
| |
| round_up = 0; |
| if (((unsigned) final_exponent) >= 3 * 256) { |
| if (final_exponent < 0) { |
| // underflow |
| if (final_exponent + 16 < 0) { |
| res = sign_x ^ sign_y; |
| __set_status_flags (pfpsf, |
| UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION); |
| if (rmode == ROUNDING_UP) |
| res |= 1; |
| BID_RETURN (res); |
| } |
| |
| uf_status = UNDERFLOW_EXCEPTION; |
| if (final_exponent == -1) { |
| __add_128_64 (PU, P, round_const_table[rmode][extra_digits]); |
| if (__unsigned_compare_ge_128 |
| (PU, power10_table_128[extra_digits + 16])) |
| uf_status = 0; |
| } |
| extra_digits -= final_exponent; |
| final_exponent = 0; |
| |
| if (extra_digits > 17) { |
| __mul_128x128_full (Q_high, Q_low, P, reciprocals10_128[16]); |
| |
| amount = recip_scale[16]; |
| __shr_128 (P, Q_high, amount); |
| |
| // get sticky bits |
| amount2 = 64 - amount; |
| remainder_h = 0; |
| remainder_h--; |
| remainder_h >>= amount2; |
| remainder_h = remainder_h & Q_high.w[0]; |
| |
| extra_digits -= 16; |
| if (remainder_h || (Q_low.w[1] > reciprocals10_128[16].w[1] |
| || (Q_low.w[1] == |
| reciprocals10_128[16].w[1] |
| && Q_low.w[0] >= |
| reciprocals10_128[16].w[0]))) { |
| round_up = 1; |
| __set_status_flags (pfpsf, |
| UNDERFLOW_EXCEPTION | |
| INEXACT_EXCEPTION); |
| P.w[0] = (P.w[0] << 3) + (P.w[0] << 1); |
| P.w[0] |= 1; |
| extra_digits++; |
| } |
| } |
| } else { |
| res = |
| fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent, |
| 1000000000000000ull, rnd_mode, |
| pfpsf); |
| BID_RETURN (res); |
| } |
| } |
| |
| |
| if (extra_digits > 0) { |
| // will divide by 10^(digits_p - 16) |
| |
| // add a constant to P, depending on rounding mode |
| // 0.5*10^(digits_p - 16) for round-to-nearest |
| __add_128_64 (P, P, round_const_table[rmode][extra_digits]); |
| |
| // get P*(2^M[extra_digits])/10^extra_digits |
| __mul_128x128_full (Q_high, Q_low, P, |
| reciprocals10_128[extra_digits]); |
| |
| // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 |
| amount = recip_scale[extra_digits]; |
| __shr_128 (C128, Q_high, amount); |
| |
| C64 = __low_64 (C128); |
| |
| #ifndef IEEE_ROUND_NEAREST_TIES_AWAY |
| #ifndef IEEE_ROUND_NEAREST |
| if (rmode == 0) //ROUNDING_TO_NEAREST |
| #endif |
| if ((C64 & 1) && !round_up) { |
| // check whether fractional part of initial_P/10^extra_digits |
| // is exactly .5 |
| // this is the same as fractional part of |
| // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero |
| |
| // get remainder |
| remainder_h = Q_high.w[0] << (64 - amount); |
| |
| // test whether fractional part is 0 |
| if (!remainder_h |
| && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] |
| || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] |
| && Q_low.w[0] < |
| reciprocals10_128[extra_digits].w[0]))) { |
| C64--; |
| } |
| } |
| #endif |
| |
| #ifdef SET_STATUS_FLAGS |
| status = INEXACT_EXCEPTION | uf_status; |
| |
| // get remainder |
| remainder_h = Q_high.w[0] << (64 - amount); |
| |
| switch (rmode) { |
| case ROUNDING_TO_NEAREST: |
| case ROUNDING_TIES_AWAY: |
| // test whether fractional part is 0 |
| if (remainder_h == 0x8000000000000000ull |
| && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] |
| || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] |
| && Q_low.w[0] < |
| reciprocals10_128[extra_digits].w[0]))) |
| status = EXACT_STATUS; |
| break; |
| case ROUNDING_DOWN: |
| case ROUNDING_TO_ZERO: |
| if (!remainder_h |
| && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] |
| || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] |
| && Q_low.w[0] < |
| reciprocals10_128[extra_digits].w[0]))) |
| status = EXACT_STATUS; |
| break; |
| default: |
| // round up |
| __add_carry_out (Stemp.w[0], CY, Q_low.w[0], |
| reciprocals10_128[extra_digits].w[0]); |
| __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], |
| reciprocals10_128[extra_digits].w[1], CY); |
| if ((remainder_h >> (64 - amount)) + carry >= |
| (((UINT64) 1) << amount)) |
| status = EXACT_STATUS; |
| } |
| |
| __set_status_flags (pfpsf, status); |
| #endif |
| |
| // convert to BID and return |
| res = |
| fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent, C64, |
| rmode, pfpsf); |
| BID_RETURN (res); |
| } |
| // go to convert_format and exit |
| C64 = __low_64 (P); |
| res = |
| get_BID64 (sign_x ^ sign_y, |
| exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64, |
| rmode, pfpsf); |
| BID_RETURN (res); |
| } |
| } |