/*

 * // 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::{dd_fmla, dyad_fmla, f_fmla};

use crate::double_double::DoubleDouble;

use crate::shared_eval::poly_dd_3;

pub(crate) static ATAN_CIRCLE: [[u16; 3]; 31] = [

    [419, 81, 0],

    [500, 81, 0],

    [582, 163, 0],

    [745, 163, 0],

    [908, 326, 0],

    [1234, 326, 0],

    [1559, 651, 0],

    [2210, 650, 1],

    [2860, 1299, 3],

    [4156, 1293, 4],

    [5444, 2569, 24],

    [7989, 2520, 32],

    [10476, 4917, 168],

    [15224, 4576, 200],

    [19601, 8341, 838],

    [27105, 6648, 731],

    [33036, 10210, 1998],

    [41266, 6292, 1117],

    [46469, 7926, 2048],

    [52375, 4038, 849],

    [55587, 4591, 1291],

    [58906, 2172, 479],

    [60612, 2390, 688],

    [62325, 1107, 247],

    [63192, 1207, 349],

    [64056, 556, 124],

    [64491, 605, 175],

    [64923, 278, 62],

    [65141, 303, 88],

    [65358, 139, 31],

    [65467, 151, 44],

];

pub(crate) static ATAN_REDUCE: [(u64, u64); 129] = [

    (0x0000000000000000, 0x0000000000000000),

    (0x3f89224e047e368e, 0x3c1a3ca6c727c59d),

    (0x3f992346247a91f0, 0x3bf138b0ef96a186),

    (0x3fa2dbaae9a05db0, 0x3c436e7f8a3f5e42),

    (0x3fa927278a3b1162, 0xbbfac986efb92662),

    (0x3faf7495ea3f3783, 0x3c406ec8011ee816),

    (0x3fb2e239ccff3831, 0xbc5858437d431332),

    (0x3fb60b9f7597fdec, 0xbc3cebd13eb7c513),

    (0x3fb936bb8c5b2da2, 0xbc5840cac0d81db5),

    (0x3fbc63ce377fc802, 0x3c5400b0fdaa109e),

    (0x3fbf93183a8db9e9, 0x3c40e04e06c86e72),

    (0x3fc1626d85a91e70, 0x3c4f7ad829163ca7),

    (0x3fc2fcac73a60640, 0xbc52680735ce2cd8),

    (0x3fc4986a74cf4e57, 0xbc690559690b42e4),

    (0x3fc635c990ce0d36, 0x3c591d29110b41aa),

    (0x3fc7d4ec54fb5968, 0xbc4ea90e27182780),

    (0x3fc975f5e0553158, 0xbc2dc82ac14e3e1c),

    (0x3fcb1909efd8b762, 0xbc573a10fd13daaf),

    (0x3fccbe4ceb4b4cf2, 0xbc63a7ffbeabda0b),

    (0x3fce65e3f27c9f2a, 0xbc6db6627a24d523),

    (0x3fd007fa758626ae, 0xbc645f97dd3099f6),

    (0x3fd0de53475f3b3c, 0xbc66293f68741816),

    (0x3fd1b6103d3597e9, 0xbc6ab240d40633e9),

    (0x3fd28f459ecad74d, 0xbc2de34d14e832e0),

    (0x3fd36a08355c63dc, 0x3c6af540d9fb4926),

    (0x3fd4466d542bac92, 0x3c6da60fdbc82ac4),

    (0x3fd5248ae1701b17, 0xbc792a601170138a),

    (0x3fd604775fbb27df, 0xbc67f1fca1d5d15b),

    (0x3fd6e649f7d78649, 0xbc64e223ea716c7b),

    (0x3fd7ca1a832d0f84, 0x3c7b24c824ac51fc),

    (0x3fd8b00196b3d022, 0x3c64314cd132ba43),

    (0x3fd998188e816bf0, 0xbc711f1e0817879a),

    (0x3fda827999fcef32, 0xbc6c3dea4dbad538),

    (0x3fdb6f3fc8c61e5b, 0x3c660d1b780ee3eb),

    (0x3fdc5e87185e67b6, 0xbc4ab5edb7dfa545),

    (0x3fdd506c82a2c800, 0xbc68e1437048b5bd),

    (0x3fde450e0d273e7a, 0xbc706951c97b050f),

    (0x3fdf3c8ad985d9ee, 0xbc414af9522ab518),

    (0x3fe01b819b5a7cf7, 0xbc7aba0d7d97d1f2),

    (0x3fe09a4c59bd0d4d, 0x3c4095bc4ebc2c42),

    (0x3fe11ab7190834ec, 0x3c8798826fa27774),

    (0x3fe19cd3fe8e405d, 0x3c8008f6258fc98f),

    (0x3fe220b5ef047825, 0xbc5462af7ceb7de6),

    (0x3fe2a6709a74f289, 0xbc71184dfd78b472),

    (0x3fe32e1889047ffd, 0x3c79141876dc40c5),

    (0x3fe3b7c3289ed6f3, 0x3c8481c20189726c),

    (0x3fe44386db9ce5db, 0x3c82e851bd025441),

    (0x3fe4d17b087b265d, 0x3c713ada9b8bc419),

    (0x3fe561b82ab7f990, 0xbc805b4c3c4cbee8),

    (0x3fe5f457e4f4812e, 0xbc85619249bd96f1),

    (0x3fe6897514751db6, 0xbc6b0a0fbcafc671),

    (0x3fe7212be621be6d, 0xbc819ff2dc66da45),

    (0x3fe7bb99ed2990cf, 0x3c81320449592d92),

    (0x3fe858de3b716571, 0xbc81fddcd2f3da8e),

    (0x3fe8f9197bf85eeb, 0x3c6d44a42e35cc97),

    (0x3fe99c6e0f634394, 0xbc7585a178b4a18d),

    (0x3fea43002ae42850, 0x3c6f95a531b3a970),

    (0x3feaecf5f9ba35a6, 0xbc396c2d43ca3392),

    (0x3feb9a77c18c1af2, 0xbc6a5bed94b05def),

    (0x3fec4bb009e77983, 0x3c454509d2bff511),

    (0x3fed00cbc7384d2e, 0xbc6b4c867cef300c),

    (0x3fedb9fa89953fcf, 0xbc1ddfac663d6bc6),

    (0x3fee776eafc91706, 0xbc7a510683ff7cb6),

    (0x3fef395d9f0e3c92, 0x3c44fdcd8e4e8710),

    (0x3ff0000000000000, 0x0000000000000000),

    (0x3ff065c900aaf2d8, 0xbc8deec7fc9042ad),

    (0x3ff0ce29d0883c99, 0xbc8395ae45e0657d),

    (0x3ff139447e6a86ee, 0x3c8332cf301a97f3),

    (0x3ff1a73d55278c4b, 0xbc86cc8c4b78213b),

    (0x3ff2183b0c4573ff, 0x3c870a90841da57a),

    (0x3ff28c66fdaf8f09, 0xbc6ba39bad450ee0),

    (0x3ff303ed61109e20, 0xbc88692946d9f93c),

    (0x3ff37efd8d87607e, 0x3c63b711bf765b58),

    (0x3ff3fdca42847507, 0x3c7c21387985b081),

    (0x3ff48089f8bf42cc, 0xbc87ddb19d3d0efc),

    (0x3ff507773c537ead, 0xbc7f5e354cf971f3),

    (0x3ff592d11142fa55, 0xbc700f0ad675330d),

    (0x3ff622db63c8ecc2, 0xbc82c93f50ab2c0e),

    (0x3ff6b7df86265200, 0x3c7bec391adc37d5),

    (0x3ff7522cbdd428a8, 0xbc69686ddc9ffcf5),

    (0x3ff7f218e25a7461, 0xbc78d16529514246),

    (0x3ff89801106cc709, 0xbc8092f51e9c2803),

    (0x3ff9444a7462122a, 0xbc807c06755404c4),

    (0x3ff9f7632fa9e871, 0x3c802e0d43abc92b),

    (0x3ffab1c35d8a74ea, 0x3c5d0184e48af6f7),

    (0x3ffb73ee3c3ef16a, 0x3c773be957380bc2),

    (0x3ffc3e738086bc0f, 0xbc702b6e26c84462),

    (0x3ffd11f0dae40609, 0x3c525c4f3ffa6e1f),

    (0x3ffdef13b73c1406, 0xbc5e302db3c6823f),

    (0x3ffed69b4153a45d, 0x3c73207830326c0e),

    (0x3fffc95abad6cf4a, 0xbc66308cee7927bf),

    (0x4000641e192ceab3, 0xbc70147ebf0df4c5),

    (0x4000ea21d716fbf7, 0xbc7168533cc41d8b),

    (0x40017749711a6679, 0xbc652a0b0333e9c5),

    (0x40020c36c6a7f38e, 0x3c68659eece35395),

    (0x4002a99f50fd4f4f, 0x3c820fcad18cb36f),

    (0x4003504f333f9de6, 0xbc752afdbd5a8c74),

    (0x4004012ce2586a17, 0xbc79747a792907d7),

    (0x4004bd3d87fe0650, 0x3c790c59393b52c8),

    (0x400585aa4e1530fa, 0x3c7af6934f13a3a8),

    (0x40065bc6cc825147, 0xbc48534dcab5ad3e),

    (0x40074118e4b6a7c8, 0xbc7555aa8bfca9a1),

    (0x400837626d70fdb8, 0xbc556b3fee9ca72b),

    (0x400940ad30abc792, 0x3c54b3fdd4fdc06c),

    (0x400a5f59e90600dd, 0x3c6285d367c55ddc),

    (0x400b9633283b6d14, 0xbc48712976f17a16),

    (0x400ce885653127e7, 0xbc3abe8ab65d49fc),

    (0x400e5a3de972a377, 0x3c5cd9be81ad764b),

    (0x400ff01305ecd8dc, 0x3c4742c2922656fa),

    (0x4010d7dc7cff4c9e, 0xbc77c842978bee09),

    (0x4011d0143e71565f, 0x3c67bc7dea7c3c03),

    (0x4012e4ff1626b949, 0x3c4aefbe25b404e9),

    (0x40141bfee2424771, 0xbc34bcfaaa95cb2c),

    (0x40157be4eaa5e11b, 0x3c50fe741e4ec679),

    (0x40170d751908c1b1, 0x3c5fe74a5b0ec709),

    (0x4018dc25c117782b, 0x3c50ca1c19f710ef),

    (0x401af73f4ca3310f, 0x3c52867b40ba77d6),

    (0x401d7398d15e70db, 0x3c60fd4e0d4b1547),

    (0x4020372fb36b87e2, 0x3c5c16c9ecc1621d),

    (0x402208dbdae055ef, 0x3c56b81a36e75e8c),

    (0x40244e6c595afdcc, 0xbc57c22045771848),

    (0x4027398c57f3f1ad, 0x3c5970503be105c0),

    (0x402b1d03c03d2f7f, 0xbc3f299d010aead2),

    (0x403046e9fe60a77e, 0x3c5d2b61deff33ec),

    (0x40345affed201b55, 0x3bf0e84d9567203a),

    (0x403b267195b1ffae, 0xbbfad44b44b92653),

    (0x40445e2455e4aaa7, 0xbc3296d577b5e21d),

    (0x40545eed6854ce99, 0x3c02db53886013ca),

    (0x0000000000000000, 0x0000000000000000),

];

#[cold]

fn atan_refine(x: f64, a: f64) -> f64 {

    const CH: [(u64, u64); 3] = [

        (0xbfd5555555555555, 0xbc75555555555555),

        (0x3fc999999999999a, 0xbc6999999999bcb8),

        (0xbfc2492492492492, 0xbc6249242093c016),

    ];

    const CL: [u64; 4] = [

        0x3fbc71c71c71c71c,

        0xbfb745d1745d1265,

        0x3fb3b13b115bcbc4,

        0xbfb1107c41ad3253,

    ];

    let phi = f_fmla(a.abs(), f64::from_bits(0x40545f306dc9c883), 256.5).to_bits();

    let i = ((phi >> (52 - 8)) & 0xff) as i64;

    let h: DoubleDouble = if i == 128 {

        DoubleDouble::from_quick_recip(-x)

    } else {

        let ta = f64::copysign(f64::from_bits(ATAN_REDUCE[i as usize].0), x);

        let dzt = DoubleDouble::from_exact_mult(x, ta);

        let zmta = x - ta;

        let v = 1.0 + dzt.hi;

        let d = 1.0 - v;

        let ev = (d + dzt.hi) - ((d + v) - 1.0) + dzt.lo;

        let r = 1.0 / v;

        let rl = f_fmla(-ev, r, dd_fmla(r, -v, 1.0)) * r;

        DoubleDouble::quick_mult_f64(DoubleDouble::new(rl, r), zmta)

    };

    let h2 = DoubleDouble::quick_mult(h, h);

    let h3 = DoubleDouble::quick_mult(h, h2);

    let h4 = h2.hi * h2.hi;

    let zw0 = f_fmla(h2.hi, f64::from_bits(CL[3]), f64::from_bits(CL[2]));

    let zw1 = f_fmla(h2.hi, f64::from_bits(CL[1]), f64::from_bits(CL[0]));

    let zfl = f_fmla(h2.hi, zw1, h4 * zw0);

    let mut f = poly_dd_3(h2, CH, zfl);

    f = DoubleDouble::quick_mult(h3, f);

    let (ah, mut az);

    if i == 0 {

        ah = h.hi;

        az = f;

    } else {

        let mut df = 0.;

        if i < 128 {

            df = f64::copysign(1.0, x) * f64::from_bits(ATAN_REDUCE[i as usize].1);

        }

        let id = f64::copysign(i as f64, x);

        ah = f64::from_bits(0x3f8921fb54442d00) * id;

        az = DoubleDouble::new(

            f64::from_bits(0xb97fc8f8cbb5bf80) * id,

            f64::from_bits(0x3c88469898cc5180) * id,

        );

        az = DoubleDouble::add(az, DoubleDouble::new(0., df));

        az = DoubleDouble::add(az, h);

        az = DoubleDouble::add(az, f);

    }

    let v0 = DoubleDouble::from_exact_add(ah, az.hi);

    let v1 = DoubleDouble::from_exact_add(v0.lo, az.lo);

    v1.hi + v0.hi

}

/// fma - fma

/// dd_fma - DD fma fallback

/// dyad_fma - mandatory fma fallback

#[inline(always)]

fn atan_gen_impl<

    Q: Fn(f64, f64, f64) -> f64,

    F: Fn(f64, f64, f64) -> f64,

    D: Fn(f64, f64, f64) -> f64,

>(

    x: f64,

    fma: Q,

    dd_fma: F,

    dyad_fma: D,

) -> f64 {

    const CH: [u64; 4] = [

        0x3ff0000000000000,

        0xbfd555555555552b,

        0x3fc9999999069c20,

        0xbfc248d2c8444ac6,

    ];

    let t = x.to_bits();

    let at = t & 0x7fffffffffffffff; // at encodes |x|

    let mut i = (at.wrapping_shr(51) as i64).wrapping_sub(2030i64); // -2030 <= i <= 2065

    if at < 0x3f7b21c475e6362au64 {

        // |x| < 0x1.b21c475e6362ap-8

        if at == 0 {

            return x;

        } // atan(+/-0) = +/-0

        const CH2: [u64; 4] = [

            0xbfd5555555555555,

            0x3fc99999999998c1,

            0xbfc249249176aec0,

            0x3fbc711fd121ae80,

        ];

        if at < 0x3e40000000000000u64 {

            // |x| < 0x1p-27

            /* underflow when 0 < |x| < 2^-1022 or when |x| = 2^-1022

            and rounding towards zero. */

            return dyad_fma(f64::from_bits(0xbc90000000000000), x, x);

        }

        let x2 = x * x;

        let x3 = x * x2;

        let x4 = x2 * x2;

        let w0 = fma(x2, f64::from_bits(CH2[3]), f64::from_bits(CH2[2]));

        let w1 = fma(x2, f64::from_bits(CH2[1]), f64::from_bits(CH2[0]));

        let f = x3 * fma(x4, w0, w1);

        let ub = fma(f, f64::from_bits(0x3cd2000000000000), f) + x;

        let lb = fma(-f, f64::from_bits(0x3cc4000000000000), f) + x;

        // Ziv's accuracy test

        if ub != lb {

            return atan_refine(x, ub);

        }

        return ub;

    }

    let h;

    let ah;

    let mut al;

    if at > 0x4062ded8e34a9035u64 {

        // |x| > 0x4062ded8e34a9035

        ah = f64::copysign(f64::from_bits(0x3ff921fb54442d18), x);

        al = f64::copysign(f64::from_bits(0x3c91a62633145c07), x);

        if at >= 0x434d02967c31cdb5u64 {

            // |x| >= 1.63312e+16

            if at > (0x7ffu64 << 52) {

                return x + x;

            } // NaN

            return ah + al;

        }

        h = -1.0 / x;

    } else {

        // now 0.006624 <= |x| <= 150.964 thus 1<=i<=30

        let u: u64 = t & 0x0007ffffffffffff;

        let ut: u64 = u.wrapping_shr(51 - 16);

        let ut2: u64 = (ut * ut).wrapping_shr(16);

        let lc = ATAN_CIRCLE[i as usize];

        i = (((lc[0] as u64)

            .wrapping_shl(16)

            .wrapping_add(ut * lc[1] as u64)

            .wrapping_sub(ut2 * lc[2] as u64))

            >> (16 + 9)) as i64;

        let la = ATAN_REDUCE[i as usize];

        let ta = f64::copysign(1.0, x) * f64::from_bits(la.0);

        let id = f64::copysign(1.0, x) * i as f64;

        al = dd_fma(

            f64::copysign(1.0, x),

            f64::from_bits(la.1),

            f64::from_bits(0x3c88469898cc5170) * id,

        );

        h = (x - ta) / fma(x, ta, 1.0);

        ah = f64::from_bits(0x3f8921fb54442d00) * id;

    }

    let h2 = h * h;

    let h4 = h2 * h2;

    let f0 = fma(h2, f64::from_bits(CH[3]), f64::from_bits(CH[2]));

    let f1 = fma(h2, f64::from_bits(CH[1]), f64::from_bits(CH[0]));

    let f = fma(h4, f0, f1);

    al = dd_fma(h, f, al);

    let ub = fma(h, f64::from_bits(0x3ccf800000000000), al) + ah;

    let lb = fma(-h, f64::from_bits(0x3ccf800000000000), al) + ah;

    // Ziv's accuracy test

    if lb != ub {

        return atan_refine(x, ub);

    }

    ub

}

