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(always)]
136#[allow(unused)]
137pub(crate) fn tanpif_eval_fma(y: f64, k: i64) -> TanfEval {
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 as u64) & 63) as usize]);
145
146    // tan(x*pi/32) generated by Sollya:
147    // d = [0, 0.5];
148    // f_tan = tan(y*pi/32)/y;
149    // Q = fpminimax(f_tan, [|0, 2, 4, 6|], [|D...|], d, relative, floating);
150    // See ./notes/tanpif.sollya
151    use crate::polyeval::d_polyeval4;
152    let tan_y = d_polyeval4(
153        y_sqr,
154        f64::from_bits(0x3fb921fb54442cef),
155        f64::from_bits(0x3f34abbce63a363e),
156        f64::from_bits(0x3eb466baced705e8),
157        f64::from_bits(0x3e346a33cde88184),
158    ) * y;
159    TanfEval { tan_y, tan_k }
160}
161
162#[inline(always)]
163pub(crate) fn tanf_eval(x: f64, x_abs: u32) -> TanfEval {
164    let (y, k) = ArgumentReducer { x, x_abs }.reduce();
165    let y_sqr = y * y;
166
167    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
168    // So k is an integer and -0.5 <= y <= 0.5.
169
170    // picking minus sin and cos according to quadrant
171    let tan_k = f64::from_bits(TAN_K_PI_OVER32[(k & 63) as usize]);
172
173    // tan(x*pi/32) generated by Sollya:
174    // d = [0, 0.5];
175    // f_tan = tan(x*pi/32)/x;
176    // Q = fpminimax(f_tan, [|0, 2, 4, 6|], [|D...|], d);
177    // See ./notes/tanf.sollya
178    let tan_y = f_polyeval4(
179        y_sqr,
180        f64::from_bits(0x3fb921fb54442cef),
181        f64::from_bits(0x3f34abbce63a363e),
182        f64::from_bits(0x3eb466baced705e8),
183        f64::from_bits(0x3e346a33cde88184),
184    ) * y;
185    TanfEval { tan_y, tan_k }
186}
187
188#[inline(always)]
189#[allow(unused)]
190pub(crate) fn tanf_eval_fma(x: f64, x_abs: u32) -> TanfEval {
191    let (y, k) = ArgumentReducer { x, x_abs }.reduce_fma();
192    let y_sqr = y * y;
193
194    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
195    // So k is an integer and -0.5 <= y <= 0.5.
196
197    // picking minus sin and cos according to quadrant
198    let tan_k = f64::from_bits(TAN_K_PI_OVER32[(k & 63) as usize]);
199
200    // tan(x*pi/32) generated by Sollya:
201    // d = [0, 0.5];
202    // f_tan = tan(x*pi/32)/x;
203    // Q = fpminimax(f_tan, [|0, 2, 4, 6|], [|D...|], d);
204    // See ./notes/tanf.sollya
205    use crate::polyeval::d_polyeval4;
206    let tan_y = d_polyeval4(
207        y_sqr,
208        f64::from_bits(0x3fb921fb54442cef),
209        f64::from_bits(0x3f34abbce63a363e),
210        f64::from_bits(0x3eb466baced705e8),
211        f64::from_bits(0x3e346a33cde88184),
212    ) * y;
213    TanfEval { tan_y, tan_k }
214}
pxfm/tangent/evalf.rs

pxfm/tangent/
evalf.rs
