[PATCH 11/17] math: Use expm1f from CORE-MATH
DJ Delorie
dj@redhat.com
Wed Oct 30 19:27:24 GMT 2024
LGTM
Reviewed-by: DJ Delorie <dj@redhat.com>
Adhemerval Zanella <adhemerval.zanella@linaro.org> writes:
> diff --git a/SHARED-FILES b/SHARED-FILES
> +sysdeps/ieee754/flt-32/s_expm1f.c
> + (file src/binary32/expm1/expm1f.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
> Function: Imaginary part of "csin":
> +float: 1
> ldouble: 1
I assume this is because csin uses expm1 internally.
> Function: Real part of "csin_downward":
> @@ -1163,7 +1164,6 @@ float: 1
>
> Function: "expm1":
> double: 1
> -float: 1
> ldouble: 2
Ok.
> 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
Ok.
> diff --git a/sysdeps/arc/nofpu/libm-test-ulps b/sysdeps/arc/nofpu/libm-test-ulps
>
> Function: "expm1":
> double: 1
> -float: 1
>
I assume the rest were :0 anyway... ok.
> diff --git a/sysdeps/arm/libm-test-ulps b/sysdeps/arm/libm-test-ulps
Ok.
> 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
Ok.
> diff --git a/sysdeps/i386/fpu/s_expm1f.S b/sysdeps/i386/fpu/s_expm1f.S
> deleted file mode 100644
Ok.
> diff --git a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps
> index 9afb63bc0f..6d9d49d153 100644
> --- a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps
> +++ b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps
> @@ -1237,19 +1237,16 @@ ldouble: 3
>
> Function: "expm1_downward":
> double: 1
> -float: 1
> float128: 2
> ldouble: 4
>
> Function: "expm1_towardzero":
> double: 1
> -float: 1
> float128: 4
> ldouble: 4
>
> Function: "expm1_upward":
> double: 1
> -float: 1
> float128: 3
> ldouble: 4
Ok.
> diff --git a/sysdeps/ieee754/flt-32/s_expm1f.c b/sysdeps/ieee754/flt-32/s_expm1f.c
> index 35f7b5214a..edd7c9acf8 100644
> --- a/sysdeps/ieee754/flt-32/s_expm1f.c
> +++ b/sysdeps/ieee754/flt-32/s_expm1f.c
> @@ -1,132 +1,124 @@
> -/* s_expm1f.c -- float version of s_expm1.c.
> - */
> +/* Correctly-rounded natural exponent function biased by 1 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.
> +
> +This file is part of the CORE-MATH project
> +project (file src/binary32/expm1/expm1f.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.
> +*/
>
> -#include <errno.h>
> -#include <float.h>
> #include <math.h>
> -#include <math-barriers.h>
> -#include <math_private.h>
> #include <math-underflow.h>
> #include <libm-alias-float.h>
> -
> -static const float huge = 1.0e+30;
> -static const float tiny = 1.0e-30;
Ok.
> -static const float
> -one = 1.0,
> -o_threshold = 8.8721679688e+01,/* 0x42b17180 */
> -ln2_hi = 6.9313812256e-01,/* 0x3f317180 */
> -ln2_lo = 9.0580006145e-06,/* 0x3717f7d1 */
> -invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */
> - /* scaled coefficients related to expm1 */
> -Q1 = -3.3333335072e-02, /* 0xbd088889 */
> -Q2 = 1.5873016091e-03, /* 0x3ad00d01 */
> -Q3 = -7.9365076090e-05, /* 0xb8a670cd */
> -Q4 = 4.0082177293e-06, /* 0x36867e54 */
> -Q5 = -2.0109921195e-07; /* 0xb457edbb */
> +#include "math_config.h"
Ok.
> float
> -__expm1f(float x)
> +__expm1f (float x)
> {
Ok.
> - float y,hi,lo,c,t,e,hxs,hfx,r1;
> - int32_t k,xsb;
> - uint32_t hx;
> -
> - GET_FLOAT_WORD(hx,x);
> - xsb = hx&0x80000000; /* sign bit of x */
> - if(xsb==0) y=x; else y= -x; /* y = |x| */
> - hx &= 0x7fffffff; /* high word of |x| */
> -
> - /* filter out huge and non-finite argument */
> - if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */
> - if(hx >= 0x42b17218) { /* if |x|>=88.721... */
> - if(hx>0x7f800000)
> - return x+x; /* NaN */
> - if(hx==0x7f800000)
> - return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
> - if(x > o_threshold) {
> - __set_errno (ERANGE);
> - return huge*huge; /* overflow */
> - }
> - }
> - if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */
> - math_force_eval(x+tiny);/* raise inexact */
> - return tiny-one; /* return -1 */
> - }
> - }
> -
> - /* argument reduction */
> - if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
> - if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */
> - if(xsb==0)
> - {hi = x - ln2_hi; lo = ln2_lo; k = 1;}
> - else
> - {hi = x + ln2_hi; lo = -ln2_lo; k = -1;}
> - } else {
> - k = invln2*x+((xsb==0)?(float)0.5:(float)-0.5);
> - t = k;
> - hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
> - lo = t*ln2_lo;
> - }
> - x = hi - lo;
> - c = (hi-x)-lo;
> - }
> - else if(hx < 0x33000000) { /* when |x|<2**-25, return x */
> - math_check_force_underflow (x);
> - t = huge+x; /* return x with inexact flags when x!=0 */
> - return x - (t-(huge+x));
Ok.
> + static const double c[] =
> + {
> + 1, 0x1.62e42fef4c4e7p-6, 0x1.ebfd1b232f475p-13, 0x1.c6b19384ecd93p-20
> + };
> + static const double ch[] =
> + {
> + 0x1.62e42fefa39efp-6, 0x1.ebfbdff82c58fp-13, 0x1.c6b08d702e0edp-20,
> + 0x1.3b2ab6fb92e5ep-27, 0x1.5d886e6d54203p-35, 0x1.430976b8ce6efp-43
> + };
> + static const double td[] =
> + {
> + 0x1p+0, 0x1.059b0d3158574p+0, 0x1.0b5586cf9890fp+0,
> + 0x1.11301d0125b51p+0, 0x1.172b83c7d517bp+0, 0x1.1d4873168b9aap+0,
> + 0x1.2387a6e756238p+0, 0x1.29e9df51fdee1p+0, 0x1.306fe0a31b715p+0,
> + 0x1.371a7373aa9cbp+0, 0x1.3dea64c123422p+0, 0x1.44e086061892dp+0,
> + 0x1.4bfdad5362a27p+0, 0x1.5342b569d4f82p+0, 0x1.5ab07dd485429p+0,
> + 0x1.6247eb03a5585p+0, 0x1.6a09e667f3bcdp+0, 0x1.71f75e8ec5f74p+0,
> + 0x1.7a11473eb0187p+0, 0x1.82589994cce13p+0, 0x1.8ace5422aa0dbp+0,
> + 0x1.93737b0cdc5e5p+0, 0x1.9c49182a3f09p+0, 0x1.a5503b23e255dp+0,
> + 0x1.ae89f995ad3adp+0, 0x1.b7f76f2fb5e47p+0, 0x1.c199bdd85529cp+0,
> + 0x1.cb720dcef9069p+0, 0x1.d5818dcfba487p+0, 0x1.dfc97337b9b5fp+0,
> + 0x1.ea4afa2a490dap+0, 0x1.f50765b6e454p+0
> + };
> + const double iln2 = 0x1.71547652b82fep+5;
> + const double big = 0x1.8p52;
> + double z = x;
> + uint32_t ux = asuint (x);
> + uint32_t ax = ux << 1;
> + if (__glibc_likely (ax < 0x7c400000u))
> + { /* |x| < 0.15625 */
> + if (__glibc_unlikely (ax < 0x676a09e8u))
> + { /* |x| < 0x1.6a09e8p-24 */
> + if (__glibc_unlikely (ax == 0x0u))
> + return x; /* x = +-0 */
> + return fmaf (fabsf (x), 0x1p-25f, x);
> }
Ok.
> - else k = 0;
> -
> - /* x is now in primary range */
> - hfx = (float)0.5*x;
> - hxs = x*hfx;
> - r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
> - t = (float)3.0-r1*hfx;
> - e = hxs*((r1-t)/((float)6.0 - x*t));
> - if(k==0) return x - (x*e-hxs); /* c is 0 */
> - else {
> - e = (x*(e-c)-c);
> - e -= hxs;
> - if(k== -1) return (float)0.5*(x-e)-(float)0.5;
> - if(k==1) {
> - if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5));
> - else return one+(float)2.0*(x-e);
> - }
> - if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
> - int32_t i;
> - y = one-(e-x);
> - GET_FLOAT_WORD(i,y);
> - SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */
> - return y-one;
> - }
> - t = one;
> - if(k<23) {
> - int32_t i;
> - SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
> - y = t-(e-x);
> - GET_FLOAT_WORD(i,y);
> - SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */
> - } else {
> - int32_t i;
> - SET_FLOAT_WORD(t,((0x7f-k)<<23)); /* 2^-k */
> - y = x-(e+t);
> - y += one;
> - GET_FLOAT_WORD(i,y);
> - SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */
> - }
Ok.
> + static const double b[] =
> + {
> + 0x1.fffffffffffc2p-2, 0x1.55555555555fep-3, 0x1.555555559767fp-5,
> + 0x1.1111111098dc1p-7, 0x1.6c16bca988aa9p-10, 0x1.a01a07658483fp-13,
> + 0x1.a05b04d2c3503p-16, 0x1.71de3a960b5e3p-19
> + };
> + double z2 = z * z, z4 = z2 * z2;
> + double r = z + z2
> + * ((b[0] + z * b[1]) + z2 * (b[2] + z * b[3])
> + + z4 * ((b[4] + z * b[5]) + z2 * (b[6] + z * b[7])));
> + return r;
> + }
> + if (__glibc_unlikely (ax >= 0x8562e430u))
> + { /* |x| > 88.72 */
> + if (ax > (0xffu << 24))
> + return x + x; /* nan */
> + if (__glibc_unlikely (ux >> 31))
> + { /* x < 0 */
> + if (ax == (0xffu << 24))
> + return -1.0f;
> + return -1.0f + 0x1p-26f;
> }
> - return y;
> + if (ax == (0xffu << 24))
> + return INFINITY;
> + return __math_oflowf (0);
> + }
Ok.
> + double a = iln2 * z;
> + double ia = roundeven (a);
> + double h = a - ia;
> + double h2 = h * h;
> + uint64_t u = asuint64 (ia + big);
> + double c2 = c[2] + h * c[3], c0 = c[0] + h * c[1];
> + const uint64_t *tdl = (uint64_t *) ((void *) td);
> + double sv = asdouble (tdl[u & 0x1f] + ((u >> 5) << 52));
> + double r = (c0 + h2 * c2) * sv - 1.0;
> + float ub = r, lb = r - sv * 0x1.3b3p-33;
> + if (__glibc_unlikely (ub != lb))
> + {
> + if (__glibc_unlikely (ux > 0xc18aa123u)) /* x < -17.32 */
> + return -1.0f + 0x1p-26f;
> + const double iln2h = 0x1.7154765p+5;
> + const double iln2l = 0x1.5c17f0bbbe88p-26;
> + double s = sv;
> + h = (iln2h * z - ia) + iln2l * z;
> + h2 = h * h;
> + double w = s * h;
> + r = (s - 1) + w
> + * ((ch[0] + h * ch[1])
> + + h2 * ((ch[2] + h * ch[3]) + h2 * (ch[4] + h * ch[5])));
> + ub = r;
> + }
> + return ub;
> }
> libm_alias_float (__expm1, expm1)
Ok.
> diff --git a/sysdeps/loongarch/lp64/libm-test-ulps b/sysdeps/loongarch/lp64/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.
More information about the Libc-alpha
mailing list