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::bessel::i0::bessel_rsqrt_hard;
use crate::bessel::i0_exp;
use crate::bessel::i1::i1_0_to_7p75;
use crate::common::f_fmla;
use crate::double_double::DoubleDouble;
use crate::dyadic_float::{DyadicFloat128, DyadicSign};

/// Modified exponentially scaled Bessel of the first kind of order 1
///
/// Computes exp(-|x|)*I1(x)
pub fn f_i1e(x: f64) -> f64 {
    let ux = x.to_bits().wrapping_shl(1);

    if ux >= 0x7ffu64 << 53 || ux <= 0x7960000000000000u64 {
        // |x| <= f64::EPSILON, |x| == inf, x == NaN
        if ux <= 0x760af31dc4611874u64 {
            // |x| <= 2.2204460492503131e-24
            return x * 0.5;
        }
        if ux <= 0x7960000000000000u64 {
            // |x| <= f64::EPSILON
            // Power series of I1(x)*exp(-|x|) ~ x/2 - x^2/2 + O(x^3)
            return f_fmla(x, -x * 0.5, x * 0.5);
        }
        if x.is_infinite() {
            return 0.;
        }
        return x + f64::NAN; // x == NaN
    }

    let xb = x.to_bits() & 0x7fff_ffff_ffff_ffff;

    static SIGN: [f64; 2] = [1., -1.];

    let sign_scale = SIGN[x.is_sign_negative() as usize];

    if xb < 0x401f000000000000u64 {
        // |x| <= 7.75
        let v_exp = i0_exp(-f64::from_bits(xb));
        let vi1 = i1_0_to_7p75(f64::from_bits(xb));
        let r = DoubleDouble::quick_mult(vi1, v_exp);
        return f64::copysign(r.to_f64(), sign_scale);
    }

    i1e_asympt(f64::from_bits(xb), sign_scale)
}

