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