[PATCH v2 18/25] math: Use asinf from CORE-MATH
Adhemerval Zanella Netto
adhemerval.zanella@linaro.org
Wed Dec 18 13:41:19 GMT 2024
On 07/12/24 01:37, DJ Delorie wrote:
>
> Adhemerval Zanella <adhemerval.zanella@linaro.org> writes:
>> The CORE-MATH implementation is correctly rounded (for any rounding mode)
>> and shows slight better performance to the generic asinf.
>
> A few comments about moving common functions into a header file, but
> LGTM.
> Reviewed-by: DJ Delorie <dj@redhat.com>
>
>> diff --git a/SHARED-FILES b/SHARED-FILES
>> index c29c1c8106..6a6ff76b31 100644
>> --- a/SHARED-FILES
>> +++ b/SHARED-FILES
>> @@ -302,3 +302,7 @@ sysdeps/ieee754/flt-32/e_acoshf.c:
>> (src/binary32/acosh/acoshf.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/e_asinf.c:
>> + (src/binary32/asin/asinf.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
>
> Ok.
>
>> diff --git a/sysdeps/i386/fpu/e_asinf.S b/sysdeps/i386/fpu/e_asinf.S
>> deleted file mode 100644
>
> Ok.
>
>> 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/e_asinf.c b/sysdeps/ieee754/flt-32/e_asinf.c
>> index d9f446d1b1..944bf6f5ce 100644
>> --- a/sysdeps/ieee754/flt-32/e_asinf.c
>> +++ b/sysdeps/ieee754/flt-32/e_asinf.c
>> @@ -1,105 +1,131 @@
>> -/* e_asinf.c -- float version of e_asin.c.
>> - */
>> +/* Correctly-rounded arc-sine 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) 2023-2024 Alexei Sibidanov.
>>
>> -/*
>> - Modifications for single precision expansion are
>> - Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
>> - and are incorporated herein by permission of the author. The author
>> - reserves the right to distribute this material elsewhere under different
>> - copying permissions. These modifications are distributed here under
>> - the following terms:
>> +The original version of this file was copied from the CORE-MATH
>> +project (file src/binary32/asin/asinf.c, revision bc385c2).
>>
>> - This library is free software; you can redistribute it and/or
>> - modify it under the terms of the GNU Lesser General Public
>> - License as published by the Free Software Foundation; either
>> - version 2.1 of the License, or (at your option) any later version.
>> +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:
>>
>> - This library is distributed in the hope that it will be useful,
>> - but WITHOUT ANY WARRANTY; without even the implied warranty of
>> - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
>> - Lesser General Public License for more details.
>> +The above copyright notice and this permission notice shall be included in all
>> +copies or substantial portions of the Software.
>>
>> - You should have received a copy of the GNU Lesser General Public
>> - License along with this library; if not, see
>> - <https://www.gnu.org/licenses/>. */
>> +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.
>> +*/
>>
>> -#if defined(LIBM_SCCS) && !defined(lint)
>> -static char rcsid[] = "$NetBSD: e_asinf.c,v 1.5 1995/05/12 04:57:25 jtc Exp $";
>> -#endif
>> -
>> -#include <float.h>
>> -#include <math.h>
>> -#include <math_private.h>
>> -#include <math-underflow.h>
>> +#include <stdint.h>
>> +#include <errno.h>
>> #include <libm-alias-finite.h>
>> +#include "math_config.h"
>
> Ok.
>
>> -static const float
>> -one = 1.0000000000e+00, /* 0x3F800000 */
>> -huge = 1.000e+30,
>> -
>> -pio2_hi = 1.57079637050628662109375f,
>> -pio2_lo = -4.37113900018624283e-8f,
>> -pio4_hi = 0.785398185253143310546875f,
>
> Ok.
>
>> +static __attribute__ ((noinline)) float
>> +as_special (float x)
>> +{
>> + uint32_t ax = asuint (x) << 1;
>> + if (ax > (0xffu << 24))
>> + return x + x; /* nan */
>> + return __math_invalidf (0.0);
>> +}
>
> Ok.
>
>> -/* asin x = x + x^3 p(x^2)
>> - -0.5 <= x <= 0.5;
>> - Peak relative error 4.8e-9 */
>> -p0 = 1.666675248e-1f,
>> -p1 = 7.495297643e-2f,
>> -p2 = 4.547037598e-2f,
>> -p3 = 2.417951451e-2f,
>> -p4 = 4.216630880e-2f;
>
> Ok.
>
>> +static double
>> +poly12 (double z, const double *c)
>> +{
>> + double z2 = z * z;
>> + double 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];
>> + double c8 = c[8] + z * c[9];
>> + double c10 = c[10] + z * c[11];
>> + c0 += c2 * z2;
>> + c4 += c6 * z2;
>> + c8 += z2 * c10;
>> + c0 += z4 * (c4 + z4 * c8);
>> + return c0;
>> +}
>
> Ok. This is the kind of thing that belongs in a header, though, but I
> won't advocate for it if core-math doesn't do it first ;-)
I will propose this to CORE-MATH, although its idea is to provide self-contained
implementation with as minimal as possible dependencies (such as multiple files).
>
>> -float __ieee754_asinf(float x)
>> +float
>> +__ieee754_asinf (float x)
>
> Ok.
>
>> {
>> - float t,w,p,q,c,r,s;
>> - int32_t hx,ix;
>> - GET_FLOAT_WORD(hx,x);
>> - ix = hx&0x7fffffff;
>> - if(ix==0x3f800000) {
>> - /* asin(1)=+-pi/2 with inexact */
>> - return x*pio2_hi+x*pio2_lo;
>> - } else if(ix> 0x3f800000) { /* |x|>= 1 */
>> - return (x-x)/(x-x); /* asin(|x|>1) is NaN */
>> - } else if (ix<0x3f000000) { /* |x|<0.5 */
>> - if(ix<0x32000000) { /* if |x| < 2**-27 */
>> - math_check_force_underflow (x);
>> - if(huge+x>one) return x;/* return x with inexact if x!=0*/
>> - } else {
>> - t = x*x;
>> - w = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4))));
>> - return x+x*w;
>> - }
>> - }
>> - /* 1> |x|>= 0.5 */
>> - w = one-fabsf(x);
>> - t = w*0.5f;
>> - p = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4))));
>> - s = sqrtf(t);
>> - if(ix>=0x3F79999A) { /* if |x| > 0.975 */
>> - t = pio2_hi-(2.0f*(s+s*p)-pio2_lo);
>> - } else {
>> - int32_t iw;
>> - w = s;
>> - GET_FLOAT_WORD(iw,w);
>> - SET_FLOAT_WORD(w,iw&0xfffff000);
>> - c = (t-w*w)/(s+w);
>> - r = p;
>> - p = 2.0f*s*r-(pio2_lo-2.0f*c);
>> - q = pio4_hi-2.0f*w;
>> - t = pio4_hi-(p-q);
>> - }
>> - if(hx>0) return t; else return -t;
>
> Ok.
>
>> + const double pi2 = 0x1.921fb54442d18p+0;
>> + double xs = x;
>> + double r;
>> + uint32_t ax = asuint (x) << 1;
>> + if (__glibc_unlikely (ax > 0x7f << 24))
>> + return as_special(x);
>
> Ok.
>
>> + if (__glibc_likely (ax < 0x7ec29000u))
>> + {
>> + if (__glibc_unlikely (ax < 115 << 24))
>> + return fmaf (x, 0x1p-25, x);
>> + static const double b[] =
>> + {
>> + 0x1.0000000000005p+0, 0x1.55557aeca105dp-3, 0x1.3314ec3db7d12p-4,
>> + 0x1.775738a5a6f92p-5, 0x1.5d5f7ce1c8538p-8, 0x1.605c6d58740fp-2,
>> + -0x1.5728b732d73c6p+1, 0x1.f152170f151ebp+3, -0x1.f962ea3ca992ep+5,
>> + 0x1.71971e17375ap+7, -0x1.860512b4ba23p+8, 0x1.26a3b8d4bdb14p+9,
>> + -0x1.36f2ea5698b51p+9, 0x1.b3d722aebfa2ep+8, -0x1.6cf89703b1289p+7,
>> + 0x1.1518af6a65e2dp+5
>> + };
>> + double z = xs;
>> + double z2 = z * z;
>> + double z4 = z2 * z2;
>> + double z8 = z4 * z4;
>> + double z16 = z8 * z8;
>> + r = z * ((((b[0] + z2 * b[1]) + z4 * (b[2] + z2 * b[3]))
>> + + z8 * ((b[4] + z2 * b[5]) + z4 * (b[6] + z2 * b[7])))
>> + + z16 * (((b[8] + z2 * b[9]) + z4 * (b[10] + z2 * b[11]))
>> + + z8 * ((b[12] + z2 * b[13])
>> + + z4 * (b[14] + z2 * b[15]))));
>
> This could go in a header (poly16) along with poly12... but ok.
Not a bad idea, I will also propose this to CORE-MATH.
>
>> + float ub = r;
>> + float lb = r - z * 0x1.efa8ebp-31;
>> + if (ub == lb)
>> + return ub;
>> + }
>
> Ok.
>
>> + if (ax < (0x7eu << 24))
>> + {
>> + static const double c[] =
>> + {
>> + 0x1.555555555529cp-3, 0x1.333333337e0ddp-4, 0x1.6db6db3b4465ep-5,
>> + 0x1.f1c72e13ac306p-6, 0x1.6e89cebe06bc4p-6, 0x1.1c6dcf5289094p-6,
>> + 0x1.c6dbbcc7c6315p-7, 0x1.8f8dc2615e996p-7, 0x1.a5833b7bf15e8p-8,
>> + 0x1.43f44ace1665cp-6, -0x1.0fb17df881c73p-6, 0x1.07520c026b2d6p-5
>> + };
>> + double z = xs;
>> + double z2 = z * z;
>> + double c0 = poly12 (z2, c);
>> + r = z + (z * z2) * c0;
>> + }
>
> Ok.
>
>> + else
>> + {
>> + if (__glibc_unlikely (ax == 0x7e55688au))
>> + return copysignf (0x1.75b8a2p-1f, x) + copysignf (0x1p-26f, x);
>> + if (__glibc_unlikely (ax == 0x7e107434u))
>> + return copysignf (0x1.1f4b64p-1f, x) + copysignf (0x1p-26f, x);
>> + double bx = fabs (xs);
>> + double z = 1.0 - bx;
>> + double s = sqrt (z);
>> + static const double c[] =
>> + {
>> + 0x1.6a09e667f3bcbp+0, 0x1.e2b7dddff2db9p-4, 0x1.b27247ab42dbcp-6,
>> + 0x1.02995cc4e0744p-7, 0x1.5ffb0276ec8eap-9, 0x1.033885a928decp-10,
>> + 0x1.911f2be23f8c7p-12, 0x1.4c3c55d2437fdp-13, 0x1.af477e1d7b461p-15,
>> + 0x1.abd6bdff67dcbp-15, -0x1.1717e86d0fa28p-16, 0x1.6ff526de46023p-16
>> + };
>> + r = pi2 - s * poly12 (z, c);
>> + r = copysign (r, xs);
>> + }
>> + return r;
>> }
>> libm_alias_finite (__ieee754_asinf, __asinf)
>
> 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/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