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

pxfm/tangent/
cotpi.rs
