/*

 * // Copyright (c) Radzivon Bartoshyk 7/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, f_fmlaf, pow2if};

use crate::polyeval::f_polyeval6;

use crate::rounding::CpuRound;

use std::hint::black_box;

const TBLSIZE: usize = 64;

#[repr(align(64))]

struct Exp2Table([(u32, u32); TBLSIZE]);

#[rustfmt::skip]

static EXP2FT: Exp2Table = Exp2Table([(0x3F3504F3, 0xB2D4175E),(0x3F36FD92, 0x3268D5EF),(0x3F38FBAF, 0xB30E8719),(0x3F3AFF5B, 0x3319E7DA),(0x3F3D08A4, 0x333CD82F),(0x3F3F179A, 0x330E1902),(0x3F412C4D, 0x32CCF4D7),(0x3F4346CD, 0x328F330E),(0x3F45672A, 0xB201B5B7),(0x3F478D75, 0x32CCCE34),(0x3F49B9BE, 0x335E937C),(0x3F4BEC15, 0x2FF41909),(0x3F4E248C, 0xB21760EA),(0x3F506334, 0x3283628B),(0x3F52A81E, 0x3340F500),(0x3F54F35B, 0x331202BD),(0x3F5744FD, 0x32B66A3E),(0x3F599D16, 0x32D0D9B1),(0x3F5BFBB8, 0x332ED93F),(0x3F5E60F5, 0x3350A709),(0x3F60CCDF, 0x32025744),(0x3F633F89, 0xB33A7C4D),(0x3F65B907, 0x321DA4E9),(0x3F68396A, 0xB2FF36A7),(0x3F6AC0C7, 0x3217E40E),(0x3F6D4F30, 0xB2400CBB),(0x3F6FE4BA, 0x331A2ACC),(0x3F728177, 0xB2B7D3E5),(0x3F75257D, 0xB1FED2BE),(0x3F77D0DF, 0xB32B73BA),(0x3F7A83B3, 0x32579081),(0x3F7D3E0C, 0xB19726B5),(0x3F800000, 0x00000000),(0x3F8164D2, 0x320C09FB),(0x3F82CD87, 0x3391E031),(0x3F843A29, 0x33287EEF),(0x3F85AAC3, 0xB38F6665),(0x3F871F62, 0x339004AB),(0x3F88980F, 0x33AC4561),(0x3F8A14D5, 0xB39CDAEA),(0x3F8B95C2, 0x32949D5C),(0x3F8D1ADF, 0xB36F79FA),(0x3F8EA43A, 0x33971DC2),(0x3F9031DC, 0xB32BD022),(0x3F91C3D3, 0xB3928952),(0x3F935A2B, 0xB2EBFECF),(0x3F94F4F0, 0x3357B8BB),(0x3F96942D, 0xB307353B),(0x3F9837F0, 0xB345DFE9),(0x3F99E046, 0x3382A804),(0x3F9B8D3A, 0x3326993E),(0x3F9D3EDA, 0x3350A029),(0x3F9EF532, 0xB3605F62),(0x3FA0B051, 0xB210909B),(0x3FA27043, 0xB0DDC369),(0x3FA43516, 0x33385844),(0x3FA5FED7, 0x33400757),(0x3FA7CD94, 0x3325446E),(0x3FA9A15B, 0x33237A50),(0x3FAB7A3A, 0x33201CA4),(0x3FAD583F, 0x32394687),(0x3FAF3B79, 0x332E1225),(0x3FB123F6, 0x33838969),(0x3FB311C4, 0xB219F2BA)]);

/**

Generated by SageMath:

```python

print("[")

for k in range(64):

    k = RealField(150)(2)**(RealField(150)(k) / RealField(150)(64))

    print(double_to_hex(k) + ",")

print("];")

```

**/

