/*

 * // Copyright (c) Radzivon Bartoshyk 8/2025. All rights reserved.

 * //

 * // Redistribution and use in source and binary forms, with or without modification,

 * // are permitted provided that the following conditions are met:

 * //

 * // 1.  Redistributions of source code must retain the above copyright notice, this

 * // list of conditions and the following disclaimer.

 * //

 * // 2.  Redistributions in binary form must reproduce the above copyright notice,

 * // this list of conditions and the following disclaimer in the documentation

 * // and/or other materials provided with the distribution.

 * //

 * // 3.  Neither the name of the copyright holder nor the names of its

 * // contributors may be used to endorse or promote products derived from

 * // this software without specific prior written permission.

 * //

 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE

 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE

 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL

 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR

 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER

 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,

 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE

 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 */

use crate::common::f_fmla;

use crate::sin_cosf::ArgumentReducerPi;

use crate::tangent::evalf::tanpif_eval;

/// Computes 1/tan(PI*x)

///

/// Max found ULP 0.5

#[inline]

pub fn f_cotpif(x: f32) -> f32 {

    let ix = x.to_bits();

    let e = ix & (0xff << 23);

    if e > (150 << 23) {

        // |x| > 2^23

        if e == (0xff << 23) {

            // x = nan or inf

            if (ix.wrapping_shl(9)) == 0 {

                // x = inf

                return f32::NAN;

            }

            return x + x; // x = nan

        }

        return f32::INFINITY;

    }

    let argument_reduction = ArgumentReducerPi { x: x as f64 };

    let (y, k) = argument_reduction.reduce();

    if y == 0.0 {

        let km = (k.abs() & 31) as i32; // k mod 32

        match km {

            0 => return f32::copysign(f32::INFINITY, x), // cotpi(n) = ∞

            16 => return 0.0f32.copysign(x),             // cotpi(n+0.5) = 0

            8 => return f32::copysign(1.0, x),           // cotpi(n+0.25) = 1

            24 => return -f32::copysign(1.0, x),         // cotpi(n+0.75) = -1

            _ => {}

        }

    }

    let ax = ix & 0x7fff_ffff;

    if ax < 0x3bc49ba6u32 {

        // taylor series for cot(PI*x) where |x| < 0.006

        let dx = x as f64;

        let dx_sqr = dx * dx;

        // cot(PI*x) ~ 1/(PI*x) - PI*x/3 - PI^3*x^3/45 + O(x^5)

        const ONE_OVER_PI: f64 = f64::from_bits(0x3fd45f306dc9c883);

        let r = f_fmla(

            dx_sqr,

            f64::from_bits(0xbfe60c8539c1dc14),

            f64::from_bits(0xbff0c152382d7366),

        );

        let rcp = 1. / dx;

        return f_fmla(rcp, ONE_OVER_PI, r * dx) as f32;

    }

    // tanpif_eval returns:

    // - rs.tan_y = tan(pi/32 * y)          -> tangent of the remainder

    // - rs.tan_k = tan(pi/32 * k)          -> tan of the main angle multiple

    let rs = tanpif_eval(y, k);

    // Then computing tan through identities

    // num = tan(k*pi/32) + tan(y*pi/32)

    let num = rs.tan_y + rs.tan_k;

    // den = 1 - tan(k*pi/32) * tan(y*pi/32)

    let den = f_fmla(rs.tan_y, -rs.tan_k, 1.);

    // cotangent is tangent in inverse order

    (den / num) as f32

}

#[inline]

pub(crate) fn cotpif_core(x: f64) -> f64 {

    let argument_reduction = ArgumentReducerPi { x };

    let (y, k) = argument_reduction.reduce();

    // tanpif_eval returns:

    // - rs.tan_y = tan(pi/32 * y)          -> tangent of the remainder

    // - rs.tan_k = tan(pi/32 * k)          -> tan of the main angle multiple

    let rs = tanpif_eval(y, k);

    // Then computing tan through identities

    // num = tan(k*pi/32) + tan(y*pi/32)

    let num = rs.tan_y + rs.tan_k;

    // den = 1 - tan(k*pi/32) * tan(y*pi/32)

    let den = f_fmla(rs.tan_y, -rs.tan_k, 1.);

    // cotangent is tangent in inverse order

    den / num

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_cotpif() {

        assert_eq!(f_cotpif(0.00046277765), 687.82416);

        assert_eq!(f_cotpif(2.3588752e-6), 134941.39);

        assert_eq!(f_cotpif(10775313000000000000000000000000.), f32::INFINITY);

        assert_eq!(f_cotpif(5.5625), -0.19891237);

        assert_eq!(f_cotpif(-29.75), 1.0);

        assert_eq!(f_cotpif(-21.5625), 0.19891237);

        assert_eq!(f_cotpif(-15.611655), 0.3659073);

        assert_eq!(f_cotpif(115.30706), 0.693186);

        assert_eq!(f_cotpif(0.), f32::INFINITY);

        assert!(f_cotpif(f32::INFINITY).is_nan());

        assert!(f_cotpif(f32::NAN).is_nan());

    }

}