Hi! This patch on top of http://sources.redhat.com/ml/libc-alpha/2010-10/msg00022.html implements extended long double and IEEE quad long double fmal and fixes some bugs in fma, in particular fma (DBL_MAX, DBL_MAX, -__builtin_inf ()) should be __builtin_inf () rather than NaN, as the computation is supposed to happen with infinite precision before rounding and thus x * y will be still finite, albeit very large. Also, unfortunately it seems two fetestexcept calls for adjust == -1 are unavoidable as shown by one of the added testcases - if result is slightly above DBL_MIN, we compute (a1 + u.d) * 0x1p106 in round to nearest, but if the inexact bit from first addition hasn't been ored into u.d, it is lost. Attached is also greatly improved mpfr testing proglet, now handles all of float (-DDO_FLT), double (-DDO_DBL) and long double (-DDO_LDBL, both extended and quad). Tested on x86_64-linux and i686-linux and slightly on s390x-linux, unfortunately I've discovered mpfr_set_ld is buggy there, so the tester is currently unusable there. 2010-10-15 Jakub Jelinek [BZ #3268] * math/libm-test.inc (fma_test): Some new testcases. * sysdeps/ieee754/ldbl-128/s_fmal.c: New file. * sysdeps/ieee754/ldbl-96/s_fma.c (__fma): Fix fma with finite x and y and infinite z. Do multiplication by C already in long double. * sysdeps/ieee754/ldbl-96/s_fmal.c: New file. * sysdeps/ieee754/dbl-64/s_fma.c (__fma): Fix fma with finite x and y and infinite z. Do bitwise or of inexact bit into u.d. * sysdeps/ieee754/ldbl-64-128/s_fmal.c: New file. * sysdeps/i386/fpu/s_fmaf.S: Removed. * sysdeps/i386/fpu/s_fma.S: Removed. * sysdeps/i386/fpu/s_fmal.S: Removed. --- libc/math/libm-test.inc.jj 2010-10-14 20:21:24.000000000 +0200 +++ libc/math/libm-test.inc 2010-10-15 17:51:16.000000000 +0200 @@ -2787,8 +2787,24 @@ fma_test (void) TEST_fff_f (fma, minus_infty, plus_infty, plus_infty, nan_value, INVALID_EXCEPTION); TEST_fff_f (fma, plus_infty, minus_infty, plus_infty, nan_value, INVALID_EXCEPTION); TEST_fff_f (fma, minus_infty, minus_infty, minus_infty, nan_value, INVALID_EXCEPTION); + TEST_fff_f (fma, plus_infty, 3.5L, minus_infty, nan_value, INVALID_EXCEPTION); + TEST_fff_f (fma, minus_infty, -7.5L, minus_infty, nan_value, INVALID_EXCEPTION); + TEST_fff_f (fma, -13.5L, plus_infty, plus_infty, nan_value, INVALID_EXCEPTION); + TEST_fff_f (fma, minus_infty, 7.5L, plus_infty, nan_value, INVALID_EXCEPTION); TEST_fff_f (fma, 1.25L, 0.75L, 0.0625L, 1.0L); + + FLOAT fltmax = CHOOSE (LDBL_MAX, DBL_MAX, FLT_MAX, + LDBL_MAX, DBL_MAX, FLT_MAX); + TEST_fff_f (fma, -fltmax, -fltmax, minus_infty, minus_infty); + TEST_fff_f (fma, fltmax / 2, fltmax / 2, minus_infty, minus_infty); + TEST_fff_f (fma, -fltmax, fltmax, plus_infty, plus_infty); + TEST_fff_f (fma, fltmax / 2, -fltmax / 4, plus_infty, plus_infty); + TEST_fff_f (fma, plus_infty, 4, plus_infty, plus_infty); + TEST_fff_f (fma, 2, minus_infty, minus_infty, minus_infty); + TEST_fff_f (fma, minus_infty, minus_infty, plus_infty, plus_infty); + TEST_fff_f (fma, plus_infty, minus_infty, minus_infty, minus_infty); + #if defined (TEST_FLOAT) && FLT_MANT_DIG == 24 TEST_fff_f (fma, 0x1.7ff8p+13, 0x1.000002p+0, 0x1.ffffp-24, 0x1.7ff802p+13); TEST_fff_f (fma, 0x1.fffp+0, 0x1.00001p+0, -0x1.fffp+0, 0x1.fffp-20); @@ -2818,6 +2834,15 @@ fma_test (void) TEST_fff_f (fma, -0x1.19cab66d73e17p-959, 0x1.c7108a8c5ff51p-107, -0x0.80b0ad65d9b64p-1022, -0x0.80b0ad65d9d59p-1022); TEST_fff_f (fma, -0x1.d2eaed6e8e9d3p-979, -0x1.4e066c62ac9ddp-63, -0x0.9245e6b003454p-1022, -0x0.9245c09c5fb5dp-1022); TEST_fff_f (fma, 0x1.153d650bb9f06p-907, 0x1.2d01230d48407p-125, -0x0.b278d5acfc3cp-1022, -0x0.b22757123bbe9p-1022); + TEST_fff_f (fma, -0x1.fffffffffffffp-711, 0x1.fffffffffffffp-275, 0x1.fffffe00007ffp-983, 0x1.7ffffe00007ffp-983); +#endif +#if defined (TEST_LDOUBLE) && LDBL_MANT_DIG == 64 + TEST_fff_f (fma, -0x8.03fcp+3696L, 0xf.fffffffffffffffp-6140L, 0x8.3ffffffffffffffp-2450L, -0x8.01ecp-2440L); + TEST_fff_f (fma, 0x9.fcp+2033L, -0x8.000e1f000ff800fp-3613L, -0xf.fffffffffffc0ffp-1579L, -0xd.fc119fb093ed092p-1577L); + TEST_fff_f (fma, 0xc.7fc000003ffffffp-1194L, 0x8.1e0003fffffffffp+15327L, -0x8.fffep+14072L, 0xc.ae9f164020effffp+14136L); + TEST_fff_f (fma, -0x8.0001fc000000003p+1798L, 0xcp-2230L, 0x8.f7e000000000007p-468L, -0xc.0002f9ffee10404p-429L); + TEST_fff_f (fma, 0xc.0000000000007ffp+10130L, -0x8.000000000000001p+4430L, 0xc.07000000001ffffp+14513L, -0xb.fffffffffffd7e4p+14563L); + TEST_fff_f (fma, 0xb.ffffp-4777L, 0x8.000000fffffffffp-11612L, -0x0.3800fff8p-16385L, 0x5.c7fe80c7ffeffffp-16385L); #endif END (fma); --- libc/sysdeps/ieee754/ldbl-128/s_fmal.c.jj 2010-10-15 18:03:31.000000000 +0200 +++ libc/sysdeps/ieee754/ldbl-128/s_fmal.c 2010-10-15 18:08:43.000000000 +0200 @@ -0,0 +1,221 @@ +/* Compute x * y + z as ternary operation. + Copyright (C) 2010 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Jakub Jelinek , 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, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include +#include +#include +#include + +/* This implementation uses rounding to odd to avoid problems with + double rounding. See a paper by Boldo and Melquiond: + http://www.lri.fr/~melquion/doc/08-tc.pdf */ + +long double +__fmal (long double x, long double y, long double z) +{ + union ieee854_long_double u, v, w; + int adjust = 0; + u.d = x; + v.d = y; + w.d = z; + if (__builtin_expect (u.ieee.exponent + v.ieee.exponent + >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS + - LDBL_MANT_DIG, 0) + || __builtin_expect (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0) + || __builtin_expect (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0) + || __builtin_expect (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0) + || __builtin_expect (u.ieee.exponent + v.ieee.exponent + <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG, 0)) + { + /* If z is Inf, but x and y are finite, the result should be + z rather than NaN. */ + if (w.ieee.exponent == 0x7fff + && u.ieee.exponent != 0x7fff + && v.ieee.exponent != 0x7fff) + return (z + x) + y; + /* If x or y or z is Inf/NaN, or if fma will certainly overflow, + or if x * y is less than half of LDBL_DENORM_MIN, + compute as x * y + z. */ + if (u.ieee.exponent == 0x7fff + || v.ieee.exponent == 0x7fff + || w.ieee.exponent == 0x7fff + || u.ieee.exponent + v.ieee.exponent + > 0x7fff + IEEE854_LONG_DOUBLE_BIAS + || u.ieee.exponent + v.ieee.exponent + < IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG - 2) + return x * y + z; + if (u.ieee.exponent + v.ieee.exponent + >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG) + { + /* Compute 1p-113 times smaller result and multiply + at the end. */ + if (u.ieee.exponent > v.ieee.exponent) + u.ieee.exponent -= LDBL_MANT_DIG; + else + v.ieee.exponent -= LDBL_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 > LDBL_MANT_DIG) + w.ieee.exponent -= LDBL_MANT_DIG; + adjust = 1; + } + else if (w.ieee.exponent >= 0x7fff - LDBL_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 > LDBL_MANT_DIG) + u.ieee.exponent -= LDBL_MANT_DIG; + } + else if (v.ieee.exponent > LDBL_MANT_DIG) + v.ieee.exponent -= LDBL_MANT_DIG; + w.ieee.exponent -= LDBL_MANT_DIG; + adjust = 1; + } + else if (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG) + { + u.ieee.exponent -= LDBL_MANT_DIG; + if (v.ieee.exponent) + v.ieee.exponent += LDBL_MANT_DIG; + else + v.d *= 0x1p113L; + } + else if (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG) + { + v.ieee.exponent -= LDBL_MANT_DIG; + if (u.ieee.exponent) + u.ieee.exponent += LDBL_MANT_DIG; + else + u.d *= 0x1p113L; + } + else /* if (u.ieee.exponent + v.ieee.exponent + <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */ + { + if (u.ieee.exponent > v.ieee.exponent) + u.ieee.exponent += 2 * LDBL_MANT_DIG; + else + v.ieee.exponent += 2 * LDBL_MANT_DIG; + if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 4) + { + if (w.ieee.exponent) + w.ieee.exponent += 2 * LDBL_MANT_DIG; + else + w.d *= 0x1p226L; + 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 ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) + long double x1 = x * C; + long double y1 = y * C; + long double m1 = x * y; + x1 = (x - x1) + x1; + y1 = (y - y1) + y1; + long double x2 = x - x1; + long double y2 = y - y1; + long double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2; + + /* Addition a1 + a2 = z + m1 using Knuth's algorithm. */ + long double a1 = z + m1; + long double t1 = a1 - z; + long double t2 = a1 - t1; + t1 = m1 - t1; + t2 = z - t2; + long double a2 = t1 + t2; + + fenv_t env; + feholdexcept (&env); + fesetround (FE_TOWARDZERO); + /* Perform m2 + a2 addition with round to odd. */ + u.d = a2 + m2; + + if (__builtin_expect (adjust == 0, 1)) + { + if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff) + u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0; + feupdateenv (&env); + /* Result is a1 + u.d. */ + return a1 + u.d; + } + else if (__builtin_expect (adjust > 0, 1)) + { + if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff) + u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0; + feupdateenv (&env); + /* Result is a1 + u.d, scaled up. */ + return (a1 + u.d) * 0x1p113L; + } + else + { + if ((u.ieee.mantissa3 & 1) == 0) + u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0; + v.d = a1 + u.d; + int j = fetestexcept (FE_INEXACT) != 0; + feupdateenv (&env); + /* Ensure the following computations are performed in default rounding + mode instead of just reusing the round to zero computation. */ + asm volatile ("" : "=m" (u) : "m" (u)); + /* If a1 + u.d is exact, the only rounding happens during + scaling down. */ + if (j == 0) + return v.d * 0x1p-226L; + /* If result rounded to zero is not subnormal, no double + rounding will occur. */ + if (v.ieee.exponent > 226) + return (a1 + u.d) * 0x1p-226L; + /* If v.d * 0x1p-226L with round to zero is a subnormal above + or equal to LDBL_MIN / 2, then v.d * 0x1p-226L shifts mantissa + down just by 1 bit, which means v.ieee.mantissa3 |= j would + change the round bit, not sticky or guard bit. + v.d * 0x1p-226L never normalizes by shifting up, + so round bit plus sticky bit should be already enough + for proper rounding. */ + if (v.ieee.exponent == 226) + { + /* v.ieee.mantissa3 & 2 is LSB bit of the result before rounding, + v.ieee.mantissa3 & 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.mantissa3 & 3) == 1) + { + v.d *= 0x1p-226L; + if (v.ieee.negative) + return v.d - 0x1p-16493L /* __LDBL_DENORM_MIN__ */; + else + return v.d + 0x1p-16493L /* __LDBL_DENORM_MIN__ */; + } + else + return v.d * 0x1p-226L; + } + v.ieee.mantissa3 |= j; + return v.d * 0x1p-226L; + } +} +weak_alias (__fmal, fmal) --- libc/sysdeps/ieee754/ldbl-96/s_fma.c.jj 2010-10-14 09:08:58.000000000 +0200 +++ libc/sysdeps/ieee754/ldbl-96/s_fma.c 2010-10-15 13:41:15.000000000 +0200 @@ -30,11 +30,20 @@ double __fma (double x, double y, double z) { + if (__builtin_expect (isinf (z), 0)) + { + /* If z is Inf, but x and y are finite, the result should be + z rather than NaN. */ + if (finite (x) && finite (y)) + return (z + x) + y; + return (x * y) + z; + } + /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ #define C ((1ULL << (LDBL_MANT_DIG + 1) / 2) + 1) - long double x1 = x * C; - long double y1 = y * C; - long double m1 = x * y; + long double x1 = (long double) x * C; + long double y1 = (long double) y * C; + long double m1 = (long double) x * y; x1 = (x - x1) + x1; y1 = (y - y1) + y1; long double x2 = x - x1; --- libc/sysdeps/ieee754/ldbl-96/s_fmal.c.jj 2010-10-14 22:39:23.000000000 +0200 +++ libc/sysdeps/ieee754/ldbl-96/s_fmal.c 2010-10-15 17:51:46.000000000 +0200 @@ -0,0 +1,221 @@ +/* Compute x * y + z as ternary operation. + Copyright (C) 2010 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Jakub Jelinek , 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, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include +#include +#include +#include + +/* This implementation uses rounding to odd to avoid problems with + double rounding. See a paper by Boldo and Melquiond: + http://www.lri.fr/~melquion/doc/08-tc.pdf */ + +long double +__fmal (long double x, long double y, long double z) +{ + union ieee854_long_double u, v, w; + int adjust = 0; + u.d = x; + v.d = y; + w.d = z; + if (__builtin_expect (u.ieee.exponent + v.ieee.exponent + >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS + - LDBL_MANT_DIG, 0) + || __builtin_expect (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0) + || __builtin_expect (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0) + || __builtin_expect (w.ieee.exponent >= 0x7fff - LDBL_MANT_DIG, 0) + || __builtin_expect (u.ieee.exponent + v.ieee.exponent + <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG, 0)) + { + /* If z is Inf, but x and y are finite, the result should be + z rather than NaN. */ + if (w.ieee.exponent == 0x7fff + && u.ieee.exponent != 0x7fff + && v.ieee.exponent != 0x7fff) + return (z + x) + y; + /* If x or y or z is Inf/NaN, or if fma will certainly overflow, + or if x * y is less than half of LDBL_DENORM_MIN, + compute as x * y + z. */ + if (u.ieee.exponent == 0x7fff + || v.ieee.exponent == 0x7fff + || w.ieee.exponent == 0x7fff + || u.ieee.exponent + v.ieee.exponent + > 0x7fff + IEEE854_LONG_DOUBLE_BIAS + || u.ieee.exponent + v.ieee.exponent + < IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG - 2) + return x * y + z; + if (u.ieee.exponent + v.ieee.exponent + >= 0x7fff + IEEE854_LONG_DOUBLE_BIAS - LDBL_MANT_DIG) + { + /* Compute 1p-64 times smaller result and multiply + at the end. */ + if (u.ieee.exponent > v.ieee.exponent) + u.ieee.exponent -= LDBL_MANT_DIG; + else + v.ieee.exponent -= LDBL_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 > LDBL_MANT_DIG) + w.ieee.exponent -= LDBL_MANT_DIG; + adjust = 1; + } + else if (w.ieee.exponent >= 0x7fff - LDBL_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 > LDBL_MANT_DIG) + u.ieee.exponent -= LDBL_MANT_DIG; + } + else if (v.ieee.exponent > LDBL_MANT_DIG) + v.ieee.exponent -= LDBL_MANT_DIG; + w.ieee.exponent -= LDBL_MANT_DIG; + adjust = 1; + } + else if (u.ieee.exponent >= 0x7fff - LDBL_MANT_DIG) + { + u.ieee.exponent -= LDBL_MANT_DIG; + if (v.ieee.exponent) + v.ieee.exponent += LDBL_MANT_DIG; + else + v.d *= 0x1p64L; + } + else if (v.ieee.exponent >= 0x7fff - LDBL_MANT_DIG) + { + v.ieee.exponent -= LDBL_MANT_DIG; + if (u.ieee.exponent) + u.ieee.exponent += LDBL_MANT_DIG; + else + u.d *= 0x1p64L; + } + else /* if (u.ieee.exponent + v.ieee.exponent + <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) */ + { + if (u.ieee.exponent > v.ieee.exponent) + u.ieee.exponent += 2 * LDBL_MANT_DIG; + else + v.ieee.exponent += 2 * LDBL_MANT_DIG; + if (w.ieee.exponent <= 4 * LDBL_MANT_DIG + 4) + { + if (w.ieee.exponent) + w.ieee.exponent += 2 * LDBL_MANT_DIG; + else + w.d *= 0x1p128L; + 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 ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) + long double x1 = x * C; + long double y1 = y * C; + long double m1 = x * y; + x1 = (x - x1) + x1; + y1 = (y - y1) + y1; + long double x2 = x - x1; + long double y2 = y - y1; + long double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2; + + /* Addition a1 + a2 = z + m1 using Knuth's algorithm. */ + long double a1 = z + m1; + long double t1 = a1 - z; + long double t2 = a1 - t1; + t1 = m1 - t1; + t2 = z - t2; + long double a2 = t1 + t2; + + fenv_t env; + feholdexcept (&env); + fesetround (FE_TOWARDZERO); + /* Perform m2 + a2 addition with round to odd. */ + u.d = a2 + m2; + + if (__builtin_expect (adjust == 0, 1)) + { + if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff) + u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0; + feupdateenv (&env); + /* Result is a1 + u.d. */ + return a1 + u.d; + } + else if (__builtin_expect (adjust > 0, 1)) + { + if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff) + u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0; + feupdateenv (&env); + /* Result is a1 + u.d, scaled up. */ + return (a1 + u.d) * 0x1p64L; + } + else + { + if ((u.ieee.mantissa1 & 1) == 0) + u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0; + v.d = a1 + u.d; + int j = fetestexcept (FE_INEXACT) != 0; + feupdateenv (&env); + /* Ensure the following computations are performed in default rounding + mode instead of just reusing the round to zero computation. */ + asm volatile ("" : "=m" (u) : "m" (u)); + /* If a1 + u.d is exact, the only rounding happens during + scaling down. */ + if (j == 0) + return v.d * 0x1p-128L; + /* If result rounded to zero is not subnormal, no double + rounding will occur. */ + if (v.ieee.exponent > 128) + return (a1 + u.d) * 0x1p-128L; + /* If v.d * 0x1p-128L with round to zero is a subnormal above + or equal to LDBL_MIN / 2, then v.d * 0x1p-128L 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-128L never normalizes by shifting up, + so round bit plus sticky bit should be already enough + for proper rounding. */ + if (v.ieee.exponent == 128) + { + /* 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-128L; + if (v.ieee.negative) + return v.d - 0x1p-16445L /* __LDBL_DENORM_MIN__ */; + else + return v.d + 0x1p-16445L /* __LDBL_DENORM_MIN__ */; + } + else + return v.d * 0x1p-128L; + } + v.ieee.mantissa1 |= j; + return v.d * 0x1p-128L; + } +} +weak_alias (__fmal, fmal) --- libc/sysdeps/ieee754/dbl-64/s_fma.c.jj 2010-10-14 20:26:05.000000000 +0200 +++ libc/sysdeps/ieee754/dbl-64/s_fma.c 2010-10-15 17:56:30.000000000 +0200 @@ -43,6 +43,12 @@ __fma (double x, double y, double z) || __builtin_expect (u.ieee.exponent + v.ieee.exponent <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG, 0)) { + /* 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) + return (z + x) + y; /* 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. */ @@ -165,6 +171,8 @@ __fma (double x, double y, double z) } else { + if ((u.ieee.mantissa1 & 1) == 0) + u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0; v.d = a1 + u.d; int j = fetestexcept (FE_INEXACT) != 0; feupdateenv (&env); --- libc/sysdeps/ieee754/ldbl-64-128/s_fmal.c.jj 2010-10-15 18:05:04.000000000 +0200 +++ libc/sysdeps/ieee754/ldbl-64-128/s_fmal.c 2010-10-15 18:05:17.000000000 +0200 @@ -0,0 +1,5 @@ +#include +#undef weak_alias +#define weak_alias(n,a) +#include +long_double_symbol (libm, __fmal, fmal); --- libc/sysdeps/i386/fpu/s_fmaf.S.jj 2009-05-16 19:23:38.000000000 +0200 +++ libc/sysdeps/i386/fpu/s_fmaf.S 2010-09-14 11:29:01.831482015 +0200 @@ -1,31 +0,0 @@ -/* Compute (X * Y) + Z as ternary operation. - Copyright (C) 1997 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper , 1997. - - 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, write to the Free - Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA - 02111-1307 USA. */ - -#include - - .text -ENTRY(__fmaf) - flds 4(%esp) // x - fmuls 8(%esp) // x * y - flds 12(%esp) // z : x * y - faddp // (x * y) + z - ret -END(__fmaf) -weak_alias (__fmaf, fmaf) --- libc/sysdeps/i386/fpu/s_fma.S.jj 2009-05-16 19:23:38.000000000 +0200 +++ libc/sysdeps/i386/fpu/s_fma.S 2010-09-14 11:29:01.831482015 +0200 @@ -1,31 +0,0 @@ -/* Compute (X * Y) + Z as ternary operation. - Copyright (C) 1997, 1998 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper , 1997. - - 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, write to the Free - Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA - 02111-1307 USA. */ - -#include - - .text -ENTRY(__fma) - fldl 4(%esp) // x - fmull 12(%esp) // x * y - fldl 20(%esp) // z : x * y - faddp // (x * y) + z - ret -END(__fma) -weak_alias (__fma, fma) --- libc/sysdeps/i386/fpu/s_fmal.S.jj 2009-05-16 19:23:38.000000000 +0200 +++ libc/sysdeps/i386/fpu/s_fmal.S 2010-09-14 11:29:01.831482015 +0200 @@ -1,32 +0,0 @@ -/* Compute (X * Y) + Z as ternary operation. - Copyright (C) 1997 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper , 1997. - - 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, write to the Free - Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA - 02111-1307 USA. */ - -#include - - .text -ENTRY(__fmal) - fldt 4(%esp) // x - fldt 16(%esp) // x : y - fmulp // x * y - fldt 28(%esp) // z : x * y - faddp // (x * y) + z - ret -END(__fmal) -weak_alias (__fmal, fmal) Jakub