pxfm/bessel/
i2f.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 8/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::bessel::j0f::j1f_rsqrt;
30use crate::exponents::core_expf;
31use crate::polyeval::{f_estrin_polyeval8, f_estrin_polyeval9};
32
33/// Modified Bessel of the first kind of order 2
34///
35/// ULP 0.5
36pub fn f_i2f(x: f32) -> f32 {
37    let ux = x.to_bits().wrapping_shl(1);
38    if ux >= 0xffu32 << 24 || ux == 0 {
39        // |x| == 0, |x| == inf, |x| == NaN
40        if ux == 0 {
41            // |x| == 0
42            return 0.;
43        }
44        if x.is_infinite() {
45            return f32::INFINITY;
46        }
47        return x + f32::NAN; // x == NaN
48    }
49
50    let xb = x.to_bits() & 0x7fff_ffff;
51
52    if xb >= 0x42b7d875u32 {
53        // |x| >= 91.92277 it's infinity
54        return f32::INFINITY;
55    }
56
57    if xb <= 0x40f80000u32 {
58        // |x| <= 7.75
59        if xb <= 0x34000000u32 {
60            // |x| <= f32::EPSILON
61            let dx = x as f64;
62            const R: f64 = 1. / 8.;
63            return (dx * dx * R) as f32;
64        }
65        return i2f_small(f32::from_bits(xb));
66    }
67
68    i2f_asympt(f32::from_bits(xb))
69}
70
71/**
72Computes
73I2(x) = x^2 * R(x^2)
74
75Generated by Wolfram Mathematica:
76
77```text
78<<FunctionApproximations`
79ClearAll["Global`*"]
80f[x_]:=BesselI[2,x]/x^2
81g[z_]:=f[Sqrt[z]]
82{err,approx}=MiniMaxApproximation[g[z],{z,{0.000000000001,7.75},8,7},WorkingPrecision->75]
83poly=Numerator[approx][[1]];
84coeffs=CoefficientList[poly,z];
85TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
86poly=Denominator[approx][[1]];
87coeffs=CoefficientList[poly,z];
88TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
89```
90**/
91#[inline]
92fn i2f_small(x: f32) -> f32 {
93    let dx = x as f64;
94    let x_sqr = dx * dx;
95
96    let p_num = f_estrin_polyeval9(
97        x_sqr,
98        f64::from_bits(0x3fc0000000000000),
99        f64::from_bits(0x3f831469a38d72c7),
100        f64::from_bits(0x3f2f453dd3dd98f4),
101        f64::from_bits(0x3ec8af52ee8fce9b),
102        f64::from_bits(0x3e5589f2cb4e0ec9),
103        f64::from_bits(0x3dd60fa268a4206d),
104        f64::from_bits(0x3d4ab3091ee18d6b),
105        f64::from_bits(0x3cb1efec43b15186),
106        f64::from_bits(0x3c050992c6e9e63f),
107    );
108    let p_den = f_estrin_polyeval8(
109        x_sqr,
110        f64::from_bits(0x3ff0000000000000),
111        f64::from_bits(0xbf82075d8e3f1476),
112        f64::from_bits(0x3f03ef86564a284b),
113        f64::from_bits(0xbe7c498fab4a57d8),
114        f64::from_bits(0x3dec162ca0f68486),
115        f64::from_bits(0xbd53bb6398461540),
116        f64::from_bits(0x3cb265215261e64a),
117        f64::from_bits(0xbc01cf52cc350e81),
118    );
119    let p = p_num / p_den;
120    (p * x_sqr) as f32
121}
122
123/**
124Asymptotic expansion for I2.
125
126Computes:
127sqrt(x) * exp(-x) * I2(x) = Pn(1/x)/Qn(1/x)
128hence:
129I2(x) = Pn(1/x)/Qm(1/x)*exp(x)/sqrt(x)
130
131Generated by Mathematica:
132```text
133<<FunctionApproximations`
134ClearAll["Global`*"]
135f[x_]:=Sqrt[x] Exp[-x] BesselI[2,x]
136g[z_]:=f[1/z]
137{err,approx}=MiniMaxApproximation[g[z],{z,{1/92.3,1/7.5},8,8},WorkingPrecision->70]
138poly=Numerator[approx][[1]];
139coeffs=CoefficientList[poly,z];
140TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
141poly=Denominator[approx][[1]];
142coeffs=CoefficientList[poly,z];
143TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
144```
145**/
146#[inline]
147fn i2f_asympt(x: f32) -> f32 {
148    let dx = x as f64;
149    let recip = 1. / dx;
150    let p_num = f_estrin_polyeval9(
151        recip,
152        f64::from_bits(0x3fd9884533d45f46),
153        f64::from_bits(0xc02b979526807e1e),
154        f64::from_bits(0x406b1dd3e795bbed),
155        f64::from_bits(0xc09e43629031ec91),
156        f64::from_bits(0x40c48c03a39aec1d),
157        f64::from_bits(0xc0e0f022ccb8807a),
158        f64::from_bits(0x40f0302eeb22a776),
159        f64::from_bits(0xc0f02b01549d38b8),
160        f64::from_bits(0x40dad4e70f2bc264),
161    );
162    let p_den = f_estrin_polyeval9(
163        recip,
164        f64::from_bits(0x3ff0000000000000),
165        f64::from_bits(0xc0405a71a88b191c),
166        f64::from_bits(0x407e19f7d247d098),
167        f64::from_bits(0xc0aeaac6e0ca17fe),
168        f64::from_bits(0x40d2301702f40a98),
169        f64::from_bits(0xc0e7e6c6c01841b3),
170        f64::from_bits(0x40ed67317e9e46cc),
171        f64::from_bits(0xc0d13786aa1ef416),
172        f64::from_bits(0xc0a6c9cfe579ae22),
173    );
174    let z = p_num / p_den;
175
176    let e = core_expf(x);
177    let r_sqrt = j1f_rsqrt(dx);
178    (z * r_sqrt * e) as f32
179}
180
181#[cfg(test)]
182mod tests {
183    use super::*;
184    #[test]
185    fn test_i2f() {
186        assert_eq!(f_i2f(0.), 0.);
187        assert_eq!(f_i2f(f32::INFINITY), f32::INFINITY);
188        assert_eq!(f_i2f(f32::NEG_INFINITY), f32::INFINITY);
189        assert_eq!(f_i2f(1.), 0.13574767);
190        assert_eq!(f_i2f(-1.), 0.13574767);
191        assert_eq!(f_i2f(9.432), 1314.6553);
192        assert_eq!(f_i2f(-9.432), 1314.6553);
193    }
194}
pxfm/bessel/i2f.rs

pxfm/bessel/
i2f.rs
