1*f3087befSAndrew Turner// polynomial for approximating tanpi/f(x) 2*f3087befSAndrew Turner// 3*f3087befSAndrew Turner// Copyright (c) 2024, Arm Limited. 4*f3087befSAndrew Turner// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 5*f3087befSAndrew Turner 6*f3087befSAndrew Turner// 0 for tanpi/f [0,0.25], 1 for tanpi/f [0.25,1] 7*f3087befSAndrew Turnermethod = 0; 8*f3087befSAndrew Turnerdtype = double; 9*f3087befSAndrew Turner 10*f3087befSAndrew Turnerif (dtype == single) then { 11*f3087befSAndrew Turner if (method == 0) then { deg = 5; } 12*f3087befSAndrew Turner else if (method == 1) then { deg = 3; }; 13*f3087befSAndrew Turner} else if (dtype == double) then { 14*f3087befSAndrew Turner if (method == 0) then { deg = 13; } 15*f3087befSAndrew Turner else if (method == 1) then { deg = 8; }; 16*f3087befSAndrew Turner}; 17*f3087befSAndrew Turner 18*f3087befSAndrew Turnera = 0x1.0p-126; 19*f3087befSAndrew Turnerb = 1/4; 20*f3087befSAndrew Turner 21*f3087befSAndrew Turnerif (method == 0) then { 22*f3087befSAndrew Turner g = tan(pi * x); 23*f3087befSAndrew Turner F = proc(P) { return pi * x + x^3 * P(x^2); }; 24*f3087befSAndrew Turner f = (g(sqrt(x)) - pi * sqrt(x))/(x^(3/2)); 25*f3087befSAndrew Turner} else if (method == 1) then { 26*f3087befSAndrew Turner g = 1/tan(pi * x); 27*f3087befSAndrew Turner F = proc(P) { return 1/(pi * x) + x * P(x^2); }; 28*f3087befSAndrew Turner f = (g(sqrt(x)) / sqrt(x)) - 1/(pi * x); 29*f3087befSAndrew Turner}; 30*f3087befSAndrew Turner 31*f3087befSAndrew Turnerpoly = fpminimax(f, deg, [|dtype ...|], [a*a;b*b]); 32*f3087befSAndrew Turner 33*f3087befSAndrew Turner// 34*f3087befSAndrew Turner// Display coefficients in Sollya 35*f3087befSAndrew Turner// 36*f3087befSAndrew Turnerdisplay = hexadecimal!; 37*f3087befSAndrew Turnerif (dtype==double) then { prec = 53!; } 38*f3087befSAndrew Turnerelse if (dtype==single) then { prec = 23!; }; 39*f3087befSAndrew Turnerprint("_coeffs :_ hex"); 40*f3087befSAndrew Turnerfor i from 0 to deg do coeff(poly, i); 41*f3087befSAndrew Turner 42*f3087befSAndrew Turner// Compute errors 43*f3087befSAndrew Turner//display = hexadecimal!; 44*f3087befSAndrew Turnerd_rel_err = dirtyinfnorm(1-F(poly)/g(x), [a;b]); 45*f3087befSAndrew Turnerd_abs_err = dirtyinfnorm(g(x)-F(poly), [a;b]); 46*f3087befSAndrew Turnerprint("dirty rel error:", d_rel_err); 47*f3087befSAndrew Turnerprint("dirty abs error:", d_abs_err); 48*f3087befSAndrew Turnerprint("in [",a,b,"]"); 49