pxfm/bessel/

i1ef.rs

/*
 * // 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::j0f::j1f_rsqrt;
use crate::common::f_fmla;
use crate::exponents::core_expf;
use crate::polyeval::{f_estrin_polyeval7, f_estrin_polyeval9, f_polyeval10};

/// Modified exponentially scaled Bessel of the first kind of order 1
///
/// Computes exp(-|x|)*I1(x)
///
/// Max ULP 0.5
pub fn f_i1ef(x: f32) -> f32 {
    let ux = x.to_bits().wrapping_shl(1);
    if ux >= 0xffu32 << 24 || ux == 0 {
        // |x| == 0, |x| == inf, x == NaN
        if ux == 0 {
            // |x| == 0
            return 0.;
        }
        if x.is_infinite() {
            return if x.is_sign_positive() { 0. } else { -0. };
        }
        return x + f32::NAN; // |x| == NaN
    }

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

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

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

    if xb <= 0x40f80000u32 {
        // |x| <= 7.75
        let core_expf = core_expf(-f32::from_bits(xb));
        if xb <= 0x34000000u32 {
            // |x| <= f32::EPSILON
            // taylor series for I1(x)/Exp(x) ~ x/2 - x^2/2 + O(x^3)
            #[cfg(any(
                all(
                    any(target_arch = "x86", target_arch = "x86_64"),
                    target_feature = "fma"
                ),
                all(target_arch = "aarch64", target_feature = "neon")
            ))]
            {
                use crate::common::f_fmlaf;
                let half_x = 0.5 * x;
                return f_fmlaf(x, -half_x, half_x);
            }
            #[cfg(not(any(
                all(
                    any(target_arch = "x86", target_arch = "x86_64"),
                    target_feature = "fma"
                ),
                all(target_arch = "aarch64", target_feature = "neon")
            )))]
            {
                let dx = x as f64;
                let half_x = 0.5 * dx;
                return f_fmla(dx, -half_x, half_x) as f32;
            }
        }
        return i1ef_small(f32::from_bits(xb), sign_scale, core_expf) as f32;
    }

    i1ef_asympt(f32::from_bits(xb), sign_scale)
}

/**
Computes
I1(x) = x/2 * (1 + 1 * (x/2)^2 + (x/2)^4 * P((x/2)^2))

Generated by Woflram Mathematica:

```text
<<FunctionApproximations`
ClearAll["Global`*"]
f[x_]:=(BesselI[1,x]*2/x-1-1/2(x/2)^2)/(x/2)^4
g[z_]:=f[2 Sqrt[z]]
{err, approx}=MiniMaxApproximation[g[z],{z,{0.000000001,7.75},6,6},WorkingPrecision->60]
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 i1ef_small(x: f32, sign_scale: f64, core_expf: f64) -> f64 {
    let dx = x as f64;
    let x_over_two = dx * 0.5;
    let x_over_two_sqr = x_over_two * x_over_two;
    let x_over_two_p4 = x_over_two_sqr * x_over_two_sqr;

    let p_num = f_estrin_polyeval7(
        x_over_two_sqr,
        f64::from_bits(0x3fb5555555555555),
        f64::from_bits(0x3f706cdccca396c4),
        f64::from_bits(0x3f23f9e12bdbba92),
        f64::from_bits(0x3ec8e39208e926b2),
        f64::from_bits(0x3e62e53b433c42ff),
        f64::from_bits(0x3def7cb16d10fb46),
        f64::from_bits(0x3d6747cd73d9d783),
    );
    let p_den = f_estrin_polyeval7(
        x_over_two_sqr,
        f64::from_bits(0x3ff0000000000000),
        f64::from_bits(0xbfa2075f77b54885),
        f64::from_bits(0x3f438c6d797c29f5),
        f64::from_bits(0xbeda57e2a258c6da),
        f64::from_bits(0x3e677e777c569432),
        f64::from_bits(0xbdea9212a96babc1),
        f64::from_bits(0x3d5e183186d5d782),
    );
    let p = p_num / p_den;

    let p1 = f_fmla(0.5, x_over_two_sqr, 1.);
    let p2 = f_fmla(x_over_two_p4, p, p1);
    p2 * x_over_two * sign_scale * core_expf
}

/**
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,{2^-33,1/7.75},9,8},WorkingPrecision->60]
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 i1ef_asympt(x: f32, sign_scale: f64) -> f32 {
    let dx = x as f64;
    let recip = 1. / dx;
    let p_num = f_polyeval10(
        recip,
        f64::from_bits(0x3fd9884533d43652),
        f64::from_bits(0xc030686a3694d13c),
        f64::from_bits(0x407344697f45c2ee),
        f64::from_bits(0xc0aa037ee36a8967),
        f64::from_bits(0x40d5b2eab8cf5b17),
        f64::from_bits(0xc0f65addf81dbee8),
        f64::from_bits(0x410afc22ec1f9b8b),
        f64::from_bits(0xc110821dd0fc12b4),
        f64::from_bits(0x40feb3452c93aada),
        f64::from_bits(0xc0c6d04e8c5d02f3),
    );
    let p_den = f_estrin_polyeval9(
        recip,
        f64::from_bits(0x3ff0000000000000),
        f64::from_bits(0xc04460707a9ceed4),
        f64::from_bits(0x4087ac89fcf51e9b),
        f64::from_bits(0xc0bf830689f31b42),
        f64::from_bits(0x40e9c281c367fab2),
        f64::from_bits(0xc109b59ade76eb8c),
        f64::from_bits(0x411d553a9f5673c5),
        f64::from_bits(0xc11f9dbe0665523b),
        f64::from_bits(0x4103b62a329b60d7),
    );
    let z = p_num / p_den;
    let r_sqrt = j1f_rsqrt(dx);
    (z * r_sqrt * sign_scale) as f32
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_i1ef() {
        assert!(f_i1ef(f32::NAN).is_nan());
        assert_eq!(f_i1ef(f32::INFINITY), 0.0);
        assert_eq!(f_i1ef(f32::NEG_INFINITY), 0.0);
        assert_eq!(f_i1ef(0.), 0.);
        assert_eq!(f_i1ef(1.), 0.20791042);
        assert_eq!(f_i1ef(-1.), -0.20791042);
        assert_eq!(f_i1ef(9.), 0.12722498);
        assert_eq!(f_i1ef(-9.), -0.12722498);
        assert_eq!(f_i1ef(0.000000000543453), 2.717265e-10);
        assert_eq!(f_i1ef(-0.000000000543453), -2.717265e-10);
    }
}