pxfm 0.1.28

/*
 * // 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::*;
use crate::double_double::DoubleDouble;
use crate::sin::{
    cos_accurate_near_zero, range_reduction_small, sin_accurate_near_zero, sincos_eval,
};
use crate::sin_helper::sincos_eval_dd;
use crate::sin_table::SIN_K_PI_OVER_128;
use crate::sincos_reduce::LargeArgumentReduction;
use std::hint::black_box;

/// Sine and cosine for double precision
///
/// ULP 0.5
pub fn f_sincos(x: f64) -> (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();

    // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA)
    if x_e < E_BIAS + 16 {
        // |x| < 2^-26
        let ax = x.to_bits() & 0x7fff_ffff_ffff_ffff;
        if ax <= 0x3fa921fbd34a9550 {
            // |x| <= 0.0490874
            if x_e < E_BIAS - 27 {
                // Signed zeros.
                if x == 0.0 {
                    return (x, 1.0);
                }
                // For |x| < 2^-26, |sin(x) - x| < ulp(x)/2.
                let s_sin = dyad_fmla(x, f64::from_bits(0xbc90000000000000), x);
                let s_cos = black_box(1.0) - min_normal_f64();
                return (s_sin, s_cos);
            }

            // Polynomial for sin(x)/x
            // Generated by Sollya:
            // d = [0, 0.0490874];
            // f_sin = sin(y)/y;
            // Q = fpminimax(f_sin, [|0, 2, 4, 6, 8|], [|1, D...|], d, relative, floating);
            const SIN_C: [u64; 4] = [
                0xbfc5555555555555,
                0x3f8111111110e45a,
                0xbf2a019ffd7fdaaf,
                0x3ec71819b9bf01ef,
            ];

            let x2 = x * x;
            let x4 = x2 * x2;

            let p01 = f_fmla(x2, f64::from_bits(SIN_C[1]), f64::from_bits(SIN_C[0]));
            let p23 = f_fmla(x2, f64::from_bits(SIN_C[3]), f64::from_bits(SIN_C[2]));
            let w0 = f_fmla(x4, p23, p01);
            let w1 = x2 * w0 * x;
            let r_sin = DoubleDouble::from_exact_add(x, w1);

            // Polynomial for cos(x)
            // Generated by Sollya:
            // d = [0, 0.0490874];
            // f_cos = cos(y);
            // Q = fpminimax(f_cos, [|0, 2, 4, 6, 8|], [|1, D...|], d, relative, floating);
            const COS_C: [u64; 4] = [
                0xbfe0000000000000,
                0x3fa55555555554a4,
                0xbf56c16c1619b84a,
                0x3efa013d3d01cf7f,
            ];
            let p01 = f_fmla(x2, f64::from_bits(COS_C[1]), f64::from_bits(COS_C[0]));
            let p23 = f_fmla(x2, f64::from_bits(COS_C[3]), f64::from_bits(COS_C[2]));
            let w0 = f_fmla(x4, p23, p01);
            let w1 = x2 * w0;
            let r_cos = DoubleDouble::from_exact_add(1., w1);

            let err = f_fmla(
                x2,
                f64::from_bits(0x3cb0000000000000), // 2^-52
                f64::from_bits(0x3be0000000000000), // 2^-65
            );

            let sin_ub = r_sin.hi + (r_sin.lo + err);
            let sin_lb = r_sin.hi + (r_sin.lo - err);
            let sin_x = if sin_ub == sin_lb {
                sin_ub
            } else {
                sin_accurate_near_zero(x)
            };

            let cos_ub = r_cos.hi + (r_cos.lo + err);
            let cos_lb = r_cos.hi + (r_cos.lo - err);
            let cos_x = if cos_ub == cos_lb {
                cos_ub
            } else {
                cos_accurate_near_zero(x)
            };
            return (sin_x, cos_x);
        } else {
            // // Small range reduction.
            (y, k) = range_reduction_small(x);
        }
    } else {
        // Inf or NaN
        if x_e > 2 * E_BIAS {
            // sin(+-Inf) = NaN
            return (x + f64::NAN, x + f64::NAN);
        }

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

    let r_sincos = sincos_eval(y);
    let (sin_y, cos_y) = (r_sincos.v_sin, r_sincos.v_cos);

    // Fast look up version, but needs 256-entry table.
    // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
    let sk = SIN_K_PI_OVER_128[(k & 255) as usize];
    let ck = SIN_K_PI_OVER_128[((k.wrapping_add(64)) & 255) as usize];
    let sin_k = DoubleDouble::from_bit_pair(sk);
    let cos_k = DoubleDouble::from_bit_pair(ck);

    let msin_k = -sin_k;

    // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128).
    // So k is an integer and -pi / 256 <= y <= pi / 256.
    // Then sin(x) = sin((k * pi/128 + y)
    //             = sin(y) * cos(k*pi/128) + cos(y) * sin(k*pi/128)
    let sin_k_cos_y = DoubleDouble::quick_mult(cos_y, sin_k);
    let cos_k_sin_y = DoubleDouble::quick_mult(sin_y, cos_k);
    //      cos(x) = cos((k * pi/128 + y)
    //             = cos(y) * cos(k*pi/128) - sin(y) * sin(k*pi/128)
    let cos_k_cos_y = DoubleDouble::quick_mult(cos_y, cos_k);
    let msin_k_sin_y = DoubleDouble::quick_mult(sin_y, msin_k);

    // cos_k_sin_y is always >> sin_k_cos_y
    let mut sin_dd = DoubleDouble::from_exact_add(sin_k_cos_y.hi, cos_k_sin_y.hi);
    // cos_k_cos_y is always >> msin_k_sin_y
    let mut cos_dd = DoubleDouble::from_exact_add(cos_k_cos_y.hi, msin_k_sin_y.hi);
    sin_dd.lo += sin_k_cos_y.lo + cos_k_sin_y.lo;
    cos_dd.lo += msin_k_sin_y.lo + cos_k_cos_y.lo;

    let sin_lp = sin_dd.lo + r_sincos.err;
    let sin_lm = sin_dd.lo - r_sincos.err;
    let cos_lp = cos_dd.lo + r_sincos.err;
    let cos_lm = cos_dd.lo - r_sincos.err;

    let sin_upper = sin_dd.hi + sin_lp;
    let sin_lower = sin_dd.hi + sin_lm;
    let cos_upper = cos_dd.hi + cos_lp;
    let cos_lower = cos_dd.hi + cos_lm;

    // Ziv's rounding test.
    if sin_upper == sin_lower && cos_upper == cos_lower {
        return (sin_upper, cos_upper);
    }

    sincos_hard(y, sin_k, cos_k, sin_upper, sin_lower, cos_upper, cos_lower)
}

#[cold]
#[inline(never)]
#[allow(clippy::too_many_arguments)]
fn sincos_hard(
    y: DoubleDouble,
    sin_k: DoubleDouble,
    cos_k: DoubleDouble,
    sin_upper: f64,
    sin_lower: f64,
    cos_upper: f64,
    cos_lower: f64,
) -> (f64, f64) {
    let r_sincos = sincos_eval_dd(y);

    let msin_k = -sin_k;

    let sin_x = if sin_upper == sin_lower {
        sin_upper
    } else {
        // sin(x) = sin((k * pi/128 + u)
        //        = sin(u) * cos(k*pi/128) + cos(u) * sin(k*pi/128)

        DoubleDouble::mul_add(sin_k, r_sincos.v_cos, cos_k * r_sincos.v_sin).to_f64()
    };

    let cos_x = if cos_upper == cos_lower {
        cos_upper
    } else {
        // cos(x) = cos((k * pi/128 + u)
        //        = cos(u) * cos(k*pi/128) - sin(u) * sin(k*pi/128)
        DoubleDouble::mul_add(cos_k, r_sincos.v_cos, msin_k * r_sincos.v_sin).to_f64()
    };
    (sin_x, cos_x)
}

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

    #[test]
    fn f_sincos_test() {
        let subnormal = f_sincos(0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000015708065354637772);
        assert_eq!(subnormal.0, 1.5708065354637772e-307);
        assert_eq!(subnormal.1, 1.0);
        let zx_0 = f_sincos(0.0);
        assert_eq!(zx_0.0, 0.0);
        assert_eq!(zx_0.1, 1.0);
        let zx_1 = f_sincos(1.0);
        assert_eq!(zx_1.0, 0.8414709848078965);
        assert_eq!(zx_1.1, 0.5403023058681398);
        let zx_0_p5 = f_sincos(-0.5);
        assert_eq!(zx_0_p5.0, -0.479425538604203);
        assert_eq!(zx_0_p5.1, 0.8775825618903728);

        let zx_1 = f_sincos(0.002341235432);
        assert_eq!(zx_1.0, 0.0023412332931324344);
        assert_eq!(zx_1.1, 0.9999972593095778);

        let zx_1 = f_sincos(0.0198676543432);
        assert_eq!(zx_1.0, 0.019866347330026367);
        assert_eq!(zx_1.1, 0.9998026446473137);
    }
}