imrcpan/freebsd
1
0 Fork 0
mirror of https://github.com/freebsd/freebsd-src.git synced 2026-10-28 11:25:21 +01:00
Code Issues Packages Projects Releases Wiki Activity

msun: add asinpi, acospi, and atanpi

Browse Source 
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
This commit is contained in:
Steve Kargl authored and Robert Clausecker committed 2026-08-02 13:17:45 +02:00
1 parent 00a79975c0
commit ae417b3194
15 files changed
+1735 -15
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

lib/msun/Makefile
+9 -5
View File
@@ -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 \
lib/msun/Symbol.map
+9
View File
@@ -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;
lib/msun/ld128/s_asinpil.c
+173
View File
@@ -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);
}
lib/msun/ld128/s_atanpil.c
+157
View File
@@ -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);
}
lib/msun/ld80/s_asinpil.c
+191
View File
@@ -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);
}
lib/msun/ld80/s_atanpil.c
+175
View File
@@ -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);
}
lib/msun/man/acospi.3
+77
View File
@@ -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)
.\" 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.
.\"
.Dd June 5, 2026
.Dt ACOSPI 3
.Os
.Sh NAME
.Nm acospi ,
.Nm acospif ,
.Nm acospil
.Nd half-cycle arc cosine functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft double
.Fn acospi "double x"
.Ft float
.Fn acospif "float x"
.Ft long double
.Fn acospil "long double x"
.Sh DESCRIPTION
The
.Fn acospi ,
.Fn acospif ,
and
.Fn acospil
functions compute the principal value of the arc cosine of
.Fa x ,
divided by \*(Pi, thus measuring the angle in half-cycles.
.Sh RETURN VALUES
The
.Fn acospi ,
.Fn acospif ,
and
.Fn acospil
functions return values in the range
.Bk -words
.Bq -1, +1 .
.Ek
If
.Pf \&| Ns Fa x Ns \&| > 1 ,
.Fn acospi x
raise an invalid exception and returns an \*(Na.
.Sh SEE ALSO
.Xr asinpi 3 ,
.Xr atanpi 3 ,
.Xr acos 3 ,
.Xr cospi 3 ,
.Sh STANDARDS
The
.Fn acospi ,
.Fn acospif ,
and
.Fn acospil
functions conform to
.St -isoC-2023 .
lib/msun/man/asinpi.3
+77
View File
@@ -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)
.\" 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.
.\"
.Dd June 5, 2026
.Dt ASINPI 3
.Os
.Sh NAME
.Nm asinpi ,
.Nm asinpif ,
.Nm asinpil
.Nd half-cycle arc sine functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft double
.Fn asinpi "double x"
.Ft float
.Fn asinpif "float x"
.Ft long double
.Fn asinpil "long double x"
.Sh DESCRIPTION
The
.Fn asinpi ,
.Fn asinpif ,
and
.Fn asinpil
functions compute the principal value of the arc sine of
.Fa x ,
divided by \*(Pi, thus measuring the angle in half-cycles.
.Sh RETURN VALUES
The
.Fn asinpi ,
.Fn asinpif ,
and
.Fn asinpil
functions return values in the range
.Bk -words
.Bq -1, +1 .
.Ek
If
.Pf \&| Ns Fa x Ns \&| > 1 ,
.Fn asinpi x
raises an invalid exception and returns an \*(Na.
.Sh SEE ALSO
.Xr acospi 3 ,
.Xr atanpi 3 ,
.Xr asin 3 ,
.Xr sinpi 3 ,
.Sh STANDARDS
The
.Fn asinpi ,
.Fn asinpif ,
and
.Fn asinpil
functions conform to
.St -isoC-2023 .
lib/msun/man/atanpi.3
+73
View File
@@ -0,0 +1,73 @@
.\" 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)
.\" 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.
.\"
.Dd June 5, 2026
.Dt ATANPI 3
.Os
.Sh NAME
.Nm atanpi ,
.Nm atanpif ,
.Nm atanpil
.Nd half-cycle arc tangent functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft double
.Fn atanpi "double x"
.Ft float
.Fn atanpif "float x"
.Ft long double
.Fn atanpil "long double x"
.Sh DESCRIPTION
The
.Fn atanpi ,
.Fn atanpif ,
and
.Fn atanpil
functions compute the principal value of the arc tangent of
.Fa x ,
divided by \*(Pi, thus measuring the angle in half-cycles.
.Sh RETURN VALUES
The
.Fn atanpi ,
.Fn atanpif ,
and
.Fn atanpil
functions return values in the range
.Bk -words
.Bq -1/2 , +1/2 .
.Ek
.Sh SEE ALSO
.Xr acospi 3 ,
.Xr asinpi 3 ,
.Xr atan 3 ,
.Xr tanpi 3 ,
.Sh STANDARDS
The
.Fn atanpi ,
.Fn atanpif ,
and
.Fn atanpil
functions conform to
.St -isoC-2023 .
lib/msun/src/math.h
+9
View File
@@ -511,6 +511,15 @@ void sincosl(long double, long double *, long double *);
#endif /* __BSD_VISIBLE */
#if __ISO_C_VISIBLE >= 2023
double acospi(double);
float acospif(float);
long double acospil(long double);
double asinpi(double);
float asinpif(float);
long double asinpil(long double);
double atanpi(double);
float atanpif(float);
long double atanpil(long double);
double cospi(double);
float cospif(float);
long double cospil(long double);
lib/msun/src/math_private.h
+28 -10
View File
@@ -482,7 +482,7 @@ do { \
*/
#define _SPLIT(x, xh, xl) \
do { \
typeof(x) __t1; \
typeof(xl) __t1; \
__t1 = (x) * _CC; \
xh = __t1 + ((x) - __t1); \
xl = (x) - xh; \
@@ -505,8 +505,8 @@ do { \
*/
#define _SLOW2SUM(x, y, hi, lo) \
do { \
volatile typeof(x) __t1; \
typeof(x) __t2; \
volatile typeof(lo) __t1; \
typeof(lo) __t2; \
hi = (x) + (y); \
__t1 = hi - (y); \
__t2 = hi - __t1; \
@@ -520,7 +520,7 @@ do { \
*/
#define _XADD(xh, xl, yh, yl, zh, zl) \
do { \
typeof(xh) __s1, __s2, __s3, __s4, __s5, __s6; \
typeof(zl) __s1, __s2, __s3, __s4, __s5, __s6; \
_SLOW2SUM(xh, yh, __s1, __s2); \
_SLOW2SUM(xl, yl, __s3, __s4); \
_FAST2SUM(__s1, __s2 + __s3, __s5, __s6); \
@@ -533,8 +533,7 @@ do { \
*/
#define _MUL(x, y, r1 ,r2) \
do { \
typeof(x) __xh, __xl, __yh, __yl; \
typeof(x) __t1; \
typeof(r2) __xh, __xl, __yh, __yl, __t1; \
_SPLIT(x, __xh, __xl); \
_SPLIT(y, __yh, __yl); \
r1 = (x) * (y); \
@@ -547,12 +546,31 @@ do { \
* and low parts via the _SPLIT macro. x and y are multiplied to give z
* as high and low parts.
*/
#define _XMUL(xh, xl, yh, yl, ph, pl) \
do { \
_MUL(xh, yh, ph, pl); \
pl += xl * yl + xl * yh + xh * yl; \
#define _XMUL(xh, xl, yh, yl, ph, pl) \
do { \
_MUL(xh, yh, ph, pl); \
pl += xl * yl + xl * yh + xh * yl; \
} while(0)
/*
* Use the final Newton-Raphson iteration for a / sqrt(a) from Karp and
* Markstein to get nearly 2*p bits of precision split into zh and zl
* for sqrt(a).
*/
#define _SQRT(x, zh, zl) \
do { \
typeof(zl) s, ph, pl, xh, xl, yh, yl; \
_SPLIT(x, xh, xl); \
s = 1 / _ROOT(x); \
yh = s * xh; \
yl = s * xl; \
_XMUL(yh, yl, yh, yl, ph, pl); \
_XADD(xh, xl, -ph, -pl, zh, zl); \
s /= 2; \
_MUL(s, zh, ph, pl); \
pl += s * zl; \
_XADD(yh, yl, ph, pl, zh, zl); \
} while(0)
/*
* Common routine to process the arguments to nan(), nanf(), and nanl().
lib/msun/src/s_asinpi.c
+243
View File
@@ -0,0 +1,243 @@
/*-
* 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.
*/
/**
* asinpi(x) = asin(x) / pi Eq. (1)
*
* The rational approximation for asinpi(x) has the following form:
*
* x R(x^2)
* asinpi(x) = ---- + x * x^2 * ------------ Eq. (2)
* pi 1 + S(x^2)
*
* with x^2 = x * x. Define r(x^2) = x^2 * [R / (1 + S)], one then has
*
* asinpi(x) = x * [1 / pi + r(x^2)] Eq. (3)
*
* For |x| << 1, asinpi(x) = x / pi. That is, the 2nd term in the righthand
* side of Eq. (2) can be neglected for |x| < 0x1p{-N/2} where {N/2} is half
* the precision (e.g., N = 53, {N/2} = 26). It is noted that for
* |x| < 0x1p{emin+m} with m chosen through testing, x / pi is approaching
* or is subnormal. To compute the result, x is scaled by 0x1p{N+1}, x / pi
* is computed, and finally rescaled by 0x1p{-(N+1)}.
*
* In the domain, 0x1p{-N/2} <= |x| < 0.5, the approximation becomes
*
* asinpi(x) = x * (lo + r(x^2) + hi)
*
* where lo and hi are full-precision low and high parts of 1 / pi.
*
* In the interval [0.5,1), the following relationship
*
* asin(x) = pi / 2 - 2 * asin(t)
*
* with t = [(1 - x) / 2]^{1/2} is used to rewrite Eq. (1). Thus,
*
* asinpi(x) = 1 / 2 - 2 * t * [1 / pi + r(t^2)] Eq. (4)
*
* Note, special cases:
*
* asinpi(+-0) = +-0, exactly.
* asinpi(+-1) = +-1/2, exactly.
* asinpi(x) = nan for |x| > 1
* asinpi(nan) = nan
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
#define _CC (0x1p27 + 1)
#define _ROOT sqrt
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 double
invpihi = 3.1830988618379069e-01,
invpilo = -1.9678676675182486e-17;
/*
* R(x^2)
* __r(x^2) = x^2 * ------------
* 1 + S(x^2)
*
* Prior to the leading multiplication by x^2, the rational approximation
* has an absolute minimax error less than 8.57e-21 over the [0x1p-40,0.5]
* domain (or log2(error) = -66.7).
*/
static inline double
__r(double xs)
{
static const double
R0 = 5.3051647697298449e-02,
R1 = -1.2219903601836109e-01,
R2 = 9.7236612309627199e-02,
R3 = -3.0778625727037261e-02,
R4 = 3.1527637063244254e-03,
R5 = -1.9159514282614908e-05,
S1 = -2.7533975629862262e+00,
S2 = 2.8040387218379421e+00,
S3 = -1.2867553139513013e+00,
S4 = 2.5507476275412666e-01,
S5 = -1.6150977787265989e-02;
double r, s;
r = R0 + (R1 + (R2 + (R3 + (R4 + R5 * xs) * xs) * xs) * xs) * xs;
s = 1 + (S1 + (S2 + (S3 + (S4 + S5 * xs) * xs) * xs) * xs) * xs;
return (xs * (r / s));
}
#include <stdio.h>
double
asinpi(double x)
{
double ax, hi, lo, xh, xl, y, zh, zl;
uint32_t hx, ix, lx;
EXTRACT_WORDS(hx, lx, x);
ix = hx & 0x7fffffff;
if (ix > 0x3ff00000) /* |x| > 1 */
return ((x - x) / (x - x));
INSERT_WORDS(ax, ix, lx);
if (ix <= 0x3fe00000) { /* |x| <= 0.5 */
if (ix < 0x3e400000) { /* |x| < 0x1p-27 */
if (ix < 0x00800000) { /* |x| < 0x1p-1015 */
if ((ix | lx) == 0)
return (x);
/* Scale for near subnormal. */
ax *= 0x1p54;
_XMUL(ax, 0, invpihi, invpilo, hi, lo);
y = (hi + lo) * 0x1p-54;
} 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 < 0x3ff00000) { /* |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 ((hx & 0x80000000) ? -y : y);
}
#if LDBL_MANT_DIG == 53
__weak_reference(asinpi, asinpil);
#endif
/*
* acospi(x) = acos(x) / pi
*
* The implementation uses two identities:
*
* acos(x) = pi / 2 - asin(x) Eq. (5)
* acos(-|x|) = pi - acos(|x|) Eq. (6)
*
* and the definitions for asinpi(x) above. Conversion of Eq. (5)
* with the aid of Eq. (3) leads to the form:
*
* acospi(x) = 1 / 2 - x * [1 / pi + r(x^2)] Eq. (7)
*
* where 0 <= |x| < 0.5. There are two thresholds. For |x| < 0x1p{-N}
* acospi(x) = 1/2 - tiny, which raises FE_INEXACT while preventing spurious
* underflow. For |x| < 0x1p{-M}, acospi(x) = 1/2 - x / pi where M is
* a sloppy threshold determined from testing. For 0.5 <= |x| < 1, there
* are two approximations:
*
* acospi(x) = 2 * t * [1 / pi + r(t^2)] Eq. (8)
*
* for 0.5 <= x < 1. When -1 < x <= -0.5, the relevant expression is
*
* acospi(x) = 1 - 2 * t * [1 / pi + r(t^2)] Eq. (9)
*
* Note, special cases:
*
* acospi(+-0) = 1/2, exactly
* acospi(1) = 0, exactly
* acospi(-1) = 1, exactly
* asinpi(x) = nan for |x| > 1
* asinpi(nan) = nan
*/
double
acospi(double x)
{
double ax, hi, lo, xh, xl, y, zh, zl;
uint32_t hx, ix, lx;
EXTRACT_WORDS(hx, lx, x);
ix = hx & 0x7fffffff;
if (ix > 0x3ff00000) /* |x| > 1 */
return ((x - x) / (x - x));
if (ix <= 0x3fe00000) { /* |x| <= 0.5 */
if (ix < 0x3eb00000) { /* |x| < 0x1p-20 */
y = ((ix | lx) == 0) ? half : ((ix < 0x3ca00000) ?
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 < 0x3ff00000) { /* |x| < 1 */
INSERT_WORDS(ax, ix, lx);
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 (hx & 0x80000000)
_XADD(one, 0, -hi, -lo, y, ax);
else
y = hi + lo;
} else /* |x| == 1 */
y = hx & 0x80000000 ? 1 : 0;
return (y);
}
#if LDBL_MANT_DIG == 53
__weak_reference(acospi, acospil);
#endif
lib/msun/src/s_asinpif.c
+157
View File
@@ -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_asinpi.c for implementation details.
*/
#include "math.h"
#include "math_private.h"
#define _CC (0x1p12F + 1)
#define _ROOT sqrtf
volatile static const float tiny = 1.e-30;
static const float half = 0.5f, one = 1.f;
/* Full precision high and low parts of 1 / pi. */
static const float
invpihi = 3.18309873e-01f,
invpilo = 1.28412765e-08f;
/*
* Prior to the leading multiplication by x^2, the rational approximation
* has an absolute minimax error less than 2.8e-10 over the [0x1p-12,0.5]
* domain (or log2(error) = -31.7).
*/
static inline float
__r(float xs)
{
static const float
R0 = 5.30516468e-02f,
R1 = -3.80413085e-02f,
R2 = 1.74116367e-03f,
S1 = -1.16706085e+00f,
S2 = 2.90115148e-01f;
float r, s;
r = R0 + (R1 + R2 * xs) * xs;
s = 1 + (S1 + S2 * xs) * xs;
return (xs * (r / s));
}
float
asinpif(float x)
{
float ax, hi, lo, xh, xl, y, zh, zl;
uint32_t hx, ix;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
/* |x| > 1 */
if (ix > 0x3f800000)
return ((x - x) / (x - x));
SET_FLOAT_WORD(ax, ix);
if (ix <= 0x3f000000) { /* |x| <= 0.5 */
if (ix < 0x39800000) { /* |x| < 0x1p-12 */
if (ix < 0x03800000) { /* |x| < 0x1p-120 */
if (ix == 0)
return (x);
/* Scale for near subnormal. */
ax *= 0x1p25f;
_XMUL(ax, 0, invpihi, invpilo, hi, lo);
y = (hi + lo) * 0x1p-25f;
} 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 < 0x3f800000) { /* |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 ((hx & 0x80000000) ? -y : y);
}
/*
* See src/s_asinpi.c for implementation details.
*/
float
acospif(float x)
{
float ax, hi, lo, xh, xl, y, zh, zl;
uint32_t hx, ix;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
/* |x| > 1 */
if (ix > 0x3f800000)
return ((x - x) / (x - x));
if (ix <= 0x3f000000) { /* |x| <= 0.5 */
if (ix <= 0x3a800000) { /* |x| <= 0x1p-10 */
y = (ix == 0) ? half : ((ix < 0x33800000) ?
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 < 0x3f800000) { /* |x| < 1 */
SET_FLOAT_WORD(ax, ix);
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 (hx & 0x80000000)
_XADD(one, 0, -hi, -lo, y, ax);
else
y = hi + lo;
} else /* |x| == 1 */
y = hx & 0x80000000 ? 1 : 0;
return (y);
}
lib/msun/src/s_atanpi.c
+215
View File
@@ -0,0 +1,215 @@
/*-
* 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.
*/
/**
* atanpi(x) = atan(x) / pi Eq. (1)
*
* Note, special cases:
*
* atanpi(+-0) = +-0, exactly.
* atanpi(+-inf) = +-1/2, exactly.
* atanpi(nan) = nan
*
* Reflection symmetry atanpi(-|x|) = - atanpi(|x|) allows the
* implementation to be defined for x >= 0.
*
* A rational approximation for atanpi(x) has the following form:
*
* x R(x^2)
* atanpi(x) = ---- + x * x^2 * ------------ Eq. (2)
* pi 1 + S(x^2)
*
* with x^2 = x * x. Define r(x^2) = x^2 * [R / (1 + S)], one then has
*
* atanpi(x) = x * [1 / pi + r(x^2)] Eq. (3)
*
* In addition, for some subdomains of x, the addition formula is used.
*
* atanpi(x) = atanpi(v) + atanpi[(x - v) / (1 + v * x)] Eq. (4)
*
* In the interval [0,0x1p{-N/2}) with N the precision of the floating
* point type, Eq. (2) can be reduced to
*
* atanpi(x) = x / pi. Eq. (5)
*
* However, the division by pi (or more appropriately multiplication]
* by the reciprocal) causes issues with |x| < 0x1p{emin+m} with 'm'
* determined from testing. The result of Eq. (5) approaches or is a
* subnormal. Here, x is scaled by 0x1p{N+1}, Eq. (4) is evaluated, and
* then the result is scaled by 0x1p{-(N+1)}.
*
* In the interval [0xp{-N/2}, 0.5], Eq. (2) is evaluated where the
* rational approximation has be found by a minimax procedure.
*
* In the interval [0.5,0.75), the addition formula gives
*
* atanpi(x) = atanpi(x0) + atanpi[(x - x0)/(1 + x0 * x)] Eq. (6)
*
* with x0 = 5/8 chosen at the center of the interval.
*
* In the interval [0.75,1), the addition formula gives
*
* atanpi(x) = atanpi(x1) + atanpi[(x - x1)/(1 + x1 * x)] Eq. (7)
*
* with x1 = 7/8 chosen at the center of the interval.
*
* In the interval [1,2), the addition formula gives
*
* atanpi(x) = atanpi(x2) + atanpi[(x - x2)/(1 + x2 * x)] Eq. (8)
*
* with x2 = 1.5 chosen at the center of the interval.
*
* Finally, in the interval [2,inf) the identity
*
* atanpi(x) = 1/2 - atanpi(1 / x) Eq. (9)
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
#define _CC (0x1p27 + 1)
#define _ROOT sqrt
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 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) */
/*
* R(x^2)
* __r(x^2) = x^2 * ------------
* 1 + S(x^2)
*
* Prior to the leading multiplication by x^2, the rational approximation
* has an absolute minimax error less than 6.24e-19 over the [0x1p-40,0.5]
* domain (or log2(error) = -63.8).
*/
static inline double
__r(double xs)
{
static const double
R0 = -1.0610329539459690e-01,
R1 = -2.0683077993309035e-01,
R2 = -1.3099673469163398e-01,
R3 = -2.9655125284635996e-02,
R4 = -1.7208096636878276e-03,
S1 = 2.5493341763221311e+00,
S2 = 2.3356442152763948e+00,
S3 = 9.2164104384874268e-01,
S4 = 1.4526328350834117e-01,
S5 = 6.2132943401189099e-03;
double r, s;
r = R0 + (R1 + (R2 + (R3 + R4 * xs) * xs) * xs) * xs;
s = 1 + (S1 + (S2 + (S3 + (S4 + S5 * xs) * xs) * xs) * xs) * xs;
return (xs * (r / s));
}
double
atanpi(double x)
{
double ax, hi, lo, xh, xl, y, zh, zl;
uint32_t hx, ix, lx;
EXTRACT_WORDS(hx, lx, x);
ix = hx & 0x7fffffff;
/* x = +-inf, nan */
if (ix >= 0x7ff00000) {
if (ix > 0x7ff00000)
return (x + x);
return ((hx & 0x80000000) ? -half : half);
}
INSERT_WORDS(ax, ix, lx);
if (ix <= 0x3fe00000) { /* |x| <= 0.5 */
if (ix < 0x3e400000) { /* |x| < 0x1p-27 */
if (ix < 0x00800000) { /* |x| < 0x1p-1015 */
if ((ix | lx) == 0)
return (x);
/* Scale for near subnormal. */
ax *= 0x1p54;
_XMUL(ax, 0, invpihi, invpilo, hi, lo);
y = (hi + lo) * 0x1p-54;
} 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 < 0x3ff00000) { /* |x| < 1 */
if (ix < 0x3fe80000) { /* |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 (ix < 0x40000000) { /* |x| < 2 */
if (ix == 0x3ff00000 && lx == 0)
return ((hx & 0x80000000) ? -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 ((hx & 0x80000000) ? -y : y);
}
#if LDBL_MANT_DIG == 53
__weak_reference(atanpi, atanpil);
#endif
lib/msun/src/s_atanpif.c
+142
View File
@@ -0,0 +1,142 @@
/*-
* 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 implementation details.
*/
#include "math.h"
#include "math_private.h"
#define _CC (0x1p12F + 1)
#define _ROOT sqrtf
volatile static const float tiny = 1.e-30;
static const float half = 0.5f, one = 1.f, qrtr = 0.25f;
static const float x0 = 0.625, x1 = 0.875, x2 = 1.5;
/* Full precision high and low parts. */
static const float
invpihi = 3.18309873e-01f, /* 1 / pi */
invpilo = 1.28412765e-08f, /* 1 / pi */
a0hi = 1.77807689e-01f, /* atanpi(x0) */
a0lo = -4.22372093e-09f, /* atanpi(x0) */
a1hi = 2.28810698e-01f, /* atanpi(x1) */
a1lo = -2.42890708e-09f, /* atanpi(x1) */
a2hi = 3.12832952e-01f, /* atanpi(x2) */
a2lo = 6.64328592e-09f; /* atanpi(x2) */
/*
* Prior to the leading multiplication by x^2, the rational approximation
* has an absolute minimax error less than 1.59e-10 over the [0x1p-12,0.5]
* domain (or log2(error) = -32.5).
*/
static inline float
__r(float xs)
{
static const float
R0 = -1.06103294e-01f,
R1 = -6.81197494e-02f,
R2 = -1.61480496e-03f,
S1 = 1.24201322e+00f,
S2 = 3.31868112e-01f;
float r, s;
r = R0 + (R1 + R2 * xs) * xs;
s = 1 + (S1 + S2 * xs) * xs;
return (xs * (r / s));
}
float
atanpif(float x)
{
float ax, hi, lo, xh, xl, y, zh, zl;
uint32_t hx, ix;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
/* x = +-inf, nan */
if (ix >= 0x7f800000) {
if (ix > 0x7f800000)
return (x + x);
return ((hx & 0x80000000) ? -half : half);
}
SET_FLOAT_WORD(ax, ix);
if (ix <= 0x3f000000) { /* |x| <= 0.5 */
if (ix < 0x39000000) { /* |x| < 0x1p-13 */
if (ix < 0x03800000) { /* |x| < 0x1p-120 */
if (ix == 0)
return (x);
/* Scale for near subnormal. */
ax *= 0x1p25f;
_XMUL(ax, 0, invpihi, invpilo, hi, lo);
y = (hi + lo) * 0x1p-25f;
} 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 < 0x3f800000) { /* |x| < 1 */
if (ix < 0x3f400000) { /* |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| < 1 */
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 < 0x40000000) { /* |x| < 2 */
if (ix == 0x3f800000)
return ((hx & 0x80000000) ? -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 ((hx & 0x80000000) ? -y : y);
}