[PATCH 11/17] math: Use expm1f from CORE-MATH
Adhemerval Zanella Netto
adhemerval.zanella@linaro.org
Thu Oct 31 14:35:41 GMT 2024
On 30/10/24 16:27, DJ Delorie wrote:
>
> 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.
Yes, math/s_csin_template.c calls __ieee754_sinhf which in turn __expm1. For
'0xcp-4 + 0x1.4p+0 * I' input, __expm1 is called with input of '0x1.4p+0'.
The correctly rounded value, which CORE-MATH returns, is 0x1.3ec38ep+1 while
current glibc one returns 0x1.3ec39p+1.
>
>> 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