From bec749fda1cbc1934f7e58dd2763603f4f207f26 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Mon, 1 Oct 2012 08:30:06 +0000 Subject: Fix sign of inexact zero return from fma (bug 14645). --- sysdeps/ieee754/dbl-64/s_fma.c | 5 +++++ sysdeps/ieee754/ldbl-128/s_fmal.c | 5 +++++ sysdeps/ieee754/ldbl-96/s_fmal.c | 5 +++++ 3 files changed, 15 insertions(+) (limited to 'sysdeps') diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c index c9809fb180..5e21461a4b 100644 --- a/sysdeps/ieee754/dbl-64/s_fma.c +++ b/sysdeps/ieee754/dbl-64/s_fma.c @@ -49,6 +49,11 @@ __fma (double x, double y, double z) && u.ieee.exponent != 0x7ff && v.ieee.exponent != 0x7ff) return (z + x) + y; + /* If z is zero and x are y are nonzero, compute the result + as x * y to avoid the wrong sign of a zero result if x * y + underflows to 0. */ + if (z == 0 && x != 0 && y != 0) + return 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. */ diff --git a/sysdeps/ieee754/ldbl-128/s_fmal.c b/sysdeps/ieee754/ldbl-128/s_fmal.c index df68ade435..46b3d81ce5 100644 --- a/sysdeps/ieee754/ldbl-128/s_fmal.c +++ b/sysdeps/ieee754/ldbl-128/s_fmal.c @@ -50,6 +50,11 @@ __fmal (long double x, long double y, long double z) && u.ieee.exponent != 0x7fff && v.ieee.exponent != 0x7fff) return (z + x) + y; + /* If z is zero and x are y are nonzero, compute the result + as x * y to avoid the wrong sign of a zero result if x * y + underflows to 0. */ + if (z == 0 && x != 0 && y != 0) + return 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. */ diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c index c27b0bd852..d125124286 100644 --- a/sysdeps/ieee754/ldbl-96/s_fmal.c +++ b/sysdeps/ieee754/ldbl-96/s_fmal.c @@ -50,6 +50,11 @@ __fmal (long double x, long double y, long double z) && u.ieee.exponent != 0x7fff && v.ieee.exponent != 0x7fff) return (z + x) + y; + /* If z is zero and x are y are nonzero, compute the result + as x * y to avoid the wrong sign of a zero result if x * y + underflows to 0. */ + if (z == 0 && x != 0 && y != 0) + return 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. */ -- cgit 1.4.1