From 5b5b04d6282df0364424c6f2c0462e5c1a4394b0 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Sat, 3 Nov 2012 19:48:53 +0000 Subject: Make fma use of Dekker and Knuth algorithms use round-to-nearest (bug 14796). --- sysdeps/generic/math_private.h | 16 ++++++++++++++++ sysdeps/i386/fpu/fclrexcpt.c | 3 ++- sysdeps/i386/fpu/fenv_private.h | 23 +++++++++++++++++++++++ sysdeps/ieee754/dbl-64/s_fma.c | 18 ++++++++++++++++-- sysdeps/ieee754/ldbl-128/s_fmal.c | 18 ++++++++++++++++-- sysdeps/ieee754/ldbl-96/s_fma.c | 18 ++++++++++++++++-- sysdeps/ieee754/ldbl-96/s_fmal.c | 18 ++++++++++++++++-- sysdeps/powerpc/fpu/fclrexcpt.c | 3 ++- sysdeps/s390/fpu/fclrexcpt.c | 3 ++- sysdeps/sh/sh4/fpu/fclrexcpt.c | 1 + sysdeps/sparc/fpu/fclrexcpt.c | 3 ++- sysdeps/sparc/fpu/fenv_private.h | 12 ++++++++++++ sysdeps/x86_64/fpu/fclrexcpt.c | 3 ++- 13 files changed, 126 insertions(+), 13 deletions(-) (limited to 'sysdeps') diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h index b375bc0c56..7661788e6d 100644 --- a/sysdeps/generic/math_private.h +++ b/sysdeps/generic/math_private.h @@ -401,6 +401,22 @@ default_libc_feholdexcept (fenv_t *e) # define libc_feholdexceptl default_libc_feholdexcept #endif +static __always_inline void +default_libc_fesetround (int r) +{ + (void) fesetround (r); +} + +#ifndef libc_fesetround +# define libc_fesetround default_libc_fesetround +#endif +#ifndef libc_fesetroundf +# define libc_fesetroundf default_libc_fesetround +#endif +#ifndef libc_fesetroundl +# define libc_fesetroundl default_libc_fesetround +#endif + static __always_inline void default_libc_feholdexcept_setround (fenv_t *e, int r) { diff --git a/sysdeps/i386/fpu/fclrexcpt.c b/sysdeps/i386/fpu/fclrexcpt.c index f24d07f3d1..f28ef6b3fa 100644 --- a/sysdeps/i386/fpu/fclrexcpt.c +++ b/sysdeps/i386/fpu/fclrexcpt.c @@ -1,5 +1,5 @@ /* Clear given exceptions in current floating-point environment. - Copyright (C) 1997,99,2000, 2001, 2003, 2004 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper , 1997. @@ -65,4 +65,5 @@ strong_alias (__feclearexcept, __old_feclearexcept) compat_symbol (libm, __old_feclearexcept, feclearexcept, GLIBC_2_1); #endif +libm_hidden_ver (__feclearexcept, feclearexcept) versioned_symbol (libm, __feclearexcept, feclearexcept, GLIBC_2_2); diff --git a/sysdeps/i386/fpu/fenv_private.h b/sysdeps/i386/fpu/fenv_private.h index f33f57c39c..03f4c97a9c 100644 --- a/sysdeps/i386/fpu/fenv_private.h +++ b/sysdeps/i386/fpu/fenv_private.h @@ -76,6 +76,24 @@ libc_feholdexcept_387 (fenv_t *e) "st(4)", "st(5)", "st(6)", "st(7)"); } +static __always_inline void +libc_fesetround_sse (int r) +{ + unsigned int mxcsr; + asm (STMXCSR " %0" : "=m" (*&mxcsr)); + mxcsr = (mxcsr & ~0x6000) | (r << 3); + asm volatile (LDMXCSR " %0" : : "m" (*&mxcsr)); +} + +static __always_inline void +libc_fesetround_387 (int r) +{ + fpu_control_t cw; + _FPU_GETCW (cw); + cw = (cw & ~0xc00) | r; + _FPU_SETCW (cw); +} + static __always_inline void libc_feholdexcept_setround_sse (fenv_t *e, int r) { @@ -247,6 +265,7 @@ libc_feresetround_387 (fenv_t *e) #ifdef __SSE_MATH__ # define libc_feholdexceptf libc_feholdexcept_sse +# define libc_fesetroundf libc_fesetround_sse # define libc_feholdexcept_setroundf libc_feholdexcept_setround_sse # define libc_fetestexceptf libc_fetestexcept_sse # define libc_fesetenvf libc_fesetenv_sse @@ -256,6 +275,7 @@ libc_feresetround_387 (fenv_t *e) # define libc_feresetroundf libc_feresetround_sse #else # define libc_feholdexceptf libc_feholdexcept_387 +# define libc_fesetroundf libc_fesetround_387 # define libc_feholdexcept_setroundf libc_feholdexcept_setround_387 # define libc_fetestexceptf libc_fetestexcept_387 # define libc_fesetenvf libc_fesetenv_387 @@ -267,6 +287,7 @@ libc_feresetround_387 (fenv_t *e) #ifdef __SSE2_MATH__ # define libc_feholdexcept libc_feholdexcept_sse +# define libc_fesetround libc_fesetround_sse # define libc_feholdexcept_setround libc_feholdexcept_setround_sse # define libc_fetestexcept libc_fetestexcept_sse # define libc_fesetenv libc_fesetenv_sse @@ -276,6 +297,7 @@ libc_feresetround_387 (fenv_t *e) # define libc_feresetround libc_feresetround_sse #else # define libc_feholdexcept libc_feholdexcept_387 +# define libc_fesetround libc_fesetround_387 # define libc_feholdexcept_setround libc_feholdexcept_setround_387 # define libc_fetestexcept libc_fetestexcept_387 # define libc_fesetenv libc_fesetenv_387 @@ -286,6 +308,7 @@ libc_feresetround_387 (fenv_t *e) #endif /* __SSE2_MATH__ */ #define libc_feholdexceptl libc_feholdexcept_387 +#define libc_fesetroundl libc_fesetround_387 #define libc_feholdexcept_setroundl libc_feholdexcept_setround_387 #define libc_fetestexceptl libc_fetestexcept_387 #define libc_fesetenvl libc_fesetenv_387 diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c index 68c63a1987..07fe715617 100644 --- a/sysdeps/ieee754/dbl-64/s_fma.c +++ b/sysdeps/ieee754/dbl-64/s_fma.c @@ -167,6 +167,9 @@ __fma (double x, double y, double z) if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) return x * y + z; + fenv_t env; + libc_feholdexcept_setround (&env, FE_TONEAREST); + /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ #define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1) double x1 = x * C; @@ -185,9 +188,20 @@ __fma (double x, double y, double z) t1 = m1 - t1; t2 = z - t2; double a2 = t1 + t2; + feclearexcept (FE_INEXACT); - fenv_t env; - libc_feholdexcept_setround (&env, FE_TOWARDZERO); + /* If the result is an exact zero, ensure it has the correct + sign. */ + if (a1 == 0 && m2 == 0) + { + libc_feupdateenv (&env); + /* Ensure that round-to-nearest value of z + m1 is not + reused. */ + asm volatile ("" : "=m" (z) : "m" (z)); + return z + m1; + } + + libc_fesetround (FE_TOWARDZERO); /* Perform m2 + a2 addition with round to odd. */ u.d = a2 + m2; diff --git a/sysdeps/ieee754/ldbl-128/s_fmal.c b/sysdeps/ieee754/ldbl-128/s_fmal.c index c6a3d71f1c..d576403b25 100644 --- a/sysdeps/ieee754/ldbl-128/s_fmal.c +++ b/sysdeps/ieee754/ldbl-128/s_fmal.c @@ -170,6 +170,10 @@ __fmal (long double x, long double y, long double z) if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) return x * y + z; + fenv_t env; + feholdexcept (&env); + fesetround (FE_TONEAREST); + /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ #define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) long double x1 = x * C; @@ -188,9 +192,19 @@ __fmal (long double x, long double y, long double z) t1 = m1 - t1; t2 = z - t2; long double a2 = t1 + t2; + feclearexcept (FE_INEXACT); + + /* If the result is an exact zero, ensure it has the correct + sign. */ + if (a1 == 0 && m2 == 0) + { + feupdateenv (&env); + /* Ensure that round-to-nearest value of z + m1 is not + reused. */ + asm volatile ("" : "=m" (z) : "m" (z)); + return z + m1; + } - fenv_t env; - feholdexcept (&env); fesetround (FE_TOWARDZERO); /* Perform m2 + a2 addition with round to odd. */ u.d = a2 + m2; diff --git a/sysdeps/ieee754/ldbl-96/s_fma.c b/sysdeps/ieee754/ldbl-96/s_fma.c index 001d8063d4..bf2d794c9b 100644 --- a/sysdeps/ieee754/ldbl-96/s_fma.c +++ b/sysdeps/ieee754/ldbl-96/s_fma.c @@ -42,6 +42,10 @@ __fma (double x, double y, double z) if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) return x * y + z; + fenv_t env; + feholdexcept (&env); + fesetround (FE_TONEAREST); + /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ #define C ((1ULL << (LDBL_MANT_DIG + 1) / 2) + 1) long double x1 = (long double) x * C; @@ -60,9 +64,19 @@ __fma (double x, double y, double z) t1 = m1 - t1; t2 = z - t2; long double a2 = t1 + t2; + feclearexcept (FE_INEXACT); + + /* If the result is an exact zero, ensure it has the correct + sign. */ + if (a1 == 0 && m2 == 0) + { + feupdateenv (&env); + /* Ensure that round-to-nearest value of z + m1 is not + reused. */ + asm volatile ("" : "=m" (z) : "m" (z)); + return z + m1; + } - fenv_t env; - feholdexcept (&env); fesetround (FE_TOWARDZERO); /* Perform m2 + a2 addition with round to odd. */ a2 = a2 + m2; diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c index b5fdcba6bc..32e71a18ba 100644 --- a/sysdeps/ieee754/ldbl-96/s_fmal.c +++ b/sysdeps/ieee754/ldbl-96/s_fmal.c @@ -168,6 +168,10 @@ __fmal (long double x, long double y, long double z) if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0)) return x * y + z; + fenv_t env; + feholdexcept (&env); + fesetround (FE_TONEAREST); + /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */ #define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) long double x1 = x * C; @@ -186,9 +190,19 @@ __fmal (long double x, long double y, long double z) t1 = m1 - t1; t2 = z - t2; long double a2 = t1 + t2; + feclearexcept (FE_INEXACT); + + /* If the result is an exact zero, ensure it has the correct + sign. */ + if (a1 == 0 && m2 == 0) + { + feupdateenv (&env); + /* Ensure that round-to-nearest value of z + m1 is not + reused. */ + asm volatile ("" : "=m" (z) : "m" (z)); + return z + m1; + } - fenv_t env; - feholdexcept (&env); fesetround (FE_TOWARDZERO); /* Perform m2 + a2 addition with round to odd. */ u.d = a2 + m2; diff --git a/sysdeps/powerpc/fpu/fclrexcpt.c b/sysdeps/powerpc/fpu/fclrexcpt.c index 87376bfe97..e99cc58ace 100644 --- a/sysdeps/powerpc/fpu/fclrexcpt.c +++ b/sysdeps/powerpc/fpu/fclrexcpt.c @@ -1,5 +1,5 @@ /* Clear given exceptions in current floating-point environment. - Copyright (C) 1997,99,2000,01 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -44,4 +44,5 @@ strong_alias (__feclearexcept, __old_feclearexcept) compat_symbol (libm, __old_feclearexcept, feclearexcept, GLIBC_2_1); #endif +libm_hidden_ver (__feclearexcept, feclearexcept) versioned_symbol (libm, __feclearexcept, feclearexcept, GLIBC_2_2); diff --git a/sysdeps/s390/fpu/fclrexcpt.c b/sysdeps/s390/fpu/fclrexcpt.c index 2352d74b07..3e8d9bb3aa 100644 --- a/sysdeps/s390/fpu/fclrexcpt.c +++ b/sysdeps/s390/fpu/fclrexcpt.c @@ -1,5 +1,5 @@ /* Clear given exceptions in current floating-point environment. - Copyright (C) 2000 Free Software Foundation, Inc. + Copyright (C) 2000-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -37,3 +37,4 @@ feclearexcept (int excepts) /* Success. */ return 0; } +libm_hidden_def (feclearexcept) diff --git a/sysdeps/sh/sh4/fpu/fclrexcpt.c b/sysdeps/sh/sh4/fpu/fclrexcpt.c index b4b2ead02c..2d31fa6516 100644 --- a/sysdeps/sh/sh4/fpu/fclrexcpt.c +++ b/sysdeps/sh/sh4/fpu/fclrexcpt.c @@ -39,3 +39,4 @@ feclearexcept (int excepts) return 0; } +libm_hidden_def (feclearexcept) diff --git a/sysdeps/sparc/fpu/fclrexcpt.c b/sysdeps/sparc/fpu/fclrexcpt.c index 82b6ac4e2e..50ec4f4ff3 100644 --- a/sysdeps/sparc/fpu/fclrexcpt.c +++ b/sysdeps/sparc/fpu/fclrexcpt.c @@ -1,5 +1,5 @@ /* Clear given exceptions in current floating-point environment. - Copyright (C) 1997, 1999, 2000 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -39,4 +39,5 @@ strong_alias (__feclearexcept, __old_feclearexcept) compat_symbol (libm, __old_feclearexcept, feclearexcept, GLIBC_2_1); #endif +libm_hidden_ver (__feclearexcept, feclearexcept) versioned_symbol (libm, __feclearexcept, feclearexcept, GLIBC_2_2); diff --git a/sysdeps/sparc/fpu/fenv_private.h b/sysdeps/sparc/fpu/fenv_private.h index a6e8e95a55..8690879ddc 100644 --- a/sysdeps/sparc/fpu/fenv_private.h +++ b/sysdeps/sparc/fpu/fenv_private.h @@ -13,6 +13,15 @@ libc_feholdexcept (fenv_t *e) __fenv_ldfsr(etmp); } +static __always_inline void +libc_fesetround (int r) +{ + fenv_t etmp; + __fenv_stfsr(etmp); + etmp = (etmp & ~__FE_ROUND_MASK) | (r); + __fenv_ldfsr(etmp); +} + static __always_inline void libc_feholdexcept_setround (fenv_t *e, int r) { @@ -79,6 +88,7 @@ libc_feresetround (fenv_t *e) } #define libc_feholdexceptf libc_feholdexcept +#define libc_fesetroundf libc_fesetround #define libc_feholdexcept_setroundf libc_feholdexcept_setround #define libc_fetestexceptf libc_fetestexcept #define libc_fesetenvf libc_fesetenv @@ -87,6 +97,7 @@ libc_feresetround (fenv_t *e) #define libc_feholdsetroundf libc_feholdsetround #define libc_feresetroundf libc_feresetround #define libc_feholdexcept libc_feholdexcept +#define libc_fesetround libc_fesetround #define libc_feholdexcept_setround libc_feholdexcept_setround #define libc_fetestexcept libc_fetestexcept #define libc_fesetenv libc_fesetenv @@ -95,6 +106,7 @@ libc_feresetround (fenv_t *e) #define libc_feholdsetround libc_feholdsetround #define libc_feresetround libc_feresetround #define libc_feholdexceptl libc_feholdexcept +#define libc_fesetroundl libc_fesetround #define libc_feholdexcept_setroundl libc_feholdexcept_setround #define libc_fetestexceptl libc_fetestexcept #define libc_fesetenvl libc_fesetenv diff --git a/sysdeps/x86_64/fpu/fclrexcpt.c b/sysdeps/x86_64/fpu/fclrexcpt.c index 5ef316212a..2d40313394 100644 --- a/sysdeps/x86_64/fpu/fclrexcpt.c +++ b/sysdeps/x86_64/fpu/fclrexcpt.c @@ -1,5 +1,5 @@ /* Clear given exceptions in current floating-point environment. - Copyright (C) 2001 Free Software Foundation, Inc. + Copyright (C) 2001-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -49,3 +49,4 @@ feclearexcept (int excepts) /* Success. */ return 0; } +libm_hidden_def (feclearexcept) -- cgit 1.4.1