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

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)
+.\" 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 .
diff --git a/lib/msun/man/asinpi.3 b/lib/msun/man/asinpi.3
new file mode 100644
index 000000000000..465a924b3428
--- /dev/null
+++ b/lib/msun/man/asinpi.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)
+.\" 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 .
diff --git a/lib/msun/man/atanpi.3 b/lib/msun/man/atanpi.3
new file mode 100644
index 000000000000..740f1da78841
--- /dev/null
+++ b/lib/msun/man/atanpi.3
@@ -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 .
diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index e5c2ccb786a4..3ef629cda6e2 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -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);
diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h
index e9055a9a8c61..0134a044dae4 100644
--- a/lib/msun/src/math_private.h
+++ b/lib/msun/src/math_private.h
@@ -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().
diff --git a/lib/msun/src/s_asinpi.c b/lib/msun/src/s_asinpi.c
new file mode 100644
index 000000000000..b51704f36521
--- /dev/null
+++ b/lib/msun/src/s_asinpi.c
@@ -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
diff --git a/lib/msun/src/s_asinpif.c b/lib/msun/src/s_asinpif.c
new file mode 100644
index 000000000000..ac414400564e
--- /dev/null
+++ b/lib/msun/src/s_asinpif.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_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);
+}
diff --git a/lib/msun/src/s_atanpi.c b/lib/msun/src/s_atanpi.c
new file mode 100644
index 000000000000..088d8d54dc9a
--- /dev/null
+++ b/lib/msun/src/s_atanpi.c
@@ -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
diff --git a/lib/msun/src/s_atanpif.c b/lib/msun/src/s_atanpif.c
new file mode 100644
index 000000000000..81791f78c38e
--- /dev/null
+++ b/lib/msun/src/s_atanpif.c
@@ -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);
+}