pxfm/tangent/
cotpi.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 8/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, reduce_pi_64};
31use crate::tangent::tanpi::tanpi_eval;
32use crate::tangent::tanpi_table::TANPI_K_PI_OVER_64;
33
34#[cold]
35fn cotpi_hard(x: f64, tan_k: DoubleDouble) -> DoubleDouble {
36    const C: [(u64, u64); 6] = [
37        (0x3ca1a62632712fc8, 0x400921fb54442d18),
38        (0xbcc052338fbb4528, 0x4024abbce625be53),
39        (0x3ced42454c5f85b3, 0x404466bc6775aad9),
40        (0xbd00c7d6a971a560, 0x40645fff9b4b244d),
41        (0x3d205970eff53274, 0x40845f46e96c3a0b),
42        (0xbd3589489ad24fc4, 0x40a4630551cd123d),
43    ];
44    let x2 = DoubleDouble::from_exact_mult(x, x);
45    let mut tan_y = DoubleDouble::quick_mul_add(
46        x2,
47        DoubleDouble::from_bit_pair(C[5]),
48        DoubleDouble::from_bit_pair(C[4]),
49    );
50    tan_y = DoubleDouble::quick_mul_add(x2, tan_y, DoubleDouble::from_bit_pair(C[3]));
51    tan_y = DoubleDouble::quick_mul_add(x2, tan_y, DoubleDouble::from_bit_pair(C[2]));
52    tan_y = DoubleDouble::quick_mul_add(x2, tan_y, DoubleDouble::from_bit_pair(C[1]));
53    tan_y = DoubleDouble::quick_mul_add(x2, tan_y, DoubleDouble::from_bit_pair(C[0]));
54    tan_y = DoubleDouble::quick_mult_f64(tan_y, x);
55
56    // num = tan(y*pi/64) + tan(k*pi/64)
57    let num = DoubleDouble::full_dd_add(tan_y, tan_k);
58    // den = 1 - tan(y*pi/64)*tan(k*pi/64)
59    let den = DoubleDouble::mul_add_f64(tan_y, -tan_k, 1.);
60    // cot = den / num
61    DoubleDouble::div(den, num)
62}
63
64#[inline(always)]
65fn cotpi_gen_impl<B: SinCosPiBackend>(x: f64, backend: B) -> f64 {
66    if x == 0. {
67        return if x.is_sign_negative() {
68            f64::NEG_INFINITY
69        } else {
70            f64::INFINITY
71        };
72    }
73    let ax = x.to_bits() & 0x7fff_ffff_ffff_ffff;
74    if ax >= (0x7ffu64 << 52) {
75        // NaN, Inf
76        if ax > (0x7ffu64 << 52) {
77            return x + x;
78        } // NaN
79        return f64::NAN; // x=Inf
80    }
81    let e: i32 = (ax >> 52) as i32 - 1023;
82    if e > 0 {
83        if e >= 52 {
84            // when |x| > 2^53 it's always an integer
85            return f64::copysign(f64::INFINITY, x);
86        }
87        // |x| > 1 and |x| < 2^53
88        let m = (ax & ((1u64 << 52) - 1)) | (1u64 << 52); // mantissa with hidden 1
89        let shift = 52 - e;
90
91        let frac = m & ((1u64 << shift) - 1);
92        if frac == (1u64 << (shift - 1)) {
93            // |x| is always integer.5 means it's inf
94            return f64::copysign(0., x);
95        }
96    }
97
98    if ax <= 0x3cb0000000000000 {
99        // for tiny x ( |x| < f64::EPSILON ) just small taylor expansion
100        // cot(PI*x) ~ 1/(PI*x) + O(x^3)
101        const ONE_OVER_PI: DoubleDouble =
102            DoubleDouble::from_bit_pair((0xbc76b01ec5417056, 0x3fd45f306dc9c883));
103        if ax <= 0x3ca0000000000000 {
104            // |x| <= 2^-53, renormalize value
105            let e: i32 = (ax >> 52) as i32;
106            let sc = f64::from_bits((2045i64 - e as i64).wrapping_shl(52) as u64);
107            let dx = x * sc;
108            let q0 = backend.quick_mult(ONE_OVER_PI, DoubleDouble::from_quick_recip(dx));
109            let r = q0.to_f64() * sc;
110            return r;
111        }
112        let q0 = backend.quick_mult(ONE_OVER_PI, DoubleDouble::from_quick_recip(x));
113        let r = q0.to_f64();
114        return r;
115    }
116
117    // argument reduction
118    let (y, k) = backend.arg_reduce_pi_64(x);
119
120    if y == 0.0 {
121        let km = (k.abs() & 63) as i32; // k mod 64
122
123        match km {
124            0 => return f64::copysign(f64::INFINITY, x), // cotpi(n) = 0
125            32 => return f64::copysign(0., x),           // cotpi(n+0.5) = ±∞
126            16 => return f64::copysign(1.0, x),          // cotpi(n+0.25) = ±1
127            48 => return -f64::copysign(1.0, x),         // cotpi(n+0.75) = ∓1
128            _ => {}
129        }
130    }
131
132    let tan_k = DoubleDouble::from_bit_pair(TANPI_K_PI_OVER_64[((k as u64) & 127) as usize]);
133
134    // Computes tan(pi*x) through identities
135    // 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))
136    let tan_y = tanpi_eval(y, &backend);
137    // num = tan(y*pi/64) + tan(k*pi/64)
138    let num = DoubleDouble::add(tan_k, tan_y);
139    // den = 1 - tan(y*pi/64)*tan(k*pi/64)
140    let den = backend.mul_add_f64(tan_y, -tan_k, 1.);
141    // cot = den / num
142    let tan_value = backend.div(den, num);
143
144    let err = backend.fma(
145        tan_value.hi,
146        f64::from_bits(0x3bf0000000000000), // 2^-64
147        f64::from_bits(0x3b60000000000000), // 2^-73
148    );
149    let ub = tan_value.hi + (tan_value.lo + err);
150    let lb = tan_value.hi + (tan_value.lo - err);
151    if ub == lb {
152        return tan_value.to_f64();
153    }
154    cotpi_hard(y, tan_k).to_f64()
155}
156
157#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
158#[target_feature(enable = "avx", enable = "fma")]
159unsafe fn cotpi_fma_impl(x: f64) -> f64 {
160    use crate::sincospi::FmaSinCosPiBackend;
161    cotpi_gen_impl(x, FmaSinCosPiBackend {})
162}
163
164/// Computes cotangent 1/tan(PI*x)
165///
166/// ulp 0.5
167pub fn f_cotpi(x: f64) -> f64 {
168    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
169    {
170        cotpi_gen_impl(x, GenSinCosPiBackend {})
171    }
172    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
173    {
174        use std::sync::OnceLock;
175        static EXECUTOR: OnceLock<unsafe fn(f64) -> f64> = OnceLock::new();
176        let q = EXECUTOR.get_or_init(|| {
177            if std::arch::is_x86_feature_detected!("avx")
178                && std::arch::is_x86_feature_detected!("fma")
179            {
180                cotpi_fma_impl
181            } else {
182                fn def_cotpi(x: f64) -> f64 {
183                    cotpi_gen_impl(x, GenSinCosPiBackend {})
184                }
185                def_cotpi
186            }
187        });
188        unsafe { q(x) }
189    }
190}
191
192#[inline]
193pub(crate) fn cotpi_core(x: f64) -> DoubleDouble {
194    // argument reduction
195    let (y, k) = reduce_pi_64(x);
196
197    let tan_k = DoubleDouble::from_bit_pair(TANPI_K_PI_OVER_64[((k as u64) & 127) as usize]);
198
199    // Computes tan(pi*x) through identities.
200    // 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))
201    let tan_y = tanpi_eval(y, &GenSinCosPiBackend {});
202    // num = tan(y*pi/64) + tan(k*pi/64)
203    let num = DoubleDouble::add(tan_k, tan_y);
204    // den = 1 - tan(y*pi/64)*tan(k*pi/64)
205    let den = DoubleDouble::mul_add_f64(tan_y, -tan_k, 1.);
206    // cot = den / num
207    DoubleDouble::div(den, num)
208}
209
210#[cfg(test)]
211mod tests {
212    use super::*;
213
214    #[test]
215    fn test_cotpi() {
216        assert_eq!(f_cotpi(3.382112265043898e-306), 9.411570676518013e304);
217        assert_eq!(f_cotpi(0.0431431231), 7.332763436038805);
218        assert_eq!(f_cotpi(-0.0431431231), -7.332763436038805);
219        assert_eq!(f_cotpi(0.52324), -0.07314061937774036);
220        assert!(f_cotpi(f64::INFINITY).is_nan());
221        assert!(f_cotpi(f64::NAN).is_nan());
222        assert!(f_cotpi(f64::NEG_INFINITY).is_nan());
223    }
224}
pxfm/tangent/cotpi.rs

pxfm/tangent/
cotpi.rs
