diff options
author | Elliott Hughes <enh@google.com> | 2023-01-25 16:31:42 +0000 |
---|---|---|
committer | Elliott Hughes <enh@google.com> | 2023-01-25 18:11:56 +0000 |
commit | 62662f115a17a5548402fcd02f3b7e1b9b38e087 (patch) | |
tree | 945eb1832687ad5eb97c433013617dbe107718e5 /pl/math/tools | |
parent | 04d0bc466ce39c702459d77a15754a19533eaec9 (diff) | |
parent | 56e3bf05c19c4e28e1f5edd9093c712f16c5c32a (diff) | |
download | arm-optimized-routines-62662f115a17a5548402fcd02f3b7e1b9b38e087.tar.gz |
Upgrade ARM-software/optimized-routines to v23.01
This project was upgraded with external_updater.
Usage: tools/external_updater/updater.sh update arm-optimized-routines
For more info, check https://cs.android.com/android/platform/superproject/+/master:tools/external_updater/README.md
Test: TreeHugger
Change-Id: I03de0ecb0ea89c2a77b1c411685af09fd480b63e
Diffstat (limited to 'pl/math/tools')
-rw-r--r-- | pl/math/tools/asinh.sollya | 28 | ||||
-rw-r--r-- | pl/math/tools/asinhf.sollya | 29 | ||||
-rw-r--r-- | pl/math/tools/atan.sollya | 23 | ||||
-rw-r--r-- | pl/math/tools/atanf.sollya | 20 | ||||
-rw-r--r-- | pl/math/tools/cbrt.sollya | 20 | ||||
-rw-r--r-- | pl/math/tools/cbrtf.sollya | 20 | ||||
-rw-r--r-- | pl/math/tools/erfc.sollya | 23 | ||||
-rw-r--r-- | pl/math/tools/erfcf.sollya | 31 | ||||
-rw-r--r-- | pl/math/tools/expm1.sollya | 21 | ||||
-rw-r--r-- | pl/math/tools/expm1f.sollya | 21 | ||||
-rw-r--r-- | pl/math/tools/log10.sollya | 44 | ||||
-rw-r--r-- | pl/math/tools/log10f.sollya | 37 | ||||
-rw-r--r-- | pl/math/tools/log1p.sollya | 30 | ||||
-rw-r--r-- | pl/math/tools/log1pf.sollya | 21 | ||||
-rw-r--r-- | pl/math/tools/tan.sollya | 20 | ||||
-rw-r--r-- | pl/math/tools/tanf.sollya | 78 | ||||
-rw-r--r-- | pl/math/tools/v_erf.sollya | 20 | ||||
-rw-r--r-- | pl/math/tools/v_erfc.sollya | 46 | ||||
-rw-r--r-- | pl/math/tools/v_log10.sollya | 38 | ||||
-rw-r--r-- | pl/math/tools/v_log10f.sollya | 45 | ||||
-rw-r--r-- | pl/math/tools/v_log2f.sollya | 38 |
21 files changed, 653 insertions, 0 deletions
diff --git a/pl/math/tools/asinh.sollya b/pl/math/tools/asinh.sollya new file mode 100644 index 0000000..663ee92 --- /dev/null +++ b/pl/math/tools/asinh.sollya @@ -0,0 +1,28 @@ +// polynomial for approximating asinh(x) +// +// Copyright (c) 2022-2023, 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); diff --git a/pl/math/tools/asinhf.sollya b/pl/math/tools/asinhf.sollya new file mode 100644 index 0000000..ab115b5 --- /dev/null +++ b/pl/math/tools/asinhf.sollya @@ -0,0 +1,29 @@ +// polynomial for approximating asinh(x) +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 9; + +a = 0x1.0p-12; +b = 1.0; + +f = proc(y) { + return asinh(x); +}; + +approx = proc(poly, d) { + return remez(1 - poly(x)/f(x), deg-d, [a;b], x^d/f(x), 1e-10); +}; + +poly = x; +for i from 2 to deg do { + p = roundcoefficients(approx(poly,i), [|SG ...|]); + poly = poly + x^i*coeff(p,0); +}; + +display = hexadecimal; +print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30)); +print("in [",a,b,"]"); +print("coeffs:"); +for i from 2 to deg do coeff(poly,i); diff --git a/pl/math/tools/atan.sollya b/pl/math/tools/atan.sollya new file mode 100644 index 0000000..ad4f33b --- /dev/null +++ b/pl/math/tools/atan.sollya @@ -0,0 +1,23 @@ +// polynomial for approximating atan(x) and atan2(y, x) +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +// atan is odd, so approximate with an odd polynomial: +// x + ax^3 + bx^5 + cx^7 + ... +// We generate a, b, c, ... such that we can approximate atan(x) by: +// x + x^3 * (a + bx^2 + cx^4 + ...) + +// Assemble monomials +deg = 20; +mons = [|1,...,deg|]; +for i from 0 to deg-1 do mons[i] = mons[i] * 2 + 1; + +a = 0x1.0p-1022; +b = 1; + +poly = fpminimax(atan(x)-x, mons, [|double ...|], [a;b]); + +display = hexadecimal; +print("coeffs:"); +for i from 0 to deg-1 do coeff(poly,mons[i]); diff --git a/pl/math/tools/atanf.sollya b/pl/math/tools/atanf.sollya new file mode 100644 index 0000000..ed88d0b --- /dev/null +++ b/pl/math/tools/atanf.sollya @@ -0,0 +1,20 @@ +// polynomial for approximating atanf(x) +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +// Generate list of monomials: +// Taylor series of atan is of the form x + ax^3 + bx^5 + cx^7 + ... +// So generate a, b, c, ... such that we can approximate atan(x) by: +// x + x^3 * (a + bx^2 + cx^4 + ...) + +deg = 7; + +a = 1.1754943508222875e-38; +b = 1; + +poly = fpminimax((atan(sqrt(x))-sqrt(x))/x^(3/2), deg, [|single ...|], [a;b]); + +display = hexadecimal; +print("coeffs:"); +for i from 0 to deg do coeff(poly,i); diff --git a/pl/math/tools/cbrt.sollya b/pl/math/tools/cbrt.sollya new file mode 100644 index 0000000..1d43dc7 --- /dev/null +++ b/pl/math/tools/cbrt.sollya @@ -0,0 +1,20 @@ +// polynomial for approximating cbrt(x) in double precision +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 3; + +a = 0.5; +b = 1; + + +f = x^(1/3); + +poly = fpminimax(f, deg, [|double ...|], [a;b]); + +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), D, RN); diff --git a/pl/math/tools/cbrtf.sollya b/pl/math/tools/cbrtf.sollya new file mode 100644 index 0000000..4e0cc69 --- /dev/null +++ b/pl/math/tools/cbrtf.sollya @@ -0,0 +1,20 @@ +// polynomial for approximating cbrt(x) in single precision +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 3; + +a = 0.5; +b = 1; + + +f = x^(1/3); + +poly = fpminimax(f, deg, [|single ...|], [a;b]); + +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), SG, RN); diff --git a/pl/math/tools/erfc.sollya b/pl/math/tools/erfc.sollya new file mode 100644 index 0000000..8c40b4b --- /dev/null +++ b/pl/math/tools/erfc.sollya @@ -0,0 +1,23 @@ +// polynomial for approximating erfc(x)*exp(x*x) +// +// Copyright (c) 2022-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; + +f = proc(y) { + t = y + a; + return erfc(t) * exp(t*t); +}; + +poly = remez(f(x), deg, [0;b-a], 1, 1e-16); + +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); diff --git a/pl/math/tools/erfcf.sollya b/pl/math/tools/erfcf.sollya new file mode 100644 index 0000000..69c6836 --- /dev/null +++ b/pl/math/tools/erfcf.sollya @@ -0,0 +1,31 @@ +// polynomial for approximating erfc(x)*exp(x*x) +// +// Copyright (c) 2022-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); +}; + +poly = 0; +for i from 0 to deg do { + p = roundcoefficients(approx(poly,i), [|D ...|]); + poly = poly + x^i*coeff(p,0); + print(i); +}; + +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); diff --git a/pl/math/tools/expm1.sollya b/pl/math/tools/expm1.sollya new file mode 100644 index 0000000..7b6f324 --- /dev/null +++ b/pl/math/tools/expm1.sollya @@ -0,0 +1,21 @@ +// polynomial for approximating exp(x)-1 in double precision +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 12; + +a = -log(2)/2; +b = log(2)/2; + +f = proc(y) { + return exp(y)-1; +}; + +poly = fpminimax(f(x), deg, [|double ...|], [a;b]); + +display = hexadecimal; +print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30)); +print("in [",a,b,"]"); +print("coeffs:"); +for i from 2 to deg do round(coeff(poly,i), D, RN); diff --git a/pl/math/tools/expm1f.sollya b/pl/math/tools/expm1f.sollya new file mode 100644 index 0000000..efdf1bd --- /dev/null +++ b/pl/math/tools/expm1f.sollya @@ -0,0 +1,21 @@ +// polynomial for approximating exp(x)-1 in single precision +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 5; + +a = -log(2)/2; +b = log(2)/2; + +f = proc(y) { + return exp(y)-1; +}; + +poly = fpminimax(f(x), deg, [|single ...|], [a;b]); + +display = hexadecimal; +print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30)); +print("in [",a,b,"]"); +print("coeffs:"); +for i from 2 to deg do round(coeff(poly,i), SG, RN); diff --git a/pl/math/tools/log10.sollya b/pl/math/tools/log10.sollya new file mode 100644 index 0000000..85d1d15 --- /dev/null +++ b/pl/math/tools/log10.sollya @@ -0,0 +1,44 @@ +// polynomial for approximating log10(1+x) +// +// Copyright (c) 2019-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 6; // poly degree +// |log10(1+x)| > 0x1p-5 outside the interval +a = -0x1.p-5; +b = 0x1.p-5; + +ln10 = evaluate(log(10),0); +invln10hi = double(1/ln10 + 0x1p21) - 0x1p21; // round away last 21 bits +invln10lo = double(1/ln10 - invln10hi); + +// find log10(1+x)/x polynomial with minimal relative error +// (minimal relative error polynomial for log10(1+x) is the same * x) +deg = deg-1; // because of /x + +// f = log(1+x)/x; using taylor series +f = 0; +for i from 0 to 60 do { f = f + (-x)^i/(i+1); }; +f = f/ln10; + +// 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); +}; + +// first coeff is fixed, iteratively find optimal double prec coeffs +poly = invln10hi + invln10lo; +for i from 1 to deg do { + p = roundcoefficients(approx(poly,i), [|D ...|]); + poly = poly + x^i*coeff(p,0); +}; +display = hexadecimal; +print("invln10hi:", invln10hi); +print("invln10lo:", invln10lo); +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); + +display = decimal; +print("in [",a,b,"]"); diff --git a/pl/math/tools/log10f.sollya b/pl/math/tools/log10f.sollya new file mode 100644 index 0000000..94bf32f --- /dev/null +++ b/pl/math/tools/log10f.sollya @@ -0,0 +1,37 @@ +// polynomial for approximating log10f(1+x) +// +// Copyright (c) 2019-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +// Computation of log10f(1+x) will be carried out in double precision + +deg = 4; // poly degree +// [OFF; 2*OFF] is divided in 2^4 intervals with OFF~0.7 +a = -0.04375; +b = 0.04375; + +// find log(1+x)/x polynomial with minimal relative error +// (minimal relative error polynomial for log(1+x) is the same * x) +deg = deg-1; // because of /x + +// f = log(1+x)/x; using taylor series +f = 0; +for i from 0 to 60 do { f = f + (-x)^i/(i+1); }; + +// 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); +}; + +// first coeff is fixed, iteratively find optimal double prec coeffs +poly = 1; +for i from 1 to deg do { + p = roundcoefficients(approx(poly,i), [|D ...|]); + poly = poly + x^i*coeff(p,0); +}; + +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 double(coeff(poly,i)); diff --git a/pl/math/tools/log1p.sollya b/pl/math/tools/log1p.sollya new file mode 100644 index 0000000..598a36a --- /dev/null +++ b/pl/math/tools/log1p.sollya @@ -0,0 +1,30 @@ +// polynomial for approximating log(1+x) in double precision +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 20; + +a = sqrt(2)/2-1; +b = sqrt(2)-1; + +f = proc(y) { + return log(1+y); +}; + +approx = proc(poly, d) { + return remez(1 - poly(x)/f(x), deg-d, [a;b], x^d/f(x), 1e-10); +}; + +poly = x; +for i from 2 to deg do { + p = roundcoefficients(approx(poly,i), [|D ...|]); + poly = poly + x^i*coeff(p,0); +}; + + +print("coeffs:"); +display = hexadecimal; +for i from 2 to deg do coeff(poly,i); +print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30)); +print("in [",a,b,"]"); diff --git a/pl/math/tools/log1pf.sollya b/pl/math/tools/log1pf.sollya new file mode 100644 index 0000000..cc1db10 --- /dev/null +++ b/pl/math/tools/log1pf.sollya @@ -0,0 +1,21 @@ +// polynomial for approximating log(1+x) in single precision +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 10; + +a = -0.25; +b = 0.5; + +f = proc(y) { + return log(1+y); +}; + +poly = fpminimax(f(x), deg, [|single ...|], [a;b]); + +display = hexadecimal; +print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30)); +print("in [",a,b,"]"); +print("coeffs:"); +for i from 2 to deg do round(coeff(poly,i), SG, RN); diff --git a/pl/math/tools/tan.sollya b/pl/math/tools/tan.sollya new file mode 100644 index 0000000..bb0bb28 --- /dev/null +++ b/pl/math/tools/tan.sollya @@ -0,0 +1,20 @@ +// polynomial for approximating double precision tan(x) +// +// Copyright (c) 2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 8; + +// interval bounds +a = 0x1.0p-126; +b = pi / 8; + +display = hexadecimal; + +f = (tan(sqrt(x))-sqrt(x))/x^(3/2); +poly = fpminimax(f, deg, [|double ...|], [a*a;b*b]); + +//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); diff --git a/pl/math/tools/tanf.sollya b/pl/math/tools/tanf.sollya new file mode 100644 index 0000000..f4b49b4 --- /dev/null +++ b/pl/math/tools/tanf.sollya @@ -0,0 +1,78 @@ +// polynomial for approximating single precision tan(x) +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +dtype = single; + +mthd = 0; // approximate tan +deg = 5; // poly degree + +// // Uncomment for cotan +// mthd = 1; // approximate cotan +// deg = 3; // poly degree + +// interval bounds +a = 0x1.0p-126; +b = pi / 4; + +print("Print some useful constants"); +display = hexadecimal!; +if (dtype==double) then { prec = 53!; } +else if (dtype==single) then { prec = 23!; }; + +print("pi/4"); +pi/4; + +// Setup precisions (display and computation) +display = decimal!; +prec=128!; +save_prec=prec; + +// +// Select function to approximate with Sollya +// +if(mthd==0) then { + s = "x + x^3 * P(x^2)"; + g = tan(x); + F = proc(P) { return x + x^3 * P(x^2); }; + f = (g(sqrt(x))-sqrt(x))/(x*sqrt(x)); + init_poly = 0; + // Display info + print("Approximate g(x) =", g, "as F(x)=", s, "."); + poly = fpminimax(f, deg, [|dtype ...|], [a*a;b*b]); +} +else if (mthd==1) then { + s = "1/x + x * P(x^2)"; + g = 1 / tan(x); + F = proc(P) { return 1/x + x * P(x^2); }; + f = (g(sqrt(x))-1/sqrt(x))/(sqrt(x)); + init_poly = 0; + deg_init_poly = -1; // a value such that we actually start by building constant coefficient + // Display info + print("Approximate g(x) =", g, "as F(x)=", s, "."); + // Fpminimax used to minimise absolute error + approx_fpminimax = proc(func, poly, d) { + return fpminimax(func - poly / x^-(deg-d), 0, [|dtype|], [a;b], absolute, floating); + }; + // Optimise all coefficients at once + poly = fpminimax(f, [|0,...,deg|], [|dtype ...|], [a;b], absolute, floating); +}; + + +// +// Display coefficients in Sollya +// +display = hexadecimal!; +if (dtype==double) then { prec = 53!; } +else if (dtype==single) then { prec = 23!; }; +print("_coeffs :_ hex"); +for i from 0 to deg do coeff(poly, i); + +// Compute errors +display = hexadecimal!; +d_rel_err = dirtyinfnorm(1-F(poly)/g(x), [a;b]); +d_abs_err = dirtyinfnorm(g(x)-F(poly), [a;b]); +print("dirty rel error:", d_rel_err); +print("dirty abs error:", d_abs_err); +print("in [",a,b,"]"); diff --git a/pl/math/tools/v_erf.sollya b/pl/math/tools/v_erf.sollya new file mode 100644 index 0000000..394ba37 --- /dev/null +++ b/pl/math/tools/v_erf.sollya @@ -0,0 +1,20 @@ +// polynomial for approximating erf(x). +// To generate coefficients for interval i (0 to 47) do: +// $ sollya v_erf.sollya $i +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +scale = 1/8; +deg = 9; + +itv = parse(__argv[0]); +if (itv == 0) then { a = 0x1p-1022; } +else { a = itv * scale; }; + +prec=256; + +poly = fpminimax(erf(scale*x+a), deg, [|D ...|], [0; 1]); + +display = hexadecimal; +for i from 0 to deg do coeff(poly, i);
\ No newline at end of file diff --git a/pl/math/tools/v_erfc.sollya b/pl/math/tools/v_erfc.sollya new file mode 100644 index 0000000..3b03ba0 --- /dev/null +++ b/pl/math/tools/v_erfc.sollya @@ -0,0 +1,46 @@ +// polynomial for approximating erfc(x)*exp(x*x) +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 12; // poly degree + +itv = parse(__argv[0]); + +bounds = [|3.725290298461914e-9, + 0.18920711500272103, + 0.41421356237309515, + 0.681792830507429, + 1, + 1.378414230005442, + 1.8284271247461903, + 2.363585661014858, + 3, + 3.756828460010884, + 4.656854249492381, + 5.727171322029716, + 7, + 8.513656920021768, + 10.313708498984761, + 12.454342644059432, + 15, + 18.027313840043536, + 21.627416997969522, + 25.908685288118864, + 31|]; + +a = bounds[itv]; +b = bounds[itv + 1]; + +f = proc(y) { + t = y + a; + return erfc(t) * exp(t*t); +}; + +poly = fpminimax(f(x), deg, [|double ...|], [0;b-a]); + +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); diff --git a/pl/math/tools/v_log10.sollya b/pl/math/tools/v_log10.sollya new file mode 100644 index 0000000..e2df436 --- /dev/null +++ b/pl/math/tools/v_log10.sollya @@ -0,0 +1,38 @@ +// polynomial used for __v_log10(x) +// +// Copyright (c) 2019-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 6; // poly degree +a = -0x1.fc1p-9; +b = 0x1.009p-8; + +// find log(1+x)/x polynomial with minimal relative error +// (minimal relative error polynomial for log(1+x) is the same * x) +deg = deg-1; // because of /x + +// f = log(1+x)/x; using taylor series +f = 0; +for i from 0 to 60 do { f = f + (-x)^i/(i+1); }; + +// 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); +}; + +// first coeff is fixed, iteratively find optimal double prec coeffs +poly = 1; +for i from 1 to deg do { + p = roundcoefficients(approx(poly,i), [|D ...|]); + poly = poly + x^i*coeff(p,0); +}; + +// scale coefficients by 1/ln(10) +ln10 = evaluate(log(10),0); +poly = poly/ln10; + +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 double(coeff(poly,i)); diff --git a/pl/math/tools/v_log10f.sollya b/pl/math/tools/v_log10f.sollya new file mode 100644 index 0000000..396d5a9 --- /dev/null +++ b/pl/math/tools/v_log10f.sollya @@ -0,0 +1,45 @@ +// polynomial for approximating v_log10f(1+x) +// +// Copyright (c) 2019-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 9; // poly degree +// |log10(1+x)| > 0x1p-4 outside the interval +a = -1/3; +b = 1/3; + +display = hexadecimal; +print("log10(2) = ", single(log10(2))); + +ln10 = evaluate(log(10),0); +invln10 = single(1/ln10); + +// find log10(1+x)/x polynomial with minimal relative error +// (minimal relative error polynomial for log10(1+x) is the same * x) +deg = deg-1; // because of /x + +// f = log(1+x)/x; using taylor series +f = 0; +for i from 0 to 60 do { f = f + (-x)^i/(i+1); }; +f = f/ln10; + +// 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); +}; + +// first coeff is fixed, iteratively find optimal double prec coeffs +poly = invln10; +for i from 1 to deg do { + p = roundcoefficients(approx(poly,i), [|SG ...|]); + poly = poly + x^i*coeff(p,0); +}; +display = hexadecimal; +print("invln10:", invln10); +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 single(coeff(poly,i)); + +display = decimal; +print("in [",a,b,"]"); diff --git a/pl/math/tools/v_log2f.sollya b/pl/math/tools/v_log2f.sollya new file mode 100644 index 0000000..99e050c --- /dev/null +++ b/pl/math/tools/v_log2f.sollya @@ -0,0 +1,38 @@ +// polynomial used for __v_log2f(x) +// +// Copyright (c) 2022-2023, Arm Limited. +// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + +deg = 9; // poly degree +a = -1/3; +b = 1/3; + +ln2 = evaluate(log(2),0); +invln2 = single(1/ln2); + +// find log2(1+x)/x polynomial with minimal relative error +// (minimal relative error polynomial for log2(1+x) is the same * x) +deg = deg-1; // because of /x + +// f = log2(1+x)/x; using taylor series +f = 0; +for i from 0 to 60 do { f = f + (-x)^i/(i+1); }; +f = f * invln2; + +// 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); +}; + +// first coeff is fixed, iteratively find optimal double prec coeffs +poly = invln2; +for i from 1 to deg do { + p = roundcoefficients(approx(poly,i), [|SG ...|]); + poly = poly + x^i*coeff(p,0); +}; + +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); |