/*

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

// Polynomials approximating x/erf(x) on ( k/8, (k + 1)/8 ) generated by Sollya and SageMath

// ```text

// def build_sollya_script(idx):

//     return f"""

// d = [{idx}/8, {idx + 1}/8];

//

// f = x/erf(x);

// Q = fpminimax(f, [|0, 2, 4, 6, 8, 10, 12, 14|], [|D...|], d);

//

// for i from 0 to degree(Q) by 2 do {{

//     write(coeff(Q, i)) >> "coefficients.txt";

//     write("\\n") >> "coefficients.txt";

// }};

// """

//

// def load_coefficients(filename):

//     with open(filename, "r") as f:

//         return [RealField(500)(line.strip()) for line in f if line.strip()]

//

// def call_sollya_on_interval(idx):

//     sollya_script = build_sollya_script(idx)

//     with open("tmp_interval.sollya", "w") as f:

//         f.write(sollya_script)

//     import subprocess

//     if os.path.exists("coefficients.txt"):

//         os.remove("coefficients.txt")

//     try:

//         result = subprocess.run(

//             ["sollya", "tmp_interval.sollya"],

//             check=True,

//             capture_output=True,

//             text=True

//         )

//     except subprocess.CalledProcessError as e:

//         return

//

// def print_coeffs(poly):

//     print("[")

//     for i in range(len(poly)):

//         coeff = poly[i]

//         print(double_to_hex(coeff), ",")

//     print("],")

//

// print(f"static COEFFS: [[u64; 8]; 32] = [")

// for i in range(0, 32):

//     call_sollya_on_interval(i)

//     coeffs = load_coefficients(f"coefficients.txt")

//     print_coeffs(coeffs)

// print("];")

// ```

static COEFFS: [[u64; 8]; 32] = [

    [

        0x3fec5bf891b4ef6b,

        0x3fd2e7fb0bcdee7f,

        0x3f842aa5641a200a,

        0xbf752081ae81d16e,

        0x3f2e1a191fb85592,

        0xbf203a20ad500043,

        0x3f861a864f719e76,

        0xbfc79f68bad20bd1,

    ],

    [

        0x3fec5bf891b4ef6b,

        0x3fd2e7fb0bcdf45b,

        0x3f842aa561f35512,

        0xbf75207c8167ac1d,

        0x3f2db4b119b4ce20,

        0x3f20fa28737c4219,

        0xbef38e74cca2219a,

        0xbec5d70713fc621e,

    ],

    [

        0x3fec5bf891b4ef30,

        0x3fd2e7fb0bce1c0f,

        0x3f842aa56138541f,

        0xbf75207c6197eb7c,

        0x3f2db4799120e074,

        0x3f20fc28d915a6e9,

        0xbef3ea5b479dc053,

        0xbebbffe6df8ec372,

    ],

    [

        0x3fec5bf891b4bf18,

        0x3fd2e7fb0bde166f,

        0x3f842aa53c721766,

        0xbf7520796733bbec,

        0x3f2db21eebf4144f,

        0x3f210545cc78d0f0,

        0xbef48ad7e4aa2d10,

        0xbeb24a043ad31907,

    ],

    [

        0x3fec5bf891ab16e9,

        0x3fd2e7fb0dc7b919,

        0x3f842aa29d8381e7,

        0xbf7520592585d601,

        0x3f2da30df1566e43,

        0x3f212780ff325aa6,

        0xbef5e98ea9819e42,

        0xbe9849d52099dcb9,

    ],

    [

        0x3fec5bf890ddfa8d,

        0x3fd2e7fb28aab312,

        0x3f842a8a461f0eb7,

        0xbf751f93b2d27114,

        0x3f2d66789eed5f95,

        0x3f21818ff1832f50,

        0xbef84264724049ef,

        0x3e9df12b02e82a5a,

    ],

    [

        0x3fec5bf887f64fa4,

        0x3fd2e7fbfcc05f75,

        0x3f842a02323e2099,

        0xbf751c86d291ced6,

        0x3f2cbd5653cde433,

        0x3f223299b32b8583,

        0xbefb7fc6e286cd94,

        0x3eb49676cb3da393,

    ],

    [

        0x3fec5bf84f8e2488,

        0x3fd2e7ffe83d2974,

        0x3f842821c5cc659c,

        0xbf7514805a6196e3,

        0x3f2b723680f64bb5,

        0x3f233416dcfcd366,

        0xbefefe55300afaa7,

        0x3ebf0c475fb71e7a,

    ],

    [

        0x3fec5bf7999e6afe,

        0x3fd2e809c6d4caa7,

        0x3f84247256be4a56,

        0xbf750838db0c0cf5,

        0x3f29e7e867267388,

        0x3f24226adee5ce74,

        0xbf00c0830af2bf01,

        0x3ec26fb6b18e628b,

    ],

    [

        0x3fec5bf801fc5ad5,

        0x3fd2e80618e8941e,

        0x3f84254c04b0b234,

        0xbf7509d7cf351201,

        0x3f2a01829944820c,

        0x3f241d7bb0c7e2de,

        0xbf00c2d844916d26,

        0x3ec2817d39abc26b,

    ],

    [

        0x3fec5c0938a12f13,

        0x3fd2e7706c510d79,

        0x3f8448392db86aae,

        0xbf75526e9c6046f0,

        0x3f2facd0bc0e7862,

        0x3f21fc4093e1e6b7,

        0xbefdf54af68ba968,

        0x3ebfe348fc246c15,

    ],

    [

        0x3fec5c6dcdadc5d8,

        0x3fd2e495072afff3,

        0x3f84d6f390564d4d,

        0xbf764a7e85749c85,

        0x3f37effb62caee80,

        0x3f19cb39bc236ae6,

        0xbef6d7035785e8f3,

        0x3eb755aa2e58fc52,

    ],

    [

        0x3fec5dd74381acff,

        0x3fd2dbe68140f116,

        0x3f86459e1acfda0f,

        0xbf7865203923a03d,

        0x3f43665053a48049,

        0x3f0409e353b761ea,

        0xbeeb0b00f567c9f8,

        0x3eabc33000611b25,

    ],

    [

        0x3fec6175431226d1,

        0x3fd2c8dcbb0babcc,

        0x3f88f5bfd61e5d2e,

        0xbf7bc60de8dff620,

        0x3f4d9b7076c7767c,

        0xbf0106584fac3547,

        0xbed0a56cd1030deb,

        0x3e970ee11e7beb48,

    ],

    [

        0x3fec68445d99a8e9,

        0x3fd2a9d608dbfea2,

        0x3f8cc072ddf22cb6,

        0xbf7fe5f2efdc5f5c,

        0x3f5431d1deff38bc,

        0xbf197220e4a1dda8,

        0x3ec9e2469e6c1c67,

        0x3e4be72535d53d7b,

    ],

    [

        0x3fec713c415bf088,

        0x3fd28610e83aa38c,

        0x3f9049ee1942f46b,

        0xbf81c513d165d6fd,

        0x3f585bc13e0fcaba,

        0xbf22715362e30768,

        0x3ede6bfa3c69e8e3,

        0xbe852cd85f8dea5b,

    ],

    [

        0x3fec770e08b47107,

        0x3fd2716324b22047,

        0x3f91460d403e6b9c,

        0xbf829ab46375f10d,

        0x3f5a0e7f00c76fb5,

        0xbf2484890f2d7eeb,

        0x3ee207b21bbd8496,

        0xbe8bbee036671b6a,

    ],

    [

        0x3fec6f4a2d01088d,

        0x3fd288e494bc89b7,

        0x3f905203788a2821,

        0xbf81eab2727ce365,

        0x3f58ddba75a3c100,

        0xbf2347c9a317a175,

        0x3ee099c93ce5f44f,

        0xbe88e9f9c064f833,

    ],

    [

        0x3fec4c9bbce50c7d,

        0x3fd2e8175b0e1837,

        0x3f89a2d1518c7a4c,

        0xbf7f3fa91859127e,

        0x3f55431c495b1077,

        0xbf1fc1af665bb1f8,

        0x3eda0f1d735195cb,

        0xbe827b8d6fa224ed,

    ],

    [

        0x3fec03cce39d7213,

        0x3fd39c2316e290bf,

        0x3f7b674438899313,

        0xbf783644c88c71fb,

        0x3f5047a3da485180,

        0xbf1748b54f823d57,

        0x3ed20c86d3302f22,

        0xbe77f94cafe045a8,

    ],

    [

        0x3feb913f0adf6c4b,

        0x3fd49c4cedae09fc,

        0xbf4a6dea9778f474,

        0xbf7006dc4e6c8125,

        0x3f461483c254fa5f,

        0xbf0e75052760bf18,

        0x3ec65425869bc096,

        0xbe6bc2df9fbc0f82,

    ],

    [

        0x3feafbeda3b7d400,

        0x3fd5cb900ee1fb5e,

        0xbf8228d16e87de3d,

        0xbf6011d44e155bf5,

        0x3f3993b736442257,

        0xbf017c7ee5efa6ad,

        0x3eb886e337d2e3c2,

        0xbe5cba4b79e90043,

    ],

    [

        0x3fea54849d309eba,

        0x3fd701afa55e3d21,

        0xbf90c72bb2e2799f,

        0xbf33c92573294e34,

        0x3f265284f7a6d53a,

        0xbef09f09298ed1e8,

        0x3ea7153a46cb2e27,

        0xbe49ef6ec79265fd,

    ],

    [

        0x3fe9b128df667870,

        0x3fd816d295a867cb,

        0xbf9713f11ea84a26,

        0x3f4edcbdd63903bb,

        0x3ef44f54fc7a6024,

        0xbed45da547d2fcb8,

        0x3e9049754d57a9a7,

        0xbe32aba05ca26a69,

    ],

    [

        0x3fe927f49edf4ace,

        0x3fd8ecd207c6a7d1,

        0xbf9b8cd215124008,

        0x3f5cbab209dd389d,

        0xbf12a8920ea6230f,

        0x3eb442dfce60b0e2,

        0x3e52494e415c7728,

        0xbe09a1b1bbb9cee4,

    ],

    [

        0x3fe8ca3d7437d06f,

        0x3fd973c08b6d33fb,

        0xbf9e272ca1fccc06,

        0x3f61efd00e2016b6,

        0xbf1e6dab18e9d45a,

        0x3ed0b446e3469be1,

        0xbe7503c584488bed,

        0x3e069968660290a4,

    ],

    [

        0x3fe8a1f4b154f663,

        0x3fd9a9a8b81692d4,

        0xbf9f1e9312dd4501,

        0x3f632b4d20599404,

        0xbf2119c1b5e43b24,

        0x3ed42b9874284d56,

        0xbe7c17cc1eef4b9d,

        0x3e117f0a9057a689,

    ],

    [

        0x3fe8b15bfcf78f33,

        0x3fd99720c884ab33,

        0xbf9ed2265979f5a6,

        0x3f62d3c30432692b,

        0xbf20a17346c37362,

        0x3ed36538f2d21c31,

        0xbe7aac6bb10f8b90,

        0x3e1061d3a1737044,

    ],

    [

        0x3fe8f479e98cb825,

        0x3fd94ab3f8d0c80c,

        0xbf9da7afe85abf94,

        0x3f618fe28f71a3d4,

        0xbf1df723b2a63e38,

        0x3ed0d190252a7f7c,

        0xbe7631fdd49272b0,

        0x3e0a17567cab4a94,

    ],

    [

        0x3fe9636d647b61c0,

        0x3fd8d4aaba0e0212,

        0xbf9bf904574e56ea,

        0x3f5fb68684d8555d,

        0xbf19d06f9cf17bbf,

        0x3ecb92b9f0b8acf3,

        0xbe7145bde0c499ae,

        0x3e033cf1cb08ce4c,

    ],

    [

        0x3fe9f4c3301b6d33,

        0x3fd844100b4598b3,

        0xbf9a0b94e19be990,

        0x3f5c0ed55c70532f,

        0xbf15a786c9e62b23,

        0x3ec5e3f05a04f5c6,

        0xbe69ea9db2e37883,

        0x3dfb3e5ad2cd0fb2,

    ],

    [

        0x3fea9f469c75536c,

        0x3fd7a51b3d9eda10,

        0xbf980f63a2cb486c,

        0x3f5887f72a9f07e0,

        0xbf11e4d454f2f994,

        0x3ec113d0aed8bdef,

        0xbe6311f84083acf4,

        0x3df2e4dc2e50e3fa,

    ],

];

/// Computes 1/erf(x)

///

/// Max ulp 0.5

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

    let x_u = x.to_bits();

    let x_abs = x_u & 0x7fff_ffffu32;

    if x == 0. {

        return if x.is_sign_negative() {

            f32::NEG_INFINITY

        } else {

            f32::INFINITY

        };

    }

    if x_abs >= 0x4080_0000u32 {

        static ONE: [f32; 2] = [1.0, -1.0];

        static SMALL: [f32; 2] = [f32::from_bits(0xb3000000), f32::from_bits(0x33000000)];

        let sign = x.is_sign_negative() as usize;

        if x_abs >= 0x7f80_0000u32 {

            return if x_abs > 0x7f80_0000 { x } else { ONE[sign] };

        }

        return ONE[sign] + SMALL[sign];

    }

    // Polynomial approximation see [COEFFS] for generation:

    //   1/erf(x) ~ (c0 + c1 * x^2 + c2 * x^4 + ... + c7 * x^14) / x

    let xd = x as f64;

    let xsq = xd * xd;

    const EIGHT: u32 = 3 << 23;

    let idx = unsafe { f32::from_bits(x_abs.wrapping_add(EIGHT)).to_int_unchecked::<usize>() };

    let c = COEFFS[idx];

    let x4 = xsq * xsq;

    let c0 = f_fmla(xsq, f64::from_bits(c[1]), f64::from_bits(c[0]));

    let c1 = f_fmla(xsq, f64::from_bits(c[3]), f64::from_bits(c[2]));

    let c2 = f_fmla(xsq, f64::from_bits(c[5]), f64::from_bits(c[4]));

    let c3 = f_fmla(xsq, f64::from_bits(c[7]), f64::from_bits(c[6]));

    let x8 = x4 * x4;

    let p0 = f_fmla(x4, c1, c0);

    let p1 = f_fmla(x4, c3, c2);

    ((f_fmla(x8, p1, p0)) / xd) as f32

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn f_erff_test() {

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

        assert_eq!(f_rerff(0.0), f32::INFINITY);

        assert_eq!(f_rerff(-0.0), f32::NEG_INFINITY);

        assert_eq!(f_rerff(0.015255669), 58.096153);

        assert_eq!(f_rerff(1.0), 1.1866608);

        assert_eq!(f_rerff(0.5), 1.9212301);

    }

}