pub(crate) static EXP2F_TABLE: [u64; 64] = [

    0x3ff0000000000000,

    0x3ff02c9a3e778061,

    0x3ff059b0d3158574,

    0x3ff0874518759bc8,

    0x3ff0b5586cf9890f,

    0x3ff0e3ec32d3d1a2,

    0x3ff11301d0125b51,

    0x3ff1429aaea92de0,

    0x3ff172b83c7d517b,

    0x3ff1a35beb6fcb75,

    0x3ff1d4873168b9aa,

    0x3ff2063b88628cd6,

    0x3ff2387a6e756238,

    0x3ff26b4565e27cdd,

    0x3ff29e9df51fdee1,

    0x3ff2d285a6e4030b,

    0x3ff306fe0a31b715,

    0x3ff33c08b26416ff,

    0x3ff371a7373aa9cb,

    0x3ff3a7db34e59ff7,

    0x3ff3dea64c123422,

    0x3ff4160a21f72e2a,

    0x3ff44e086061892d,

    0x3ff486a2b5c13cd0,

    0x3ff4bfdad5362a27,

    0x3ff4f9b2769d2ca7,

    0x3ff5342b569d4f82,

    0x3ff56f4736b527da,

    0x3ff5ab07dd485429,

    0x3ff5e76f15ad2148,

    0x3ff6247eb03a5585,

    0x3ff6623882552225,

    0x3ff6a09e667f3bcd,

    0x3ff6dfb23c651a2f,

    0x3ff71f75e8ec5f74,

    0x3ff75feb564267c9,

    0x3ff7a11473eb0187,

    0x3ff7e2f336cf4e62,

    0x3ff82589994cce13,

    0x3ff868d99b4492ed,

    0x3ff8ace5422aa0db,

    0x3ff8f1ae99157736,

    0x3ff93737b0cdc5e5,

    0x3ff97d829fde4e50,

    0x3ff9c49182a3f090,

    0x3ffa0c667b5de565,

    0x3ffa5503b23e255d,

    0x3ffa9e6b5579fdbf,

    0x3ffae89f995ad3ad,

    0x3ffb33a2b84f15fb,

    0x3ffb7f76f2fb5e47,

    0x3ffbcc1e904bc1d2,

    0x3ffc199bdd85529c,

    0x3ffc67f12e57d14b,

    0x3ffcb720dcef9069,

    0x3ffd072d4a07897c,

    0x3ffd5818dcfba487,

    0x3ffda9e603db3285,

    0x3ffdfc97337b9b5f,

    0x3ffe502ee78b3ff6,

    0x3ffea4afa2a490da,

    0x3ffefa1bee615a27,

    0x3fff50765b6e4540,

    0x3fffa7c1819e90d8,

];

/* ULP 0.508 method

  let redux = f32::from_bits(0x4b400000) / TBLSIZE as f32;

  let ui = f32::to_bits(d + redux);

  let mut i0 = ui;

  i0 = i0.wrapping_add(TBLSIZE as u32 / 2);

  let k = i0 / TBLSIZE as u32;

  i0 &= TBLSIZE as u32 - 1;

  let mut uf = f32::from_bits(ui);

  uf -= redux;

  let item = EXP2FT.0[i0 as usize];

  let z0: f32 = f32::from_bits(item.0);

  let z1: f32 = f32::from_bits(item.1);

  let f: f32 = d - uf - z1;

  let mut u = 0.055504108664458832;

  u = f_fmlaf(u, f, 0.24022650695908768);

  u = f_fmlaf(u, f, 0.69314718055994973);

  u *= f;

  let i2 = pow2if(k as i32);

  f_fmlaf(u, z0, z0) * i2

*/

/// Computing exp2f

///

/// ULP 0.4999994

