pxfm/tangent/

tanpif.rs

/*
 * // Copyright (c) Radzivon Bartoshyk 6/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;

#[inline(always)]
fn tanpif_gen_impl(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::copysign(0.0, x);
    }

    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 0.0f32.copysign(x),               // tanpi(n) = 0
            16 => return f32::copysign(f32::INFINITY, x), // tanpi(n+0.5) = ±∞
            8 => return f32::copysign(1.0, x),            // tanpi(n+0.25) = ±1
            24 => return -f32::copysign(1.0, x),          // tanpi(n+0.75) = ∓1
            _ => {}
        }
    }

    let ax = ix & 0x7fff_ffff;
    if ax < 0x38d1b717u32 {
        // taylor series for tan(PI*x) where |x| < 0.0001
        let dx = x as f64;
        let dx_sqr = dx * dx;
        // tan(PI*x) ~ PI*x + PI^3*x^3/3 + O(x^5)
        let r = f_fmla(
            dx_sqr,
            f64::from_bits(0x4024abbce625be53),
            f64::from_bits(0x400921fb54442d18),
        );
        return (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.);
    (num / den) as f32
}

#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
#[target_feature(enable = "avx", enable = "fma")]
unsafe fn tanpif_fma_impl(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::copysign(0.0, x);
    }

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

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

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

        match km {
            0 => return 0.0f32.copysign(x),               // tanpi(n) = 0
            16 => return f32::copysign(f32::INFINITY, x), // tanpi(n+0.5) = ±∞
            8 => return f32::copysign(1.0, x),            // tanpi(n+0.25) = ±1
            24 => return -f32::copysign(1.0, x),          // tanpi(n+0.75) = ∓1
            _ => {}
        }
    }

    let ax = ix & 0x7fff_ffff;
    if ax < 0x38d1b717u32 {
        // taylor series for tan(PI*x) where |x| < 0.0001
        let dx = x as f64;
        let dx_sqr = dx * dx;
        // tan(PI*x) ~ PI*x + PI^3*x^3/3 + O(x^5)
        let r = f64::mul_add(
            dx_sqr,
            f64::from_bits(0x4024abbce625be53),
            f64::from_bits(0x400921fb54442d18),
        );
        return (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
    use crate::tangent::evalf::tanpif_eval_fma;
    let rs = tanpif_eval_fma(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 = f64::mul_add(rs.tan_y, -rs.tan_k, 1.);
    (num / den) as f32
}

/// Computes tan(PI*x)
///
/// Max found ULP 0.5
#[inline]
pub fn f_tanpif(x: f32) -> f32 {
    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
    {
        tanpif_gen_impl(x)
    }
    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
    {
        use std::sync::OnceLock;
        static EXECUTOR: OnceLock<unsafe fn(f32) -> f32> = OnceLock::new();
        let q = EXECUTOR.get_or_init(|| {
            if std::arch::is_x86_feature_detected!("avx")
                && std::arch::is_x86_feature_detected!("fma")
            {
                tanpif_fma_impl
            } else {
                tanpif_gen_impl
            }
        });
        unsafe { q(x) }
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_tanpif() {
        assert_eq!(f_tanpif(3.666738e-5), 0.00011519398);
        assert_eq!(f_tanpif(1.0355987e-25), 3.2534293e-25);
        assert_eq!(f_tanpif(5.5625), -5.0273395);
        assert_eq!(f_tanpif(-29.75), 1.0);
        assert_eq!(f_tanpif(-21.5625), 5.0273395);
        assert_eq!(f_tanpif(-15.611655), 2.7329326);
        assert_eq!(f_tanpif(115.30706), 1.4426143);
        assert!(f_tanpif(f32::INFINITY).is_nan());
        assert!(f_tanpif(f32::NAN).is_nan());
    }
}