xref: /freebsd-src/contrib/arm-optimized-routines/math/tools/log2_abs.sollya (revision 072a4ba82a01476eaee33781ccd241033eefcf0b)
131914882SAlex Richardson// polynomial for approximating log2(1+x)
231914882SAlex Richardson//
331914882SAlex Richardson// Copyright (c) 2019, Arm Limited.
4*072a4ba8SAndrew Turner// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
531914882SAlex Richardson
631914882SAlex Richardsondeg = 7; // poly degree
731914882SAlex Richardson// interval ~= 1/(2*N), where N is the table entries
831914882SAlex Richardsona= -0x1.f45p-8;
931914882SAlex Richardsonb=  0x1.f45p-8;
1031914882SAlex Richardson
1131914882SAlex Richardsonln2 = evaluate(log(2),0);
1231914882SAlex Richardsoninvln2hi = double(1/ln2 + 0x1p21) - 0x1p21; // round away last 21 bits
1331914882SAlex Richardsoninvln2lo = double(1/ln2 - invln2hi);
1431914882SAlex Richardson
1531914882SAlex Richardson// find log2(1+x) polynomial with minimal absolute error
1631914882SAlex Richardsonf = log(1+x)/ln2;
1731914882SAlex Richardson
1831914882SAlex Richardson// return p that minimizes |f(x) - poly(x) - x^d*p(x)|
1931914882SAlex Richardsonapprox = proc(poly,d) {
2031914882SAlex Richardson  return remez(f(x) - poly(x), deg-d, [a;b], x^d, 1e-10);
2131914882SAlex Richardson};
2231914882SAlex Richardson
2331914882SAlex Richardson// first coeff is fixed, iteratively find optimal double prec coeffs
2431914882SAlex Richardsonpoly = x*(invln2lo + invln2hi);
2531914882SAlex Richardsonfor i from 2 to deg do {
2631914882SAlex Richardson  p = roundcoefficients(approx(poly,i), [|D ...|]);
2731914882SAlex Richardson  poly = poly + x^i*coeff(p,0);
2831914882SAlex Richardson};
2931914882SAlex Richardson
3031914882SAlex Richardsondisplay = hexadecimal;
3131914882SAlex Richardsonprint("invln2hi:", invln2hi);
3231914882SAlex Richardsonprint("invln2lo:", invln2lo);
3331914882SAlex Richardsonprint("abs error:", accurateinfnorm(f(x)-poly(x), [a;b], 30));
3431914882SAlex Richardson//// relative error computation fails if f(0)==0
3531914882SAlex Richardson//// g = f(x)/x = log2(1+x)/x; using taylor series
3631914882SAlex Richardson//g = 0;
3731914882SAlex Richardson//for i from 0 to 60 do { g = g + (-x)^i/(i+1)/ln2; };
3831914882SAlex Richardson//print("rel error:", accurateinfnorm(1-(poly(x)/x)/g(x), [a;b], 30));
3931914882SAlex Richardsonprint("in [",a,b,"]");
4031914882SAlex Richardsonprint("coeffs:");
4131914882SAlex Richardsonfor i from 0 to deg do coeff(poly,i);
42