10 files changed, 324 insertions, 41 deletions
diff --git a/pl/math/tools/asin.sollya b/pl/math/tools/asin.sollya
new file mode 100644
index 000000000000..8ef861d0898b
--- /dev/null
+++ b/pl/math/tools/asin.sollya
@@ -0,0 +1,29 @@
+// polynomial for approximating asin(x)
+//
+// Copyright (c) 2023, Arm Limited.
+// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+
+f = asin(x);
+dtype = double;
+
+prec=256;
+
+a = 0x1p-106;
+b = 0.25;
+
+deg = 11;
+
+backward = proc(poly, d) {
+  return d + d ^ 3 * poly(d * d);
+};
+
+forward = proc(f, d) {
+  return (f(sqrt(d))-sqrt(d))/(d*sqrt(d));
+};
+
+poly = fpminimax(forward(f, x), [|0,...,deg|], [|dtype ...|], [a;b], relative, floating);
+
+display = hexadecimal!;
+print("rel error:", dirtyinfnorm(1-backward(poly, x)/f(x), [a;b]));
+print("in [", a, b, "]");
+for i from 0 to deg do print(coeff(poly, i));
diff --git a/pl/math/tools/asinf.sollya b/pl/math/tools/asinf.sollya
new file mode 100644
index 000000000000..5b627e546c73
--- /dev/null
+++ b/pl/math/tools/asinf.sollya
@@ -0,0 +1,36 @@
+// polynomial for approximating asinf(x)
+//
+// Copyright (c) 2023, Arm Limited.
+// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+
+f = asin(x);
+dtype = single;
+
+a = 0x1p-24;
+b = 0.25;
+
+deg = 4;
+
+backward = proc(poly, d) {
+  return d + d ^ 3 * poly(d * d);
+};
+
+forward = proc(f, d) {
+  return (f(sqrt(d))-sqrt(d))/(d*sqrt(d));
+};
+
+approx = proc(poly, d) {
+  return remez(1 - poly(x) / forward(f, x), deg - d, [a;b], x^d/forward(f, x), 1e-16);
+};
+
+poly = 0;
+for i from 0 to deg do {
+  i;
+  p = roundcoefficients(approx(poly,i), [|dtype ...|]);
+  poly = poly + x^i*coeff(p,0);
+};
+
+display = hexadecimal!;
+print("rel error:", accurateinfnorm(1-backward(poly, x)/f(x), [a;b], 30));
+print("in [", a, b, "]");
+for i from 0 to deg do print(coeff(poly, i));
diff --git a/pl/math/tools/erf.sollya b/pl/math/tools/erf.sollya
new file mode 100644
index 000000000000..b2fc559b511e
--- /dev/null
+++ b/pl/math/tools/erf.sollya
@@ -0,0 +1,25 @@
+// tables and constants for approximating erf(x).
+//
+// Copyright (c) 2023, Arm Limited.
+// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+
+display = hexadecimal;
+prec=128;
+
+// Tables
+print("{ i, r, erf(r), 2/sqrt(pi) * exp(-r^2)}");
+for i from 0 to 768 do {
+  r = i / 128;
+  t0 = double(erf(r));
+  t1 = double(2/sqrt(pi) * exp(-r * r));
+  print("{ " @ i @ ",\t" @ r @ ",\t" @ t0 @ ",\t" @ t1 @ " },");
+};
+
+// Constants
+double(1/3);
+double(1/10);
+double(2/15);
+double(2/9);
+double(2/45);
+double(2/sqrt(pi));
+
diff --git a/pl/math/tools/erfc.sollya b/pl/math/tools/erfc.sollya
index 8c40b4b5db6b..1e2791291ebb 100644
--- a/pl/math/tools/erfc.sollya
+++ b/pl/math/tools/erfc.sollya
@@ -1,23 +1,51 @@
-// polynomial for approximating erfc(x)*exp(x*x)
+// tables and constants for approximating erfc(x).
 //
