pxfm 0.1.28

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

#[cold]
#[inline(never)]
fn cosm1_accurate(y: DoubleDouble, sin_k: DoubleDouble, cos_k: DoubleDouble) -> f64 {
    let r_sincos = sincos_eval_dd(y);

    // 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(r_sincos.v_cos, sin_k);
    let cos_k_sin_y = DoubleDouble::quick_mult(r_sincos.v_sin, cos_k);

    let mut rr = DoubleDouble::full_dd_add(sin_k_cos_y, cos_k_sin_y);

    // Computing cos(x) - 1 as follows:
    // cos(x) - 1 = -2*sin^2(x/2)
    rr = DoubleDouble::from_exact_add(rr.hi, rr.lo);
    rr = DoubleDouble::quick_mult(rr, rr);
    rr = DoubleDouble::quick_mult_f64(rr, -2.);

    rr.to_f64()
}

#[cold]
fn cosm1_tiny_hard(x: f64) -> f64 {
    // Generated by Sollya:
    // d = [2^-27, 2^-7];
    // f_cosm1 = cos(x) - 1;
    // Q = fpminimax(f_cosm1, [|2,4,6,8|], [|0, 107...|], d);
    // See ./notes/cosm1_hard.sollya
    const C: [(u64, u64); 3] = [
        (0x3c453997dc8ae20d, 0x3fa5555555555555),
        (0x3bf6100c76a1827a, 0xbf56c16c16c15749),
        (0x3b918f45acdd1fb2, 0x3efa019ddf5a583a),
    ];
    let x2 = DoubleDouble::from_exact_mult(x, x);
    let mut p = DoubleDouble::mul_add(
        x2,
        DoubleDouble::from_bit_pair(C[2]),
        DoubleDouble::from_bit_pair(C[1]),
    );
    p = DoubleDouble::mul_add(x2, p, DoubleDouble::from_bit_pair(C[0]));
    p = DoubleDouble::mul_add_f64(x2, p, f64::from_bits(0xbfe0000000000000));
    p = DoubleDouble::quick_mult(p, x2);
    p.to_f64()
}

/// Computes cos(x) - 1
pub fn f_cosm1(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();

    // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA)
    if x_e < E_BIAS + 16 {
        // |x| < 2^-7
        if x_e < E_BIAS - 7 {
            // |x| < 2^-26
            if x_e < E_BIAS - 27 {
                // Signed zeros.
                if x == 0.0 {
                    return 0.0;
                }
                // Taylor expansion for small cos(x) - 1 ~ -x^2/2 + x^4/24 + O(x^6)
                let x_sqr = x * x;
                const A0: f64 = -1. / 2.;
                const A1: f64 = 1. / 24.;
                let r0 = f_fmla(x_sqr, A1, A0);
                return r0 * x_sqr;
            }

            // Generated by Sollya:
            // d = [2^-27, 2^-7];
            // f_cosm1 = (cos(x) - 1);
            // Q = fpminimax(f_cosm1, [|2,4,6,8|], [|0, D...|], d);
            // See ./notes/cosm1.sollya

            let x2 = DoubleDouble::from_exact_mult(x, x);
            let p = f_polyeval4(
                x2.hi,
                f64::from_bits(0xbfe0000000000000),
                f64::from_bits(0x3fa5555555555555),
                f64::from_bits(0xbf56c16c16b9c2b7),
                f64::from_bits(0x3efa014d03f38855),
            );

            let r = DoubleDouble::quick_mult_f64(x2, p);

            let eps = x * f_fmla(
                x2.hi,
                f64::from_bits(0x3d00000000000000), // 2^-47
                f64::from_bits(0x3be0000000000000), // 2^-65
            );

            let ub = r.hi + (r.lo + eps);
            let lb = r.hi + (r.lo - eps);
            if ub == lb {
                return r.to_f64();
            }
            return cosm1_tiny_hard(x);
        } else {
            // // Small range reduction.
            (y, k) = range_reduction_small(x * 0.5);
        }
    } else {
        // Inf or NaN
        if x_e > 2 * E_BIAS {
            // cos(+-Inf) = NaN
            return x + f64::NAN;
        }

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

    // Computing cos(x) - 1 as follows:
    // cos(x) - 1 = -2*sin^2(x/2)

    let r_sincos = sincos_eval(y);

    // 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 sin_k_cos_y = DoubleDouble::quick_mult(r_sincos.v_cos, sin_k);
    let cos_k_sin_y = DoubleDouble::quick_mult(r_sincos.v_sin, cos_k);

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

    rr = DoubleDouble::from_exact_add(rr.hi, rr.lo);
    rr = DoubleDouble::quick_mult(rr, rr);
    rr = DoubleDouble::quick_mult_f64(rr, -2.);

    let rlp = rr.lo + r_sincos.err;
    let rlm = rr.lo - r_sincos.err;

    let r_upper = rr.hi + rlp; // (rr.lo + ERR);
    let r_lower = rr.hi + rlm; // (rr.lo - ERR);

    // Ziv's accuracy test
    if r_upper == r_lower {
        return rr.to_f64();
    }

    cosm1_accurate(y, sin_k, cos_k)
}

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

    #[test]
    fn f_cosm1f_test() {
        assert_eq!(f_cosm1(0.0017700195313803402), -0.000001566484161754997);
        assert_eq!(
            f_cosm1(0.0000000011641532182693484),
            -0.0000000000000000006776263578034406
        );
        assert_eq!(f_cosm1(0.006164513528517324), -0.000019000553351160402);
        assert_eq!(f_cosm1(6.2831853071795862), -2.999519565323715e-32);
        assert_eq!(f_cosm1(0.00015928394), -1.2685686744140693e-8);
        assert_eq!(f_cosm1(0.0), 0.0);
        assert_eq!(f_cosm1(0.0), 0.0);
        assert_eq!(f_cosm1(std::f64::consts::PI), -2.);
        assert_eq!(f_cosm1(0.5), -0.12241743810962728);
        assert_eq!(f_cosm1(0.7), -0.23515781271551153);
        assert_eq!(f_cosm1(1.7), -1.1288444942955247);
        assert!(f_cosm1(f64::INFINITY).is_nan());
        assert!(f_cosm1(f64::NEG_INFINITY).is_nan());
        assert!(f_cosm1(f64::NAN).is_nan());
        assert_eq!(f_cosm1(0.0002480338), -3.0760382813519806e-8);
    }
}