pxfm/bessel/
i0ef.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 7/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::common::f_fmla;
31use crate::exponents::core_expf;
32use crate::polyeval::{
33    f_estrin_polyeval5, f_estrin_polyeval7, f_estrin_polyeval8, f_polyeval6, f_polyeval10,
34};
35
36/// Modified exponentially scaled Bessel of the first kind of order 0
37///
38/// Computes exp(-|x|)*I0(x)
39///
40/// Max ULP 0.5
41pub fn f_i0ef(x: f32) -> f32 {
42    let ux = x.to_bits().wrapping_shl(1);
43    if ux >= 0xffu32 << 24 || ux == 0 {
44        // |x| == 0, |x| == inf, |x| == NaN
45        if ux == 0 {
46            // |x| == 0
47            return 1.;
48        }
49        if x.is_infinite() {
50            return 0.;
51        }
52        return x + f32::NAN; // x == NaN
53    }
54
55    let xb = x.to_bits() & 0x7fff_ffff;
56
57    if xb <= 0x40f00000u32 {
58        // |x| <= 7.5
59        let core_expf = core_expf(-f32::from_bits(xb));
60        if xb < 0x3f800000u32 {
61            if xb <= 0x34000000u32 {
62                // |x| <= f32::EPSILON
63                // taylor series for I0(x) * exp(-x) ~ 1 - x + O(x^2)
64                return 1. - x;
65            }
66            // |x| < 1
67            return i0f_small(f32::from_bits(xb), core_expf);
68        } else if xb <= 0x40600000u32 {
69            // |x| <= 3.5
70            return i0ef_1_to_3p5(f32::from_bits(xb), core_expf);
71        } else if xb <= 0x40c00000u32 {
72            // |x| <= 6
73            return i0f_3p5_to_6(f32::from_bits(xb), core_expf);
74        }
75        return i0f_6_to_7p5(f32::from_bits(xb), core_expf);
76    }
77
78    i0ef_asympt(f32::from_bits(xb))
79}
80
81/**
82How polynomial is obtained described at [i0f_1_to_7p5].
83
84Computes I0(x) as follows:
85I0(x) = 1 + (x/2)^2 * P(x)
86
87This method valid only [0;1]
88
89Generated by Wolfram Mathematica:
90```text
91<<FunctionApproximations`
92ClearAll["Global`*"]
93f[x_]:=(BesselI[0,x]-1)/(x/2)^2
94g[z_]:=f[2 Sqrt[z]]
95{err, approx}=MiniMaxApproximation[g[z],{z,{0.0000001,1},6,0},WorkingPrecision->60]
96poly=Numerator[approx][[1]];
97coeffs=CoefficientList[poly,z];
98TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
99```
100**/
101#[inline]
102pub(crate) fn i0f_small(x: f32, v_exp: f64) -> f32 {
103    let dx = x as f64;
104    const C: f64 = 1. / 4.;
105    let eval_x = dx * dx * C;
106
107    let p = f_estrin_polyeval7(
108        eval_x,
109        f64::from_bits(0x3ff000000000013a),
110        f64::from_bits(0x3fcffffffffc20b6),
111        f64::from_bits(0x3f9c71c71e6cd6a2),
112        f64::from_bits(0x3f5c71c65b0af15f),
113        f64::from_bits(0x3f1234796fceb081),
114        f64::from_bits(0x3ec0280faf31678c),
115        f64::from_bits(0x3e664fd494223545),
116    );
117    (f_fmla(p, eval_x, 1.) * v_exp) as f32
118}
119
120/**
121Computes I0.
122
123/// Valid only on interval [1;3.5]
124
125as rational approximation I0 = 1 + (x/2)^2 * Pn((x/2)^2)/Qm((x/2)^2))
126
127Generated by Wolram Mathematica:
128```python
129<<FunctionApproximations`
130ClearAll["Global`*"]
131f[x_]:=(BesselI[0,x]-1)/(x/2)^2
132g[z_]:=f[2 Sqrt[z]]
133{err, approx}=MiniMaxApproximation[g[z],{z,{1,3.5},5,4},WorkingPrecision->60]
134poly=Numerator[approx][[1]];
135coeffs=CoefficientList[poly,z];
136TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
137poly=Denominator[approx][[1]];
138coeffs=CoefficientList[poly,z];
139TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
140```
141**/
142#[inline]
143fn i0ef_1_to_3p5(x: f32, v_exp: f64) -> f32 {
144    let dx = x as f64;
145    const C: f64 = 1. / 4.;
146    let eval_x = dx * dx * C;
147
148    let p_num = f_polyeval6(
149        eval_x,
150        f64::from_bits(0x3feffffffffffb69),
151        f64::from_bits(0x3fc9ed7bd9dc97a7),
152        f64::from_bits(0x3f915c14693c842e),
153        f64::from_bits(0x3f45c6dc6a719e42),
154        f64::from_bits(0x3eeacb79eba725f7),
155        f64::from_bits(0x3e7b51e2acfc4355),
156    );
157    let p_den = f_estrin_polyeval5(
158        eval_x,
159        f64::from_bits(0x3ff0000000000000),
160        f64::from_bits(0xbfa84a10988f28eb),
161        f64::from_bits(0x3f50f5599197a4be),
162        f64::from_bits(0xbeea420cf9b13b1b),
163        f64::from_bits(0x3e735d0c1eb6ed7d),
164    );
165
166    (f_fmla(p_num / p_den, eval_x, 1.) * v_exp) as f32
167}
168
169// Valid only on interval [6;7]
170// Generated by Wolfram Mathematica:
171// <<FunctionApproximations`
172// ClearAll["Global`*"]
173// f[x_]:=(BesselI[0,x]-1)/(x/2)^2
174// g[z_]:=f[2 Sqrt[z]]
175// {err, approx}=MiniMaxApproximation[g[z],{z,{6,7},7,6},WorkingPrecision->60]
176// poly=Numerator[approx][[1]];
177// coeffs=CoefficientList[poly,z];
178// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
179// poly=Denominator[approx][[1]];
180// coeffs=CoefficientList[poly,z];
181// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
182#[inline]
183fn i0f_6_to_7p5(x: f32, v_exp: f64) -> f32 {
184    let dx = x as f64;
185    const C: f64 = 1. / 4.;
186    let eval_x = dx * dx * C;
187
188    let p_num = f_estrin_polyeval8(
189        eval_x,
190        f64::from_bits(0x3fefffffffffff7d),
191        f64::from_bits(0x3fcb373b00569ccf),
192        f64::from_bits(0x3f939069c3363b81),
193        f64::from_bits(0x3f4c2095c90c66b3),
194        f64::from_bits(0x3ef6713f648413db),
195        f64::from_bits(0x3e947efa2f9936b4),
196        f64::from_bits(0x3e2486a182f49420),
197        f64::from_bits(0x3da213034a33de33),
198    );
199    let p_den = f_estrin_polyeval7(
200        eval_x,
201        f64::from_bits(0x3ff0000000000000),
202        f64::from_bits(0xbfa32313fea59d9e),
203        f64::from_bits(0x3f460594c2ec6706),
204        f64::from_bits(0xbedf725fb714690f),
205        f64::from_bits(0x3e6d9cb39b19555c),
206        f64::from_bits(0xbdf1900e3abcb7a6),
207        f64::from_bits(0x3d64a21a2ea78ef6),
208    );
209
210    (f_fmla(p_num / p_den, eval_x, 1.) * v_exp) as f32
211}
212
213// Valid only on interval [3.5;6]
214// Generated in Wolfram Mathematica:
215// <<FunctionApproximations`
216// ClearAll["Global`*"]
217// f[x_]:=(BesselI[0,x]-1)/(x/2)^2
218// g[z_]:=f[2 Sqrt[z]]
219// {err, approx}=MiniMaxApproximation[g[z],{z,{3.5,6},5,5},WorkingPrecision->60]
220// poly=Numerator[approx][[1]];
221// coeffs=CoefficientList[poly,z];
222// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
223// poly=Denominator[approx][[1]];
224// coeffs=CoefficientList[poly,z];
225// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
226#[inline]
227fn i0f_3p5_to_6(x: f32, v_exp: f64) -> f32 {
228    let dx = x as f64;
229    const C: f64 = 1. / 4.;
230    let eval_x = dx * dx * C;
231
232    let p_num = f_polyeval6(
233        eval_x,
234        f64::from_bits(0x3feffffffffd9550),
235        f64::from_bits(0x3fc97e18ee033fb4),
236        f64::from_bits(0x3f90b3199079bce1),
237        f64::from_bits(0x3f442c300a425372),
238        f64::from_bits(0x3ee7831030ae18ca),
239        f64::from_bits(0x3e76387d67354932),
240    );
241    let p_den = f_polyeval6(
242        eval_x,
243        f64::from_bits(0x3ff0000000000000),
244        f64::from_bits(0xbfaa079c484e406a),
245        f64::from_bits(0x3f5452098f1556fb),
246        f64::from_bits(0xbef33efb4a8128ac),
247        f64::from_bits(0x3e865996e19448ca),
248        f64::from_bits(0xbe09acbb64533c3e),
249    );
250
251    (f_fmla(p_num / p_den, eval_x, 1.) * v_exp) as f32
252}
253
254/**
255Asymptotic expansion for I0.
256
257Computes:
258sqrt(x) * exp(-x) * I0(x) = Pn(1/x)/Qn(1/x)
259hence:
260I0(x)exp(-x) = Pn(1/x)/Qm(1/x)/sqrt(x)
261
262Generated by Mathematica:
263```text
264<<FunctionApproximations`
265ClearAll["Global`*"]
266f[x_]:=Sqrt[x] Exp[-x] BesselI[0,x]
267g[z_]:=f[1/z]
268{err,approx}=MiniMaxApproximation[g[z],{z,{2^-33,1/7.5},9,9},WorkingPrecision->70]
269num=Numerator[approx][[1]];
270den=Denominator[approx][[1]];
271poly=num;
272coeffs=CoefficientList[poly,z];
273TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
274poly=den;
275coeffs=CoefficientList[poly,z];
276TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
277```
278**/
279#[inline]
280fn i0ef_asympt(x: f32) -> f32 {
281    let dx = x as f64;
282    let recip = 1. / dx;
283    let p_num = f_polyeval10(
284        recip,
285        f64::from_bits(0x3fd9884533d4364f),
286        f64::from_bits(0xc02ed6c9269921a7),
287        f64::from_bits(0x4070ee77ffed64a5),
288        f64::from_bits(0xc0a4ffd558b06889),
289        f64::from_bits(0x40cf2633e2840f6f),
290        f64::from_bits(0xc0ea813a9ba42b84),
291        f64::from_bits(0x40f569bf5d63eb8c),
292        f64::from_bits(0xc0b3138874cdd180),
293        f64::from_bits(0xc0fa3152ed485937),
294        f64::from_bits(0x40ddaccbed454f47),
295    );
296    let p_den = f_polyeval10(
297        recip,
298        f64::from_bits(0x3ff0000000000000),
299        f64::from_bits(0xc0436352c350b88c),
300        f64::from_bits(0x40855eaa17b05edd),
301        f64::from_bits(0xc0baa46f155bd266),
302        f64::from_bits(0x40e3e9fd90a2e695),
303        f64::from_bits(0xc1012dc621dfc1e8),
304        f64::from_bits(0x410cafeea713e8ce),
305        f64::from_bits(0xc0e0a3ee0077d7f7),
306        f64::from_bits(0xc110bcced6a39e9e),
307        f64::from_bits(0x40f9a1e4a91be4d6),
308    );
309    let z = p_num / p_den;
310    let r_sqrt = j1f_rsqrt(dx);
311    (z * r_sqrt) as f32
312}
313
314#[cfg(test)]
315mod tests {
316    use super::*;
317
318    #[test]
319    fn test_i0f() {
320        assert!(f_i0ef(f32::NAN).is_nan());
321        assert_eq!(f_i0ef(f32::NEG_INFINITY), 0.);
322        assert_eq!(f_i0ef(f32::INFINITY), 0.);
323        assert_eq!(f_i0ef(1.), 0.4657596);
324        assert_eq!(f_i0ef(5.), 0.1835408);
325        assert_eq!(f_i0ef(16.), 0.100544125);
326        assert_eq!(f_i0ef(32.), 0.070804186);
327        assert_eq!(f_i0ef(92.0), 0.04164947);
328        assert_eq!(f_i0ef(0.), 1.0);
329        assert_eq!(f_i0ef(28.), 0.075736605);
330        assert_eq!(f_i0ef(-28.), 0.075736605);
331        assert_eq!(f_i0ef(-32.), 0.070804186);
332        assert_eq!(f_i0ef(-92.0), 0.04164947);
333        assert_eq!(f_i0ef(-0.), 1.0);
334    }
335}
pxfm/bessel/i0ef.rs

pxfm/bessel/
i0ef.rs
