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