pxfm/tangent/
evalf.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 9/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::polyeval::f_polyeval4;
30use crate::sin_cosf::ArgumentReducer;
31
32// Generated in SageMath:
33// print("[")
34// for k in range(64):
35//     k = RealField(150)(k) * RealField(150).pi() / RealField(150)(32)
36//     print(double_to_hex(k.tan()) + ",")
37// print("];")
38pub(crate) static TAN_K_PI_OVER32: [u64; 64] = [
39    0x0000000000000000,
40    0x3fb936bb8c5b2da2,
41    0x3fc975f5e0553158,
42    0x3fd36a08355c63dc,
43    0x3fda827999fcef32,
44    0x3fe11ab7190834ec,
45    0x3fe561b82ab7f990,
46    0x3fea43002ae42850,
47    0x3ff0000000000000,
48    0x3ff37efd8d87607e,
49    0x3ff7f218e25a7461,
50    0x3ffdef13b73c1406,
51    0x4003504f333f9de6,
52    0x400a5f59e90600dd,
53    0x40141bfee2424771,
54    0x40244e6c595afdcc,
55    0xc950457bf6be49c7,
56    0xc0244e6c595afdcc,
57    0xc0141bfee2424771,
58    0xc00a5f59e90600dd,
59    0xc003504f333f9de6,
60    0xbffdef13b73c1406,
61    0xbff7f218e25a7461,
62    0xbff37efd8d87607e,
63    0xbff0000000000000,
64    0xbfea43002ae42850,
65    0xbfe561b82ab7f990,
66    0xbfe11ab7190834ec,
67    0xbfda827999fcef32,
68    0xbfd36a08355c63dc,
69    0xbfc975f5e0553158,
70    0xbfb936bb8c5b2da2,
71    0x369f77598338bfdf,
72    0x3fb936bb8c5b2da2,
73    0x3fc975f5e0553158,
74    0x3fd36a08355c63dc,
75    0x3fda827999fcef32,
76    0x3fe11ab7190834ec,
77    0x3fe561b82ab7f990,
78    0x3fea43002ae42850,
79    0x3ff0000000000000,
80    0x3ff37efd8d87607e,
81    0x3ff7f218e25a7461,
82    0x3ffdef13b73c1406,
83    0x4003504f333f9de6,
84    0x400a5f59e90600dd,
85    0x40141bfee2424771,
86    0x40244e6c595afdcc,
87    0xc935b1fa9e530d0a,
88    0xc0244e6c595afdcc,
89    0xc0141bfee2424771,
90    0xc00a5f59e90600dd,
91    0xc003504f333f9de6,
92    0xbffdef13b73c1406,
93    0xbff7f218e25a7461,
94    0xbff37efd8d87607e,
95    0xbff0000000000000,
96    0xbfea43002ae42850,
97    0xbfe561b82ab7f990,
98    0xbfe11ab7190834ec,
99    0xbfda827999fcef32,
100    0xbfd36a08355c63dc,
101    0xbfc975f5e0553158,
102    0xbfb936bb8c5b2da2,
103];
104
105pub(crate) struct TanfEval {
106    pub(crate) tan_k: f64,
107    pub(crate) tan_y: f64,
108}
109
110#[inline]
111pub(crate) fn tanpif_eval(y: f64, k: i64) -> TanfEval {
112    let y_sqr = y * y;
113
114    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
115    // So k is an integer and -0.5 <= y <= 0.5.
116
117    // picking minus sin and cos according to quadrant
118    let tan_k = f64::from_bits(TAN_K_PI_OVER32[((k as u64) & 63) as usize]);
119
120    // tan(x*pi/32) generated by Sollya:
121    // d = [0, 0.5];
122    // f_tan = tan(y*pi/32)/y;
123    // Q = fpminimax(f_tan, [|0, 2, 4, 6|], [|D...|], d, relative, floating);
124    // See ./notes/tanpif.sollya
125    let tan_y = f_polyeval4(
126        y_sqr,
127        f64::from_bits(0x3fb921fb54442cef),
128        f64::from_bits(0x3f34abbce63a363e),
129        f64::from_bits(0x3eb466baced705e8),
130        f64::from_bits(0x3e346a33cde88184),
131    ) * y;
132    TanfEval { tan_y, tan_k }
133}
134
135#[inline]
136pub(crate) fn tanf_eval(x: f64, x_abs: u32) -> TanfEval {
137    let (y, k) = ArgumentReducer { x, x_abs }.reduce();
138    let y_sqr = y * y;
139
140    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
141    // So k is an integer and -0.5 <= y <= 0.5.
142
143    // picking minus sin and cos according to quadrant
144    let tan_k = f64::from_bits(TAN_K_PI_OVER32[(k & 63) as usize]);
145
146    // tan(x*pi/32) generated by Sollya:
147    // d = [0, 0.5];
148    // f_tan = tan(x*pi/32)/x;
149    // Q = fpminimax(f_tan, [|0, 2, 4, 6|], [|D...|], d);
150    // See ./notes/tanf.sollya
151    let tan_y = f_polyeval4(
152        y_sqr,
153        f64::from_bits(0x3fb921fb54442cef),
154        f64::from_bits(0x3f34abbce63a363e),
155        f64::from_bits(0x3eb466baced705e8),
156        f64::from_bits(0x3e346a33cde88184),
157    ) * y;
158    TanfEval { tan_y, tan_k }
159}
pxfm/tangent/evalf.rs

pxfm/tangent/
evalf.rs
