pxfm 0.1.29

/*
 * // 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::bits::EXP_MASK;
use crate::common::{dyad_fmla, f_fmla};
use crate::double_double::DoubleDouble;
use crate::sin::range_reduction_small;
use crate::sincos_reduce::LargeArgumentReduction;
use crate::tangent::tanpi_table::TAN_K_PI_OVER_128;

#[inline]
pub(crate) fn tan_eval(u: DoubleDouble) -> (DoubleDouble, f64) {
    // Evaluate tan(y) = tan(x - k * (pi/128))
    // We use the degree-9 Taylor approximation:
    //   tan(y) ~ P(y) = y + y^3/3 + 2*y^5/15 + 17*y^7/315 + 62*y^9/2835
    // Then the error is bounded by:
    //   |tan(y) - P(y)| < 2^-6 * |y|^11 < 2^-6 * 2^-66 = 2^-72.
    // For y ~ u_hi + u_lo, fully expanding the polynomial and drop any terms
    // < ulp(u_hi^3) gives us:
    //   P(y) = y + y^3/3 + 2*y^5/15 + 17*y^7/315 + 62*y^9/2835 = ...
    // ~ u_hi + u_hi^3 * (1/3 + u_hi^2 * (2/15 + u_hi^2 * (17/315 +
    //                                                     + u_hi^2 * 62/2835))) +
    //        + u_lo (1 + u_hi^2 * (1 + u_hi^2 * 2/3))
    let u_hi_sq = u.hi * u.hi; // Error < ulp(u_hi^2) < 2^(-6 - 52) = 2^-58.
    // p1 ~ 17/315 + u_hi^2 62 / 2835.
    let p1 = f_fmla(
        u_hi_sq,
        f64::from_bits(0x3f9664f4882c10fa),
        f64::from_bits(0x3faba1ba1ba1ba1c),
    );
    // p2 ~ 1/3 + u_hi^2 2 / 15.
    let p2 = f_fmla(
        u_hi_sq,
        f64::from_bits(0x3fc1111111111111),
        f64::from_bits(0x3fd5555555555555),
    );
    // q1 ~ 1 + u_hi^2 * 2/3.
    let q1 = f_fmla(u_hi_sq, f64::from_bits(0x3fe5555555555555), 1.0);
    let u_hi_3 = u_hi_sq * u.hi;
    let u_hi_4 = u_hi_sq * u_hi_sq;
    // p3 ~ 1/3 + u_hi^2 * (2/15 + u_hi^2 * (17/315 + u_hi^2 * 62/2835))
    let p3 = f_fmla(u_hi_4, p1, p2);
    // q2 ~ 1 + u_hi^2 * (1 + u_hi^2 * 2/3)
    let q2 = f_fmla(u_hi_sq, q1, 1.0);
    let tan_lo = f_fmla(u_hi_3, p3, u.lo * q2);
    // Overall, |tan(y) - (u_hi + tan_lo)| < ulp(u_hi^3) <= 2^-71.
    // And the relative errors is:
    // |(tan(y) - (u_hi + tan_lo)) / tan(y) | <= 2*ulp(u_hi^2) < 2^-64
    let err = f_fmla(
        u_hi_3.abs(),
        f64::from_bits(0x3cc0000000000000),
        f64::from_bits(0x3990000000000000),
    );
    (DoubleDouble::from_exact_add(u.hi, tan_lo), err)
}

#[inline]
pub(crate) fn tan_eval_dd(x: DoubleDouble) -> DoubleDouble {
    // Generated by Sollya:
    // d = [0, pi/128];
    // f_tan = tan(x)/x;
    // Q = fpminimax(f_tan, [|0, 2, 4, 6, 8, 10, 12|], [|107...|], d);
    const C: [(u64, u64); 7] = [
        (0x0000000000000000, 0x3ff0000000000000),
        (0x3c75555549d35a4b, 0x3fd5555555555555),
        (0x3c413dd6ea580288, 0x3fc1111111111111),
        (0x3c4e100b946f0c89, 0x3faba1ba1ba1b9fe),
        (0xbc3c180dfac2b8c3, 0x3f9664f488307189),
        (0x3c1fd691c256a14a, 0x3f8226e300c1988e),
        (0x3bec7b08c90fdfc0, 0x3f6d739baeacd204),
    ];
    let x2 = DoubleDouble::quick_mult(x, x);
    let mut p = DoubleDouble::quick_mul_add(
        x2,
        DoubleDouble::from_bit_pair(C[6]),
        DoubleDouble::from_bit_pair(C[5]),
    );
    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[4]));
    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[3]));
    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[2]));
    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[1]));
    p = DoubleDouble::quick_mul_add_f64(x2, p, f64::from_bits(0x3ff0000000000000));
    let r = DoubleDouble::quick_mult(p, x);
    DoubleDouble::from_exact_add(r.hi, r.lo)
}

#[cold]
fn tan_accurate(y: DoubleDouble, tan_k: DoubleDouble) -> f64 {
    // Computes tan(x) through identities
    // tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a)tan(b)) = (tan(y) + tan(k*pi/128)) / (1 - tan(y)*tan(k*pi/128))
    let tan_y = tan_eval_dd(y);

    // num = tan(y) + tan(k*pi/64)
    let num_dd = DoubleDouble::full_dd_add(tan_y, tan_k);
    // den = 1 - tan(y)*tan(k*pi/64)
    let den_dd = DoubleDouble::mul_add_f64(tan_y, -tan_k, 1.);

    let tan_x = DoubleDouble::div(num_dd, den_dd);
    tan_x.to_f64()
}

#[cold]
fn tan_near_zero_hard(x: f64) -> DoubleDouble {
    // Generated by Sollya:
    // d = [0, 0.03125];
    // f_tan = tan(x)/x;
    // Q = fpminimax(f_tan, [|0, 2, 4, 6, 8, 10, 12|], [|1, 107...|], d);
    const C: [(u64, u64); 7] = [
        (0x0000000000000000, 0x3ff0000000000000),
        (0x3c7555548455da94, 0x3fd5555555555555),
        (0x3c4306a6358cc810, 0x3fc1111111111111),
        (0x3c2ca9cd025ea98c, 0x3faba1ba1ba1b952),
        (0x3c3cb012803c55a5, 0x3f9664f4883eb962),
        (0x3c28afc1adb8f202, 0x3f8226e276097461),
        (0xbbdf8f911392f348, 0x3f6d7791601277ea),
    ];
    let x2 = DoubleDouble::from_exact_mult(x, x);
    let mut p = DoubleDouble::quick_mul_add(
        x2,
        DoubleDouble::from_bit_pair(C[6]),
        DoubleDouble::from_bit_pair(C[5]),
    );
    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[4]));
    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[3]));
    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[2]));
    p = DoubleDouble::quick_mul_add(x2, p, DoubleDouble::from_bit_pair(C[1]));
    p = DoubleDouble::quick_mul_add_f64(x2, p, f64::from_bits(0x3ff0000000000000));
    DoubleDouble::quick_mult_f64(p, x)
}

/// Tangent in double precision
///
/// ULP 0.5
pub fn f_tan(x: f64) -> f64 {
    let x_e = (x.to_bits() >> 52) & 0x7ff;
    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;

    let y: DoubleDouble;
    let k;

    let mut argument_reduction = LargeArgumentReduction::default();

    if x_e < E_BIAS + 16 {
        // |x| < 2^16
        if x_e < E_BIAS - 5 {
            // |x| < 2^-5
            if x_e < E_BIAS - 27 {
                // |x| < 2^-27, |tan(x) - x| < ulp(x)/2.
                if x == 0.0 {
                    // Signed zeros.
                    return x + x;
                }
                return dyad_fmla(x, f64::from_bits(0x3c90000000000000), x);
            }
            let x2 = x * x;
            let x4 = x2 * x2;
            // Generated by Sollya:
            // d = [0, 0.03125];
            // f_tan = tan(x)/x;
            // Q = fpminimax(f_tan, [|0, 2, 4, 6, 8|], [|D...|], d);
            const C: [u64; 4] = [
                0x3fd555555555554b,
                0x3fc1111111142d5b,
                0x3faba1b9860ed843,
                0x3f966a802210f5bb,
            ];
            let p01 = f_fmla(x2, f64::from_bits(C[1]), f64::from_bits(C[0]));
            let p23 = f_fmla(x2, f64::from_bits(C[3]), f64::from_bits(C[2]));
            let w0 = f_fmla(x4, p23, p01);
            let w1 = x2 * w0 * x;
            let r = DoubleDouble::from_exact_add(x, w1);
            let err = f_fmla(
                x2,
                f64::from_bits(0x3cb0000000000000), // 2^-52
                f64::from_bits(0x3be0000000000000), // 2^-65
            );
            let ub = r.hi + (r.lo + err);
            let lb = r.hi + (r.lo - err);
            if ub == lb {
                return ub;
            }
            return tan_near_zero_hard(x).to_f64();
        } else {
            // Small range reduction.
            (y, k) = range_reduction_small(x);
        }
    } else {
        // Inf or NaN
        if x_e > 2 * E_BIAS {
            if x.is_nan() {
                return f64::NAN;
            }
            // tan(+-Inf) = NaN
            return x + f64::NAN;
        }

        // Large range reduction.
        (k, y) = argument_reduction.reduce(x);
    }

    let (tan_y, err) = tan_eval(y);

    // Computes tan(x) through identities.
    // tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a)tan(b)) = (tan(y) + tan(k*pi/128)) / (1 - tan(y)*tan(k*pi/128))
    let tan_k = DoubleDouble::from_bit_pair(TAN_K_PI_OVER_128[(k & 255) as usize]);

    // num = tan(y) + tan(k*pi/64)
    let num_dd = DoubleDouble::add(tan_y, tan_k);
    // den = 1 - tan(y)*tan(k*pi/64)
    let den_dd = DoubleDouble::mul_add_f64(tan_y, -tan_k, 1.);

    let tan_x = DoubleDouble::div(num_dd, den_dd);

    // Simple error bound: |1 / den_dd| < 2^(1 + floor(-log2(den_dd)))).
    let den_inv = ((E_BIAS + 1) << (52 + 1)) - (den_dd.hi.to_bits() & EXP_MASK);
    // For tan_x = (num_dd + err) / (den_dd + err), the error is bounded by:
    //   | tan_x - num_dd / den_dd |  <= err * ( 1 + | tan_x * den_dd | ).
    let tan_err = err * f_fmla(f64::from_bits(den_inv), tan_x.hi.abs(), 1.0);

    let err_higher = tan_x.lo + tan_err;
    let err_lower = tan_x.lo - tan_err;

    let tan_upper = tan_x.hi + err_higher;
    let tan_lower = tan_x.hi + err_lower;

    // Ziv_s rounding test.
    if tan_upper == tan_lower {
        return tan_upper;
    }

    tan_accurate(y, tan_k)
}

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

    #[test]
    fn tan_test() {
        assert_eq!(f_tan(0.0003242153424), 0.0003242153537600293);
        assert_eq!(f_tan(-0.3242153424), -0.3360742441422686);
        assert_eq!(f_tan(0.3242153424), 0.3360742441422686);
        assert_eq!(f_tan(-97301943219974110.), 0.000000000000000481917514979303);
        assert_eq!(f_tan(0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000007291122019556397),
                   0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000007291122019556397);
        assert_eq!(f_tan(0.0), 0.0);
        assert_eq!(f_tan(0.4324324324), 0.4615683710506729);
        assert_eq!(f_tan(1.0), 1.5574077246549023);
        assert_eq!(f_tan(-0.5), -0.5463024898437905);
        assert!(f_tan(f64::INFINITY).is_nan());
        assert!(f_tan(f64::NEG_INFINITY).is_nan());
        assert!(f_tan(f64::NAN).is_nan());
    }
}