URL
https://opencores.org/ocsvn/openrisc_2011-10-31/openrisc_2011-10-31/trunk
Subversion Repositories openrisc_2011-10-31
[/] [openrisc/] [tags/] [gnu-src/] [gcc-4.5.1/] [gcc-4.5.1-or32-1.0rc1/] [libgcc/] [config/] [libbid/] [bid64_fma.c] - Rev 338
Compare with Previous | Blame | View Log
/* 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 <http://www.gnu.org/licenses/>. */ /***************************************************************************** * BID64 fma ***************************************************************************** * * Algorithm description: * * if multiplication is guranteed exact (short coefficients) * call the unpacked arg. equivalent of bid64_add(x*y, z) * else * get full coefficient_x*coefficient_y product * call subroutine to perform addition of 64-bit argument * to 128-bit product * ****************************************************************************/ #include "bid_inline_add.h" #if DECIMAL_CALL_BY_REFERENCE extern void bid64_mul (UINT64 * pres, UINT64 * px, UINT64 * py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM); #else extern UINT64 bid64_mul (UINT64 x, UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM); #endif #if DECIMAL_CALL_BY_REFERENCE void bid64_fma (UINT64 * pres, UINT64 * px, UINT64 * py, UINT64 * pz _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { UINT64 x, y, z; #else UINT64 bid64_fma (UINT64 x, UINT64 y, UINT64 z _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { #endif UINT128 P, PU, CT, CZ; UINT64 sign_x, sign_y, coefficient_x, coefficient_y, sign_z, coefficient_z; UINT64 C64, remainder_y, res; UINT64 CYh, CY0L, T, valid_x, valid_y, valid_z; int_double tempx, tempy; int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy, bin_expon_product, rmode; int digits_p, bp, final_exponent, exponent_z, digits_z, ez, ey, scale_z, uf_status; #if DECIMAL_CALL_BY_REFERENCE #if !DECIMAL_GLOBAL_ROUNDING _IDEC_round rnd_mode = *prnd_mode; #endif x = *px; y = *py; z = *pz; #endif valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x); valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y); valid_z = unpack_BID64 (&sign_z, &exponent_z, &coefficient_z, z); // unpack arguments, check for NaN, Infinity, or 0 if (!valid_x || !valid_y || !valid_z) { if ((y & MASK_NAN) == MASK_NAN) { // y is NAN // if x = {0, f, inf, NaN}, y = NaN, z = {0, f, inf, NaN} then res = Q (y) // check first for non-canonical NaN payload y = y & 0xfe03ffffffffffffull; // clear G6-G12 if ((y & 0x0003ffffffffffffull) > 999999999999999ull) { y = y & 0xfe00000000000000ull; // clear G6-G12 and the payload bits } if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN // set invalid flag *pfpsf |= INVALID_EXCEPTION; // return quiet (y) res = y & 0xfdffffffffffffffull; } else { // y is QNaN // return y res = y; // if z = SNaN or x = SNaN signal invalid exception if ((z & MASK_SNAN) == MASK_SNAN || (x & MASK_SNAN) == MASK_SNAN) { // set invalid flag *pfpsf |= INVALID_EXCEPTION; } } BID_RETURN (res) } else if ((z & MASK_NAN) == MASK_NAN) { // z is NAN // if x = {0, f, inf, NaN}, y = {0, f, inf}, z = NaN then res = Q (z) // check first for non-canonical NaN payload z = z & 0xfe03ffffffffffffull; // clear G6-G12 if ((z & 0x0003ffffffffffffull) > 999999999999999ull) { z = z & 0xfe00000000000000ull; // clear G6-G12 and the payload bits } if ((z & MASK_SNAN) == MASK_SNAN) { // z is SNAN // set invalid flag *pfpsf |= INVALID_EXCEPTION; // return quiet (z) res = z & 0xfdffffffffffffffull; } else { // z is QNaN // return z res = z; // if x = SNaN signal invalid exception if ((x & MASK_SNAN) == MASK_SNAN) { // set invalid flag *pfpsf |= INVALID_EXCEPTION; } } BID_RETURN (res) } else if ((x & MASK_NAN) == MASK_NAN) { // x is NAN // if x = NaN, y = {0, f, inf}, z = {0, f, inf} then res = Q (x) // check first for non-canonical NaN payload x = x & 0xfe03ffffffffffffull; // clear G6-G12 if ((x & 0x0003ffffffffffffull) > 999999999999999ull) { x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits } if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNAN // set invalid flag *pfpsf |= INVALID_EXCEPTION; // return quiet (x) res = x & 0xfdffffffffffffffull; } else { // x is QNaN // return x res = x; // clear out G[6]-G[16] } BID_RETURN (res) } if (!valid_x) { // x is Inf. or 0 // x is Infinity? if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) { // check if y is 0 if (!coefficient_y) { // y==0, return NaN #ifdef SET_STATUS_FLAGS if ((z & 0x7e00000000000000ull) != 0x7c00000000000000ull) __set_status_flags (pfpsf, INVALID_EXCEPTION); #endif BID_RETURN (0x7c00000000000000ull); } // test if z is Inf of oposite sign if (((z & 0x7c00000000000000ull) == 0x7800000000000000ull) && (((x ^ y) ^ z) & 0x8000000000000000ull)) { // return NaN #ifdef SET_STATUS_FLAGS __set_status_flags (pfpsf, INVALID_EXCEPTION); #endif BID_RETURN (0x7c00000000000000ull); } // otherwise return +/-Inf BID_RETURN (((x ^ y) & 0x8000000000000000ull) | 0x7800000000000000ull); } // x is 0 if (((y & 0x7800000000000000ull) != 0x7800000000000000ull) && ((z & 0x7800000000000000ull) != 0x7800000000000000ull)) { if (coefficient_z) { exponent_y = exponent_x - DECIMAL_EXPONENT_BIAS + exponent_y; sign_z = z & 0x8000000000000000ull; if (exponent_y >= exponent_z) BID_RETURN (z); res = add_zero64 (exponent_y, sign_z, exponent_z, coefficient_z, &rnd_mode, pfpsf); BID_RETURN (res); } } } if (!valid_y) { // y is Inf. or 0 // y is Infinity? if ((y & 0x7800000000000000ull) == 0x7800000000000000ull) { // check if x is 0 if (!coefficient_x) { // y==0, return NaN #ifdef SET_STATUS_FLAGS __set_status_flags (pfpsf, INVALID_EXCEPTION); #endif BID_RETURN (0x7c00000000000000ull); } // test if z is Inf of oposite sign if (((z & 0x7c00000000000000ull) == 0x7800000000000000ull) && (((x ^ y) ^ z) & 0x8000000000000000ull)) { #ifdef SET_STATUS_FLAGS __set_status_flags (pfpsf, INVALID_EXCEPTION); #endif // return NaN BID_RETURN (0x7c00000000000000ull); } // otherwise return +/-Inf BID_RETURN (((x ^ y) & 0x8000000000000000ull) | 0x7800000000000000ull); } // y is 0 if (((z & 0x7800000000000000ull) != 0x7800000000000000ull)) { if (coefficient_z) { exponent_y += exponent_x - DECIMAL_EXPONENT_BIAS; sign_z = z & 0x8000000000000000ull; if (exponent_y >= exponent_z) BID_RETURN (z); res = add_zero64 (exponent_y, sign_z, exponent_z, coefficient_z, &rnd_mode, pfpsf); BID_RETURN (res); } } } if (!valid_z) { // y is Inf. or 0 // test if y is NaN/Inf if ((z & 0x7800000000000000ull) == 0x7800000000000000ull) { BID_RETURN (coefficient_z & QUIET_MASK64); } // z is 0, return x*y if ((!coefficient_x) || (!coefficient_y)) { //0+/-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; if (exponent_x <= exponent_z) res = ((UINT64) exponent_x) << 53; else res = ((UINT64) exponent_z) << 53; if ((sign_x ^ sign_y) == sign_z) res |= sign_z; #ifndef IEEE_ROUND_NEAREST_TIES_AWAY #ifndef IEEE_ROUND_NEAREST else if (rnd_mode == ROUNDING_DOWN) res |= 0x8000000000000000ull; #endif #endif BID_RETURN (res); } } } /* get binary coefficients of x and y */ //--- 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; final_exponent = exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS; if ((final_exponent > 0) || (!coefficient_z)) { res = get_add64 (sign_x ^ sign_y, final_exponent, C64, sign_z, exponent_z, coefficient_z, rnd_mode, pfpsf); BID_RETURN (res); } else { P.w[0] = C64; P.w[1] = 0; extra_digits = 0; } } else { if (!coefficient_z) { #if DECIMAL_CALL_BY_REFERENCE bid64_mul (&res, px, py _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); #else res = bid64_mul (x, y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); #endif BID_RETURN (res); } // 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; } if (((unsigned) final_exponent) >= 3 * 256) { if (final_exponent < 0) { //--- get number of bits in the coefficients of z --- tempx.d = (double) coefficient_z; bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; // get number of decimal digits in the coeff_x digits_z = estimate_decimal_digits[bin_expon_cx]; if (coefficient_z >= power10_table_128[digits_z].w[0]) digits_z++; // underflow if ((final_exponent + 16 < 0) || (exponent_z + digits_z > 33 + final_exponent)) { res = BID_normalize (sign_z, exponent_z, coefficient_z, sign_x ^ sign_y, 1, rnd_mode, pfpsf); BID_RETURN (res); } ez = exponent_z + digits_z - 16; if (ez < 0) ez = 0; scale_z = exponent_z - ez; coefficient_z *= power10_table_128[scale_z].w[0]; ey = final_exponent - extra_digits; extra_digits = ez - ey; if (extra_digits > 33) { res = BID_normalize (sign_z, exponent_z, coefficient_z, sign_x ^ sign_y, 1, rnd_mode, pfpsf); BID_RETURN (res); } //else // extra_digits<=32 if (extra_digits > 17) { CYh = __truncate (P, 16); // get remainder T = power10_table_128[16].w[0]; __mul_64x64_to_64 (CY0L, CYh, T); remainder_y = P.w[0] - CY0L; extra_digits -= 16; P.w[0] = CYh; P.w[1] = 0; } else remainder_y = 0; // align coeff_x, CYh __mul_64x64_to_128 (CZ, coefficient_z, power10_table_128[extra_digits].w[0]); if (sign_z == (sign_y ^ sign_x)) { __add_128_128 (CT, CZ, P); if (__unsigned_compare_ge_128 (CT, power10_table_128[16 + extra_digits])) { extra_digits++; ez++; } } else { if (remainder_y && (__unsigned_compare_ge_128 (CZ, P))) { P.w[0]++; if (!P.w[0]) P.w[1]++; } __sub_128_128 (CT, CZ, P); if (((SINT64) CT.w[1]) < 0) { sign_z = sign_y ^ sign_x; CT.w[0] = 0 - CT.w[0]; CT.w[1] = 0 - CT.w[1]; if (CT.w[0]) CT.w[1]--; } else if(!(CT.w[1]|CT.w[0])) sign_z = (rnd_mode!=ROUNDING_DOWN)? 0: 0x8000000000000000ull; if (ez && (__unsigned_compare_gt_128 (power10_table_128[15 + extra_digits], CT))) { extra_digits--; ez--; } } #ifdef SET_STATUS_FLAGS uf_status = 0; if ((!ez) && __unsigned_compare_gt_128 (power10_table_128 [extra_digits + 15], CT)) { rmode = rnd_mode; if (sign_z && (unsigned) (rmode - 1) < 2) rmode = 3 - rmode; //__add_128_64(PU, CT, round_const_table[rmode][extra_digits]); PU = power10_table_128[extra_digits + 15]; PU.w[0]--; if (__unsigned_compare_gt_128 (PU, CT) || (rmode == ROUNDING_DOWN) || (rmode == ROUNDING_TO_ZERO)) uf_status = UNDERFLOW_EXCEPTION; else if (extra_digits < 2) { if ((rmode == ROUNDING_UP)) { if (!extra_digits) uf_status = UNDERFLOW_EXCEPTION; else { if (remainder_y && (sign_z != (sign_y ^ sign_x))) remainder_y = power10_table_128[16].w[0] - remainder_y; if (power10_table_128[15].w[0] > remainder_y) uf_status = UNDERFLOW_EXCEPTION; } } else // RN or RN_away { if (remainder_y && (sign_z != (sign_y ^ sign_x))) remainder_y = power10_table_128[16].w[0] - remainder_y; if (!extra_digits) { remainder_y += round_const_table[rmode][15]; if (remainder_y < power10_table_128[16].w[0]) uf_status = UNDERFLOW_EXCEPTION; } else { if (remainder_y < round_const_table[rmode][16]) uf_status = UNDERFLOW_EXCEPTION; } } //__set_status_flags (pfpsf, uf_status); } } #endif res = __bid_full_round64_remainder (sign_z, ez - extra_digits, CT, extra_digits, remainder_y, rnd_mode, pfpsf, uf_status); BID_RETURN (res); } else { if ((sign_z == (sign_x ^ sign_y)) || (final_exponent > 3 * 256 + 15)) { res = fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent, 1000000000000000ull, rnd_mode, pfpsf); BID_RETURN (res); } } } if (extra_digits > 0) { res = get_add128 (sign_z, exponent_z, coefficient_z, sign_x ^ sign_y, final_exponent, P, extra_digits, rnd_mode, pfpsf); BID_RETURN (res); } // go to convert_format and exit else { C64 = __low_64 (P); res = get_add64 (sign_x ^ sign_y, exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64, sign_z, exponent_z, coefficient_z, rnd_mode, pfpsf); BID_RETURN (res); } }