Pierre Blanchard | e9e23e5 | 2022-04-11 14:42:01 +0100 | [diff] [blame] | 1 | // polynomial for approximating log10f(1+x) |
| 2 | // |
Joe Ramsay | f0f80b8 | 2023-01-06 09:10:57 +0000 | [diff] [blame^] | 3 | // Copyright (c) 2019-2023, Arm Limited. |
Pierre Blanchard | e9e23e5 | 2022-04-11 14:42:01 +0100 | [diff] [blame] | 4 | // SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception |
| 5 | |
Pierre Blanchard | 38fb9e7 | 2022-04-19 17:54:25 +0100 | [diff] [blame] | 6 | // Computation of log10f(1+x) will be carried out in double precision |
Pierre Blanchard | e9e23e5 | 2022-04-11 14:42:01 +0100 | [diff] [blame] | 7 | |
Pierre Blanchard | 38fb9e7 | 2022-04-19 17:54:25 +0100 | [diff] [blame] | 8 | deg = 4; // poly degree |
| 9 | // [OFF; 2*OFF] is divided in 2^4 intervals with OFF~0.7 |
| 10 | a = -0.04375; |
| 11 | b = 0.04375; |
Pierre Blanchard | e9e23e5 | 2022-04-11 14:42:01 +0100 | [diff] [blame] | 12 | |
Pierre Blanchard | 38fb9e7 | 2022-04-19 17:54:25 +0100 | [diff] [blame] | 13 | // find log(1+x)/x polynomial with minimal relative error |
| 14 | // (minimal relative error polynomial for log(1+x) is the same * x) |
Pierre Blanchard | e9e23e5 | 2022-04-11 14:42:01 +0100 | [diff] [blame] | 15 | deg = deg-1; // because of /x |
| 16 | |
| 17 | // f = log(1+x)/x; using taylor series |
| 18 | f = 0; |
| 19 | for i from 0 to 60 do { f = f + (-x)^i/(i+1); }; |
Pierre Blanchard | e9e23e5 | 2022-04-11 14:42:01 +0100 | [diff] [blame] | 20 | |
| 21 | // return p that minimizes |f(x) - poly(x) - x^d*p(x)|/|f(x)| |
| 22 | approx = proc(poly,d) { |
| 23 | return remez(1 - poly(x)/f(x), deg-d, [a;b], x^d/f(x), 1e-10); |
| 24 | }; |
| 25 | |
| 26 | // first coeff is fixed, iteratively find optimal double prec coeffs |
Pierre Blanchard | 38fb9e7 | 2022-04-19 17:54:25 +0100 | [diff] [blame] | 27 | poly = 1; |
Pierre Blanchard | e9e23e5 | 2022-04-11 14:42:01 +0100 | [diff] [blame] | 28 | for i from 1 to deg do { |
Pierre Blanchard | 38fb9e7 | 2022-04-19 17:54:25 +0100 | [diff] [blame] | 29 | p = roundcoefficients(approx(poly,i), [|D ...|]); |
Pierre Blanchard | e9e23e5 | 2022-04-11 14:42:01 +0100 | [diff] [blame] | 30 | poly = poly + x^i*coeff(p,0); |
| 31 | }; |
Pierre Blanchard | 38fb9e7 | 2022-04-19 17:54:25 +0100 | [diff] [blame] | 32 | |
Pierre Blanchard | e9e23e5 | 2022-04-11 14:42:01 +0100 | [diff] [blame] | 33 | display = hexadecimal; |
Pierre Blanchard | e9e23e5 | 2022-04-11 14:42:01 +0100 | [diff] [blame] | 34 | print("rel error:", accurateinfnorm(1-poly(x)/f(x), [a;b], 30)); |
| 35 | print("in [",a,b,"]"); |
| 36 | print("coeffs:"); |
Pierre Blanchard | 38fb9e7 | 2022-04-19 17:54:25 +0100 | [diff] [blame] | 37 | for i from 0 to deg do double(coeff(poly,i)); |