pl/math: Add scalar asinh

The new routine uses a similar approach to asinhf, using a polynomial
only in the region where either returning x or calculating the result
directly is not sufficiently precise.

Worst-case error is about 2 ULP, close to 1. There are 4
intervals with slightly different error behaviour, as follows:

Interval              Worst-case accuracy (ulp)
|x| < 2^-26           0.0
|x| < 1               1.5
|x| < ~sqrt(DBL_MAX)  2.0
|x| < infinity        1.0

log has been copied from the main math directory so that it can be
used in asinh. The only modifications to the relevant files are to
remove aliases and rename log itself to an internal 'helper' name.

1 files changed, 28 insertions, 0 deletions

diff --git a/pl/math/tools/asinh.sollya b/pl/math/tools/asinh.sollya
new file mode 100644
index 0000000..6ff217f
--- /dev/null
+++ b/pl/math/tools/asinh.sollya
@@ -0,0 +1,28 @@
+// polynomial for approximating asinh(x)
+//
+// Copyright (c) 2022, Arm Limited.
+// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+
+// Polynomial is used in [2^-26, 1]. However it is least accurate close to 1, so
+// we use 2^-6 as the lower bound for coeff generation, which yields sufficiently
+// accurate results in [2^-26, 2^-6].
+a = 0x1p-6;
+b = 1.0;
+
+f = (asinh(sqrt(x)) - sqrt(x))/x^(3/2);
+
+approx = proc(poly, d) {
+  return remez(1 - poly(x)/f(x), deg-d, [a;b], x^d/f(x), 1e-10);
+};
+
+poly = 0;
+for i from 0 to deg do {
+  i;
+  p = roundcoefficients(approx(poly,i), [|D ...|]);
+  poly = poly + x^i*coeff(p,0);
+};
+
+
+display = hexadecimal;
+print("coeffs:");
+for i from 0 to deg do coeff(poly,i);