Date: Sun, 02 Aug 2026 12:10:51 +0000 From: Robert Clausecker <fuz@FreeBSD.org> To: src-committers@FreeBSD.org, dev-commits-src-all@FreeBSD.org, dev-commits-src-main@FreeBSD.org Cc: Steve Kargl <kargl@FreeBSD.org> Subject: git: ae417b3194e7 - main - msun: add asinpi, acospi, and atanpi Message-ID: <6a6f33cb.1e802.25c5cbb1@gitrepo.freebsd.org>
index | next in thread | raw e-mail
The branch main has been updated by fuz: URL: https://cgit.FreeBSD.org/src/commit/?id=ae417b3194e76ce26065dc20281493ee83619879 commit ae417b3194e76ce26065dc20281493ee83619879 Author: Steve Kargl <kargl@FreeBSD.org> AuthorDate: 2026-07-29 11:45:42 +0000 Commit: Robert Clausecker <fuz@FreeBSD.org> CommitDate: 2026-08-02 11:17:45 +0000 msun: add asinpi, acospi, and atanpi This commit implements the inverse half-cycle trigonometric functions: asinpi(x) = asin(x) / pi Eq. (1) acospi(x) = acos(x) / pi atanpi(x) = atan(x) / pi Implemention details are contained in src/s_asinpi.c and src/a_atanpi.c, where the details for acospi(x) appear in the former. ************* CAVEAT EMPTOR: The ld128 code has been only compiled. It has not been tested for correctness due to lack of hardware. ************* Code compiled on AMD Ryzen 7 7700X system run FreeBSD 16.0-CURRENT main-n284956-de9fe28ab847. Exhaustive testing of acospif(x), asinpif(x), and atanpif(x) on the indicated intervals yields % ./tlibm acospi -fPE -x 0x1p-120 -X 1 Interval tested for acospif: [7.52316e-37,1] ulp <= 0.5: 99.627% 1002878299 | 99.627% 1002878299 0.5 < ulp <= 0.6: 0.277% 2789599 | 99.904% 1005667898 0.6 < ulp <= 0.7: 0.096% 965062 | 100.000% 1006632960 Max ulp: 0.63661975 at 5.96046412e-08 0x1.fffffep-25 % ./tlibm asinpi -fPED -x 0x1p-120f -X 1.f Interval tested for asinpif: [7.52316e-37,1] ulp <= 0.5: 99.851% 1005129353 | 99.851% 1005129353 0.5 < ulp <= 0.6: 0.149% 1501097 | 100.000% 1006630450 0.6 < ulp <= 0.7: 0.000% 2510 | 100.000% 1006632960 Max ulp: 0.68957579 at 5.04878759e-01 0x1.027f78p-1 % ./tlibm atanpi -fPE -x 0x1p-120 -X max > zatanpif.txt & Interval tested for atanpif: [7.52316e-37,3.40282e+38] ulp <= 0.5: 99.865% 2077574602 | 99.865% 2077574602 0.5 < ulp <= 0.6: 0.131% 2735011 | 99.997% 2080309613 0.6 < ulp <= 0.7: 0.003% 65170 | 100.000% 2080374783 Max ulp: 0.68433094 at 5.01186252e-01 0x1.009b7cp-1 Testing the double and long double version cannot be done in an exhaustive manner. For 300 M values, uniformily distributed in the indicated interals, one finds the max ULP: Interval tested for acospi: [9.31323e-10,0.25] xm = 2.4423788416892520e-01, /* 0x3fcf432f, 0xde79920f */ libm = 4.2146222480005391e-01, /* 0x3fdaf93c, 0xb201001c */ mpfr = 4.2146222480005396e-01, /* 0x3fdaf93c, 0xb201001d */ ULP = 0.50499351466286857 Interval tested for acospi: [0.25,0.5] xm = 4.9689430915631438e-01, /* 0x3fdfcd1d, 0xc9d945c6 */ libm = 3.3447366122373884e-01, /* 0x3fd56804, 0x371513ef */ mpfr = 3.3447366122373889e-01, /* 0x3fd56804, 0x371513f0 */ ULP = 0.57195275455053829 Interval tested for acospi: [0.5,0.75] xm = 5.0238623667462079e-01, /* 0x3fe0138c, 0x4d0f4be0 */ libm = 3.3245556599062825e-01, /* 0x3fd546f3, 0xb5d36303 */ mpfr = 3.3245556599062820e-01, /* 0x3fd546f3, 0xb5d36302 */ ULP = 0.63427929243758807 Interval tested for acospi: [0.75,1] xm = 7.5853651919512177e-01, /* 0x3fe845ee, 0x60d8789f */ libm = 2.2591472240382732e-01, /* 0x3fcceac6, 0x0c3465ce */ mpfr = 2.2591472240382729e-01, /* 0x3fcceac6, 0x0c3465cd */ ULP = 0.56915750216472161 Interval tested for asinpi: [9.31323e-10,0.25] xm = 1.9502362835488171e-01, /* 0x3fc8f688, 0xc4dda0fb */ libm = 6.2478354989018887e-02, /* 0x3faffd29, 0xb6c57c61 */ mpfr = 6.2478354989018881e-02, /* 0x3faffd29, 0xb6c57c60 */ ULP = 0.52347765415885006 Interval tested for asinpi: [0.25,0.5] xm = 4.9937103583123676e-01, /* 0x3fdff5b1, 0xeeddbf62 */ libm = 1.6643553767987129e-01, /* 0x3fc54dc2, 0x7b9d15a4 */ mpfr = 1.6643553767987126e-01, /* 0x3fc54dc2, 0x7b9d15a3 */ ULP = 0.66214688371031072 Interval tested for asinpi: [0.5,0.75] xm = 5.0228515250761718e-01, /* 0x3fe012b8, 0x4fe92bbb */ libm = 1.6750722213679006e-01, /* 0x3fc570e0, 0x6c75edd5 */ mpfr = 1.6750722213679009e-01, /* 0x3fc570e0, 0x6c75edd6 */ ULP = 0.78223048105528226 Interval tested for asinpi: [0.75,1] xm = 7.5425933001419776e-01, /* 0x3fe822e4, 0x7663a4aa */ libm = 2.7200385380185182e-01, /* 0x3fd16882, 0xda1dc13b */ mpfr = 2.7200385380185188e-01, /* 0x3fd16882, 0xda1dc13c */ ULP = 0.53747973176773822 Interval tested for atanpi: [9.31323e-10,0.25] xm = 1.9666113418757322e-01, /* 0x3fc92c31, 0x29dd6d2f */ libm = 6.1810387818117797e-02, /* 0x3fafa59c, 0x7476baa5 */ mpfr = 6.1810387818117804e-02, /* 0x3fafa59c, 0x7476baa6 */ ULP = 0.54674297446584263 Interval tested for atanpi: [0.25,0.5] xm = 4.1312119637707068e-01, /* 0x3fda7093, 0xe2ee5494 */ libm = 1.2470309560460152e-01, /* 0x3fbfec8a, 0xc554ebec */ mpfr = 1.2470309560460154e-01, /* 0x3fbfec8a, 0xc554ebed */ ULP = 0.73116638175113347 Interval tested for atanpi: [0.5,0.75] xm = 5.0018949583396499e-01, /* 0x3fe0018d, 0x66cd1b82 */ libm = 1.4763186871058706e-01, /* 0x3fc2e599, 0xdffacb8f */ mpfr = 1.4763186871058709e-01, /* 0x3fc2e599, 0xdffacb90 */ ULP = 0.69192753950764663 Interval tested for atanpi: [0.75,1] xm = 7.5007880583359599e-01, /* 0x3fe800a5, 0x448f4c03 */ libm = 2.0484881828445453e-01, /* 0x3fca387c, 0x6f93f71f */ mpfr = 2.0484881828445450e-01, /* 0x3fca387c, 0x6f93f71e */ ULP = 0.65765471872064396 Interval tested for atanpi: [1,2] xm = 1.0103228000344093e+00, /* 0x3ff02a48, 0x3d88d0a2 */ libm = 2.5163447403817019e-01, /* 0x3fd01ac7, 0x7b229108 */ mpfr = 2.5163447403817013e-01, /* 0x3fd01ac7, 0x7b229107 */ ULP = 0.67409519689166042 Interval tested for atanpi: [2,4] xm = 2.0231383267437946e+00, /* 0x40002f63, 0x25a530a9 */ libm = 3.5387589538123299e-01, /* 0x3fd6a5e7, 0x156053c6 */ mpfr = 3.5387589538123293e-01, /* 0x3fd6a5e7, 0x156053c5 */ ULP = 0.69695587476021503 Interval tested for atanpi: [4,1.79769e+308] xm = 4.0000000000000000e+00, /* 0x40100000, 0x00000000 */ libm = 4.2202086962263069e-01, /* 0x3fdb0263, 0xd2508e31 */ mpfr = 4.2202086962263069e-01, /* 0x3fdb0263, 0xd2508e31 */ ULP = 0.27709400511686716 PR: 295884 MFC after: 1 month Reviewed by: fuz --- lib/msun/Makefile | 14 ++- lib/msun/Symbol.map | 9 ++ lib/msun/ld128/s_asinpil.c | 173 +++++++++++++++++++++++++++++++ lib/msun/ld128/s_atanpil.c | 157 ++++++++++++++++++++++++++++ lib/msun/ld80/s_asinpil.c | 191 ++++++++++++++++++++++++++++++++++ lib/msun/ld80/s_atanpil.c | 175 +++++++++++++++++++++++++++++++ lib/msun/man/acospi.3 | 77 ++++++++++++++ lib/msun/man/asinpi.3 | 77 ++++++++++++++ lib/msun/man/atanpi.3 | 73 +++++++++++++ lib/msun/src/math.h | 9 ++ lib/msun/src/math_private.h | 38 +++++-- lib/msun/src/s_asinpi.c | 243 ++++++++++++++++++++++++++++++++++++++++++++ lib/msun/src/s_asinpif.c | 157 ++++++++++++++++++++++++++++ lib/msun/src/s_atanpi.c | 215 +++++++++++++++++++++++++++++++++++++++ lib/msun/src/s_atanpif.c | 142 ++++++++++++++++++++++++++ 15 files changed, 1735 insertions(+), 15 deletions(-) diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 52848965ae64..745819f5268f 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -97,9 +97,8 @@ COMMON_SRCS= b_tgamma.c \ w_cabs.c w_cabsf.c w_drem.c w_dremf.c # IEEE-754 2008 and ISO/IEC TS 18661-4 half-cycle trignometric functions -COMMON_SRCS+= s_cospi.c s_cospif.c \ - s_sinpi.c s_sinpif.c \ - s_tanpi.c s_tanpif.c +COMMON_SRCS+= s_asinpi.c s_asinpif.c s_atanpi.c s_atanpif.c \ + s_cospi.c s_cospif.c s_sinpi.c s_sinpif.c s_tanpi.c s_tanpif.c # LIBCSRCDIR is defined lazily, so can't be used for exists() below LIBC_SRCTOP=${SRCTOP}/lib/libc @@ -135,7 +134,8 @@ COMMON_SRCS+= b_tgammal.c catrigl.c \ e_lgammal.c e_lgammal_r.c e_powl.c \ e_remainderl.c e_sinhl.c e_sqrtl.c \ invtrig.c k_cosl.c k_sinl.c k_tanl.c \ - s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cexpl.c \ + s_asinhl.c s_asinpil.c s_atanl.c s_atanpil.c \ + s_cbrtl.c s_ceill.c s_cexpl.c \ s_clogl.c s_cosl.c s_cospil.c s_cprojl.c \ s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \ s_fmaxl.c s_fmaximuml.c s_fmaximum_magl.c \ @@ -180,7 +180,8 @@ SRCS= ${COMMON_SRCS} ${ARCH_SRCS} INCS+= fenv.h math.h -MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \ +MAN= acos.3 acosh.3 acospi.3 asin.3 asinh.3 asinpi.3 \ + atan.3 atan2.3 atanh.3 atanpi.3 \ ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \ cimag.3 clog.3 copysign.3 cos.3 cosh.3 cospi.3 \ cpow.3 csqrt.3 erf.3 \ @@ -198,10 +199,13 @@ MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \ MLINKS+=acos.3 acosf.3 acos.3 acosl.3 MLINKS+=acosh.3 acoshf.3 acosh.3 acoshl.3 +MLINKS+=acospi.3 acospif.3 acospi.3 acospil.3 MLINKS+=asin.3 asinf.3 asin.3 asinl.3 MLINKS+=asinh.3 asinhf.3 asinh.3 asinhl.3 +MLINKS+=asinpi.3 asinpif.3 asinpi.3 asinpil.3 MLINKS+=atan.3 atanf.3 atan.3 atanl.3 MLINKS+=atanh.3 atanhf.3 atanh.3 atanhl.3 +MLINKS+=atanpi.3 atanpif.3 atanpi.3 atanpil.3 MLINKS+=atan2.3 atan2f.3 atan2.3 atan2l.3 \ atan2.3 carg.3 atan2.3 cargf.3 atan2.3 cargl.3 MLINKS+=cacos.3 cacosf.3 cacos.3 cacosl.3 \ diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map index 2484cc13013b..6f96c0d4f206 100644 --- a/lib/msun/Symbol.map +++ b/lib/msun/Symbol.map @@ -320,6 +320,15 @@ FBSD_1.7 { /* First added in 16.0-CURRENT */ FBSD_1.9 { + acospi; + acospif; + acospil; + asinpi; + asinpif; + asinpil; + atanpi; + atanpif; + atanpil; fmaximum; fmaximumf; fmaximuml; diff --git a/lib/msun/ld128/s_asinpil.c b/lib/msun/ld128/s_asinpil.c new file mode 100644 index 000000000000..d9dbee927da7 --- /dev/null +++ b/lib/msun/ld128/s_asinpil.c @@ -0,0 +1,173 @@ +/*- + * SPDX-License-Identifier: BSD-2-Clause + * + * Copyright (c) 2026 Steven G. Kargl + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice unmodified, this list of conditions, and the following + * disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/* + * src/s_asinpi.c for implemenation details. + */ + +#include "math.h" +#include "math_private.h" + +#define _CC (0x1p57L + 1) +#define _ROOT sqrtl + +volatile static const double tiny = 1.e-300; +static const double half = 0.5, one = 1.; + +/* Full precision high and low parts of 1 / pi. */ +static const long double +invpihi = 3.18309886183790671537767526745028737e-01L, +invpilo = -1.28821588763206006125693864783127482e-35L; + +/* + * Prior to the leading multiplication by x^2, the rational approximation + * has an absolute minimax error less than 2.86e-38 over the [0x1p-56,0.5] + * domain (or log2(error) = -124.7). + */ +static inline long double +__r(long double xs) +{ + static const long double + R0 = 5.30516476972984452562945877908381207e-02L, + R1 = -2.60404681812983888677486288898051068e-01L, + R2 = 5.43274368732204127849880302812989126e-01L, + R3 = -6.27476039646895838848725721148947475e-01L, + R4 = 4.37806223811161460670580937913479055e-01L, + R5 = -1.88854683672201011042310131486279147e-01L, + R6 = 4.94528645546233985859708428238164776e-02L, + R7 = -7.38228286468784921885881430313681146e-03L, + R8 = 5.47609957562826794808966842094709768e-04L, + R9 = -1.44837427671490633843418618030301533e-05L, + R10= 1.60087061374239702655150492329207605e-08L, + S1 = -5.35851261206434695164324875722635942e+00L, + S2 = 1.23839541267281630422432341920191869e+01L, + S3 = -1.61473998442126533731585929134597693e+01L, + S4 = 1.30442314991246752613457576698550221e+01L, + S5 = -6.74685395302389922545194846871599697e+00L, + S6 = 2.22958089278066716953022085768628547e+00L, + S7 = -4.55313732093612723575589438727154778e-01L, + S8 = 5.33391487207251214709627478847204628e-02L, + S9 = -3.08343599417750507732125775680767303e-03L, + S10= 6.12260957946655623349049151336307244e-05L; + + long double r, s; + r = R5 + (R6 + (R7 + (R8 + (R9 + R10 * xs) * xs) *xs) * xs) * xs; + r = R0 + (R1 + (R2 + (R3 + (R4 + r * xs) * xs) * xs) * xs) * xs; + s = S5 + (S6 + (S7 + (S8 + (S9 + S10 * xs) * xs) *xs) * xs) * xs; + s = 1 + (S1 + (S2 + (S3 + (S4 + r * xs) * xs) * xs) * xs) * xs; + return (xs * (r / s)); +} + +long double +asinpil(long double x) +{ + long double ax, hi, lo, xh, xl, y, zh, zl; + + if (isnan(x) || isinf(x)) + return ((x - x) / (x - x)); + + ax = fabsl(x); + + if (ax > 1) /* |x| > 1 */ + return ((x - x) / (x - x)); + + + if (ax <= 0.5) { /* |x| <= 0.5 */ + if (ax < 0x1p-57L) { /* |x| < 0x1p-57 */ + if (ax < 0x1p-16340L) { /* |x| < 0x1p-16340 */ + if (ax == 0) + return (x); + /* Scale for near subnormal. */ + ax *= 0x1p114; + _XMUL(ax, 0, invpihi, invpilo, hi, lo); + y = (hi + lo) * 0x1p-114; + } else { + _XMUL(ax, 0, invpihi, invpilo, hi, lo); + y = hi + lo; + } + } else { + y = __r(ax * ax); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(ax, 0, xh, xl, hi, lo); + y = hi + lo; + } + } else if (ax < 1) { /* |x| < 1 */ + y = 1 - ax; + x = __r(y / 2); + _XADD(invpihi, invpilo, x, 0, xh, xl); + _SQRT(2 * y, zh, zl); + _XMUL(xh, xl, zh, zl, hi, lo); + _XADD(half, 0, -hi, -lo, y, x); + } else /* |x| == 1 */ + y = half; + + return (x < 0 ? -y : y); +} + +/* + * See src/s_asinpi.c for implementation details. + */ + +long double +acospil(long double x) +{ + long double ax, hi, lo, xh, xl, y, zh, zl; + + if (isnan(x) || isinf(x)) + return ((x - x) / (x - x)); + + ax = fabsl(x); + + if (ax > 1) /* |x| > 1 */ + return ((x - x) / (x - x)); + + if (ax <= 0.5) { /* |x| <= 0.5 */ + if (ax < 0x1p-55L) { /* |x| < 0x1p-55 */ + y = (ax == 0) ? half : (ax < 0x1p-113 ? + half - tiny : half - x * invpihi); + } else { + y = __r(x * x); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(x, 0, xh, xl, hi, lo); + _XADD(half, 0, -hi, -lo, y, ax); + } + } else if (ax < 1) { /* |x| < 1 */ + y = 1 - ax; + ax = __r(y / 2); + _XADD(invpihi, invpilo, ax, 0, xh, xl); + _SQRT(2 * y, zh, zl); + _XMUL(xh, xl, zh, zl, hi, lo); + if (x < 0) + _XADD(one, 0, -hi, -lo, y, ax); + else + y = hi + lo; + } else /* |x| == 1 */ + y = half; + + return (y); +} diff --git a/lib/msun/ld128/s_atanpil.c b/lib/msun/ld128/s_atanpil.c new file mode 100644 index 000000000000..3f3b72785bb1 --- /dev/null +++ b/lib/msun/ld128/s_atanpil.c @@ -0,0 +1,157 @@ +/*- + * SPDX-License-Identifier: BSD-2-Clause + * + * Copyright (c) 2026 Steven G. Kargl + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice unmodified, this list of conditions, and the following + * disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/* + * src/s_atanpi.c for implemenation details. + */ + +#include "math.h" +#include "math_private.h" + +#define _CC (0x1p57L + 1) +#define _ROOT sqrtl + +volatile static const double tiny = 1.e-300; +static const double half = 0.5, one = 1., qrtr = 0.25; +static const double x0 = 0.625, x1 = 0.875, x2 = 1.5; + +/* Full precision high and low parts. */ +static const long double +invpihi = 3.18309886183790671537767526745028737e-01L, /* 1/pi */ +invpilo = -1.28821588763206006125693864783127482e-35L, /* 1/pi */ +a0hi = 1.77807684489352753115503587502248118e-01L, /* atanpi(x0) */ +a0lo = -1.04565241713355884231841538298149212e-35L, /* atanpi(x0) */ +a1hi = 2.28810695365053587806047702039884546e-01L, /* atanpi(x1) */ +a1lo = -5.01451892949247960667313017509908719e-37L, /* atanpi(x1) */ +a2hi = 3.12832958189001183813747252435221446e-01L, /* atanpi(x2) */ +a2lo = 1.35071436051692259858806390941906106e-36L; /* atanpi(x2) */ + +/* + * Prior to the leading multiplication by x^2, the rational approximation + * has an absolute minimax error less than 2.86e-38 over the [0x1p-56,0.5] + * domain (or log2(error) = -124.7). + */ +static inline long double +__r(long double xs) +{ + static const long double + R0 = 5.30516476972984452562945877908381207e-02L, + R1 = -2.60404681812983888677486288898051068e-01L, + R2 = 5.43274368732204127849880302812989126e-01L, + R3 = -6.27476039646895838848725721148947475e-01L, + R4 = 4.37806223811161460670580937913479055e-01L, + R5 = -1.88854683672201011042310131486279147e-01L, + R6 = 4.94528645546233985859708428238164776e-02L, + R7 = -7.38228286468784921885881430313681146e-03L, + R8 = 5.47609957562826794808966842094709768e-04L, + R9 = -1.44837427671490633843418618030301533e-05L, + R10= 1.60087061374239702655150492329207605e-08L, + S1 = -5.35851261206434695164324875722635942e+00L, + S2 = 1.23839541267281630422432341920191869e+01L, + S3 = -1.61473998442126533731585929134597693e+01L, + S4 = 1.30442314991246752613457576698550221e+01L, + S5 = -6.74685395302389922545194846871599697e+00L, + S6 = 2.22958089278066716953022085768628547e+00L, + S7 = -4.55313732093612723575589438727154778e-01L, + S8 = 5.33391487207251214709627478847204628e-02L, + S9 = -3.08343599417750507732125775680767303e-03L, + S10= 6.12260957946655623349049151336307244e-05L; + + long double r, s; + r = R5 + (R6 + (R7 + (R8 + (R9 + R10 * xs) * xs) *xs) * xs) * xs; + r = R0 + (R1 + (R2 + (R3 + (R4 + r * xs) * xs) * xs) * xs) * xs; + s = S5 + (S6 + (S7 + (S8 + (S9 + S10 * xs) * xs) *xs) * xs) * xs; + s = 1 + (S1 + (S2 + (S3 + (S4 + r * xs) * xs) * xs) * xs) * xs; + return (xs * (r / s)); +} + +long double +atanpil(long double x) +{ + long double ax, hi, lo, xh, xl, y, zh, zl; + + if (isnan(x) || isinf(x)) + return ((x - x) / (x - x)); + + ax = fabsl(x); + + if (ax > 1) /* |x| > 1 */ + return ((x - x) / (x - x)); + + + if (ax <= 0.5) { /* |x| <= 0.5 */ + if (ax < 0x1p-57L) { /* |x| < 0x1p-57 */ + if (ax < 0x1p-16340L) { /* |x| < 0x1p-16340 */ + if (ax == 0) + return (x); + /* Scale for near subnormal. */ + ax *= 0x1p114; + _XMUL(ax, 0, invpihi, invpilo, hi, lo); + y = (hi + lo) * 0x1p-114; + } else { + _XMUL(ax, 0, invpihi, invpilo, hi, lo); + y = hi + lo; + } + } else { + y = __r(ax * ax); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(ax, 0, xh, xl, hi, lo); + y = hi + lo; + } + } else if (ax < 1) { /* |x| < 1 */ + if (ax < 0.75) { /* |x| < 0.75 */ + x = (ax - x0) / (1 + x0 * ax); + y = __r(x * x); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(x, 0, xh, xl, hi, lo); + _XADD(a0hi, a0lo, hi, lo, y, xl); + } else { + x = (ax - x1) / (1 + x1 * ax); + y = __r(x * x); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(x, 0, xh, xl, hi, lo); + _XADD(a1hi, a1lo, hi, lo, y, xl); + } + } else if (ax < 2) { /* |x| < 2 */ + if (ax == 1) + return (x < 0 ? -qrtr : qrtr); + x = (ax - x2) / (1 + x2 * ax); + y = __r(x * x); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(x, 0, xh, xl, hi, lo); + _XADD(a2hi, a2lo, hi, lo, y, xl); + } else { /* |x| > 2 */ + x = 1 / ax; + y = __r(x * x); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(x, 0, xh, xl, hi, lo); + _XADD(half, 0, -hi, -lo, y, x); + } + + return (x < 0 ? -y : y); +} diff --git a/lib/msun/ld80/s_asinpil.c b/lib/msun/ld80/s_asinpil.c new file mode 100644 index 000000000000..a806c47bd970 --- /dev/null +++ b/lib/msun/ld80/s_asinpil.c @@ -0,0 +1,191 @@ +/*- + * SPDX-License-Identifier: BSD-2-Clause + * + * Copyright (c) 2026 Steven G. Kargl + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice unmodified, this list of conditions, and the following + * disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/* + * src/s_asinpi.c for implemenation details. + */ + +#ifdef __i386__ +#include <ieeefp.h> +#endif +#include <stdint.h> + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" + +#define _CC (0x1p32L + 1) +#define _ROOT sqrtl + +volatile static const double tiny = 1.e-300; +static const double half = 0.5, one = 1.; + +/* 1/pi split into the leading and trailing 53 bits. */ +static const double +invpihi = 3.1830988618379069e-01, +invpilo = -1.9678676675182486e-17; + +/* + * Prior to the leading multiplication by x^2, the rational approximation + * has an absolute minimax error less than 1.36e-22 over the [0x1p-32,0.5] + * domain (or log2(error) = -72.6). + */ +static inline long double +__r(long double xs) +{ + static const union IEEEl2bits + R0u = LD80C(0xd94caf3dbdb01c38, -5, 5.30516476972984452564e-02L), + R1u = LD80C(0x8d346599ebe3212a, -3, -1.37895190734499743665e-01L), + R2u = LD80C(0x852c1751d7112655, -3, 1.30051006670116064731e-01L), + R3u = LD80C(0xdb55697553742657, -5, -5.35482520546937456640e-02L), + R4u = LD80C(0x931ba448fbd6e8c3, -7, 8.97875827280130762861e-03L), + R5u = LD80C(0xdbbe8e80762881b0, -12, -4.19129108247665906395e-04L), + S1u = LD80C(0xc3272074944c6389, 1, -3.04926310906120631911e+00L), + S2u = LD80C(0xe390d58f48c1a96d, 1, 3.55571497910116337692e+00L), + S3u = LD80C(0xfccb6664d1fbfc39, 0, -1.97495727465499715865e+00L), + S4u = LD80C(0x86f4f140e3be416b, -1, 5.27175024358933283746e-01L), + S5u = LD80C(0xf218070935962bdb, -5, -5.91049456446391222176e-02L), + S6u = LD80C(0xec656d26215b624b, -10, 1.80355985054588957542e-03L); + +#define R0 (R0u.e) +#define R1 (R1u.e) +#define R2 (R2u.e) +#define R3 (R3u.e) +#define R4 (R4u.e) +#define R5 (R5u.e) +#define S1 (S1u.e) +#define S2 (S2u.e) +#define S3 (S3u.e) +#define S4 (S4u.e) +#define S5 (S5u.e) +#define S6 (S6u.e) + + long double r, s; + r = R0 + (R1 + (R2 + (R3 + (R4 + R5 * xs) * xs) * xs) * xs) * xs; + s = 1 + (S1 + (S2 + (S3 + (S4 + (S5 + S6 * xs) * xs) * xs) * + xs) * xs) * xs; + return (xs * (r / s)); +} + +#define GREATER(a) (ix == a && lx > 0x8000000000000000ull) +#define LESSEQ(a) (ix == a && lx <= 0x8000000000000000ull) + +long double +asinpil(long double x) +{ + long double ax, hi, lo, xh, xl, y, zh, zl; + uint64_t lx; + uint16_t hx, ix; + + EXTRACT_LDBL80_WORDS(hx, lx, x); + ix = hx & 0x7fff; + + if (ix >= 0x4000 || GREATER(0x3fff)) /* |x| > 1 */ + return ((x - x) / (x - x)); + + ENTERI(); + + INSERT_LDBL80_WORDS(ax, ix, lx); + + if (ix < 0x3ffe || LESSEQ(0x3ffe)) { /* |x| <= 0.5 */ + if (ix < 0x3fde) { /* |x| < 0x1p-33 */ + if (ix < 0x002b) { /* |x| < 0x1p-16340 */ + if ((ix | lx) == 0) + RETURNI(x); + /* Scale for near subnormal. */ + ax *= 0x1p65; + _XMUL(ax, 0, invpihi, invpilo, hi, lo); + y = (hi + lo) * 0x1p-65; + } else { + _XMUL(ax, 0, invpihi, invpilo, hi, lo); + y = hi + lo; + } + } else { + y = __r(ax * ax); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(ax, 0, xh, xl, hi, lo); + y = hi + lo; + } + } else if (ix < 0x3fff) { /* |x| < 1 */ + y = 1 - ax; + x = __r(y / 2); + _XADD(invpihi, invpilo, x, 0, xh, xl); /* 1 / pi + r(t^2) */ + _SQRT(2 * y, zh, zl); /* 2 * t */ + _XMUL(xh, xl, zh, zl, hi, lo); + _XADD(half, 0, -hi, -lo, y, x); + } else /* |x| == 1 */ + y = half; + + RETURNI((hx & 0x8000) ? -y : y); +} + +/* + * See src/s_asinpi.c for implementation details. + */ + +long double +acospil(long double x) +{ + long double ax, hi, lo, xh, xl, y, zh, zl; + uint64_t lx; + uint16_t hx, ix; + + EXTRACT_LDBL80_WORDS(hx, lx, x); + ix = hx & 0x7fff; + + if (ix >= 0x4000 || GREATER(0x3fff)) /* |x| > 1 */ + return ((x - x) / (x - x)); + + ENTERI(); + + if (ix < 0x3ffe || LESSEQ(0x3ffe)) { /* |x| <= 0.5 */ + if (ix < 0x3fe9) { /* |x| < 0x1p-22 */ + y = ((ix | lx) == 0) ? half : (LESSEQ(0x3fbf) ? + half - tiny : half - x * invpihi); + } else { + y = __r(x * x); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(x, 0, xh, xl, hi, lo); + _XADD(half, 0, -hi, -lo, y, ax); + } + } else if (ix < 0x3fff) { /* |x| < 1 */ + INSERT_LDBL80_WORDS(ax, ix, lx); + y = 1 - ax; + ax = __r(y / 2); + _XADD(invpihi, invpilo, ax, 0, xh, xl); /* 1 / pi + r(t^2) */ + _SQRT(2 * y, zh, zl); /* 2 * t */ + _XMUL(xh, xl, zh, zl, hi, lo); + if (hx & 0x8000) + _XADD(one, 0, -hi, -lo, y, ax); + else + y = hi + lo; + } else /* |x| == 1 */ + y = half; + + RETURNI(y); +} diff --git a/lib/msun/ld80/s_atanpil.c b/lib/msun/ld80/s_atanpil.c new file mode 100644 index 000000000000..27f5fee6dd32 --- /dev/null +++ b/lib/msun/ld80/s_atanpil.c @@ -0,0 +1,175 @@ +/*- + * SPDX-License-Identifier: BSD-2-Clause + * + * Copyright (c) 2026 Steven G. Kargl + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice unmodified, this list of conditions, and the following + * disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/* + * src/s_atanpi.c for implemenation details. + */ + +#ifdef __i386__ +#include <ieeefp.h> +#endif +#include <stdint.h> + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" + +#define _CC (0x1p32L + 1) +#define _ROOT sqrtl +#define NBIT (0x8000000000000000ull) + +volatile static const double tiny = 1.e-300; +static const double half = 0.5, one = 1., qrtr = 0.25; +static const double x0 = 0.625, x1 = 0.875, x2 = 1.5; + +/* 53-bit high and low parts. */ +static const double +invpihi = 3.1830988618379069e-01, /* 1/pi */ +invpilo = -1.9678676675182486e-17, /* 1/pi */ +a0hi = 1.7780768448935275e-01, /* atanpi(x0) */ +a0lo = 6.7223942595197191e-18, /* atanpi(x0) */ +a1hi = 2.2881069536505358e-01, /* atanpi(x1) */ +a1lo = 8.7193139538130510e-18, /* atanpi(x1) */ +a2hi = 3.1283295818900120e-01, /* atanpi(x2) */ +a2lo = -1.4076885713501453e-17; /* atanpi(x2) */ + +/* + * Prior to the leading multiplication by x^2, the rational approximation + * has an absolute minimax error less than 1.22e-23 over the [0x1p-32,0.5] + * domain (or log2(error) = -76.1). + */ +static inline long double +__r(long double xs) +{ + static const union IEEEl2bits + R0u = LD80C(0xd94caf3dbdb01c38, -4, -1.06103295394596890513e-01L), + R1u = LD80C(0x845b6d12f35f0ccb, -2, -2.58510025561632250608e-01L), + R2u = LD80C(0xe679d75bf98585ae, -3, -2.25074162472429082081e-01L), + R3u = LD80C(0xaba63bb6f14e17e4, -4, -8.38131585316324449243e-02L), + R4u = LD80C(0xca768e4a7d2fd52e, -7, -1.23573674736290235872e-02L), + R5u = LD80C(0x802e09b751d4c6bc, -11, -4.88967268965928487608e-04L), + S1u = LD80C(0xc2545ef3d335e598, 1, 3.03639959155120062705e+00L), + S2u = LD80C(0xe0ee42f67b6b24bc, 1, 3.51454233236806805243e+00L), + S3u = LD80C(0xf7200820c8d271d1, 0, 1.93066503144077077959e+00L), + S4u = LD80C(0x820ec69663c6f1f3, -1, 5.08037959781883845928e-01L), + S5u = LD80C(0xe61a8fde47e2f517, -5, 5.61776752333507006525e-02L), + S6u = LD80C(0xe260573083b4d35b, -10, 1.72711433720654801154e-03L); + +#define R0 (R0u.e) +#define R1 (R1u.e) +#define R2 (R2u.e) +#define R3 (R3u.e) +#define R4 (R4u.e) +#define R5 (R5u.e) +#define S1 (S1u.e) +#define S2 (S2u.e) +#define S3 (S3u.e) +#define S4 (S4u.e) +#define S5 (S5u.e) +#define S6 (S6u.e) + + long double r, s; + r = R0 + (R1 + (R2 + (R3 + (R4 + R5 * xs) * xs) * xs) * xs) * xs; + s = 1 + (S1 + (S2 + (S3 + (S4 + (S5 + S6 * xs) * xs) * xs) * + xs) * xs) * xs; + return (xs * (r / s)); +} + +long double +atanpil(long double x) +{ + long double ax, hi, lo, xh, xl, y, zh, zl; + uint64_t lx; + uint16_t hx, ix; + + EXTRACT_LDBL80_WORDS(hx, lx, x); + ix = hx & 0x7fff; + + /* x = +-inf, nan */ + if (ix >= 0x7fff && lx >= 0x8000000000000000ull) { + if (lx > 0x8000000000000000ull) + return (x + x); + return ((hx & 0x8000) ? -half : half); + } + + ENTERI(); + + INSERT_LDBL80_WORDS(ax, ix, lx); + + if (ix < 0x3ffe ) { /* |x| < 0.5 */ + if (ix < 0x3fde) { /* |x| < 0x1p-33 */ + if (ix < 0x002b) { /* |x| < 0x1p-16340 */ + if ((ix | lx) == 0) + RETURNI(x); + /* Scale for near subnormal. */ + ax *= 0x1p65; + _XMUL(ax, 0, invpihi, invpilo, hi, lo); + y = (hi + lo) * 0x1p-65; + } else { + _XMUL(ax, 0, invpihi, invpilo, hi, lo); + y = hi + lo; + } + } else { + y = __r(ax * ax); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(ax, 0, xh, xl, hi, lo); + y = hi + lo; + } + } else if (ax < 1) { /* |x| < 1 */ + /* |x| < 0.75 */ + if (ix == 0x3ffe && lx < 0xc000000000000000ull) { + x = (ax - x0) / (1 + x0 * ax); + y = __r(x * x); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(x, 0, xh, xl, hi, lo); + _XADD(a0hi, a0lo, hi, lo, y, xl); + } else { + x = (ax - x1) / (1 + x1 * ax); + y = __r(x * x); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(x, 0, xh, xl, hi, lo); + _XADD(a1hi, a1lo, hi, lo, y, xl); + } + } else if (ix < 0x4000) { /* |x| < 2 */ + if (ix == 0x3fff && lx == NBIT) + return ((hx & 0x8000) ? -qrtr : qrtr); + x = (ax - x2) / (1 + x2 * ax); + y = __r(x * x); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(x, 0, xh, xl, hi, lo); + _XADD(a2hi, a2lo, hi, lo, y, xl); + } else { /* |x| > 2 */ + x = 1 / ax; + y = __r(x * x); + _XADD(invpihi, invpilo, y, 0, xh, xl); + _XMUL(x, 0, xh, xl, hi, lo); + _XADD(half, 0, -hi, -lo, y, x); + } + + RETURNI((hx & 0x8000) ? -y : y); +} diff --git a/lib/msun/man/acospi.3 b/lib/msun/man/acospi.3 new file mode 100644 index 000000000000..f5d68b379022 --- /dev/null +++ b/lib/msun/man/acospi.3 @@ -0,0 +1,77 @@ +.\" Copyright (c) 2026 Steven G. Kargl. +.\" +.\" Redistribution and use in source and binary forms, with or without +.\" modification, are permitted provided that the following conditions +.\" are met: +.\" 1. Redistributions of source code must retain the above copyright +.\" notice, this list of conditions and the following disclaimer. +.\" 2. Redistributions in binary form must reproduce the above copyright +.\" notice, this list of conditions and the following disclaimer in the +.\" documentation and/or other materials provided with the distribution. +.\" +.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND +.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +.\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE +.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) *** 1101 LINES SKIPPED ***home | help
Want to link to this message? Use this
URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?6a6f33cb.1e802.25c5cbb1>