/* Copyright (C) 2007, 2009 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 . */ /***************************************************************************** * 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); } }