pxfm/tangent/
cotpif.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 */
29
30use crate::common::f_fmla;
31use crate::sin_cosf::ArgumentReducerPi;
32use crate::tangent::evalf::tanpif_eval;
33
34/// Computes 1/tan(PI*x)
35///
36/// Max found ULP 0.5
37#[inline]
38pub fn f_cotpif(x: f32) -> f32 {
39    let ix = x.to_bits();
40    let e = ix & (0xff << 23);
41    if e > (150 << 23) {
42        // |x| > 2^23
43        if e == (0xff << 23) {
44            // x = nan or inf
45            if (ix.wrapping_shl(9)) == 0 {
46                // x = inf
47                return f32::NAN;
48            }
49            return x + x; // x = nan
50        }
51        return f32::INFINITY;
52    }
53
54    let argument_reduction = ArgumentReducerPi { x: x as f64 };
55
56    let (y, k) = argument_reduction.reduce();
57
58    if y == 0.0 {
59        let km = (k.abs() & 31) as i32; // k mod 32
60
61        match km {
62            0 => return f32::copysign(f32::INFINITY, x), // cotpi(n) = ∞
63            16 => return 0.0f32.copysign(x),             // cotpi(n+0.5) = 0
64            8 => return f32::copysign(1.0, x),           // cotpi(n+0.25) = 1
65            24 => return -f32::copysign(1.0, x),         // cotpi(n+0.75) = -1
66            _ => {}
67        }
68    }
69
70    let ax = ix & 0x7fff_ffff;
71    if ax < 0x3bc49ba6u32 {
72        // taylor series for cot(PI*x) where |x| < 0.006
73        let dx = x as f64;
74        let dx_sqr = dx * dx;
75        // cot(PI*x) ~ 1/(PI*x) - PI*x/3 - PI^3*x^3/45 + O(x^5)
76        const ONE_OVER_PI: f64 = f64::from_bits(0x3fd45f306dc9c883);
77        let r = f_fmla(
78            dx_sqr,
79            f64::from_bits(0xbfe60c8539c1dc14),
80            f64::from_bits(0xbff0c152382d7366),
81        );
82        let rcp = 1. / dx;
83        return f_fmla(rcp, ONE_OVER_PI, r * dx) as f32;
84    }
85
86    // tanpif_eval returns:
87    // - rs.tan_y = tan(pi/32 * y)          -> tangent of the remainder
88    // - rs.tan_k = tan(pi/32 * k)          -> tan of the main angle multiple
89    let rs = tanpif_eval(y, k);
90
91    // Then computing tan through identities
92    // num = tan(k*pi/32) + tan(y*pi/32)
93    let num = rs.tan_y + rs.tan_k;
94    // den = 1 - tan(k*pi/32) * tan(y*pi/32)
95    let den = f_fmla(rs.tan_y, -rs.tan_k, 1.);
96    // cotangent is tangent in inverse order
97    (den / num) as f32
98}
99
100#[inline]
101pub(crate) fn cotpif_core(x: f64) -> f64 {
102    let argument_reduction = ArgumentReducerPi { x };
103
104    let (y, k) = argument_reduction.reduce();
105
106    // tanpif_eval returns:
107    // - rs.tan_y = tan(pi/32 * y)          -> tangent of the remainder
108    // - rs.tan_k = tan(pi/32 * k)          -> tan of the main angle multiple
109    let rs = tanpif_eval(y, k);
110
111    // Then computing tan through identities
112    // num = tan(k*pi/32) + tan(y*pi/32)
113    let num = rs.tan_y + rs.tan_k;
114    // den = 1 - tan(k*pi/32) * tan(y*pi/32)
115    let den = f_fmla(rs.tan_y, -rs.tan_k, 1.);
116    // cotangent is tangent in inverse order
117    den / num
118}
119
120#[cfg(test)]
121mod tests {
122    use super::*;
123
124    #[test]
125    fn test_cotpif() {
126        assert_eq!(f_cotpif(0.00046277765), 687.82416);
127        assert_eq!(f_cotpif(2.3588752e-6), 134941.39);
128        assert_eq!(f_cotpif(10775313000000000000000000000000.), f32::INFINITY);
129        assert_eq!(f_cotpif(5.5625), -0.19891237);
130        assert_eq!(f_cotpif(-29.75), 1.0);
131        assert_eq!(f_cotpif(-21.5625), 0.19891237);
132        assert_eq!(f_cotpif(-15.611655), 0.3659073);
133        assert_eq!(f_cotpif(115.30706), 0.693186);
134        assert_eq!(f_cotpif(0.), f32::INFINITY);
135        assert!(f_cotpif(f32::INFINITY).is_nan());
136        assert!(f_cotpif(f32::NAN).is_nan());
137    }
138}
pxfm/tangent/cotpif.rs

pxfm/tangent/
cotpif.rs