Unexecuted instantiation: pxfm::tangent::atan::atan_gen_impl::<<f64>::mul_add, <f64>::mul_add, <f64>::mul_add>

Unexecuted instantiation: pxfm::tangent::atan::atan_gen_impl::<pxfm::common::f_fmla, pxfm::common::dd_fmla, pxfm::common::dyad_fmla>

#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]

#[target_feature(enable = "avx", enable = "fma")]

unsafe fn atan_fma_impl(x: f64) -> f64 {

    atan_gen_impl(x, f64::mul_add, f64::mul_add, f64::mul_add)

}

/// Computes atan in double precision

///

/// ULP 0.5

pub fn f_atan(x: f64) -> f64 {

    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]

    {

        atan_gen_impl(x, f_fmla, dd_fmla, dyad_fmla)

    }

    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]

    {

        use std::sync::OnceLock;

        static EXECUTOR: OnceLock<unsafe fn(f64) -> f64> = OnceLock::new();

        let q = EXECUTOR.get_or_init(|| {

            if std::arch::is_x86_feature_detected!("avx")

                && std::arch::is_x86_feature_detected!("fma")

            {

                atan_fma_impl

            } else {

                fn def_atan(x: f64) -> f64 {

                    atan_gen_impl(x, f_fmla, dd_fmla, dyad_fmla)

                }

                def_atan

            }

        });

        unsafe { q(x) }

    }

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn atan_test() {

        assert_eq!(f_atan(7.7877585082074305E-308), 7.78775850820743e-308);

        assert_eq!(f_atan(0.0), 0.0);

        assert_eq!(f_atan(1.0), 0.7853981633974483096156608458198);

        assert_eq!(f_atan(35.9), 1.542948374599341097473183563168947);

        assert_eq!(f_atan(-35.9), -1.542948374599341097473183563168947);

        assert_eq!(f_atan(f64::INFINITY), 1.5707963267948966);

    }

}