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::{dd_fmla, f_fmla};
use std::hint::black_box;

static ERR0: [u64; 128] = [
    0x3ff0000000000000,
    0x3ff0163da9fb3335,
    0x3ff02c9a3e778061,
    0x3ff04315e86e7f85,
    0x3ff059b0d3158574,
    0x3ff0706b29ddf6de,
    0x3ff0874518759bc8,
    0x3ff09e3ecac6f383,
    0x3ff0b5586cf9890f,
    0x3ff0cc922b7247f7,
    0x3ff0e3ec32d3d1a2,
    0x3ff0fb66affed31b,
    0x3ff11301d0125b51,
    0x3ff12abdc06c31cc,
    0x3ff1429aaea92de0,
    0x3ff15a98c8a58e51,
    0x3ff172b83c7d517b,
    0x3ff18af9388c8dea,
    0x3ff1a35beb6fcb75,
    0x3ff1bbe084045cd4,
    0x3ff1d4873168b9aa,
    0x3ff1ed5022fcd91d,
    0x3ff2063b88628cd6,
    0x3ff21f49917ddc96,
    0x3ff2387a6e756238,
    0x3ff251ce4fb2a63f,
    0x3ff26b4565e27cdd,
    0x3ff284dfe1f56381,
    0x3ff29e9df51fdee1,
    0x3ff2b87fd0dad990,
    0x3ff2d285a6e4030b,
    0x3ff2ecafa93e2f56,
    0x3ff306fe0a31b715,
    0x3ff32170fc4cd831,
    0x3ff33c08b26416ff,
    0x3ff356c55f929ff1,
    0x3ff371a7373aa9cb,
    0x3ff38cae6d05d866,
    0x3ff3a7db34e59ff7,
    0x3ff3c32dc313a8e5,
    0x3ff3dea64c123422,
    0x3ff3fa4504ac801c,
    0x3ff4160a21f72e2a,
    0x3ff431f5d950a897,
    0x3ff44e086061892d,
    0x3ff46a41ed1d0057,
    0x3ff486a2b5c13cd0,
    0x3ff4a32af0d7d3de,
    0x3ff4bfdad5362a27,
    0x3ff4dcb299fddd0d,
    0x3ff4f9b2769d2ca7,
    0x3ff516daa2cf6642,
    0x3ff5342b569d4f82,
    0x3ff551a4ca5d920f,
    0x3ff56f4736b527da,
    0x3ff58d12d497c7fd,
    0x3ff5ab07dd485429,
    0x3ff5c9268a5946b7,
    0x3ff5e76f15ad2148,
    0x3ff605e1b976dc09,
    0x3ff6247eb03a5585,
    0x3ff6434634ccc320,
    0x3ff6623882552225,
    0x3ff68155d44ca973,
    0x3ff6a09e667f3bcd,
    0x3ff6c012750bdabf,
    0x3ff6dfb23c651a2f,
    0x3ff6ff7df9519484,
    0x3ff71f75e8ec5f74,
    0x3ff73f9a48a58174,
    0x3ff75feb564267c9,
    0x3ff780694fde5d3f,
    0x3ff7a11473eb0187,
    0x3ff7c1ed0130c132,
    0x3ff7e2f336cf4e62,
    0x3ff80427543e1a12,
    0x3ff82589994cce13,
    0x3ff8471a4623c7ad,
    0x3ff868d99b4492ed,
    0x3ff88ac7d98a6699,
    0x3ff8ace5422aa0db,
    0x3ff8cf3216b5448c,
    0x3ff8f1ae99157736,
    0x3ff9145b0b91ffc6,
    0x3ff93737b0cdc5e5,
    0x3ff95a44cbc8520f,
    0x3ff97d829fde4e50,
    0x3ff9a0f170ca07ba,
    0x3ff9c49182a3f090,
    0x3ff9e86319e32323,
    0x3ffa0c667b5de565,
    0x3ffa309bec4a2d33,
    0x3ffa5503b23e255d,
    0x3ffa799e1330b358,
    0x3ffa9e6b5579fdbf,
    0x3ffac36bbfd3f37a,
    0x3ffae89f995ad3ad,
    0x3ffb0e07298db666,
    0x3ffb33a2b84f15fb,
    0x3ffb59728de5593a,
    0x3ffb7f76f2fb5e47,
    0x3ffba5b030a1064a,
    0x3ffbcc1e904bc1d2,
    0x3ffbf2c25bd71e09,
    0x3ffc199bdd85529c,
    0x3ffc40ab5fffd07a,
    0x3ffc67f12e57d14b,
    0x3ffc8f6d9406e7b5,
    0x3ffcb720dcef9069,
    0x3ffcdf0b555dc3fa,
    0x3ffd072d4a07897c,
    0x3ffd2f87080d89f2,
    0x3ffd5818dcfba487,
    0x3ffd80e316c98398,
    0x3ffda9e603db3285,
    0x3ffdd321f301b460,
    0x3ffdfc97337b9b5f,
    0x3ffe264614f5a129,
    0x3ffe502ee78b3ff6,
    0x3ffe7a51fbc74c83,
    0x3ffea4afa2a490da,
    0x3ffecf482d8e67f1,
    0x3ffefa1bee615a27,
    0x3fff252b376bba97,
    0x3fff50765b6e4540,
    0x3fff7bfdad9cbe14,
    0x3fffa7c1819e90d8,
    0x3fffd3c22b8f71f1,
];