-// Copyright (c) 2022-2023, Arm Limited.
+// Copyright (c) 2023, Arm Limited.
 // SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
 
-deg = 12; // poly degree
-
-// interval bounds
-a = 0x1.60dfc14636e2ap0;
-b = 0x1.d413cccfe779ap0;
+display = hexadecimal;
+prec=128;
 
-f = proc(y) {
-  t = y + a;
-  return erfc(t) * exp(t*t);
+// Tables
+print("{ i, r, erfc(r), 2/sqrt(pi) * exp(-r^2) }");
+for i from 0 to 3787 do {
+  r = 0.0 + i / 128;
+  t0 = double(erfc(r) * 2^128);
+  t1 = double(2/sqrt(pi) * exp(-r * r) * 2^128);
+  print("{ " @ t0 @ ",\t" @ t1 @ " },");
 };
 
-poly = remez(f(x), deg, [0;b-a], 1, 1e-16);
+// Constants
+print("> 2/sqrt(pi)");
+double(2/sqrt(pi));
 
-display = hexadecimal;
-print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30));
-print("in [",a,b,"]");
-print("coeffs:");
-for i from 0 to deg do round(coeff(poly,i), 52, RN);
+print("> 1/3");
+double(1/3);
+
+print("> P5");
+double(2/15);
+double(1/10);
+double(2/9);
+double(2/45);
+
+print("> P6");
+double(1/42);
+double(1/7);
+double(2/21);
+double(4/315);
+
+print("> Q");
+double( 5.0 / 4.0);
+double( 6.0 / 5.0);
+double( 7.0 / 6.0);
+double( 8.0 / 7.0);
+double( 9.0 / 8.0);
+double(10.0 / 9.0);
+
+print("> R");
+double(-2.0 * 4.0 / (5.0 * 6.0));
+double(-2.0 * 5.0 / (6.0 * 7.0));
+double(-2.0 * 6.0 / (7.0 * 8.0));
+double(-2.0 * 7.0 / (8.0 * 9.0));
+double(-2.0 * 8.0 / (9.0 * 10.0));
+double(-2.0 * 9.0 / (10.0 * 11.0));
diff --git a/pl/math/tools/erfcf.sollya b/pl/math/tools/erfcf.sollya
index 69c683647af7..1d7fc264d99d 100644
--- a/pl/math/tools/erfcf.sollya
+++ b/pl/math/tools/erfcf.sollya
@@ -1,31 +1,22 @@
-// polynomial for approximating erfc(x)*exp(x*x)
+// tables and constants for approximating erfcf(x).
 //
-// Copyright (c) 2022-2023, Arm Limited.
+// Copyright (c) 2023, Arm Limited.
 // SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
 
-deg = 15; // poly degree
-
-// interval bounds
-a = 0x1.0p-26;
-b = 2;
-
-f = proc(y) {
-  return erfc(y) * exp(y*y);
-};
-
-approx = proc(poly, d) {
-  return remez(1 - poly(x)/f(x), deg-d, [a;b], x^d/f(x), 1e-10);
-};
+display = hexadecimal;
+prec=128;
 