/**
Asymptotic expansion for I1.

Computes:
sqrt(x) * exp(-x) * I1(x) = Pn(1/x)/Qn(1/x)
hence:
I1(x)exp(-|x|) = Pn(1/x)/Qm(1/x)/sqrt(x)

Generated by Wolfram Mathematica:
```text
<<FunctionApproximations`
ClearAll["Global`*"]
f[x_]:=Sqrt[x] Exp[-x] BesselI[1,x]
g[z_]:=f[1/z]
{err,approx}=MiniMaxApproximation[g[z],{z,{1/713.98,1/7.75},11,11},WorkingPrecision->120]
poly=Numerator[approx][[1]];
coeffs=CoefficientList[poly,z];
TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
poly=Denominator[approx][[1]];
coeffs=CoefficientList[poly,z];
TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
```
**/
#[inline]
fn i1e_asympt(x: f64, sign_scale: f64) -> f64 {
    let dx = x;
    let recip = DoubleDouble::from_quick_recip(x);
    const P: [(u64, u64); 12] = [
        (0xbc73a823f28a2f5e, 0x3fd9884533d43651),
        (0x3cc0d5bb78e674b3, 0xc0354325c8029263),
        (0x3d20e1155aaaa283, 0x4080c09b027c46a4),
        (0xbd5b90dcf81b99c1, 0xc0bfc1311090c839),
        (0xbd98f2fda9e8fa1b, 0x40f3bb81bb190ae2),
        (0xbdcec960752b60da, 0xc1207c0bbbc31cd9),
        (0x3dd3c9a299c9c41f, 0x414253e25c4584af),
        (0xbde82e7b9a3e1acc, 0xc159a656aece377c),
        (0x3e0d3d30d701a8ab, 0x416398df24c74ef2),
        (0xbdf57b85ab7006e2, 0xc151fd119be1702b),
        (0x3dd760928f4515fd, 0xc1508327e42639bc),
        (0x3dc09e71bc648589, 0x4143e4933afa621c),
    ];

    let x2 = DoubleDouble::quick_mult(recip, recip);
    let x4 = DoubleDouble::quick_mult(x2, x2);
    let x8 = DoubleDouble::quick_mult(x4, x4);

    let e0 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(P[1]),
        DoubleDouble::from_bit_pair(P[0]),
    );
    let e1 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(P[3]),
        DoubleDouble::from_bit_pair(P[2]),
    );
    let e2 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(P[5]),
        DoubleDouble::from_bit_pair(P[4]),
    );
    let e3 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(P[7]),
        DoubleDouble::from_bit_pair(P[6]),
    );
    let e4 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(P[9]),
        DoubleDouble::from_bit_pair(P[8]),
    );
    let e5 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(P[11]),
        DoubleDouble::from_bit_pair(P[10]),
    );

    let f0 = DoubleDouble::mul_add(x2, e1, e0);
    let f1 = DoubleDouble::mul_add(x2, e3, e2);
    let f2 = DoubleDouble::mul_add(x2, e5, e4);

    let g0 = DoubleDouble::mul_add(x4, f1, f0);

    let p_num = DoubleDouble::mul_add(x8, f2, g0);

    const Q: [(u64, u64); 12] = [
        (0x0000000000000000, 0x3ff0000000000000),
        (0xbcb334d5a476d9ad, 0xc04a75f94c1a0c1a),
        (0xbd324d58ed98bfae, 0x4094b00e60301c42),
        (0x3d7c8725666c4360, 0xc0d36b9d28d45928),
        (0x3d7f8457c2945822, 0x4107d6c398a174ed),
        (0x3dbc655ea216594b, 0xc1339393e6776e38),
        (0xbdebb5dffef78272, 0x415537198d23f6a1),
        (0xbdb577f8abad883e, 0xc16c6c399dcd6949),
        (0x3e14261c5362f109, 0x4173c02446576949),
        (0x3dc382ededad42c5, 0xc1547dff5770f4ec),
        (0xbe05c0f74d4c7956, 0xc165c88046952562),
        (0xbdbf9145927aa2c7, 0x414395e46bc45d5b),
    ];

    let e0 = DoubleDouble::mul_add_f64(
        recip,
        DoubleDouble::from_bit_pair(Q[1]),
        f64::from_bits(0x3ff0000000000000),
    );
    let e1 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(Q[3]),
        DoubleDouble::from_bit_pair(Q[2]),
    );
    let e2 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(Q[5]),
        DoubleDouble::from_bit_pair(Q[4]),
    );
    let e3 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(Q[7]),
        DoubleDouble::from_bit_pair(Q[6]),
    );
    let e4 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(Q[9]),
        DoubleDouble::from_bit_pair(Q[8]),
    );
    let e5 = DoubleDouble::mul_add(
        recip,
        DoubleDouble::from_bit_pair(Q[11]),
        DoubleDouble::from_bit_pair(Q[10]),
    );

    let f0 = DoubleDouble::mul_add(x2, e1, e0);
    let f1 = DoubleDouble::mul_add(x2, e3, e2);
    let f2 = DoubleDouble::mul_add(x2, e5, e4);

    let g0 = DoubleDouble::mul_add(x4, f1, f0);

    let p_den = DoubleDouble::mul_add(x8, f2, g0);

    let z = DoubleDouble::div(p_num, p_den);

    let r_sqrt = DoubleDouble::from_rsqrt_fast(dx);

    let r = z * r_sqrt;

    let err = f_fmla(
        r.hi,
        f64::from_bits(0x3c40000000000000), // 2^-59
        f64::from_bits(0x3ba0000000000000), // 2^-69
    );
    let up = r.hi + (r.lo + err);
    let lb = r.hi + (r.lo - err);
    if up == lb {
        return f64::copysign(r.to_f64(), sign_scale);
    }
    i1e_asympt_hard(x, sign_scale)
}