static ERFC_COEFFS: [[u64; 16]; 2] = [
    [
        0x3fec162355429b28,
        0x400d99999999999a,
        0x3fdda951cece2b85,
        0xbff70ef6cff4bcc4,
        0x4003d7f7b3d617de,
        0xc009d0aa47537c51,
        0x4009754ea9a3fcb1,
        0xc0027a5453fcc015,
        0x3ff1ef2e0531aeba,
        0xbfceca090f5a1c06,
        0xbfb7a3cd173a063c,
        0x3fb30fa68a68fddd,
        0x3f555ad9a326993a,
        0xbf907e7b0bb39fbf,
        0x3f52328706c0e950,
        0x3f6d6aa0b7b19cfe,
    ],
    [
        0x401137c8983f8516,
        0x400799999999999a,
        0x3fc05b53aa241333,
        0xbfca3f53872bf870,
        0x3fbde4c30742c9d5,
        0xbfacb24bfa591986,
        0x3f9666aec059ca5f,
        0xbf7a61250eb26b0b,
        0x3f52b28b7924b34d,
        0x3f041b13a9d45013,
        0xbf16dd5e8a273613,
        0x3ef09ce8ea5e8da5,
        0x3ed33923b4102981,
        0xbec1dfd161e3f984,
        0xbe8c87618fcae3b3,
        0x3e8e8a6ffa0ba2c7,
    ],
];

