/* Copyright (C) 2007 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 2, or (at your option) any later version. In addition to the permissions in the GNU General Public License, the Free Software Foundation gives you unlimited permission to link the compiled version of this file into combinations with other programs, and to distribute those combinations without any restriction coming from the use of this file. (The General Public License restrictions do apply in other respects; for example, they cover modification of the file, and distribution when not linked into a combine executable.) 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. You should have received a copy of the GNU General Public License along with GCC; see the file COPYING. If not, write to the Free Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ /***************************************************************************** * 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 = 0, exponent_y = 0, 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 (x & QUIET_MASK64); } // x is Infinity? if ((x & INFINITY_MASK64) == INFINITY_MASK64) { // check if y is 0 if (((y & SPECIAL_ENCODING_MASK64) != SPECIAL_ENCODING_MASK64) && !(y << (64 - 53))) { #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 (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 (y & QUIET_MASK64); } // y is Infinity? if ((y & INFINITY_MASK64) == INFINITY_MASK64) { // check if x is 0 if (((x & SPECIAL_ENCODING_MASK64) != SPECIAL_ENCODING_MASK64) && !(x << (64 - 53))) { __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 = __bid_estimate_decimal_digits[bp]; if (!(__unsigned_compare_gt_128 (__bid_power10_table_128[digits_p], P))) digits_p++; // if __bid_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, __bid_round_const_table[rmode][extra_digits]); if (__unsigned_compare_ge_128 (PU, __bid_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, __bid_reciprocals10_128[16]); amount = __bid_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] > __bid_reciprocals10_128[16].w[1] || (Q_low.w[1] == __bid_reciprocals10_128[16].w[1] && Q_low.w[0] >= __bid_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, __bid_round_const_table[rmode][extra_digits]); // get P*(2^M[extra_digits])/10^extra_digits __mul_128x128_full (Q_high, Q_low, P, __bid_reciprocals10_128[extra_digits]); // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 amount = __bid_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] < __bid_reciprocals10_128[extra_digits].w[1] || (Q_low.w[1] == __bid_reciprocals10_128[extra_digits].w[1] && Q_low.w[0] < __bid_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] < __bid_reciprocals10_128[extra_digits].w[1] || (Q_low.w[1] == __bid_reciprocals10_128[extra_digits].w[1] && Q_low.w[0] < __bid_reciprocals10_128[extra_digits].w[0]))) status = EXACT_STATUS; break; case ROUNDING_DOWN: case ROUNDING_TO_ZERO: if (!remainder_h && (Q_low.w[1] < __bid_reciprocals10_128[extra_digits].w[1] || (Q_low.w[1] == __bid_reciprocals10_128[extra_digits].w[1] && Q_low.w[0] < __bid_reciprocals10_128[extra_digits].w[0]))) status = EXACT_STATUS; break; default: // round up __add_carry_out (Stemp.w[0], CY, Q_low.w[0], __bid_reciprocals10_128[extra_digits].w[0]); __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], __bid_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); } }