[PATCH 14/15] math: Use tanpif from CORE-MATH
Adhemerval Zanella Netto
adhemerval.zanella@linaro.org
Tue Feb 11 14:09:15 GMT 2025
On 10/02/25 22:15, DJ Delorie wrote:
>
> I've suggested some readability changes below, but they're cosmetic.
>
> LGTM
Ack, Paul has added them on CORE-MATH and I will sync them.
>
> Reviewed-by: DJ Delorie <dj@redhat.com>
>
> Adhemerval Zanella <adhemerval.zanella@linaro.org> writes:
>> diff --git a/SHARED-FILES b/SHARED-FILES
>> (src/binary32/sinpi/sinpif.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_tanpif.c:
>> + (src/binary32/tanpi/tanpif.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/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/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/math_config.h b/sysdeps/ieee754/flt-32/math_config.h
>> index 035461199c..8d9c8ee3ad 100644
>> --- a/sysdeps/ieee754/flt-32/math_config.h
>> +++ b/sysdeps/ieee754/flt-32/math_config.h
>> @@ -84,6 +84,31 @@ roundeven_finite (double x)
>> #endif
>> }
>>
>> +#ifndef ROUNDEVENF_INTRINSICS
>> +/* When set, roundevenf_finite will route to the internal roundevenf function. */
>> +# define ROUNDEVENF_INTRINSICS 1
>> +#endif
>
> So ROUNDEVENF_INTRINSICS will only ever be 0 or 1 (default), but always defined.
>
>> +static inline float
>> +roundevenf_finite (float x)
>> +{
>> + if (!isfinite (x))
>> + __builtin_unreachable ();
>> +#if ROUNDEVENF_INTRINSICS
>> + return roundevenf (x);
>
> THe default is to just be roundevenf() then.
>
>> +#else
>> + float y = roundf (x);
>> + if (fabs (x - y) == 0.5)
>> + {
>> + union { float f; uint32_t i; } u = {y};
>> + union { float f; uint32_t i; } v = {y - copysignf (1.0, x)};
>> + if (__builtin_ctzl (v.i) > __builtin_ctzl (u.i))
>> + y = v.f;
>> + }
>> + return y;
>> +#endif
>> +}
>
> But the caller can choose this inlined version. Ok.
>
>> diff --git a/sysdeps/ieee754/flt-32/s_tanpif.c b/sysdeps/ieee754/flt-32/s_tanpif.c
>> +/* Correctly-rounded tangent of binary32 value for angles in half-revolutions
>> +
>> +Copyright (c) 2022-2025 Alexei Sibidanov.
>> +
>> +The original version of this file was copied from the CORE-MATH
>> +project (src/binary32/tanpi/tanpif.c, revision 3bbf907).
>> +
>> +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 <stdint.h>
>> +#include <errno.h>
>> +#include <libm-alias-float.h>
>> +#include "math_config.h"
>
> Ok.
>
>> +float
>> +__tanpif (float x)
>> +{
>> + uint32_t ix = asuint (x);
>> + uint32_t e = ix & (0xff << 23);
>> + if (__glibc_unlikely (e > (150 << 23)))
>
> or > 2 ^ 23, ok.
>
>> + {
>> + if (e == (0xff << 23))
>> + {
>> + if (!(ix << 9))
>> + return __math_invalidf (x);
>> + return x + x; /* nan */
>> + }
>> + return copysign (0.0f, x);
>> + }
>
> Ok.
>
>> + float x4 = 4.0f * x;
>> + float nx4 = roundevenf_finite (x4);
>> + float dx4 = x4 - nx4;
>> + float ni = roundevenf_finite (x);
>> + float zf = x - ni;
>> + if (__glibc_unlikely (dx4 == 0.0f))
>> + {
>> + int k = x4;
>> + if (k & 1)
>> + return copysignf (1.0f, zf);
>> + k &= 7;
>
> For the sake of readability, this should be 6, not 7. It makes the
> assumption of k&1==0 explicit instead of implicit.
>
>> + if (k == 0)
>> + return copysignf (0.0f, x);
>> + if (k == 4)
>> + return -copysignf (0.0f, x);
>> + __set_errno (ERANGE);
>> + if (k == 2)
>> + return 1.0f / 0.0f;
>> + if (k == 6)
>
> This last comparison is superfluous, and could be included as "/* if
> (...) */ so it's obvious that all cases were accounted for instead of
> thinking "is there a fallthrough case?"
>
>> + return -1.0f / 0.0f;
>> + }
>
> Ok.
>
>> + ix = asuint (zf);
>> + uint32_t a = ix & (~0u >> 1);
>
>> + if (__glibc_unlikely (a == 0x3e933802u))
>> + return copysignf (0x1.44cfbap+0f, zf) + copysignf (0x1p-25f, zf);
>> + if (__glibc_unlikely (a == 0x38f26685u))
>> + return copysignf (0x1.7cc304p-12, zf) + copysignf (0x1p-37f, zf);
>
> Are these special cases where the code below doesn't round correctly?
> If so, a comment explaining such is warranted.
>
>> + double z = zf, z2 = z * z;
>> +
>> + static const double cn[] = { 0x1.921fb54442d19p-1, -0x1.1f458b3e1f8d6p-2,
>> + 0x1.68a34bd0b8f6ap-6, -0x1.e4866f7a25f99p-13 };
>> + static const double cd[] = { 0x1p+0, -0x1.4b4b98d2df3a7p-1,
>> + 0x1.8e9926d2bb901p-4, -0x1.a6f77fd847eep-9 };
>> + double z4
>> + = z2 * z2,
>> + r = (z - z * z2) * ((cn[0] + z2 * cn[1]) + z4 * (cn[2] + z2 * cn[3]))
>> + / (((cd[0] + z2 * cd[1]) + z4 * (cd[2] + z2 * cd[3])) * (0.25 - z2));
>> + return r;
>> +}
>> +libm_alias_float (__tanpi, tanpi)
>
> Ok.
>
>> diff --git a/sysdeps/loongarch/lp64/libm-test-ulps b/sysdeps/loongarch/lp64/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
>
> Ok.
>
>> diff --git a/sysdeps/powerpc/fpu/math_private.h b/sysdeps/powerpc/fpu/math_private.h
>> index aace1a8708..7065d276c0 100644
>> --- a/sysdeps/powerpc/fpu/math_private.h
>> +++ b/sysdeps/powerpc/fpu/math_private.h
>> @@ -62,6 +62,7 @@ __ieee754_sqrtf128 (_Float128 __x)
>> #ifdef _ARCH_PWR6
>> /* ISA 2.03 provides frin/round() and cntlzw/ctznll(). */
>> # define ROUNDEVEN_INTRINSICS 0
>> +# define ROUNDEVENF_INTRINSICS 0
>> #endif
>
> So PowerPC doesn't have a roundevenf() intrinsic and will use the
> special code. Ok.
>
>> 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/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