/// Complementary error function
///
/// Max ULP 0.5
pub fn f_erfcf(x: f32) -> f32 {
    let ax = f32::from_bits(x.to_bits() & 0x7fff_ffff);
    let axd = ax as f64;
    let x2 = axd * axd;
    let t = x.to_bits();
    let at = t & 0x7fff_ffff;
    let sgn = t >> 31;
    let i: i64 = (at > 0x40051000) as i64;
    /* for x < -0x1.ea8f94p+1, erfc(x) rounds to 2 (to nearest) */
    if t > 0xc07547cau32 {
        // x < -0x1.ea8f94p+1
        if t >= 0xff800000u32 {
            // -Inf or NaN
            if t == 0xff800000u32 {
                return 2.0;
            } // -Inf
            return x + x; // NaN
        }
        return black_box(2.0) - black_box(f32::from_bits(0x33000000)); // rounds to 2 or nextbelow(2)
    }
    /* at is the absolute value of x
    for x >= 0x1.41bbf8p+3, erfc(x) < 2^-150, thus rounds to 0 or to 2^-149
    depending on the rounding mode */
    if at >= 0x4120ddfcu32 {
        // |x| >= 0x1.41bbf8p+3
        if at >= 0x7f800000u32 {
            // +Inf or NaN
            if at == 0x7f800000u32 {
                return 0.0;
            } // +Inf
            return x + x; // NaN
        }
        // 0x1p-149f * 0.25f rounds to 0 or 2^-149 depending on rounding
        return black_box(f32::from_bits(0x00000001)) * black_box(0.25);
    }
    if at <= 0x3db80000u32 {
        // |x| <= 0x1.7p-4
        if t == 0xb76c9f62u32 {
            // x = -0x1.d93ec4p-17
            return black_box(f32::from_bits(0x3f800085)) + black_box(f32::from_bits(0x33000000)); // exceptional case
        }
        /* for |x| <= 0x1.c5bf88p-26. erfc(x) rounds to 1 (to nearest) */
        if at <= 0x32e2dfc4u32 {
            // |x| <= 0x1.c5bf88p-26
            if at == 0 {
                return 1.0;
            }
            static D: [f32; 2] = [f32::from_bits(0xb2800000), f32::from_bits(0x33000000)];
            return 1.0 + D[sgn as usize];
        }
        /* around 0, erfc(x) behaves as 1 - (odd polynomial) */
        const C: [u64; 5] = [
            0x3ff20dd750429b6d,
            0xbfd812746b03610b,
            0x3fbce2f218831d2f,
            0xbf9b82c609607dcb,
            0x3f7553af09b8008e,
        ];

        let fw0 = f_fmla(x2, f64::from_bits(C[4]), f64::from_bits(C[3]));
        let fw1 = f_fmla(x2, fw0, f64::from_bits(C[2]));
        let fw2 = f_fmla(x2, fw1, f64::from_bits(C[1]));

        let f0 = x as f64 * f_fmla(x2, fw2, f64::from_bits(C[0]));
        return (1.0 - f0) as f32;
    }

    /* now -0x1.ea8f94p+1 <= x <= 0x1.41bbf8p+3, with |x| > 0x1.7p-4 */
    const ILN2: f64 = f64::from_bits(0x3ff71547652b82fe);
    const LN2H: f64 = f64::from_bits(0x3f762e42fefa0000);
    const LN2L: f64 = f64::from_bits(0x3d0cf79abd6f5dc8);

    let jt = dd_fmla(x2, ILN2, -(1024. + f64::from_bits(0x3f70000000000000))).to_bits();
    let j: i64 = ((jt << 12) as i64) >> 48;
    let sf = ((j >> 7) as u64)
        .wrapping_add(0x3ffu64 | (sgn as u64) << 11)
        .wrapping_shl(52);

    const CH: [u64; 4] = [
        0xbfdffffffffff333,
        0x3fc5555555556a14,
        0xbfa55556666659b4,
        0x3f81111074cc7b22,
    ];
    let d = f_fmla(LN2L, j as f64, f_fmla(LN2H, j as f64, x2));
    let d2 = d * d;
    let e0 = f64::from_bits(ERR0[(j & 127) as usize]);

    let fw0 = f_fmla(d, f64::from_bits(CH[3]), f64::from_bits(CH[2]));
    let fw1 = f_fmla(d, f64::from_bits(CH[1]), f64::from_bits(CH[0]));
    let fw2 = f_fmla(d2, fw0, fw1);

    let f = f_fmla(d2, fw2, d);

    let ct = ERFC_COEFFS[i as usize];
    let z = (axd - f64::from_bits(ct[0])) / (axd + f64::from_bits(ct[1]));
    let z2 = z * z;
    let z4 = z2 * z2;
    let z8 = z4 * z4;
    let c = &ct[3..];

    let sw0 = f_fmla(z, f64::from_bits(c[1]), f64::from_bits(c[0]));
    let sw1 = f_fmla(z, f64::from_bits(c[3]), f64::from_bits(c[2]));
    let sw2 = f_fmla(z, f64::from_bits(c[5]), f64::from_bits(c[4]));
    let sw3 = f_fmla(z, f64::from_bits(c[7]), f64::from_bits(c[6]));

    let zw0 = f_fmla(z2, sw1, sw0);
    let zw1 = f_fmla(z2, sw3, sw2);

    let sw4 = f_fmla(z, f64::from_bits(c[9]), f64::from_bits(c[8]));
    let sw5 = f_fmla(z, f64::from_bits(c[11]), f64::from_bits(c[10]));

    let zw2 = f_fmla(z4, zw1, zw0);
    let zw3 = f_fmla(z2, sw5, sw4);
    let zw4 = f_fmla(z4, f64::from_bits(c[12]), zw3);

    let mut s = f_fmla(z8, zw4, zw2);

    s = f_fmla(z, s, f64::from_bits(ct[2]));

    static OFF: [f64; 2] = [0., 2.];

    let r = (f64::from_bits(sf) * f_fmla(-f, e0, e0)) * s;
    let y = OFF[sgn as usize] + r;
    y as f32
}

#[cfg(test)]
mod tests {
    use crate::f_erfcf;

    #[test]
    fn test_erfc() {
        assert_eq!(f_erfcf(0.0), 1.0);
        assert_eq!(f_erfcf(0.5), 0.47950011);
        assert_eq!(f_erfcf(1.0), 0.1572992);
        assert!(f_erfcf(f32::NAN).is_nan());
        assert_eq!(f_erfcf(f32::INFINITY), 0.0);
        assert_eq!(f_erfcf(f32::NEG_INFINITY), 2.0);
    }
}