[PATCH 08/11] math: Use erff from CORE-MATH
DJ Delorie
dj@redhat.com
Wed Nov 20 23:10:53 GMT 2024
Adhemerval Zanella <adhemerval.zanella@linaro.org> writes:
> The CORE-MATH implementation is correctly rounded (for any rounding mode)
> and shows better performance to the generic erff.
> diff --git a/SHARED-FILES b/SHARED-FILES
> index d367f4b62f..ccc5179f80 100644
> --- a/SHARED-FILES
> +++ b/SHARED-FILES
> @@ -272,3 +272,7 @@ sysdeps/ieee754/flt-32/s_cbrtf.c
> (file src/binary32/cbrt/cbrtf.c in CORE-MATH)
> - The code was adapted to use glibc code style and internal
> functions to handle errno, overflow, and underflow.
> +sysdeps/ieee754/flt-32/s_erff.c
> + (file src/binary32/erf/erff.c in CORE-MATH)
> + - The code was adapted to use glibc code style and internal
> + functions to handle errno, overflow, and underflow.
Ok.
> diff --git a/sysdeps/aarch64/libm-test-ulps b/sysdeps/aarch64/libm-test-ulps
> diff --git a/sysdeps/alpha/fpu/libm-test-ulps b/sysdeps/alpha/fpu/libm-test-ulps
> diff --git a/sysdeps/arc/fpu/libm-test-ulps b/sysdeps/arc/fpu/libm-test-ulps
> diff --git a/sysdeps/arc/nofpu/libm-test-ulps b/sysdeps/arc/nofpu/libm-test-ulps
> diff --git a/sysdeps/arm/libm-test-ulps b/sysdeps/arm/libm-test-ulps
> diff --git a/sysdeps/csky/fpu/libm-test-ulps b/sysdeps/csky/fpu/libm-test-ulps
> diff --git a/sysdeps/csky/nofpu/libm-test-ulps b/sysdeps/csky/nofpu/libm-test-ulps
> diff --git a/sysdeps/hppa/fpu/libm-test-ulps b/sysdeps/hppa/fpu/libm-test-ulps
> diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
> diff --git a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps
Ok.
> diff --git a/sysdeps/ieee754/flt-32/s_erff.c b/sysdeps/ieee754/flt-32/s_erff.c
> index 6c541dba23..762f160e9f 100644
> --- a/sysdeps/ieee754/flt-32/s_erff.c
> +++ b/sysdeps/ieee754/flt-32/s_erff.c
> @@ -1,155 +1,256 @@
> -/* s_erff.c -- float version of s_erf.c.
> - */
> +/* Correctly-rounded error function for binary32 value.
>
> -/*
> - * ====================================================
> - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
> - *
> - * Developed at SunPro, a Sun Microsystems, Inc. business.
> - * Permission to use, copy, modify, and distribute this
> - * software is freely granted, provided that this notice
> - * is preserved.
> - * ====================================================
> - */
> +Copyright (c) 2022-2024 Alexei Sibidanov.
>
> -#if defined(LIBM_SCCS) && !defined(lint)
> -static char rcsid[] = "$NetBSD: s_erff.c,v 1.4 1995/05/10 20:47:07 jtc Exp $";
> -#endif
> +This file is part of the CORE-MATH project
> +project (file src/binary32/erf/erff.c revision bc385c2).
> +
> +Permission is hereby granted, free of charge, to any person obtaining a copy
> +of this software and associated documentation files (the "Software"), to deal
> +in the Software without restriction, including without limitation the rights
> +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
> +copies of the Software, and to permit persons to whom the Software is
> +furnished to do so, subject to the following conditions:
> +
> +The above copyright notice and this permission notice shall be included in all
> +copies or substantial portions of the Software.
> +
> +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
> +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
> +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
> +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
> +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
> +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
> +SOFTWARE.
> +*/
Ok.
> -#include <errno.h>
> -#include <float.h>
> #include <math.h>
> -#include <math-narrow-eval.h>
> -#include <math_private.h>
> -#include <math-underflow.h>
> +#include <stdint.h>
> #include <libm-alias-float.h>
> -#include <fix-int-fp-convert-zero.h>
> +#include "math_config.h"
Ok.
> -static const float
> -tiny = 1e-30,
> -one = 1.0000000000e+00, /* 0x3F800000 */
> -erx = 8.4506291151e-01, /* 0x3f58560b */
> -/*
> - * Coefficients for approximation to erf on [0,0.84375]
> - */
> -efx = 1.2837916613e-01, /* 0x3e0375d4 */
> -pp0 = 1.2837916613e-01, /* 0x3e0375d4 */
> -pp1 = -3.2504209876e-01, /* 0xbea66beb */
> -pp2 = -2.8481749818e-02, /* 0xbce9528f */
> -pp3 = -5.7702702470e-03, /* 0xbbbd1489 */
> -pp4 = -2.3763017452e-05, /* 0xb7c756b1 */
> -qq1 = 3.9791721106e-01, /* 0x3ecbbbce */
> -qq2 = 6.5022252500e-02, /* 0x3d852a63 */
> -qq3 = 5.0813062117e-03, /* 0x3ba68116 */
> -qq4 = 1.3249473704e-04, /* 0x390aee49 */
> -qq5 = -3.9602282413e-06, /* 0xb684e21a */
> -/*
> - * Coefficients for approximation to erf in [0.84375,1.25]
> - */
> -pa0 = -2.3621185683e-03, /* 0xbb1acdc6 */
> -pa1 = 4.1485610604e-01, /* 0x3ed46805 */
> -pa2 = -3.7220788002e-01, /* 0xbebe9208 */
> -pa3 = 3.1834661961e-01, /* 0x3ea2fe54 */
> -pa4 = -1.1089469492e-01, /* 0xbde31cc2 */
> -pa5 = 3.5478305072e-02, /* 0x3d1151b3 */
> -pa6 = -2.1663755178e-03, /* 0xbb0df9c0 */
> -qa1 = 1.0642088205e-01, /* 0x3dd9f331 */
> -qa2 = 5.4039794207e-01, /* 0x3f0a5785 */
> -qa3 = 7.1828655899e-02, /* 0x3d931ae7 */
> -qa4 = 1.2617121637e-01, /* 0x3e013307 */
> -qa5 = 1.3637083583e-02, /* 0x3c5f6e13 */
> -qa6 = 1.1984500103e-02, /* 0x3c445aa3 */
> -/*
> - * Coefficients for approximation to erfc in [1.25,1/0.35]
> - */
> -ra0 = -9.8649440333e-03, /* 0xbc21a093 */
> -ra1 = -6.9385856390e-01, /* 0xbf31a0b7 */
> -ra2 = -1.0558626175e+01, /* 0xc128f022 */
> -ra3 = -6.2375331879e+01, /* 0xc2798057 */
> -ra4 = -1.6239666748e+02, /* 0xc322658c */
> -ra5 = -1.8460508728e+02, /* 0xc3389ae7 */
> -ra6 = -8.1287437439e+01, /* 0xc2a2932b */
> -ra7 = -9.8143291473e+00, /* 0xc11d077e */
> -sa1 = 1.9651271820e+01, /* 0x419d35ce */
> -sa2 = 1.3765776062e+02, /* 0x4309a863 */
> -sa3 = 4.3456588745e+02, /* 0x43d9486f */
> -sa4 = 6.4538726807e+02, /* 0x442158c9 */
> -sa5 = 4.2900814819e+02, /* 0x43d6810b */
> -sa6 = 1.0863500214e+02, /* 0x42d9451f */
> -sa7 = 6.5702495575e+00, /* 0x40d23f7c */
> -sa8 = -6.0424413532e-02, /* 0xbd777f97 */
> -/*
> - * Coefficients for approximation to erfc in [1/.35,28]
> - */
> -rb0 = -9.8649431020e-03, /* 0xbc21a092 */
> -rb1 = -7.9928326607e-01, /* 0xbf4c9dd4 */
> -rb2 = -1.7757955551e+01, /* 0xc18e104b */
> -rb3 = -1.6063638306e+02, /* 0xc320a2ea */
> -rb4 = -6.3756646729e+02, /* 0xc41f6441 */
> -rb5 = -1.0250950928e+03, /* 0xc480230b */
> -rb6 = -4.8351919556e+02, /* 0xc3f1c275 */
> -sb1 = 3.0338060379e+01, /* 0x41f2b459 */
> -sb2 = 3.2579251099e+02, /* 0x43a2e571 */
> -sb3 = 1.5367296143e+03, /* 0x44c01759 */
> -sb4 = 3.1998581543e+03, /* 0x4547fdbb */
> -sb5 = 2.5530502930e+03, /* 0x451f90ce */
> -sb6 = 4.7452853394e+02, /* 0x43ed43a7 */
> -sb7 = -2.2440952301e+01; /* 0xc1b38712 */
> -
Ok.
> -float __erff(float x)
> +float
> +__erff (float x)
> {
Ok.
> - int32_t hx,ix,i;
> - float R,S,P,Q,s,y,z,r;
> - GET_FLOAT_WORD(hx,x);
> - ix = hx&0x7fffffff;
> - if(ix>=0x7f800000) { /* erf(nan)=nan */
> - i = ((uint32_t)hx>>31)<<1;
> - return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */
> - }
> -
> - if(ix < 0x3f580000) { /* |x|<0.84375 */
> - if(ix < 0x31800000) { /* |x|<2**-28 */
> - if (ix < 0x04000000)
> - {
> - /* Avoid spurious underflow. */
> - float ret = 0.0625f * (16.0f * x + (16.0f * efx) * x);
> - math_check_force_underflow (ret);
> - return ret;
> - }
> - return x + efx*x;
> - }
> - z = x*x;
> - r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
> - s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
> - y = r/s;
> - return x + x*y;
> - }
> - if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
> - s = fabsf(x)-one;
> - P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
> - Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
> - if(hx>=0) return erx + P/Q; else return -erx - P/Q;
> - }
> - if (ix >= 0x40c00000) { /* inf>|x|>=6 */
> - if(hx>=0) return one-tiny; else return tiny-one;
> - }
> - x = fabsf(x);
> - s = one/(x*x);
> - if(ix< 0x4036DB6E) { /* |x| < 1/0.35 */
> - R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
> - ra5+s*(ra6+s*ra7))))));
> - S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
> - sa5+s*(sa6+s*(sa7+s*sa8)))))));
> - } else { /* |x| >= 1/0.35 */
> - R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
> - rb5+s*rb6)))));
> - S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
> - sb5+s*(sb6+s*sb7))))));
> - }
> - GET_FLOAT_WORD(ix,x);
> - SET_FLOAT_WORD(z,ix&0xfffff000);
> - r = __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S);
> - if(hx>=0) return one-r/x; else return r/x-one;
Ok.
> + /* for 7 <= i < 63, C[i-7] is a degree-7 polynomial approximation of
> + erf(i/16+1/32+x) for -1/32 <= x <= 1/32 */
> + static const double C[56][8] = {
> + { 0x1.f86faa9428f9cp-2, 0x1.cfc41e36c7dfap-1, -0x1.b2c7dc53508b9p-2,
> + -0x1.5a9de93fa556ep-3, 0x1.731793dbb01b5p-3, 0x1.133e06426cf18p-6,
> + -0x1.a12a6289cafd8p-5, 0x1.717d6f1d6f557p-9 },
> . . .
> + { 0x1.fffffe2ba0ea5p-1, 0x1.d06ad6ecde88ep-22, -0x1.be46aa8edc9a1p-20,
> + 0x1.143860c7840b8p-18, -0x1.edaba78fb1260p-18, 0x1.52138a96ecee2p-17,
> + -0x1.6fca538c4e2eep-17, 0x1.434040640bcefp-17 },
> + { 0x1.fffffee3cc32cp-1, 0x1.1e1e857adb8ddp-22, -0x1.1769ce5f2a6e8p-20,
> + 0x1.5fe5d479b0543p-19, -0x1.405d865c94c2ap-18, 0x1.bfc94feb96afcp-18,
> + -0x1.f245d5f3e8358p-18, 0x1.c142456acf443p-18 },
> + };
Ok.
> + float ax = fabsf (x);
> + uint32_t ux = asuint (ax);
> + double s = x;
> + double z = ax;
Ok.
> + /* 0x407ad444 corresponds to x = 0x1.f5a888p+1 = 3.91921..., which is the
> + largest float such that erf(x) does not round to 1 (to nearest). */
> + if (__glibc_unlikely (ux > 0x407ad444u))
> + {
> + float os = copysignf (1.0f, x);
> + if (ux > (0xffu << 23))
> + return x + x; /* nan */
> + if (ux == (0xffu << 23))
> + return os; /* +-inf */
> + return os - 0x1p-25f * os;
> + }
Ok.
> + double v = floor (16.0 * z);
> + uint32_t i = 16.0f * ax;
> + /* 0x3ee00000 corresponds to x = 0.4375, for smaller x we have i < 7. */
> + if (__glibc_unlikely (ux < 0x3ee00000u))
> + {
> + static const double c[] =
> + {
> + 0x1.20dd750429b6dp+0, -0x1.812746b0375fbp-2,
> + 0x1.ce2f219fd6f45p-4, -0x1.b82ce2cbf0838p-6,
> + 0x1.565bb655adb85p-8, -0x1.c025bfc879c94p-11,
> + 0x1.f81718f61309cp-14, -0x1.cc67bd88f5867p-17
> + };
> + double z2 = s * s, z4 = z2 * z2, z8 = z4 * z4;
> + double c0 = c[0] + z2 * c[1];
> + double c2 = c[2] + z2 * c[3];
> + double c4 = c[4] + z2 * c[5];
> + double c6 = c[6] + z2 * c[7];
> + c0 += z4 * c2;
> + c4 += z4 * c6;
> + c0 += z8 * c4;
> + return s * c0;
> + }
Ok.
> + z = (z - 0.03125) - 0.0625 * v;
> + const double *c = C[i - 7];
> + double z2 = z * z, z4 = z2 * z2;
> + double c0 = c[0] + z * c[1];
> + double c2 = c[2] + z * c[3];
> + double c4 = c[4] + z * c[5];
> + double c6 = c[6] + z * c[7];
> + c0 += z2 * c2;
> + c4 += z2 * c6;
> + c0 += z4 * c4;
> + return copysign (c0, s);
> }
> libm_alias_float (__erf, erf)
> -
Ok.
> diff --git a/sysdeps/loongarch/lp64/libm-test-ulps b/sysdeps/loongarch/lp64/libm-test-ulps
> diff --git a/sysdeps/m68k/m680x0/fpu/libm-test-ulps b/sysdeps/m68k/m680x0/fpu/libm-test-ulps
> diff --git a/sysdeps/microblaze/libm-test-ulps b/sysdeps/microblaze/libm-test-ulps
> diff --git a/sysdeps/mips/mips32/libm-test-ulps b/sysdeps/mips/mips32/libm-test-ulps
> diff --git a/sysdeps/mips/mips64/libm-test-ulps b/sysdeps/mips/mips64/libm-test-ulps
> diff --git a/sysdeps/nios2/libm-test-ulps b/sysdeps/nios2/libm-test-ulps
> diff --git a/sysdeps/or1k/fpu/libm-test-ulps b/sysdeps/or1k/fpu/libm-test-ulps
> diff --git a/sysdeps/or1k/nofpu/libm-test-ulps b/sysdeps/or1k/nofpu/libm-test-ulps
> diff --git a/sysdeps/powerpc/fpu/libm-test-ulps b/sysdeps/powerpc/fpu/libm-test-ulps
> diff --git a/sysdeps/powerpc/nofpu/libm-test-ulps b/sysdeps/powerpc/nofpu/libm-test-ulps
> diff --git a/sysdeps/riscv/nofpu/libm-test-ulps b/sysdeps/riscv/nofpu/libm-test-ulps
> diff --git a/sysdeps/riscv/rvd/libm-test-ulps b/sysdeps/riscv/rvd/libm-test-ulps
> diff --git a/sysdeps/s390/fpu/libm-test-ulps b/sysdeps/s390/fpu/libm-test-ulps
> diff --git a/sysdeps/sh/libm-test-ulps b/sysdeps/sh/libm-test-ulps
> diff --git a/sysdeps/sparc/fpu/libm-test-ulps b/sysdeps/sparc/fpu/libm-test-ulps
> diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
Ok.
LGTM
Reviewed-by: DJ Delorie <dj@redhat.com>
More information about the Libc-alpha
mailing list