-poly = 0;
-for i from 0 to deg do {
-  p = roundcoefficients(approx(poly,i), [|D ...|]);
-  poly = poly + x^i*coeff(p,0);
-  print(i);
+// Tables
+print("{ i, r, erfc(r), 2/sqrt(pi) * exp(-r^2) }");
+for i from 0 to 644 do {
+  r = 0.0 + i / 64;
+  t0 = single(erfc(r) * 2^47);
+  t1 = single(2/sqrt(pi) * exp(-r * r) * 2^47);
+  print("{ " @ t0 @ ",\t" @ t1 @ " },");
 };
 
-display = hexadecimal;
-print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30));
-print("in [",a,b,"]");
-print("coeffs:");
-for i from 0 to deg do coeff(poly,i);
+// Constants
+single(1/3);
+single(2/15);
+single(1/10);
+single(2/sqrt(pi));
diff --git a/pl/math/tools/erff.sollya b/pl/math/tools/erff.sollya
new file mode 100644
index 000000000000..59b23ef021f0
--- /dev/null
+++ b/pl/math/tools/erff.sollya
@@ -0,0 +1,20 @@
+// tables and constants for approximating erff(x).
+//
+// Copyright (c) 2023, Arm Limited.
+// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+
+display = hexadecimal;
+prec=128;
+
+// Tables
+print("{ i, r, erf(r), 2/sqrt(pi) * exp(-r^2)}");
+for i from 0 to 512 do {
+  r = i / 128;
+  t0 = single(erf(r));
+  t1 = single(2/sqrt(pi) * exp(-r * r));
+  print("{ " @ i @ ",\t" @ r @ ",\t" @ t0 @ ",\t" @ t1 @ " },");
+};
+
+// Constants
+single(1/3);
+single(2/sqrt(pi));
diff --git a/pl/math/tools/exp10.sollya b/pl/math/tools/exp10.sollya
new file mode 100644
index 000000000000..9f30b4018209
--- /dev/null
+++ b/pl/math/tools/exp10.sollya
@@ -0,0 +1,55 @@
+// polynomial for approximating 10^x
+//
+// Copyright (c) 2023, Arm Limited.
+// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+
+// exp10f parameters
+deg = 5; // poly degree
+N = 1; // Neon 1, SVE 64
+b = log(2)/(2 * N * log(10)); // interval
+a = -b;
+wp = single;
+
+// exp10 parameters
+//deg = 4; // poly degree - bump to 5 for ~1 ULP
+//N = 128; // table size
+//b = log(2)/(2 * N * log(10)); // interval
+//a = -b;
+//wp = D;
+
+
+// find polynomial with minimal relative error
+
+f = 10^x;
+
+// return p that minimizes |f(x) - poly(x) - x^d*p(x)|/|f(x)|
+approx = proc(poly,d) {
+  return remez(1 - poly(x)/f(x), deg-d, [a;b], x^d/f(x), 1e-10);
+};
+// return p that minimizes |f(x) - poly(x) - x^d*p(x)|
+approx_abs = proc(poly,d) {
+  return remez(f(x) - poly(x), deg-d, [a;b], x^d, 1e-10);
+};
+
+// first coeff is fixed, iteratively find optimal double prec coeffs
+poly = 1;
+for i from 1 to deg do {
+  p = roundcoefficients(approx(poly,i), [|wp ...|]);
+//  p = roundcoefficients(approx_abs(poly,i), [|wp ...|]);
+  poly = poly + x^i*coeff(p,0);
+};
+
+display = hexadecimal;
+print("rel error:", accurateinfnorm(1-poly(x)/10^x, [a;b], 30));
+print("abs error:", accurateinfnorm(10^x-poly(x), [a;b], 30));
+print("in [",a,b,"]");
+print("coeffs:");
+for i from 0 to deg do coeff(poly,i);
+
+log10_2 = round(N * log(10) / log(2), wp, RN);
+log2_10 = log(2) / (N * log(10));
+log2_10_hi = round(log2_10, wp, RN);
+log2_10_lo = round(log2_10 - log2_10_hi, wp, RN);
+print(log10_2);
+print(log2_10_hi);
+print(log2_10_lo);
diff --git a/pl/math/tools/sincos.sollya b/pl/math/tools/sincos.sollya
new file mode 100644
index 000000000000..7d36266b446b
--- /dev/null
+++ b/pl/math/tools/sincos.sollya
@@ -0,0 +1,33 @@
+// polynomial for approximating cos(x)
+//
+// Copyright (c) 2023, Arm Limited.
+// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+
+// This script only finds the coeffs for cos - see math/aarch64/v_sin.c for sin coeffs
+
+deg = 14;   // polynomial degree
+a = -pi/4; // interval
+b = pi/4;
+
+// find even polynomial with minimal abs error compared to cos(x)
+
+f = cos(x);
+
+// return p that minimizes |f(x) - poly(x) - x^d*p(x)|
+approx = proc(poly,d) {
+  return remez(f(x)-poly(x), deg-d, [a;b], x^d, 1e-10);
+};
+
+// first coeff is fixed, iteratively find optimal double prec coeffs
+poly = 1;
+for i from 1 to deg/2 do {
+  p = roundcoefficients(approx(poly,2*i), [|double ...|]);
+  poly = poly + x^(2*i)*coeff(p,0);
+};
+
+display = hexadecimal;
+//print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30));
+//print("abs error:", accurateinfnorm(f(x)-poly(x), [a;b], 30));
+print("in [",a,b,"]");
+print("coeffs:");
+for i from 0 to deg do coeff(poly,i);
diff --git a/pl/math/tools/sincosf.sollya b/pl/math/tools/sincosf.sollya
new file mode 100644
index 000000000000..178ee83ac196
--- /dev/null
+++ b/pl/math/tools/sincosf.sollya
@@ -0,0 +1,33 @@
+// polynomial for approximating cos(x)
+//
+// Copyright (c) 2023, Arm Limited.
+// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+
+// This script only finds the coeffs for cos - see math/tools/sin.sollya for sin coeffs.
+
+deg = 8;   // polynomial degree
+a = -pi/4; // interval
+b = pi/4;
+
+// find even polynomial with minimal abs error compared to cos(x)
+
+f = cos(x);
+
+// return p that minimizes |f(x) - poly(x) - x^d*p(x)|
+approx = proc(poly,d) {
+  return remez(f(x)-poly(x), deg-d, [a;b], x^d, 1e-10);
+};
+
+// first coeff is fixed, iteratively find optimal double prec coeffs
+poly = 1;
+for i from 1 to deg/2 do {
+  p = roundcoefficients(approx(poly,2*i), [|single ...|]);
+  poly = poly + x^(2*i)*coeff(p,0);
+};
+
+display = hexadecimal;
+//print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30));
+//print("abs error:", accurateinfnorm(f(x)-poly(x), [a;b], 30));
+print("in [",a,b,"]");
+print("coeffs:");
+for i from 0 to deg do coeff(poly,i);
diff --git a/pl/math/tools/sinpi.sollya b/pl/math/tools/sinpi.sollya
new file mode 100644
index 000000000000..62cc87e7697d
--- /dev/null
+++ b/pl/math/tools/sinpi.sollya
@@ -0,0 +1,33 @@
+// polynomial for approximating sinpi(x)
+//
+// Copyright (c) 2023, Arm Limited.
+// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+
+deg = 19;  // polynomial degree
+a = -1/2; // interval
+b = 1/2;
+
+// find even polynomial with minimal abs error compared to sinpi(x)
+
+// f = sin(pi* x);
+f = pi*x;
+c = 1;
+for i from 1 to 80 do { c = 2*i*(2*i + 1)*c; f = f + (-1)^i*(pi*x)^(2*i+1)/c; };
+
+// return p that minimizes |f(x) - poly(x) - x^d*p(x)|
+approx = proc(poly,d) {
+  return remez(f(x)-poly(x), deg-d, [a;b], x^d, 1e-10);
+};
+
+// first coeff is predefine, iteratively find optimal double prec coeffs
+poly = pi*x;
+for i from 0 to (deg-1)/2 do {
+  p = roundcoefficients(approx(poly,2*i+1), [|D ...|]);
+  poly = poly + x^(2*i+1)*coeff(p,0);
+};
+
+display = hexadecimal;
+print("abs error:", accurateinfnorm(sin(pi*x)-poly(x), [a;b], 30));
+print("in [",a,b,"]");
+print("coeffs:");
+for i from 0 to deg do coeff(poly,i);