1*412f47f9SXin Li// polynomial for approximating 10^x 2*412f47f9SXin Li// 3*412f47f9SXin Li// Copyright (c) 2023, Arm Limited. 4*412f47f9SXin Li// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 5*412f47f9SXin Li 6*412f47f9SXin Li// exp10f parameters 7*412f47f9SXin Lideg = 5; // poly degree 8*412f47f9SXin LiN = 1; // Neon 1, SVE 64 9*412f47f9SXin Lib = log(2)/(2 * N * log(10)); // interval 10*412f47f9SXin Lia = -b; 11*412f47f9SXin Liwp = single; 12*412f47f9SXin Li 13*412f47f9SXin Li// exp10 parameters 14*412f47f9SXin Li//deg = 4; // poly degree - bump to 5 for ~1 ULP 15*412f47f9SXin Li//N = 128; // table size 16*412f47f9SXin Li//b = log(2)/(2 * N * log(10)); // interval 17*412f47f9SXin Li//a = -b; 18*412f47f9SXin Li//wp = D; 19*412f47f9SXin Li 20*412f47f9SXin Li 21*412f47f9SXin Li// find polynomial with minimal relative error 22*412f47f9SXin Li 23*412f47f9SXin Lif = 10^x; 24*412f47f9SXin Li 25*412f47f9SXin Li// return p that minimizes |f(x) - poly(x) - x^d*p(x)|/|f(x)| 26*412f47f9SXin Liapprox = proc(poly,d) { 27*412f47f9SXin Li return remez(1 - poly(x)/f(x), deg-d, [a;b], x^d/f(x), 1e-10); 28*412f47f9SXin Li}; 29*412f47f9SXin Li// return p that minimizes |f(x) - poly(x) - x^d*p(x)| 30*412f47f9SXin Liapprox_abs = proc(poly,d) { 31*412f47f9SXin Li return remez(f(x) - poly(x), deg-d, [a;b], x^d, 1e-10); 32*412f47f9SXin Li}; 33*412f47f9SXin Li 34*412f47f9SXin Li// first coeff is fixed, iteratively find optimal double prec coeffs 35*412f47f9SXin Lipoly = 1; 36*412f47f9SXin Lifor i from 1 to deg do { 37*412f47f9SXin Li p = roundcoefficients(approx(poly,i), [|wp ...|]); 38*412f47f9SXin Li// p = roundcoefficients(approx_abs(poly,i), [|wp ...|]); 39*412f47f9SXin Li poly = poly + x^i*coeff(p,0); 40*412f47f9SXin Li}; 41*412f47f9SXin Li 42*412f47f9SXin Lidisplay = hexadecimal; 43*412f47f9SXin Liprint("rel error:", accurateinfnorm(1-poly(x)/10^x, [a;b], 30)); 44*412f47f9SXin Liprint("abs error:", accurateinfnorm(10^x-poly(x), [a;b], 30)); 45*412f47f9SXin Liprint("in [",a,b,"]"); 46*412f47f9SXin Liprint("coeffs:"); 47*412f47f9SXin Lifor i from 0 to deg do coeff(poly,i); 48*412f47f9SXin Li 49*412f47f9SXin Lilog10_2 = round(N * log(10) / log(2), wp, RN); 50*412f47f9SXin Lilog2_10 = log(2) / (N * log(10)); 51*412f47f9SXin Lilog2_10_hi = round(log2_10, wp, RN); 52*412f47f9SXin Lilog2_10_lo = round(log2_10 - log2_10_hi, wp, RN); 53*412f47f9SXin Liprint(log10_2); 54*412f47f9SXin Liprint(log2_10_hi); 55*412f47f9SXin Liprint(log2_10_lo); 56