pxfm/tangent/
tan.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 6/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::bits::EXP_MASK;
30use crate::common::{dyad_fmla, f_fmla};
31use crate::double_double::DoubleDouble;
32use crate::sin::range_reduction_small;
33use crate::sincos_reduce::LargeArgumentReduction;
34use crate::tangent::tanpi_table::TAN_K_PI_OVER_128;
35
36#[inline]
37pub(crate) fn tan_eval(u: DoubleDouble) -> (DoubleDouble, f64) {
38    // Evaluate tan(y) = tan(x - k * (pi/128))
39    // We use the degree-9 Taylor approximation:
40    //   tan(y) ~ P(y) = y + y^3/3 + 2*y^5/15 + 17*y^7/315 + 62*y^9/2835
41    // Then the error is bounded by:
42    //   |tan(y) - P(y)| < 2^-6 * |y|^11 < 2^-6 * 2^-66 = 2^-72.
43    // For y ~ u_hi + u_lo, fully expanding the polynomial and drop any terms
44    // < ulp(u_hi^3) gives us:
45    //   P(y) = y + y^3/3 + 2*y^5/15 + 17*y^7/315 + 62*y^9/2835 = ...
46    // ~ u_hi + u_hi^3 * (1/3 + u_hi^2 * (2/15 + u_hi^2 * (17/315 +
47    //                                                     + u_hi^2 * 62/2835))) +
48    //        + u_lo (1 + u_hi^2 * (1 + u_hi^2 * 2/3))
49    let u_hi_sq = u.hi * u.hi; // Error < ulp(u_hi^2) < 2^(-6 - 52) = 2^-58.
50    // p1 ~ 17/315 + u_hi^2 62 / 2835.
51    let p1 = f_fmla(
52        u_hi_sq,
53        f64::from_bits(0x3f9664f4882c10fa),
54        f64::from_bits(0x3faba1ba1ba1ba1c),
55    );
56    // p2 ~ 1/3 + u_hi^2 2 / 15.
57    let p2 = f_fmla(
58        u_hi_sq,
59        f64::from_bits(0x3fc1111111111111),
60        f64::from_bits(0x3fd5555555555555),
61    );
62    // q1 ~ 1 + u_hi^2 * 2/3.
63    let q1 = f_fmla(u_hi_sq, f64::from_bits(0x3fe5555555555555), 1.0);
64    let u_hi_3 = u_hi_sq * u.hi;
65    let u_hi_4 = u_hi_sq * u_hi_sq;
66    // p3 ~ 1/3 + u_hi^2 * (2/15 + u_hi^2 * (17/315 + u_hi^2 * 62/2835))
67    let p3 = f_fmla(u_hi_4, p1, p2);
68    // q2 ~ 1 + u_hi^2 * (1 + u_hi^2 * 2/3)
69    let q2 = f_fmla(u_hi_sq, q1, 1.0);
70    let tan_lo = f_fmla(u_hi_3, p3, u.lo * q2);
71    // Overall, |tan(y) - (u_hi + tan_lo)| < ulp(u_hi^3) <= 2^-71.
72    // And the relative errors is:
73    // |(tan(y) - (u_hi + tan_lo)) / tan(y) | <= 2*ulp(u_hi^2) < 2^-64
74    let err = f_fmla(
75        u_hi_3.abs(),
76        f64::from_bits(0x3cc0000000000000),
77        f64::from_bits(0x3990000000000000),
78    );
79    (DoubleDouble::from_exact_add(u.hi, tan_lo), err)
80}
81
82#[inline]
83pub(crate) fn tan_eval_dd(x: DoubleDouble) -> DoubleDouble {
84    // Generated by Sollya:
85    // d = [0, pi/128];
86    // f_tan = tan(x)/x;
87    // Q = fpminimax(f_tan, [|0, 2, 4, 6, 8, 10, 12|], [|107...|], d);
88    const C: [(u64, u64); 7] = [
89        (0x0000000000000000, 0x3ff0000000000000),
90        (0x3c75555549d35a4b, 0x3fd5555555555555),
91        (0x3c413dd6ea580288, 0x3fc1111111111111),
92        (0x3c4e100b946f0c89, 0x3faba1ba1ba1b9fe),
93        (0xbc3c180dfac2b8c3, 0x3f9664f488307189),
94        (0x3c1fd691c256a14a, 0x3f8226e300c1988e),
95        (0x3bec7b08c90fdfc0, 0x3f6d739baeacd204),
96    ];
97    let x2 = DoubleDouble::quick_mult(x, x);
98    let mut p = DoubleDouble::quick_mul_add(
99        x2,
100        DoubleDouble::from_bit_pair(C[6]),
101        DoubleDouble::from_bit_pair(C[5]),
102    );
103    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[4]));
104    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[3]));
105    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[2]));
106    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[1]));
107    p = DoubleDouble::quick_mul_add_f64(x2, p, f64::from_bits(0x3ff0000000000000));
108    let r = DoubleDouble::quick_mult(p, x);
109    DoubleDouble::from_exact_add(r.hi, r.lo)
110}
111
112#[cold]
113fn tan_accurate(y: DoubleDouble, tan_k: DoubleDouble) -> f64 {
114    // Computes tan(x) through identities
115    // tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a)tan(b)) = (tan(y) + tan(k*pi/128)) / (1 - tan(y)*tan(k*pi/128))
116    let tan_y = tan_eval_dd(y);
117
118    // num = tan(y) + tan(k*pi/64)
119    let num_dd = DoubleDouble::full_dd_add(tan_y, tan_k);
120    // den = 1 - tan(y)*tan(k*pi/64)
121    let den_dd = DoubleDouble::mul_add_f64(tan_y, -tan_k, 1.);
122
123    let tan_x = DoubleDouble::div(num_dd, den_dd);
124    tan_x.to_f64()
125}
126
127#[cold]
128fn tan_near_zero_hard(x: f64) -> DoubleDouble {
129    // Generated by Sollya:
130    // d = [0, 0.03125];
131    // f_tan = tan(x)/x;
132    // Q = fpminimax(f_tan, [|0, 2, 4, 6, 8, 10, 12|], [|1, 107...|], d);
133    const C: [(u64, u64); 7] = [
134        (0x0000000000000000, 0x3ff0000000000000),
135        (0x3c7555548455da94, 0x3fd5555555555555),
136        (0x3c4306a6358cc810, 0x3fc1111111111111),
137        (0x3c2ca9cd025ea98c, 0x3faba1ba1ba1b952),
138        (0x3c3cb012803c55a5, 0x3f9664f4883eb962),
139        (0x3c28afc1adb8f202, 0x3f8226e276097461),
140        (0xbbdf8f911392f348, 0x3f6d7791601277ea),
141    ];
142    let x2 = DoubleDouble::from_exact_mult(x, x);
143    let mut p = DoubleDouble::quick_mul_add(
144        x2,
145        DoubleDouble::from_bit_pair(C[6]),
146        DoubleDouble::from_bit_pair(C[5]),
147    );
148    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[4]));
149    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[3]));
150    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[2]));
151    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[1]));
152    p = DoubleDouble::quick_mul_add_f64(x2, p, f64::from_bits(0x3ff0000000000000));
153    DoubleDouble::quick_mult_f64(p, x)
154}
155
156/// Tangent in double precision
157///
158/// ULP 0.5
159pub fn f_tan(x: f64) -> f64 {
160    let x_e = (x.to_bits() >> 52) & 0x7ff;
161    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
162
163    let y: DoubleDouble;
164    let k;
165
166    let mut argument_reduction = LargeArgumentReduction::default();
167
168    if x_e < E_BIAS + 16 {
169        // |x| < 2^16
170        if x_e < E_BIAS - 5 {
171            // |x| < 2^-5
172            if x_e < E_BIAS - 27 {
173                // |x| < 2^-27, |tan(x) - x| < ulp(x)/2.
174                if x == 0.0 {
175                    // Signed zeros.
176                    return x + x;
177                }
178                return dyad_fmla(x, f64::from_bits(0x3c90000000000000), x);
179            }
180            let x2 = x * x;
181            let x4 = x2 * x2;
182            // Generated by Sollya:
183            // d = [0, 0.03125];
184            // f_tan = tan(x)/x;
185            // Q = fpminimax(f_tan, [|0, 2, 4, 6, 8|], [|D...|], d);
186            const C: [u64; 4] = [
187                0x3fd555555555554b,
188                0x3fc1111111142d5b,
189                0x3faba1b9860ed843,
190                0x3f966a802210f5bb,
191            ];
192            let p01 = f_fmla(x2, f64::from_bits(C[1]), f64::from_bits(C[0]));
193            let p23 = f_fmla(x2, f64::from_bits(C[3]), f64::from_bits(C[2]));
194            let w0 = f_fmla(x4, p23, p01);
195            let w1 = x2 * w0 * x;
196            let r = DoubleDouble::from_exact_add(x, w1);
197            let err = f_fmla(
198                x2,
199                f64::from_bits(0x3cb0000000000000), // 2^-52
200                f64::from_bits(0x3be0000000000000), // 2^-65
201            );
202            let ub = r.hi + (r.lo + err);
203            let lb = r.hi + (r.lo - err);
204            if ub == lb {
205                return ub;
206            }
207            return tan_near_zero_hard(x).to_f64();
208        } else {
209            // Small range reduction.
210            (y, k) = range_reduction_small(x);
211        }
212    } else {
213        // Inf or NaN
214        if x_e > 2 * E_BIAS {
215            if x.is_nan() {
216                return f64::NAN;
217            }
218            // tan(+-Inf) = NaN
219            return x + f64::NAN;
220        }
221
222        // Large range reduction.
223        (k, y) = argument_reduction.reduce(x);
224    }
225
226    let (tan_y, err) = tan_eval(y);
227
228    // Computes tan(x) through identities.
229    // tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a)tan(b)) = (tan(y) + tan(k*pi/128)) / (1 - tan(y)*tan(k*pi/128))
230    let tan_k = DoubleDouble::from_bit_pair(TAN_K_PI_OVER_128[(k & 255) as usize]);
231
232    // num = tan(y) + tan(k*pi/64)
233    let num_dd = DoubleDouble::add(tan_y, tan_k);
234    // den = 1 - tan(y)*tan(k*pi/64)
235    let den_dd = DoubleDouble::mul_add_f64(tan_y, -tan_k, 1.);
236
237    let tan_x = DoubleDouble::div(num_dd, den_dd);
238
239    // Simple error bound: |1 / den_dd| < 2^(1 + floor(-log2(den_dd)))).
240    let den_inv = ((E_BIAS + 1) << (52 + 1)) - (den_dd.hi.to_bits() & EXP_MASK);
241    // For tan_x = (num_dd + err) / (den_dd + err), the error is bounded by:
242    //   | tan_x - num_dd / den_dd |  <= err * ( 1 + | tan_x * den_dd | ).
243    let tan_err = err * f_fmla(f64::from_bits(den_inv), tan_x.hi.abs(), 1.0);
244
245    let err_higher = tan_x.lo + tan_err;
246    let err_lower = tan_x.lo - tan_err;
247
248    let tan_upper = tan_x.hi + err_higher;
249    let tan_lower = tan_x.hi + err_lower;
250
251    // Ziv_s rounding test.
252    if tan_upper == tan_lower {
253        return tan_upper;
254    }
255
256    tan_accurate(y, tan_k)
257}
258
259#[cfg(test)]
260mod tests {
261    use super::*;
262
263    #[test]
264    fn tan_test() {
265        assert_eq!(f_tan(0.0003242153424), 0.0003242153537600293);
266        assert_eq!(f_tan(-0.3242153424), -0.3360742441422686);
267        assert_eq!(f_tan(0.3242153424), 0.3360742441422686);
268        assert_eq!(f_tan(-97301943219974110.), 0.000000000000000481917514979303);
269        assert_eq!(f_tan(0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000007291122019556397),
270                   0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000007291122019556397);
271        assert_eq!(f_tan(0.0), 0.0);
272        assert_eq!(f_tan(0.4324324324), 0.4615683710506729);
273        assert_eq!(f_tan(1.0), 1.5574077246549023);
274        assert_eq!(f_tan(-0.5), -0.5463024898437905);
275        assert!(f_tan(f64::INFINITY).is_nan());
276        assert!(f_tan(f64::NEG_INFINITY).is_nan());
277        assert!(f_tan(f64::NAN).is_nan());
278    }
279}
pxfm/tangent/tan.rs

pxfm/tangent/
tan.rs
