diff options
Diffstat (limited to 'gcc-4.9/libgcc/config/libbid/bid64_mul.c')
| -rw-r--r-- | gcc-4.9/libgcc/config/libbid/bid64_mul.c | 374 |
1 files changed, 374 insertions, 0 deletions
diff --git a/gcc-4.9/libgcc/config/libbid/bid64_mul.c b/gcc-4.9/libgcc/config/libbid/bid64_mul.c new file mode 100644 index 000000000..213fc68eb --- /dev/null +++ b/gcc-4.9/libgcc/config/libbid/bid64_mul.c @@ -0,0 +1,374 @@ +/* Copyright (C) 2007-2014 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); + } +} |