/**
Asymptotic expansion for I1.

Computes:
sqrt(x) * exp(-x) * I1(x) = Pn(1/x)/Qn(1/x)
hence:
I1(x)exp(-|x|) = Pn(1/x)/Qm(1/x)/sqrt(x)

Generated by Wolfram Mathematica:
```text
<<FunctionApproximations`
ClearAll["Global`*"]
f[x_]:=Sqrt[x] Exp[-x] BesselI[1,x]
g[z_]:=f[1/z]
{err,approx}=MiniMaxApproximation[g[z],{z,{1/713.98,1/7.75},15,15},WorkingPrecision->120]
poly=Numerator[approx][[1]];
coeffs=CoefficientList[poly,z];
TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
poly=Denominator[approx][[1]];
coeffs=CoefficientList[poly,z];
TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
```
**/
#[cold]
#[inline(never)]
fn i1e_asympt_hard(x: f64, sign_scale: f64) -> f64 {
    static P: [DyadicFloat128; 16] = [
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -129,
            mantissa: 0xcc42299e_a1b28468_bea7da47_28f13acc_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -124,
            mantissa: 0xda979406_3df6e66f_cf31c3f5_f194b48c_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -120,
            mantissa: 0xd60b7b96_c958929b_cabe1d8c_3d874767_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -113,
            mantissa: 0xd27aad9a_8fb38d56_46ab4510_8479306e_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -108,
            mantissa: 0xe0167305_c451bd1f_d2f17b68_5c62e2ff_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -103,
            mantissa: 0x8f6d238f_c80d8e4a_08c130f6_24e1c925_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -100,
            mantissa: 0xfe32280f_2ea99024_d9924472_92d7ac8f_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -96,
            mantissa: 0xa48815ac_d265609f_da4ace94_811390b2_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -93,
            mantissa: 0x9ededfe5_833b4cc1_731efd5c_f8729c6c_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -91,
            mantissa: 0xe5b43203_2784ae6a_f7458556_0a8308ea_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -89,
            mantissa: 0xf5df279a_3fb4ef60_8d10adee_7dd2f47b_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -87,
            mantissa: 0xbdb59963_7a757ed1_87280e0e_7f93ca2b_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -86,
            mantissa: 0xc87fdea5_53177ca8_c91de5fb_3f8f78d3_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -85,
            mantissa: 0x846d16a7_4663ef5d_ad42d599_5bc726b8_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -86,
            mantissa: 0xb3ed2055_74262d95_389f33e4_2ac3774a_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -88,
            mantissa: 0xa511aa32_c18c34e4_3d029a90_a71b7a55_u128,
        },
    ];
    static Q: [DyadicFloat128; 16] = [
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -127,
            mantissa: 0x80000000_00000000_00000000_00000000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -122,
            mantissa: 0x877b771a_ad8f5fd3_5aacf5f9_f04ee9de_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -118,
            mantissa: 0x89475ecd_9c84361e_800c8a3a_c8af23bf_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -111,
            mantissa: 0x837d1196_cf2723f1_23b54da8_225efe05_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -106,
            mantissa: 0x8ae3aecb_15355751_a9ee12e5_a4dd9dde_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -102,
            mantissa: 0xb0886afa_bc13f996_ab45d252_75c8f587_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -98,
            mantissa: 0x9b37d7cd_b114b86b_7d14a389_26599aa1_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -95,
            mantissa: 0xc716bf54_09d5dd9f_bc16679b_93aaeca4_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -92,
            mantissa: 0xbe0cd82e_c8af8371_ab028ed9_c7902dd2_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -89,
            mantissa: 0x875f8d91_8ef5d434_a39d00f9_2aed3d2a_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -87,
            mantissa: 0x8e030781_5aa4ce7f_70156b82_8b216b7c_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -86,
            mantissa: 0xd4dd2687_92646fbd_5ea2d422_da64fc0b_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -85,
            mantissa: 0xd6d72ab3_64b4a827_0499af0f_13a51a80_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -84,
            mantissa: 0x828f4e8b_728747a9_2cebe54a_810e2681_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -85,
            mantissa: 0x91570096_36a3fcfb_6b936d44_68dda1be_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -89,
            mantissa: 0xf082ad00_86024ed4_dd31613b_ec41e3f8_u128,
        },
    ];

    let recip = DyadicFloat128::accurate_reciprocal(x);

    let mut p_num = P[15];
    for i in (0..15).rev() {
        p_num = recip * p_num + P[i];
    }
    let mut p_den = Q[15];
    for i in (0..15).rev() {
        p_den = recip * p_den + Q[i];
    }
    let z = p_num * p_den.reciprocal();
    let r_sqrt = bessel_rsqrt_hard(x, recip);
    (z * r_sqrt).fast_as_f64() * sign_scale
}

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn test_fi1e() {
        assert_eq!(f_i1e(f64::EPSILON), 1.1102230246251563e-16);
        assert_eq!(f_i1e(7.750000000757874), 0.13605110007443239);
        assert_eq!(f_i1e(7.482812501363189), 0.13818116726273896);
        assert_eq!(f_i1e(-7.750000000757874), -0.13605110007443239);
        assert_eq!(f_i1e(-7.482812501363189), -0.13818116726273896);
        assert!(f_i1e(f64::NAN).is_nan());
        assert_eq!(f_i1e(f64::INFINITY), 0.);
        assert_eq!(f_i1e(f64::NEG_INFINITY), 0.);
        assert_eq!(f_i1e(0.01), 0.004950311047118276);
        assert_eq!(f_i1e(-0.01), -0.004950311047118276);
        assert_eq!(f_i1e(-9.01), -0.12716101566063667);
        assert_eq!(f_i1e(9.01), 0.12716101566063667);
        assert_eq!(f_i1e(763.), 0.014435579051182581);
        assert_eq!(i1e_asympt_hard(9.01, 1.), 0.12716101566063667);
    }
}