/*

 * // 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::exponents::exp2f::EXP2F_TABLE;

use crate::polyeval::f_polyeval3;

use crate::rounding::CpuRound;

/// Computes 2^x - 1

///

/// Max found ULP 0.5

#[inline]

pub fn f_exp2m1f(x: f32) -> f32 {

    let x_u = x.to_bits();

    let x_abs = x_u & 0x7fff_ffffu32;

    if x_abs >= 0x4300_0000u32 || x_abs <= 0x3d00_0000u32 {

        // |x| <= 2^-5

        if x_abs <= 0x3d00_0000u32 {

            // Minimax polynomial generated by Sollya with:

            // > display = hexadecimal;

            // > fpminimax((2^x - 1)/x, 5, [|D...|], [-2^-5, 2^-5]);

            const C: [u64; 6] = [

                0x3fe62e42fefa39f3,

                0x3fcebfbdff82c57b,

                0x3fac6b08d6f2d7aa,

                0x3f83b2ab6fc92f5d,

                0x3f55d897cfe27125,

                0x3f243090e61e6af1,

            ];

            let xd = x as f64;

            let xsq = xd * xd;

            let c0 = f_fmla(xd, f64::from_bits(C[1]), f64::from_bits(C[0]));

            let c1 = f_fmla(xd, f64::from_bits(C[3]), f64::from_bits(C[2]));

            let c2 = f_fmla(xd, f64::from_bits(C[5]), f64::from_bits(C[4]));

            let p = f_polyeval3(xsq, c0, c1, c2);

            return (p * xd) as f32;

        }

        // x >= 128, or x is nan

        if x.is_sign_positive() {

            // x >= 128 and 2^x - 1 rounds to +inf, or x is +inf or nan

            return x + f32::INFINITY;

        }

    }

    if x <= -25.0 {

        // 2^(-inf) - 1 = -1

        if x.is_infinite() {

            return -1.0;

        }

        // 2^nan - 1 = nan

        if x.is_nan() {

            return x;

        }

        return -1.0;

    }

    // For -25 < x < 128, to compute 2^x, we perform the following range

    // reduction: find hi, mid, lo such that:

    //   x = hi + mid + lo, in which:

    //     hi is an integer,

    //     0 <= mid * 2^5 < 32 is an integer,

    //     -2^(-6) <= lo <= 2^(-6).

    // In particular,

    //   hi + mid = round(x * 2^5) * 2^(-5).

    // Then,

    //   2^x = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo.

    // 2^mid is stored in the lookup table of 32 elements.

    // 2^lo is computed using a degree-4 minimax polynomial generated by Sollya.

    // We perform 2^hi * 2^mid by simply add hi to the exponent field of 2^mid.

    // kf = (hi + mid) * 2^5 = round(x * 2^5)

    let xd = x as f64;

    let kf = (x * 64.0).cpu_round();

    let k = unsafe { kf.to_int_unchecked::<i32>() }; // it's already not indeterminate.

    // dx = lo = x - (hi + mid) = x - kf * 2^(-6)

    let dx = f_fmla(f64::from_bits(0xbf90000000000000), kf as f64, xd);

    const TABLE_BITS: u32 = 6;

    const TABLE_MASK: u64 = (1u64 << TABLE_BITS) - 1;

    // hi = floor(kf * 2^(-5))

    // exp_hi = shift hi to the exponent field of double precision.

    let exp_hi: i64 = ((k >> TABLE_BITS) as i64).wrapping_shl(52);

    // mh = 2^hi * 2^mid

    // mh_bits = bit field of mh

    let mh_bits = (EXP2F_TABLE[((k as u64) & TABLE_MASK) as usize] as i64).wrapping_add(exp_hi);

    let mh = f64::from_bits(mh_bits as u64);

    // Degree-4 polynomial approximating (2^x - 1)/x generated by Sollya with:

    // > P = fpminimax((2^y - 1)/y, 4, [|D...|], [-1/64. 1/64]);

    // see ./notes/exp2f.sollya

    const C: [u64; 5] = [

        0x3fe62e42fefa39ef,

        0x3fcebfbdff8131c4,

        0x3fac6b08d7061695,

        0x3f83b2b1bee74b2a,

        0x3f55d88091198529,

    ];

    let dx_sq = dx * dx;

    let c1 = f_fmla(dx, f64::from_bits(C[0]), 1.0);

    let c2 = f_fmla(dx, f64::from_bits(C[2]), f64::from_bits(C[1]));

    let c3 = f_fmla(dx, f64::from_bits(C[4]), f64::from_bits(C[3]));

    let p = f_polyeval3(dx_sq, c1, c2, c3);

    // 2^x = 2^(hi + mid + lo)

    //     = 2^(hi + mid) * 2^lo

    //     ~ mh * (1 + lo * P(lo))

    //     = mh + (mh*lo) * P(lo)

    f_fmla(p, mh, -1.) as f32

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_exp2m1f() {

        assert_eq!(f_exp2m1f(0.432423), 0.34949815);

        assert_eq!(f_exp2m1f(-4.), -0.9375);

        assert_eq!(f_exp2m1f(5.43122), 42.14795);

        assert_eq!(f_exp2m1f(4.), 15.0);

        assert_eq!(f_exp2m1f(3.), 7.);

        assert_eq!(f_exp2m1f(0.1), 0.07177346);

        assert_eq!(f_exp2m1f(0.0543432432), 0.038386293);

        assert!(f_exp2m1f(f32::NAN).is_nan());

        assert_eq!(f_exp2m1f(f32::INFINITY), f32::INFINITY);

        assert_eq!(f_exp2m1f(f32::NEG_INFINITY), -1.0);

    }

}