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

pxfm/bessel/
i0f.rs
