/*

 * // 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);

    }

}