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