From 27f74049764154ec2edcba1c5e361e7e99bdcaed Mon Sep 17 00:00:00 2001
From: Taylor R Campbell <riastradh@NetBSD.org>
Date: Mon, 13 May 2024 02:21:39 +0000
Subject: [PATCH] asinf(3): Use a better polynomial approximation.
This nearly-minimax polynomial approximation is generated by Sollya
using lattice reduction to find essentially optimal coefficients
among the discrete set of IEEE 754 binary32 floating-point numbers.
The Sollya script for generating the parameters, and a bound on the
maximum polynomial approximation error, is included.
The floating-point arithmetic that we compute in C is further proven
by Gappa not to introduce too much additional error -- the new code
is analytically proven to have relative forward error bounded by
4.27294e-07.
The bound isn't very impressive -- it's slightly worse than the
actual worst-case error of the existing asinf code, 3.19666e-07.
(The actual worst-case error of the new code is about 2.01648e-07 if
internal arithmetic is in binary32, and much lower, <0.56ulp, if
internal arithmetic is in binary64.)
But this is a proof-of-concept of automatic verification integrated
into libm, where Gappa is verifying the same arithmetic expressions
that the C code uses, through judicious use of cpp. The flow is:
1. asinf.sollya generates a cutoff, coefficients, a polynomial
approximation error bound, and other parameters in asinf_param.h
(and verifies a property of the cutoff)
2. asinf_poly.h has cpp #define macros to form arithmetic expressions
using the parameters from asinf_params.h
3. asinf.gappa and asinfhi.gappa verify bounds on relative error of
the arithmetic in asinf_poly.h from the true function, in terms of
the polynomial approximation error computed by asinf.sollya.
We break the domain into two parts (three, really); asinf.gappa
verifies the approach to the domain [alpha, 1/2], and
asinfhi.gappa verifies the approach to the domain (1/2, 1), where
alpha is the cutoff below which |x| is a good approximation to
asin |x|.
4. e_asinf.c uses exactly the same arithmetic expressions that are
verified by asinf.gappa and asinfhi.gappa, but in C code.
Absolutely no magic inexplicable number constants more exciting
than 1/2 in the hand-written C code.
If you want to change the polynomial, you can edit asinf.sollya and
run `make verify' to regenerate asinf_poly.h and see if the bounds
are still satisfied. Or replace Horner's rule by another polynomial
evaluation method in asinf_poly.h, and run `make verify' to confirm
the change in arithmetic can't raise the relative error above
2^{-21}.
This experiment has some rough edges -- Gappa wasn't really designed
to work with cpp, so the conjuction operator /\ has to be written /\\
instead, and some things don't work like #@cmdlineargs. And it'd be
nice to put a better bound on asinf. But it's a proof of concept.
---
lib/libm/src/Makefile | 40 ++++++++
lib/libm/src/asinf.gappa | 83 ++++++++++++++++
lib/libm/src/asinf.sollya | 168 ++++++++++++++++++++++++++++++++
lib/libm/src/asinf_param.h | 14 +++
lib/libm/src/asinf_poly.h | 118 +++++++++++++++++++++++
lib/libm/src/asinfhi.gappa | 91 ++++++++++++++++++
lib/libm/src/e_asinf.c | 191 +++++++++++++++++++++++--------------
7 files changed, 634 insertions(+), 71 deletions(-)
create mode 100644 lib/libm/src/Makefile
create mode 100644 lib/libm/src/asinf.gappa
create mode 100644 lib/libm/src/asinf.sollya
create mode 100644 lib/libm/src/asinf_param.h
create mode 100644 lib/libm/src/asinf_poly.h
create mode 100644 lib/libm/src/asinfhi.gappa
diff --git a/lib/libm/src/Makefile b/lib/libm/src/Makefile
new file mode 100644
index 000000000000..637e67b7bc85
--- /dev/null
+++ b/lib/libm/src/Makefile
@@ -0,0 +1,40 @@
+# $NetBSD$
+#
+
+default-target: all
+default-target: .PHONY
+.PHONY:
+
+GAPPA= gappa
+SOLLYA= sollya
+
+all: .PHONY
+all: verify
+verify: .PHONY
+proof: .PHONY
+
+proof: asinf.coq
+asinf.coq: asinf.gappa
+asinf.coq: asinf_param.h
+asinf.coq: asinf_poly.h
+ ${CPP} -undef -traditional asinf.gappa \
+ | ${GAPPA} -Bcoq >${.TARGET}.tmp && \
+ mv -f ${.TARGET}.tmp ${.TARGET}
+
+verify: verify-asinf
+verify-asinf: .PHONY
+verify-asinf: asinf.gappa
+verify-asinf: asinf_param.h
+verify-asinf: asinf_poly.h
+ ${CPP} -undef -traditional asinf.gappa | ${GAPPA}
+
+verify: verify-asinfhi
+verify-asinfhi: .PHONY
+verify-asinfhi: asinfhi.gappa
+verify-asinfhi: asinf_param.h
+verify-asinfhi: asinf_poly.h
+ ${CPP} -undef -traditional asinfhi.gappa | ${GAPPA}
+
+asinf_param.h: asinf.sollya
+ ${SOLLYA} asinf.sollya ${.TARGET}.tmp && \
+ mv -f ${.TARGET}.tmp ${.TARGET}
diff --git a/lib/libm/src/asinf.gappa b/lib/libm/src/asinf.gappa
new file mode 100644
index 000000000000..82b063c6dcc5
--- /dev/null
+++ b/lib/libm/src/asinf.gappa
@@ -0,0 +1,83 @@
+/* Gappa 1.3.5 (via cpp) */
+
+/* $NetBSD$ */
+
+/*-
+ * Copyright (c) 2024 The NetBSD Foundation, Inc.
+ * 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, 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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.
+ */
+
+/*
+ * Usage: cpp -undef -traditional asinf.gappa | gappa
+ *
+ * Proves an error bound and prints it.
+ */
+
+#include "asinf_param.h"
+#include "asinf_poly.h"
+
+@RN32 = float<ieee_32,ne>;
+@RN = float<ieee_32,ne>;
+
+/*
+ * Input is always the outcome of rounding to binary32.
+ */
+x = RN32(xx);
+mx = x;
+
+/*
+ * fx := f(x) = asin(x) mathematical function to approximate
+ * x2 := x^2 squared evaluation point for asin
+ * y := x*P(x^2) polynomial approximation to asin(x)
+ * z := RN32(y) result rounded to float
+ */
+
+x2 RN= S(x);
+mx2 = S(mx);
+
+y RN= P(x, x2);
+my = P(mx, mx2);
+
+/*
+ * Ouptut is always rounded to float even if the intermediate
+ * arithmetic is done in double.
+ */
+z = RN32(y);
+mz = my;
+
+e_approx = (mz - fx)/fx; /* polynomial approximation error */
+e_round = (z - mz)/mz; /* floating-point rounding error */
+e = (z - fx)/fx; /* total error */
+
+{
+ x in [ALPHA, 1b-1]
+ /\\
+ fx <> 0 /* asin is nonzero on the domain of x */
+ /\\
+ |e_approx| <= E_APPROX
+->
+ |e| in ?
+ /\\
+ |e| <= 0x1p-22
+}
diff --git a/lib/libm/src/asinf.sollya b/lib/libm/src/asinf.sollya
new file mode 100644
index 000000000000..5b44e180fafc
--- /dev/null
+++ b/lib/libm/src/asinf.sollya
@@ -0,0 +1,168 @@
+/* Sollya 8.0 */
+
+/* $NetBSD$ */
+
+/*-
+ * Copyright (c) 2024 The NetBSD Foundation, Inc.
+ * 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, 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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.
+ */
+
+/*
+ * Usage: sollya asinf.sollya asinf.h
+ */
+
+dieonerrormode = on!;
+verbosity = 1!;
+out = __argv[0];
+write("/*\n") >out;
+write(" * DO NOT EDIT -- automatically generated by Sollya:\n") >>out;
+write(" */\n") >>out;
+write("\n") >>out;
+
+/*
+ * The function we aim to approximate.
+ */
+f = asin(x);
+
+/*
+ * Floating-point precision parameter ulp(1), the distance between 1
+ * and the next larger floating-point number.
+ */
+ulp1 = 2^-23;
+
+/*
+ * Compute a cutoff below which x is a correctly rounded approximation
+ * to asin(x).
+ *
+ * The Taylor expansion of asin x about 0 is:
+ *
+ * x + x^3/6 + x^5/40 + ...
+ *
+ * Hence, (asin x - x)/x is
+ *
+ * x^2/6 + x^4/40 + ...
+ *
+ * If 0 <= x <= 2^{-k} for some k, then the relative error of asin x
+ * from x is
+ *
+ * (asin x - x)/x = x^2/6 + x^4/40 + ...
+ * <= (2^{-k})^2/6 + (2^{-k})^4/40 + ...
+ * <= (2^{-k})^2 + (2^{-k})^4 + ...
+ * = (2^{-2k}) + (2^{-2k})^2 + ...
+ * = 1/(1 - 2^{-2k}) - 1
+ * <= 2^{1 - 2k}.
+ *
+ * In other words,
+ *
+ * asin x = x*(1 + e)
+ * for
+ * e <= 2^{1 - 2k}.
+ *
+ * Hence,
+ *
+ * x = (asin x)/(1 + e) = (asin x)*(1 + d)
+ * for
+ * d = 1/(1 + e) - 1
+ * < 2 e = 2^{2 - 2k}.
+ *
+ * So as long as 2^{2 - 2k} <= ulp(1)/2, x is the correctly rounded
+ * result for any x <= alpha := 2^{-k}.
+ *
+ * To keep it simple, we take 2^{2 - 2k} = ulp(1)/2 and derive:
+ *
+ * 2^{-k} = 2^{1 - k - 1}
+ * = 2^{1 - k}/2
+ * = (2^{2 - 2k})^{1/2}/2
+ * = sqrt(ulp(1)/2)/2.
+ *
+ * A slightly better bound is possible, RN((3/2)^{1/3}) * 2^{-11} (see
+ * Muller, Brunie, Dinechin, Jeannerod, Joldes, Lefèvre, Melquiond,
+ * Revol, and Torres, _Handbook of Floating-Point Arithmetic_, 2nd ed.,
+ * Table 10.3, p. 401) but it's not as obvious how to derive it.
+ */
+alpha = sqrt(ulp1/2)/2;
+if (abs((alpha - asin(alpha))/alpha) > ulp1/2) then {
+ print("alpha too large:", alpha);
+ 1+"error";
+};
+
+/*
+ * Compute the minimax polynomial of the form
+ *
+ * x*(1 + x^2*(P0 + x^2*(P1 + ...)))
+ *
+ * for asin on the interval [alpha, 1/2], and estimate the maximum
+ * relative error of the polynomial on that interval.
+ */
+I = [alpha, 1/2];
+print("Computing minimax polynomial for", f, "on", I);
+P = fpminimax(f, [|1,3,5,7,9,11|], [|1, single...|], I, relative);
+display = decimal!; print("P =", P);
+display = hexadecimal!; print("P =", P);
+
+print("Computing maximum approximation error");
+Pnorm = supnorm(P, f, I, relative, ulp1/2);
+print("Pnorm =", Pnorm);
+display = hexadecimal!; print("e_approx <=", sup(Pnorm));
+
+/*
+ * Additionally compute the error of
+ *
+ * round(pi/2, double, RN) - 2*x*(1 + x^2*(P0 + x^2*(P1 + ...)))
+ *
+ * as an approximation to
+ *
+ * pi/2 - 2*asin(x),
+ *
+ * on the interval [sqrt(ulp(1)/4), 1/2], which we use to compute asin on
+ * [1/2, 1] via trigonometric identities.
+ */
+print("Computing maximum approximation error on the high half");
+pio2 = round(pi/2, double, RN);
+I_hi = [sqrt(ulp1/4), 1/2];
+Pnorm_hi = supnorm(pio2 - 2*P, pi/2 - 2*f, I_hi, relative, ulp1/2);
+print("Pnorm_hi =", Pnorm_hi);
+display = hexadecimal!; print("e_approx_hi <=", sup(Pnorm));
+
+/*
+ * Fill the .h file with constant definitions.
+ */
+define = proc(n, v) {
+ display = hexadecimal!;
+ write("#define "@n@" "@v@"\n") >>out;
+};
+
+display = hexadecimal!;
+define("ULP1", ulp1);
+define("ALPHA", alpha);
+define("PIO2", pio2);
+define("E_APPROX", sup(Pnorm));
+define("E_APPROX_HI", sup(Pnorm_hi));
+
+for i from 1 to (degree(P) - 1)/2 do {
+ display = decimal!;
+ define("P"@(i - 1), coeff(P, 2*i + 1));
+};
+
+quit;
diff --git a/lib/libm/src/asinf_param.h b/lib/libm/src/asinf_param.h
new file mode 100644
index 000000000000..69efb0e0aa5f
--- /dev/null
+++ b/lib/libm/src/asinf_param.h
@@ -0,0 +1,14 @@
+/*
+ * DO NOT EDIT -- automatically generated by Sollya:
+ */
+
+#define ULP1 0x1p-23
+#define ALPHA 0x1p-13
+#define PIO2 0x1.921fb54442d18p0
+#define E_APPROX 0x1.5c3f911fb8534fe7eap-28
+#define E_APPROX_HI 0x1.49e1963dcae856469dp-27
+#define P0 0x1.5555ccp-3
+#define P1 0x1.33010ep-4
+#define P2 0x1.748f98p-5
+#define P3 0x1.8b951ep-6
+#define P4 0x1.59dbe4p-5
diff --git a/lib/libm/src/asinf_poly.h b/lib/libm/src/asinf_poly.h
new file mode 100644
index 000000000000..388c77c4a4c9
--- /dev/null
+++ b/lib/libm/src/asinf_poly.h
@@ -0,0 +1,118 @@
+/* $NetBSD$ */
+
+/*-
+ * Copyright (c) 2024 The NetBSD Foundation, Inc.
+ * 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, 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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.
+ */
+
+/*
+ * Arithmetic for asinf polynomial approximation to arcsin.
+ */
+
+#ifdef HAVE_FAST_FMA
+#define POLY(a, x, y) fmaf((x), (y), (a))
+#else
+#define POLY(a, x, y) ((a) + (x)*(y))
+#endif
+
+/*
+ * Assume alpha < u <= 1/2,
+ *
+ * P(u, u^2) = u*(1 + u*P(u^2))
+ *
+ * is a polynomial approximation to arcsin.
+ *
+ *
+ * Usage:
+ *
+ * u2 = S(u);
+ * return P(u, u2);
+ */
+#define P(u,u2) \
+ ((u)*(1 + (u2)*POLY(P0, u2,POLY(P1, u2,POLY(P2, u2,POLY(P3, u2,P4))))))
+#define S(u) ((u)*(u))
+
+/*
+ * Assume 1/2 < x < 1.
+ *
+ * By the identity
+ *
+ * sin a = cos (pi/2 - a)
+ *
+ * we can rewrite
+ *
+ * asin |x| = pi/2 - acos |x|,
+ *
+ * and solve instead for b such that
+ *
+ * |x| = cos b (0 <= b <= pi/2).
+ *
+ * By the identity
+ *
+ * cos 2t = 1 - 2 sin^2 t,
+ *
+ * we can set
+ *
+ * acos |x| = b = 2t
+ *
+ * and solve for t such that
+ *
+ * |x| = cos b
+ * = cos 2t (0 <= t <= pi/4),
+ * = 1 - 2 sin^2 t,
+ *
+ * i.e., solve for t such that
+ *
+ * sin^2 t = (1 - |x|)/2.
+ *
+ * We can find t given
+ *
+ * u^2 := sin^2 t = (1 - |x|)/2
+ *
+ * by the same minimax polynomial approximation above
+ * on
+ *
+ * asin u ~ u*(1 + P(u^2)),
+ *
+ * noting that since |x| < 1, we have 1 - |x| >= eps/2
+ * and thus u^2 = (1 - |x|)/2 >= eps/4, so
+ *
+ * alpha <= sqrt(eps/4) <= u <= 1/2.
+ *
+ * Combining it all, we get:
+ *
+ * asin |x| = pi/2 - acos |x|
+ * = pi/2 - 2 asin sqrt((1 - |x|)/2).
+ *
+ *
+ * Usage:
+ *
+ * u2 = U(x);
+ * u = R(u2);
+ * t = P(u, u2);
+ * return A(t);
+ */
+#define U(x) ((1 - (x))/2)
+#define R(u2) sqrt(u2)
+#define A(t) (PIO2 - 2*t)
diff --git a/lib/libm/src/asinfhi.gappa b/lib/libm/src/asinfhi.gappa
new file mode 100644
index 000000000000..258cce21aa1c
--- /dev/null
+++ b/lib/libm/src/asinfhi.gappa
@@ -0,0 +1,91 @@
+/* Gappa 1.3.5 (via cpp) */
+
+/* $NetBSD$ */
+
+/*-
+ * Copyright (c) 2024 The NetBSD Foundation, Inc.
+ * 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, 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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.
+ */
+
+/*
+ * Usage: cpp -undef -traditional asinfhi.gappa | gappa
+ *
+ * Proves an error bound and prints it.
+ */
+
+#include "asinf_param.h"
+#include "asinf_poly.h"
+
+@RN32 = float<ieee_32,ne>;
+@RN = float<ieee_32,ne>;
+
+/*
+ * Input is always the outcome of rounding to binary32.
+ */
+x = RN32(xx);
+mx = x;
+
+/*
+ * fx := f(x) = asin(x) = pi/2 - 2 asin(sqrt((1 - x)/2))
+ * mathematical function to approximate
+ * u2 := (1 + x)/2 squared evaluation point for asin
+ * u := sqrt(u2) evaluation point for asin
+ * t := u*P(u^2) approximation to asin(sqrt((1 - x)/2))
+ * z := pio2 - 2 u t approximation to f(x)
+ */
+
+u2 RN= U(x);
+mu2 = U(mx);
+
+u RN= R(u2);
+mu = R(mu2);
+
+t RN= P(u, u2);
+mt = P(mu, mu2);
+
+y RN= A(t);
+my = A(mt);
+
+/*
+ * Ouptut is always rounded to float even if the intermediate
+ * arithmetic is done in double.
+ */
+z = RN32(y);
+mz = my;
+
+e_approx = (mz - fx)/fx; /* polynomial approximation error */
+e_round = (z - mz)/mz; /* floating-point rounding error */
+e = (z - fx)/fx; /* total error */
+
+{
+ x in [1b-1, 0x1.fffffep-1] /* [1/2, 1 - ulp(1)/2] */
+ /\\
+ fx <> 0 /* asin is nonzero on the domain of x */
+ /\\
+ |e_approx| <= E_APPROX_HI
+->
+ |e| in ?
+ /\\
+ |e| <= 0x1p-21
+}
diff --git a/lib/libm/src/e_asinf.c b/lib/libm/src/e_asinf.c
index 713f95250f74..f0442c9d7595 100644
--- a/lib/libm/src/e_asinf.c
+++ b/lib/libm/src/e_asinf.c
@@ -1,88 +1,137 @@
-/* e_asinf.c -- float version of e_asin.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
+/* $NetBSD$ */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+/*-
+ * Copyright (c) 2024 The NetBSD Foundation, Inc.
+ * 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, 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.
*
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
+ * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. 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 FOUNDATION 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.
*/
#include <sys/cdefs.h>
-#if defined(LIBM_SCCS) && !defined(lint)
-__RCSID("$NetBSD: e_asinf.c,v 1.9 2017/05/09 02:04:38 maya Exp $");
-#endif
+__RCSID("$NetBSD$");
+
+#include "namespace.h"
-#include "math.h"
+#include <assert.h>
+#include <float.h>
+#include <math.h>
+
+#include "asinf_param.h"
+#include "asinf_poly.h"
#include "math_private.h"
-static const float
-one = 1.0000000000e+00, /* 0x3F800000 */
-huge = 1.000e+30,
-pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */
-pio2_lo = 7.5497894159e-08, /* 0x33a22168 */
-pio4_hi = 7.8539818525e-01, /* 0x3f490fdb */
- /* coefficient for R(x^2) */
-pS0 = 1.6666667163e-01, /* 0x3e2aaaab */
-pS1 = -3.2556581497e-01, /* 0xbea6b090 */
-pS2 = 2.0121252537e-01, /* 0x3e4e0aa8 */
-pS3 = -4.0055535734e-02, /* 0xbd241146 */
-pS4 = 7.9153501429e-04, /* 0x3a4f7f04 */
-pS5 = 3.4793309169e-05, /* 0x3811ef08 */
-qS1 = -2.4033949375e+00, /* 0xc019d139 */
-qS2 = 2.0209457874e+00, /* 0x4001572d */
-qS3 = -6.8828397989e-01, /* 0xbf303361 */
-qS4 = 7.7038154006e-02; /* 0x3d9dc62e */
+/*
+ * The polynomial and error bounds were derived under the assumption
+ * that ULP1 is actually ulp(1), called FLT_EPSILON in C. If ulp(1) is
+ * larger or smaller, the error analysis is wrong because rounding
+ * might incur more error or because inputs might lie within the
+ * interval (ULP1, 1).
+ */
+#if 0
+static_assert(FLT_EPSILON == ULP1,
+ "FLT_EPSILON must match ulp(1) used by Sollya");
+#endif
+/*
+ * The cutoff under which we return x as an approximation to asin(x)
+ * requires alpha <= sqrt(ulp(1)/2)/2, or alpha^2 <= ulp(1)/8.
+ */
+#if 0
+static_assert(ALPHA*ALPHA <= FLT_EPSILON/8,
+ "alpha must be below sqrt(ulp(1)/2)/2");
+#endif
+
+/*
+ * __ieee754_asinf(x)
+ *
+ * Find an approximation to t such that x = sin(t).
+ *
+ * XXX Explain why this is __ieee754_asinf and not just asinf.
+ */
float
__ieee754_asinf(float x)
{
- float t,w,p,q,c,r,s;
- int32_t hx,ix;
+ float xm = fabsf(x);
+ float ym;
- t = 0;
- GET_FLOAT_WORD(hx,x);
- ix = hx&0x7fffffff;
- if(ix==0x3f800000) {
- /* asin(1)=+-pi/2 with inexact */
- return x*pio2_hi+x*pio2_lo;
- } else if(ix> 0x3f800000) { /* |x|>= 1 */
- return (x-x)/(x-x); /* asin(|x|>1) is NaN */
- } else if (ix<0x3f000000) { /* |x|<0.5 */
- if(ix<0x32000000) { /* if |x| < 2**-27 */
- if(huge+x>one) return x;/* return x with inexact if x!=0*/
- } else {
- t = x*x;
- p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
- q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
- w = p/q;
- return x+x*w;
- }
+ /*
+ * The non-exceptional cases are |x| < 1. Check for the
+ * exceptional cases, either |x| >= 1 or x is NaN, by writing
+ * the negation !(|x| < 1) of the non-exceptional cases.
+ */
+ if (!(xm < 1)) {
+ if (xm == 1) /* x = +/-1; raise inexact, return +/-pi/2 */
+ return x*M_PI/2;
+ else /* |x| > 1; raise invalid, return NaN */
+ return (x - x)/(x - x);
}
- /* 1> |x|>= 0.5 */
- w = one-fabsf(x);
- t = w*(float)0.5;
- p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
- q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
- s = __ieee754_sqrtf(t);
- if(ix>=0x3F79999A) { /* if |x| > 0.975 */
- w = p/q;
- t = pio2_hi-((float)2.0*(s+s*w)-pio2_lo);
+
+ /*
+ * Split the non-exceptional domain into three subdomains:
+ *
+ * I: |x| <= alpha
+ * II: alpha < |x| <= 1/2
+ * III: 1/2 < |x| < 1.
+ */
+ if (xm <= ALPHA) {
+ /*
+ * alpha is the cutoff below which the relative error
+ * of |x| from asin |x| is below ulp(1)/2.
+ *
+ * Raise inexact and return x -- unconditionally; the
+ * conditional branch, which always goes the same way,
+ * is a trick to raise inexact.
+ */
+ return (__predict_true(1/FLT_EPSILON + x > 1) ? x : 0);
+ } else if (xm <= 1/2.) {
+ /*
+ * alpha < |x| <= 1/2.
+ *
+ * Compute asin(u) where u = |x| directly by minimax
+ * polynomial approximation
+ *
+ * asin u ~ u*P(u^2)
+ *
+ * on [alpha, 1/2].
+ */
+ double u = xm;
+ double u2 = S(u);
+ ym = P(u, u2);
} else {
- int32_t iw;
- w = s;
- GET_FLOAT_WORD(iw,w);
- SET_FLOAT_WORD(w,iw&0xfffff000);
- c = (t-w*w)/(s+w);
- r = p/q;
- p = (float)2.0*s*r-(pio2_lo-(float)2.0*c);
- q = pio4_hi-(float)2.0*w;
- t = pio4_hi-(p-q);
+ /*
+ * 1/2 < |x| <= 1 - ulp(1)/2 < 1.
+ *
+ * Compute pi/2 - 2 asin(u) by minimax polynomial
+ * approximation, where u^2 = (1 + |x|)/2.
+ */
+ double u2 = U(xm);
+ double u = R(u2);
+ double t = P(u, u2);
+ ym = A(t);
}
- if(hx>0) return t; else return -t;
+
+ /*
+ * Given ym = |y| = asin |x|, apply the sign of x to determine
+ * y = asin x, since asin is an odd function.
+ */
+ return copysignf(ym, x);
}