| |
| /*---------------------------------------------------------------*/ |
| /*--- begin host_generic_maddf.c ---*/ |
| /*---------------------------------------------------------------*/ |
| |
| /* |
| Compute x * y + z as ternary operation. |
| Copyright (C) 2010-2013 Free Software Foundation, Inc. |
| This file is part of the GNU C Library. |
| Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. |
| |
| The GNU C Library is free software; you can redistribute it and/or |
| modify it under the terms of the GNU Lesser General Public |
| License as published by the Free Software Foundation; either |
| version 2.1 of the License, or (at your option) any later version. |
| |
| The GNU C Library 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 |
| Lesser General Public License for more details. |
| |
| You should have received a copy of the GNU Lesser General Public |
| License along with the GNU C Library; if not, see |
| <http://www.gnu.org/licenses/>. |
| */ |
| |
| /* Generic helper functions for doing FMA, i.e. compute x * y + z |
| as ternary operation. |
| These are purely back-end entities and cannot be seen/referenced |
| from IR. */ |
| |
| #include "libvex_basictypes.h" |
| #include "host_generic_maddf.h" |
| #include "main_util.h" |
| |
| /* This implementation relies on Double being more than twice as |
| precise as Float and uses rounding to odd in order to avoid problems |
| with double rounding. |
| See a paper by Boldo and Melquiond: |
| http://www.lri.fr/~melquion/doc/08-tc.pdf */ |
| |
| #define FORCE_EVAL(X) __asm __volatile__ ("" : : "m" (X)) |
| |
| #if defined(__x86_64__) && defined(__SSE2_MATH__) |
| # define ENV_TYPE unsigned int |
| /* Save current rounding mode into ENV, hold exceptions, set rounding |
| mode to rounding toward zero. */ |
| # define ROUNDTOZERO(env) \ |
| do { \ |
| unsigned int mxcsr; \ |
| __asm __volatile__ ("stmxcsr %0" : "=m" (mxcsr)); \ |
| (env) = mxcsr; \ |
| mxcsr = (mxcsr | 0x7f80) & ~0x3f; \ |
| __asm __volatile__ ("ldmxcsr %0" : : "m" (mxcsr));\ |
| } while (0) |
| /* Restore exceptions from ENV, return if inexact exception has been raised |
| since ROUNDTOZERO. */ |
| # define RESET_TESTINEXACT(env) \ |
| ({ \ |
| unsigned int mxcsr, ret; \ |
| __asm __volatile__ ("stmxcsr %0" : "=m" (mxcsr)); \ |
| ret = (mxcsr >> 5) & 1; \ |
| mxcsr = (mxcsr & 0x3d) | (env); \ |
| __asm __volatile__ ("ldmxcsr %0" : : "m" (mxcsr));\ |
| ret; \ |
| }) |
| /* Return if inexact exception has been raised since ROUNDTOZERO. */ |
| # define TESTINEXACT() \ |
| ({ \ |
| unsigned int mxcsr; \ |
| __asm __volatile__ ("stmxcsr %0" : "=m" (mxcsr)); \ |
| (mxcsr >> 5) & 1; \ |
| }) |
| #endif |
| |
| #define DBL_MANT_DIG 53 |
| #define IEEE754_DOUBLE_BIAS 0x3ff |
| |
| union vg_ieee754_double { |
| Double d; |
| |
| /* This is the IEEE 754 double-precision format. */ |
| struct { |
| #ifdef VKI_BIG_ENDIAN |
| unsigned int negative:1; |
| unsigned int exponent:11; |
| unsigned int mantissa0:20; |
| unsigned int mantissa1:32; |
| #else |
| unsigned int mantissa1:32; |
| unsigned int mantissa0:20; |
| unsigned int exponent:11; |
| unsigned int negative:1; |
| #endif |
| } ieee; |
| }; |
| |
| void VEX_REGPARM(3) |
| h_generic_calc_MAddF32 ( /*OUT*/Float* res, |
| Float* argX, Float* argY, Float* argZ ) |
| { |
| #ifndef ENV_TYPE |
| /* Lame fallback implementation. */ |
| *res = *argX * *argY + *argZ; |
| #else |
| ENV_TYPE env; |
| /* Multiplication is always exact. */ |
| Double temp = (Double) *argX * (Double) *argY; |
| union vg_ieee754_double u; |
| |
| ROUNDTOZERO (env); |
| |
| /* Perform addition with round to odd. */ |
| u.d = temp + (Double) *argZ; |
| /* Ensure the addition is not scheduled after fetestexcept call. */ |
| FORCE_EVAL (u.d); |
| |
| /* Reset rounding mode and test for inexact simultaneously. */ |
| int j = RESET_TESTINEXACT (env); |
| |
| if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff) |
| u.ieee.mantissa1 |= j; |
| |
| /* And finally truncation with round to nearest. */ |
| *res = (Float) u.d; |
| #endif |
| } |
| |
| |
| void VEX_REGPARM(3) |
| h_generic_calc_MAddF64 ( /*OUT*/Double* res, |
| Double* argX, Double* argY, Double* argZ ) |
| { |
| #ifndef ENV_TYPE |
| /* Lame fallback implementation. */ |
| *res = *argX * *argY + *argZ; |
| #else |
| Double x = *argX, y = *argY, z = *argZ; |
| union vg_ieee754_double u, v, w; |
| int adjust = 0; |
| u.d = x; |
| v.d = y; |
| w.d = z; |
| if (UNLIKELY (u.ieee.exponent + v.ieee.exponent |
| >= 0x7ff + IEEE754_DOUBLE_BIAS - DBL_MANT_DIG) |
| || UNLIKELY (u.ieee.exponent >= 0x7ff - DBL_MANT_DIG) |
| || UNLIKELY (v.ieee.exponent >= 0x7ff - DBL_MANT_DIG) |
| || UNLIKELY (w.ieee.exponent >= 0x7ff - DBL_MANT_DIG) |
| || UNLIKELY (u.ieee.exponent + v.ieee.exponent |
| <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG)) { |
| /* If z is Inf, but x and y are finite, the result should be |
| z rather than NaN. */ |
| if (w.ieee.exponent == 0x7ff |
| && u.ieee.exponent != 0x7ff |
| && v.ieee.exponent != 0x7ff) { |
| *res = (z + x) + y; |
| return; |
| } |
| /* If x or y or z is Inf/NaN, or if fma will certainly overflow, |
| or if x * y is less than half of DBL_DENORM_MIN, |
| compute as x * y + z. */ |
| if (u.ieee.exponent == 0x7ff |
| || v.ieee.exponent == 0x7ff |
| || w.ieee.exponent == 0x7ff |
| || u.ieee.exponent + v.ieee.exponent > 0x7ff + IEEE754_DOUBLE_BIAS |
| || u.ieee.exponent + v.ieee.exponent |
| < IEEE754_DOUBLE_BIAS - DBL_MANT_DIG - 2) { |
| *res = x * y + z; |
| return; |
| } |
| if (u.ieee.exponent + v.ieee.exponent |
| >= 0x7ff + IEEE754_DOUBLE_BIAS - DBL_MANT_DIG) { |
| /* Compute 1p-53 times smaller result and multiply |
| at the end. */ |
| if (u.ieee.exponent > v.ieee.exponent) |
| u.ieee.exponent -= DBL_MANT_DIG; |
| else |
| v.ieee.exponent -= DBL_MANT_DIG; |
| /* If x + y exponent is very large and z exponent is very small, |
| it doesn't matter if we don't adjust it. */ |
| if (w.ieee.exponent > DBL_MANT_DIG) |
| w.ieee.exponent -= DBL_MANT_DIG; |
| adjust = 1; |
| } else if (w.ieee.exponent >= 0x7ff - DBL_MANT_DIG) { |
| /* Similarly. |
| If z exponent is very large and x and y exponents are |
| very small, it doesn't matter if we don't adjust it. */ |
| if (u.ieee.exponent > v.ieee.exponent) { |
| if (u.ieee.exponent > DBL_MANT_DIG) |
| u.ieee.exponent -= DBL_MANT_DIG; |
| } else if (v.ieee.exponent > DBL_MANT_DIG) |
| v.ieee.exponent -= DBL_MANT_DIG; |
| w.ieee.exponent -= DBL_MANT_DIG; |
| adjust = 1; |
| } else if (u.ieee.exponent >= 0x7ff - DBL_MANT_DIG) { |
| u.ieee.exponent -= DBL_MANT_DIG; |
| if (v.ieee.exponent) |
| v.ieee.exponent += DBL_MANT_DIG; |
| else |
| v.d *= 0x1p53; |
| } else if (v.ieee.exponent >= 0x7ff - DBL_MANT_DIG) { |
| v.ieee.exponent -= DBL_MANT_DIG; |
| if (u.ieee.exponent) |
| u.ieee.exponent += DBL_MANT_DIG; |
| else |
| u.d *= 0x1p53; |
| } else /* if (u.ieee.exponent + v.ieee.exponent |
| <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG) */ { |
| if (u.ieee.exponent > v.ieee.exponent) |
| u.ieee.exponent += 2 * DBL_MANT_DIG; |
| else |
| v.ieee.exponent += 2 * DBL_MANT_DIG; |
| if (w.ieee.exponent <= 4 * DBL_MANT_DIG + 4) { |
| if (w.ieee.exponent) |
| w.ieee.exponent += 2 * DBL_MANT_DIG; |
| else |
| w.d *= 0x1p106; |
| adjust = -1; |
| } |
| /* Otherwise x * y should just affect inexact |
| and nothing else. */ |
| } |
| x = u.d; |
| y = v.d; |
| z = w.d; |
| } |
| /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ |
| # define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1) |
| Double x1 = x * C; |
| Double y1 = y * C; |
| Double m1 = x * y; |
| x1 = (x - x1) + x1; |
| y1 = (y - y1) + y1; |
| Double x2 = x - x1; |
| Double y2 = y - y1; |
| Double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2; |
| # undef C |
| |
| /* Addition a1 + a2 = z + m1 using Knuth's algorithm. */ |
| Double a1 = z + m1; |
| Double t1 = a1 - z; |
| Double t2 = a1 - t1; |
| t1 = m1 - t1; |
| t2 = z - t2; |
| Double a2 = t1 + t2; |
| |
| ENV_TYPE env; |
| ROUNDTOZERO (env); |
| |
| /* Perform m2 + a2 addition with round to odd. */ |
| u.d = a2 + m2; |
| |
| if (UNLIKELY (adjust < 0)) { |
| if ((u.ieee.mantissa1 & 1) == 0) |
| u.ieee.mantissa1 |= TESTINEXACT (); |
| v.d = a1 + u.d; |
| /* Ensure the addition is not scheduled after fetestexcept call. */ |
| FORCE_EVAL (v.d); |
| } |
| |
| /* Reset rounding mode and test for inexact simultaneously. */ |
| int j = RESET_TESTINEXACT (env) != 0; |
| |
| if (LIKELY (adjust == 0)) { |
| if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff) |
| u.ieee.mantissa1 |= j; |
| /* Result is a1 + u.d. */ |
| *res = a1 + u.d; |
| } else if (LIKELY (adjust > 0)) { |
| if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff) |
| u.ieee.mantissa1 |= j; |
| /* Result is a1 + u.d, scaled up. */ |
| *res = (a1 + u.d) * 0x1p53; |
| } else { |
| /* If a1 + u.d is exact, the only rounding happens during |
| scaling down. */ |
| if (j == 0) { |
| *res = v.d * 0x1p-106; |
| return; |
| } |
| /* If result rounded to zero is not subnormal, no double |
| rounding will occur. */ |
| if (v.ieee.exponent > 106) { |
| *res = (a1 + u.d) * 0x1p-106; |
| return; |
| } |
| /* If v.d * 0x1p-106 with round to zero is a subnormal above |
| or equal to DBL_MIN / 2, then v.d * 0x1p-106 shifts mantissa |
| down just by 1 bit, which means v.ieee.mantissa1 |= j would |
| change the round bit, not sticky or guard bit. |
| v.d * 0x1p-106 never normalizes by shifting up, |
| so round bit plus sticky bit should be already enough |
| for proper rounding. */ |
| if (v.ieee.exponent == 106) { |
| /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding, |
| v.ieee.mantissa1 & 1 is the round bit and j is our sticky |
| bit. In round-to-nearest 001 rounds down like 00, |
| 011 rounds up, even though 01 rounds down (thus we need |
| to adjust), 101 rounds down like 10 and 111 rounds up |
| like 11. */ |
| if ((v.ieee.mantissa1 & 3) == 1) { |
| v.d *= 0x1p-106; |
| if (v.ieee.negative) |
| *res = v.d - 0x1p-1074; |
| else |
| *res = v.d + 0x1p-1074; |
| } else |
| *res = v.d * 0x1p-106; |
| return; |
| } |
| v.ieee.mantissa1 |= j; |
| *res = v.d * 0x1p-106; |
| return; |
| } |
| #endif |
| } |
| |
| /*---------------------------------------------------------------*/ |
| /*--- end host_generic_maddf.c --*/ |
| /*---------------------------------------------------------------*/ |
