pxfm/tangent/
tanpi.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 7/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::double_double::DoubleDouble;
30use crate::sincospi::{GenSinCosPiBackend, SinCosPiBackend};
31use crate::tangent::tanpi_table::TANPI_K_PI_OVER_64;
32
33#[inline]
34pub(crate) fn tanpi_eval<B: SinCosPiBackend>(x: f64, backend: &B) -> DoubleDouble {
35    let x2 = DoubleDouble::from_exact_mult(x, x);
36    // tan(pi*x) generated by Sollya:
37    // d = [0, 0.0078128];
38    // f_tan = tan(y*pi)/y;
39    // Q = fpminimax(f_tan, [|0, 2, 4, 6, 8|], [|107, D...|], d, relative, floating);
40    const C: [u64; 4] = [
41        0x4024abbce625be51,
42        0x404466bc677698e5,
43        0x40645fff70379ae3,
44        0x4084626b091b7fd0,
45    ];
46    const C0: DoubleDouble = DoubleDouble::from_bit_pair((0x3ca1a6444aa5b996, 0x400921fb54442d18));
47
48    // polyeval 4, estrin scheme
49    let u0 = backend.fma(x2.hi, f64::from_bits(C[1]), f64::from_bits(C[0]));
50    let u1 = backend.fma(x2.hi, f64::from_bits(C[3]), f64::from_bits(C[2]));
51    let tan_poly_lo = backend.fma(x2.hi * x2.hi, u1, u0);
52
53    // We're splitting polynomial in two parts, since first term dominates
54    // we compute: (a0_lo + a0_hi) * x + x * (a1 * x^2 + a2 + x^4) ...
55    let r_lo = backend.quick_mult_f64(x2, tan_poly_lo);
56    let tan_lo = backend.fma(r_lo.lo, x, r_lo.hi * x);
57    let tan_hi = backend.quick_mult_f64(C0, x);
58    DoubleDouble::full_add_f64(tan_hi, tan_lo)
59}
60
61#[cold]
62fn tanpi_hard<B: SinCosPiBackend>(x: f64, tan_k: DoubleDouble, backend: &B) -> DoubleDouble {
63    const C: [(u64, u64); 6] = [
64        (0x3ca1a62632712fc8, 0x400921fb54442d18),
65        (0xbcc052338fbb4528, 0x4024abbce625be53),
66        (0x3ced42454c5f85b3, 0x404466bc6775aad9),
67        (0xbd00c7d6a971a560, 0x40645fff9b4b244d),
68        (0x3d205970eff53274, 0x40845f46e96c3a0b),
69        (0xbd3589489ad24fc4, 0x40a4630551cd123d),
70    ];
71    let x2 = backend.exact_mult(x, x);
72    let mut tan_y = backend.quick_mul_add(
73        x2,
74        DoubleDouble::from_bit_pair(C[5]),
75        DoubleDouble::from_bit_pair(C[4]),
76    );
77    tan_y = backend.quick_mul_add(x2, tan_y, DoubleDouble::from_bit_pair(C[3]));
78    tan_y = backend.quick_mul_add(x2, tan_y, DoubleDouble::from_bit_pair(C[2]));
79    tan_y = backend.quick_mul_add(x2, tan_y, DoubleDouble::from_bit_pair(C[1]));
80    tan_y = backend.quick_mul_add(x2, tan_y, DoubleDouble::from_bit_pair(C[0]));
81    tan_y = backend.quick_mult_f64(tan_y, x);
82
83    // num = tan(y*pi/64) + tan(k*pi/64)
84    let num = DoubleDouble::full_dd_add(tan_y, tan_k);
85    // den = 1 - tan(y*pi/64)*tan(k*pi/64)
86    let den = backend.mul_add_f64(tan_y, -tan_k, 1.);
87    // tan = num / den
88    backend.div(num, den)
89}
90
91#[inline(always)]
92fn tanpi_gen_impl<B: SinCosPiBackend>(x: f64, backend: B) -> f64 {
93    if x == 0. {
94        return x;
95    }
96    let ax = x.to_bits() & 0x7fff_ffff_ffff_ffff;
97    if ax >= (0x7ffu64 << 52) {
98        // NaN, Inf
99        if ax > (0x7ffu64 << 52) {
100            return x + x;
101        } // NaN
102        return f64::NAN; // x=Inf
103    }
104    let e: i32 = (ax >> 52) as i32 - 1023;
105    if e > 0 {
106        if e >= 52 {
107            // when |x| > 2^53 it's always an integer
108            return f64::copysign(0., x);
109        }
110        // |x| > 1 and |x| < 2^53
111        let m = (ax & ((1u64 << 52) - 1)) | (1u64 << 52); // mantissa with hidden 1
112        let shift = 52 - e;
113
114        let frac = m & ((1u64 << shift) - 1);
115        if frac == (1u64 << (shift - 1)) {
116            // |x| is always integer.5 means it's inf
117            return f64::INFINITY;
118        }
119    }
120
121    if ax <= 0x3cb0000000000000 {
122        // for tiny x ( |x| < f64::EPSILON ) just small taylor expansion
123        // tan(PI*x) ~ PI*x + PI^3*x^3/3 + O(x^5)
124        const PI: DoubleDouble =
125            DoubleDouble::from_bit_pair((0x3ca1a62633145c07, 0x400921fb54442d18));
126        if ax <= 0x3ca0000000000000 {
127            // |x| <= 2^-53, renormalize value
128            let e: i32 = (ax >> 52) as i32;
129            let sc = f64::from_bits((2045i64 - e as i64).wrapping_shl(52) as u64);
130            let isc = f64::from_bits(1i64.wrapping_add(e as i64).wrapping_shl(52) as u64);
131            let dx = x * sc;
132            let q0 = backend.quick_mult_f64(PI, dx);
133            let r = q0.to_f64() * isc;
134            return r;
135        }
136        let q0 = backend.quick_mult_f64(PI, x);
137        let r = q0.to_f64();
138        return r;
139    }
140
141    // argument reduction
142    let (y, k) = backend.arg_reduce_pi_64(x);
143
144    if y == 0.0 {
145        let km = (k.abs() & 63) as i32; // k mod 64
146
147        match km {
148            0 => return f64::copysign(0f64, x),           // tanpi(n) = 0
149            32 => return f64::copysign(f64::INFINITY, x), // tanpi(n+0.5) = ±∞
150            16 => return f64::copysign(1.0, x),           // tanpi(n+0.25) = ±1
151            48 => return -f64::copysign(1.0, x),          // tanpi(n+0.75) = ∓1
152            _ => {}
153        }
154    }
155
156    let tan_k = DoubleDouble::from_bit_pair(TANPI_K_PI_OVER_64[((k as u64) & 127) as usize]);
157
158    // Computes tan(pi*x) through identities.
159    // tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a)tan(b)) = (tan(y*pi/64) + tan(k*pi/64)) / (1 - tan(y*pi/64)*tan(k*pi/64))
160    let tan_y = tanpi_eval(y, &backend);
161    // num = tan(y*pi/64) + tan(k*pi/64)
162    let num = DoubleDouble::add(tan_k, tan_y);
163    // den = 1 - tan(y*pi/64)*tan(k*pi/64)
164    let den = backend.mul_add_f64(tan_y, -tan_k, 1.);
165    // tan = num / den
166    let tan_value = backend.div(num, den);
167    let err = backend.fma(
168        tan_value.hi,
169        f64::from_bits(0x3bf0000000000000), // 2^-64
170        f64::from_bits(0x3b60000000000000), // 2^-73
171    );
172    let ub = tan_value.hi + (tan_value.lo + err);
173    let lb = tan_value.hi + (tan_value.lo - err);
174    if ub == lb {
175        return tan_value.to_f64();
176    }
177    tanpi_hard(y, tan_k, &backend).to_f64()
178}
179
180#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
181#[target_feature(enable = "avx", enable = "fma")]
182unsafe fn tanpi_fma_impl(x: f64) -> f64 {
183    use crate::sincospi::FmaSinCosPiBackend;
184    tanpi_gen_impl(x, FmaSinCosPiBackend {})
185}
186
187/// Computes tan(PI*x)
188///
189/// Max found ULP 0.5
190pub fn f_tanpi(x: f64) -> f64 {
191    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
192    {
193        tanpi_gen_impl(x, GenSinCosPiBackend {})
194    }
195    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
196    {
197        use std::sync::OnceLock;
198        static EXECUTOR: OnceLock<unsafe fn(f64) -> f64> = OnceLock::new();
199        let q = EXECUTOR.get_or_init(|| {
200            if std::arch::is_x86_feature_detected!("avx")
201                && std::arch::is_x86_feature_detected!("fma")
202            {
203                tanpi_fma_impl
204            } else {
205                fn def_tanpi(x: f64) -> f64 {
206                    tanpi_gen_impl(x, GenSinCosPiBackend {})
207                }
208                def_tanpi
209            }
210        });
211        unsafe { q(x) }
212    }
213}
214
215#[cfg(test)]
216mod tests {
217    use super::*;
218
219    #[test]
220    fn test_tanpi() {
221        assert_eq!(f_tanpi(0.4999999999119535), 3615246871.564404);
222        assert_eq!(f_tanpi(7119681148991743.0), 0.);
223        assert_eq!(f_tanpi(63.5), f64::INFINITY);
224        assert_eq!(f_tanpi(63.99935913085936), -0.0020133525045719896);
225        assert_eq!(f_tanpi(3.3821122649309461E-306), 1.0625219045122997E-305);
226        assert_eq!(f_tanpi(1.8010707049867402E-255), 5.6582304953821333E-255);
227        assert_eq!(f_tanpi(1.001000000061801), 0.0031416031832113213);
228        assert_eq!(f_tanpi(-0.5000000000000226), 14054316517702.594);
229        assert_eq!(f_tanpi(0.5000000000000001), -2867080569611329.5);
230        assert_eq!(f_tanpi(0.02131), 0.06704753721009375);
231        assert!(f_tanpi(f64::INFINITY).is_nan());
232        assert!(f_tanpi(f64::NAN).is_nan());
233        assert!(f_tanpi(f64::NEG_INFINITY).is_nan());
234    }
235}
pxfm/tangent/tanpi.rs

pxfm/tangent/
tanpi.rs
