pxfm 0.1.29

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