#[inline]

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

    let mut t = x.to_bits();

    if (t & 0xffff) == 0 {

        // x maybe integer

        let k: i32 = (((t >> 23) & 0xff) as i32).wrapping_sub(127); // 2^k <= |x| < 2^(k+1)

        if k >= 0 && k < 9 && (t << (9i32.wrapping_add(k))) == 0 {

            // x integer, with 1 <= |x| < 2^9

            let msk = (t as i32) >> 31;

            let mut m: i32 = (((t & 0x7fffff) | (1 << 23)) >> (23 - k)) as i32;

            m = (m ^ msk).wrapping_sub(msk).wrapping_add(127);

            if m > 0 && m < 255 {

                t = (m as u32).wrapping_shl(23);

                return f32::from_bits(t);

            } else if m <= 0 && m > -23 {

                t = 1i32.wrapping_shl(22i32.wrapping_add(m) as u32) as u32;

                return f32::from_bits(t);

            }

        }

    }

    let ux = t.wrapping_shl(1);

    if ux >= 0x86000000u32 || ux < 0x65000000u32 {

        // |x| >= 128 or x=nan or |x| < 0x1p-26

        if ux < 0x65000000u32 {

            return 1.0 + x;

        } // |x| < 0x1p-26

        // if x < -149 or 128 <= x is special

        if !(t >= 0xc3000000u32 && t < 0xc3150000u32) {

            if ux >= 0xffu32 << 24 {

                // x is inf or nan

                if ux > (0xffu32 << 24) {

                    return x + x;

                } // x = nan

                static IR: [f32; 2] = [f32::INFINITY, 0.];

                return IR[(t >> 31) as usize]; // x = +-inf

            }

            if t >= 0xc3150000u32 {

                // x < -149

                let z = x;

                let mut y = f_fmla(

                    z as f64 + 149.,

                    f64::from_bits(0x3690000000000000),

                    f64::from_bits(0x36a0000000000000),

                );

                y = y.max(f64::from_bits(0x3680000000000000));

                return y as f32;

            }

            // now x >= 128

            let r = black_box(f64::from_bits(0x47e0000000000000))

                * black_box(f64::from_bits(0x47e0000000000000));

            return r as f32;

        }

    }

    if ux <= 0x7a000000u32 {

        // |x| < 1/32

        // Generated by Sollya exp2 on range [-1/32;1/32]:

        // d = [-1/32, 1/32];

        // f_exp2f = (2^y - 1)/y;

        // Q = fpminimax(f_exp2f, 5, [|D...|], d, relative, floating);

        // See ./notes/exp2f_small.sollya

        const C: [u64; 6] = [

            0x3fe62e42fefa39f3,

            0x3fcebfbdff82c57b,

            0x3fac6b08d6f2d7aa,

            0x3f83b2ab6fc92f5d,

            0x3f55d897cfe27125,

            0x3f243090e61e6af1,

        ];

        let xd = x as f64;

        let p = f_polyeval6(

            xd,

            f64::from_bits(C[0]),

            f64::from_bits(C[1]),

            f64::from_bits(C[2]),

            f64::from_bits(C[3]),

            f64::from_bits(C[4]),

            f64::from_bits(C[5]),

        );

        return f_fmla(p, xd, 1.) as f32;

    }

    let x_d = x as f64;

    let kf = (x_d * 64.).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, x_d);

    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_fmla(dx_sq, c3, c2);

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

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

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

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

    f_fmla(p, dx_sq * mh, c1 * mh) as f32

}

#[inline]

pub(crate) fn dirty_exp2f(d: f32) -> f32 {

    let redux = f32::from_bits(0x4b400000) / TBLSIZE as f32;

    let ui = f32::to_bits(d + redux);

    let mut i0 = ui;

    i0 = i0.wrapping_add(TBLSIZE as u32 / 2);

    let k = i0 / TBLSIZE as u32;

    i0 &= TBLSIZE as u32 - 1;

    let mut uf = f32::from_bits(ui);

    uf -= redux;

    let item = EXP2FT.0[i0 as usize];

    let z0: f32 = f32::from_bits(item.0);

    let f: f32 = d - uf;

    let mut u = 0.24022650695908768;

    u = f_fmlaf(u, f, 0.69314718055994973);

    u *= f;

    let i2 = pow2if(k as i32);

    f_fmlaf(u, z0, z0) * i2

}

Unexecuted instantiation: pxfm::exponents::exp2f::dirty_exp2f

Unexecuted instantiation: pxfm::exponents::exp2f::dirty_exp2f

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_exp2f() {

        assert_eq!(f_exp2f(1. / 64.), 1.0108893);

        assert_eq!(f_exp2f(2.0), 4.0);

        assert_eq!(f_exp2f(3.0), 8.0);

        assert_eq!(f_exp2f(4.0), 16.0);

        assert_eq!(f_exp2f(10.0), 1024.0);

        assert_eq!(f_exp2f(-10.0), 0.0009765625);

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

        assert_eq!(f_exp2f(-0.35), 0.7845841);

        assert_eq!(f_exp2f(0.35), 1.2745606);

        assert!(f_exp2f(f32::INFINITY).is_infinite());

        assert_eq!(f_exp2f(f32::NEG_INFINITY), 0.0);

    }

    #[test]

    fn test_dirty_exp2f() {

        assert!((dirty_exp2f(0.35f32) - 0.35f32.exp2()).abs() < 1e-5);

        assert!((dirty_exp2f(-0.6f32) - (-0.6f32).exp2()).abs() < 1e-5);

    }

}