1*f3087befSAndrew Turner// polynomial for approximating single precision tan(x) 2*f3087befSAndrew Turner// 3*f3087befSAndrew Turner// Copyright (c) 2022-2024, Arm Limited. 4*f3087befSAndrew Turner// SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 5*f3087befSAndrew Turner 6*f3087befSAndrew Turnerdtype = single; 7*f3087befSAndrew Turner 8*f3087befSAndrew Turnermthd = 0; // approximate tan 9*f3087befSAndrew Turnerdeg = 5; // poly degree 10*f3087befSAndrew Turner 11*f3087befSAndrew Turner// // Uncomment for cotan 12*f3087befSAndrew Turner// mthd = 1; // approximate cotan 13*f3087befSAndrew Turner// deg = 3; // poly degree 14*f3087befSAndrew Turner 15*f3087befSAndrew Turner// interval bounds 16*f3087befSAndrew Turnera = 0x1.0p-126; 17*f3087befSAndrew Turnerb = pi / 4; 18*f3087befSAndrew Turner 19*f3087befSAndrew Turnerprint("Print some useful constants"); 20*f3087befSAndrew Turnerdisplay = hexadecimal!; 21*f3087befSAndrew Turnerif (dtype==double) then { prec = 53!; } 22*f3087befSAndrew Turnerelse if (dtype==single) then { prec = 23!; }; 23*f3087befSAndrew Turner 24*f3087befSAndrew Turnerprint("pi/4"); 25*f3087befSAndrew Turnerpi/4; 26*f3087befSAndrew Turner 27*f3087befSAndrew Turner// Setup precisions (display and computation) 28*f3087befSAndrew Turnerdisplay = decimal!; 29*f3087befSAndrew Turnerprec=128!; 30*f3087befSAndrew Turnersave_prec=prec; 31*f3087befSAndrew Turner 32*f3087befSAndrew Turner// 33*f3087befSAndrew Turner// Select function to approximate with Sollya 34*f3087befSAndrew Turner// 35*f3087befSAndrew Turnerif(mthd==0) then { 36*f3087befSAndrew Turner s = "x + x^3 * P(x^2)"; 37*f3087befSAndrew Turner g = tan(x); 38*f3087befSAndrew Turner F = proc(P) { return x + x^3 * P(x^2); }; 39*f3087befSAndrew Turner f = (g(sqrt(x))-sqrt(x))/(x*sqrt(x)); 40*f3087befSAndrew Turner init_poly = 0; 41*f3087befSAndrew Turner // Display info 42*f3087befSAndrew Turner print("Approximate g(x) =", g, "as F(x)=", s, "."); 43*f3087befSAndrew Turner poly = fpminimax(f, deg, [|dtype ...|], [a*a;b*b]); 44*f3087befSAndrew Turner} 45*f3087befSAndrew Turnerelse if (mthd==1) then { 46*f3087befSAndrew Turner s = "1/x + x * P(x^2)"; 47*f3087befSAndrew Turner g = 1 / tan(x); 48*f3087befSAndrew Turner F = proc(P) { return 1/x + x * P(x^2); }; 49*f3087befSAndrew Turner f = (g(sqrt(x))-1/sqrt(x))/(sqrt(x)); 50*f3087befSAndrew Turner init_poly = 0; 51*f3087befSAndrew Turner deg_init_poly = -1; // a value such that we actually start by building constant coefficient 52*f3087befSAndrew Turner // Display info 53*f3087befSAndrew Turner print("Approximate g(x) =", g, "as F(x)=", s, "."); 54*f3087befSAndrew Turner // Fpminimax used to minimise absolute error 55*f3087befSAndrew Turner approx_fpminimax = proc(func, poly, d) { 56*f3087befSAndrew Turner return fpminimax(func - poly / x^-(deg-d), 0, [|dtype|], [a;b], absolute, floating); 57*f3087befSAndrew Turner }; 58*f3087befSAndrew Turner // Optimise all coefficients at once 59*f3087befSAndrew Turner poly = fpminimax(f, [|0,...,deg|], [|dtype ...|], [a;b], absolute, floating); 60*f3087befSAndrew Turner}; 61*f3087befSAndrew Turner 62*f3087befSAndrew Turner 63*f3087befSAndrew Turner// 64*f3087befSAndrew Turner// Display coefficients in Sollya 65*f3087befSAndrew Turner// 66*f3087befSAndrew Turnerdisplay = hexadecimal!; 67*f3087befSAndrew Turnerif (dtype==double) then { prec = 53!; } 68*f3087befSAndrew Turnerelse if (dtype==single) then { prec = 23!; }; 69*f3087befSAndrew Turnerprint("_coeffs :_ hex"); 70*f3087befSAndrew Turnerfor i from 0 to deg do coeff(poly, i); 71*f3087befSAndrew Turner 72*f3087befSAndrew Turner// Compute errors 73*f3087befSAndrew Turnerdisplay = hexadecimal!; 74*f3087befSAndrew Turnerd_rel_err = dirtyinfnorm(1-F(poly)/g(x), [a;b]); 75*f3087befSAndrew Turnerd_abs_err = dirtyinfnorm(g(x)-F(poly), [a;b]); 76*f3087befSAndrew Turnerprint("dirty rel error:", d_rel_err); 77*f3087befSAndrew Turnerprint("dirty abs error:", d_abs_err); 78*f3087befSAndrew Turnerprint("in [",a,b,"]"); 79