From patchwork Wed May 7 17:39:44 2014 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Joseph Myers X-Patchwork-Id: 829 Return-Path: X-Original-To: siddhesh@wilcox.dreamhost.com Delivered-To: siddhesh@wilcox.dreamhost.com Received: from homiemail-mx20.g.dreamhost.com (mx2.sub5.homie.mail.dreamhost.com [208.113.200.128]) by wilcox.dreamhost.com (Postfix) with ESMTP id 53865360098 for ; Wed, 7 May 2014 10:39:57 -0700 (PDT) Received: by homiemail-mx20.g.dreamhost.com (Postfix, from userid 14307373) id 0431E41839E34; Wed, 7 May 2014 10:39:56 -0700 (PDT) X-Original-To: glibc@patchwork.siddhesh.in Delivered-To: x14307373@homiemail-mx20.g.dreamhost.com Received: from sourceware.org (server1.sourceware.org [209.132.180.131]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (No client certificate requested) by homiemail-mx20.g.dreamhost.com (Postfix) with ESMTPS id CC5A341839E28 for ; Wed, 7 May 2014 10:39:56 -0700 (PDT) DomainKey-Signature: a=rsa-sha1; c=nofws; d=sourceware.org; h=list-id :list-unsubscribe:list-subscribe:list-archive:list-post :list-help:sender:date:from:to:subject:message-id:mime-version :content-type; q=dns; s=default; b=aylEe1TpXvTWsE5GNQAMJZvV7VSru NALzfa2nMXOECO1Ohxq94MUY3w64GeDU1j+nqZnDBl5hslqUZiB2z0nBcDIkzUzm 3NOuX+Xrp9Dh6C4w89hIDOsxhyftbdMeoaKY3vvYV2DnlcWu1u1WSjGLxaDeY0tE XCGVf0FicnS/9E= DKIM-Signature: v=1; a=rsa-sha1; c=relaxed; d=sourceware.org; h=list-id :list-unsubscribe:list-subscribe:list-archive:list-post :list-help:sender:date:from:to:subject:message-id:mime-version :content-type; s=default; bh=lpTG+BtJIx4j5gnIe3Sq0xWceXk=; b=ipt /sCNd+ZalF6lNuurXEzoTOqbqTZqsbNXcYuDW6/hxmjw2+1lbddrya53EG1G3WDT PrJTfhEiOsMXfSpk6u55hM6AZziYNzmh2aOoSUuXpp3Aq6OTbd8FIyD/euAchpEY A1gr05U7/T+Y8McGe2Gqn2oWpZ/z/XwEtGvOFGsA= Received: (qmail 18764 invoked by alias); 7 May 2014 17:39:54 -0000 Mailing-List: contact libc-alpha-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: libc-alpha-owner@sourceware.org Delivered-To: mailing list libc-alpha@sourceware.org Received: (qmail 18748 invoked by uid 89); 7 May 2014 17:39:53 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-2.0 required=5.0 tests=AWL, BAYES_00 autolearn=ham version=3.3.2 X-HELO: relay1.mentorg.com Date: Wed, 7 May 2014 17:39:44 +0000 From: "Joseph S. Myers" To: Subject: Fix erf underflow handling near 0 (bug 16516) Message-ID: MIME-Version: 1.0 X-DH-Original-To: glibc@patchwork.siddhesh.in Bug 16516 reports spurious underflows from erf (for all floating-point types), when the result is close to underflowing but does not actually underflow. erf (x) is about (2/sqrt(pi))*x for x close to 0, so there are subnormal arguments for which it does not underflow. The various implementations do (x + efx*x) (for efx = 2/sqrt(pi) - 1), for greater accuracy than if just using a single multiplication by an approximation to 2/sqrt(pi) (effectively, this way there are a few more bits in the approximation to 2/sqrt(pi)). This can introduce underflows when efx*x underflows even though the final result does not, so a scaled calculation with 8*efx is done in these cases - but 8 is not a big enough scale factor to avoid all such underflows. 16 is (any underflows with a scale factor of 16 would only occur when the final result underflows), so this patch changes the code to use that factor. Rather than recomputing all the values of the efx8 variable, it is removed, leaving it to the compiler's constant folding to compute 16*efx. As such scaling can also lose underflows when the final scaling down happens to be exact, appropriate checks are added to ensure underflow exceptions occur when required in such cases. Tested x86_64 and x86; no ulps updates needed. Also spot-checked for powerpc32 and mips64 to verify the changes to the ldbl-128ibm and ldbl-128 implementations. (auto-libm-test-out diffs omitted below.) 2014-05-07 Joseph Myers [BZ #16516] * sysdeps/ieee754/dbl-64/s_erf.c (efx8): Remove variable. (__erf): Scale by 16 instead of 8 in potentially underflowing case. Ensure exception if result actually underflows. * sysdeps/ieee754/flt-32/s_erff.c (efx8): Remove variable. (__erff): Scale by 16 instead of 8 in potentially underflowing case. Ensure exception if result actually underflows. * sysdeps/ieee754/ldbl-128/s_erfl.c: Include . (efx8): Remove variable. (__erfl): Scale by 16 instead of 8 in potentially underflowing case. Ensure exception if result actually underflows. * sysdeps/ieee754/ldbl-128ibm/s_erfl.c: Include . (efx8): Remove variable. (__erfl): Scale by 16 instead of 8 in potentially underflowing case. Ensure exception if result actually underflows. * sysdeps/ieee754/ldbl-96/s_erfl.c: Include . (efx8): Remove variable. (__erfl): Scale by 16 instead of 8 in potentially underflowing case. Ensure exception if result actually underflows. * math/auto-libm-test-in: Add more tests of erf. * math/auto-libm-test-out: Regenerated. diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in index 420f657..da93a30 100644 --- a/math/auto-libm-test-in +++ b/math/auto-libm-test-in @@ -796,6 +796,10 @@ erf 4.125 erf 27.0 erf -27.0 erf -0x1.fffffffffffff8p-2 +erf 0x1.c5bf94p-127 +erf 0x3.8b7fa8p-128 +erf -0x3.8b7f12369ded8p-1024 +erf 0x3.8b7f12369ded5518p-16384 erfc 0.0 erfc -0 diff --git a/sysdeps/ieee754/dbl-64/s_erf.c b/sysdeps/ieee754/dbl-64/s_erf.c index 3f37397..ea0a734 100644 --- a/sysdeps/ieee754/dbl-64/s_erf.c +++ b/sysdeps/ieee754/dbl-64/s_erf.c @@ -128,7 +128,6 @@ static const double * Coefficients for approximation to erf on [0,0.84375] */ efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */ - efx8 = 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */ pp[] = { 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */ -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */ -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */ @@ -211,7 +210,16 @@ __erf (double x) if (ix < 0x3e300000) /* |x|<2**-28 */ { if (ix < 0x00800000) - return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ + { + /* Avoid spurious underflow. */ + double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); + if (fabs (ret) < DBL_MIN) + { + double force_underflow = ret * ret; + math_force_eval (force_underflow); + } + return ret; + } return x + efx * x; } z = x * x; diff --git a/sysdeps/ieee754/flt-32/s_erff.c b/sysdeps/ieee754/flt-32/s_erff.c index 7c09589..d23cff9 100644 --- a/sysdeps/ieee754/flt-32/s_erff.c +++ b/sysdeps/ieee754/flt-32/s_erff.c @@ -33,7 +33,6 @@ erx = 8.4506291151e-01, /* 0x3f58560b */ * Coefficients for approximation to erf on [0,0.84375] */ efx = 1.2837916613e-01, /* 0x3e0375d4 */ -efx8= 1.0270333290e+00, /* 0x3f8375d4 */ pp0 = 1.2837916613e-01, /* 0x3e0375d4 */ pp1 = -3.2504209876e-01, /* 0xbea66beb */ pp2 = -2.8481749818e-02, /* 0xbce9528f */ @@ -111,8 +110,16 @@ float __erff(float x) if(ix < 0x3f580000) { /* |x|<0.84375 */ if(ix < 0x31800000) { /* |x|<2**-28 */ if (ix < 0x04000000) - /*avoid underflow */ - return (float)0.125*((float)8.0*x+efx8*x); + { + /* Avoid spurious underflow. */ + float ret = 0.0625f * (16.0f * x + (16.0f * efx) * x); + if (fabsf (ret) < FLT_MIN) + { + float force_underflow = ret * ret; + math_force_eval (force_underflow); + } + return ret; + } return x + efx*x; } z = x*x; diff --git a/sysdeps/ieee754/ldbl-128/s_erfl.c b/sysdeps/ieee754/ldbl-128/s_erfl.c index ef65ed8..35ca8c1 100644 --- a/sysdeps/ieee754/ldbl-128/s_erfl.c +++ b/sysdeps/ieee754/ldbl-128/s_erfl.c @@ -97,6 +97,7 @@ */ #include +#include #include #include @@ -143,9 +144,7 @@ tiny = 1e-4931L, one = 1.0L, two = 2.0L, /* 2/sqrt(pi) - 1 */ - efx = 1.2837916709551257389615890312154517168810E-1L, - /* 8 * (2/sqrt(pi) - 1) */ - efx8 = 1.0270333367641005911692712249723613735048E0L; + efx = 1.2837916709551257389615890312154517168810E-1L; /* erf(x) = x + x R(x^2) @@ -782,7 +781,16 @@ __erfl (long double x) if (ix < 0x3fc60000) /* |x|<2**-57 */ { if (ix < 0x00080000) - return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ + { + /* Avoid spurious underflow. */ + long double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); + if (fabsl (ret) < LDBL_MIN) + { + long double force_underflow = ret * ret; + math_force_eval (force_underflow); + } + return ret; + } return x + efx * x; } y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c index 95dc415..f55e8b7 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_erfl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_erfl.c @@ -102,6 +102,7 @@ */ #include +#include #include #include #include @@ -149,9 +150,7 @@ tiny = 1e-300L, one = 1.0L, two = 2.0L, /* 2/sqrt(pi) - 1 */ - efx = 1.2837916709551257389615890312154517168810E-1L, - /* 8 * (2/sqrt(pi) - 1) */ - efx8 = 1.0270333367641005911692712249723613735048E0L; + efx = 1.2837916709551257389615890312154517168810E-1L; /* erf(x) = x + x R(x^2) @@ -801,10 +800,17 @@ __erfl (long double x) if (ix < 0x00800000) { /* erf (-0) = -0. Unfortunately, for IBM extended double - 0.125 * (8.0 * x + efx8 * x) for x = -0 evaluates to 0. */ + 0.0625 * (16.0 * x + (16.0 * efx) * x) for x = -0 + evaluates to 0. */ if (x == 0) return x; - return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ + long double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); + if (fabsl (ret) < LDBL_MIN) + { + long double force_underflow = ret * ret; + math_force_eval (force_underflow); + } + return ret; } return x + efx * x; } diff --git a/sysdeps/ieee754/ldbl-96/s_erfl.c b/sysdeps/ieee754/ldbl-96/s_erfl.c index 47e4b9e..c27de81 100644 --- a/sysdeps/ieee754/ldbl-96/s_erfl.c +++ b/sysdeps/ieee754/ldbl-96/s_erfl.c @@ -105,6 +105,7 @@ #include +#include #include #include @@ -120,8 +121,6 @@ tiny = 1e-4931L, */ /* 2/sqrt(pi) - 1 */ efx = 1.2837916709551257389615890312154517168810E-1L, - /* 8 * (2/sqrt(pi) - 1) */ - efx8 = 1.0270333367641005911692712249723613735048E0L, pp[6] = { 1.122751350964552113068262337278335028553E6L, @@ -272,7 +271,16 @@ __erfl (long double x) if (ix < 0x3fde8000) /* |x|<2**-33 */ { if (ix < 0x00080000) - return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ + { + /* Avoid spurious underflow. */ + long double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); + if (fabsl (ret) < LDBL_MIN) + { + long double force_underflow = ret * ret; + math_force_eval (force_underflow); + } + return ret; + } return x + efx * x; } z = x * x;