pxfm/bessel/
i1f.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::{f_estrin_polyeval7, f_estrin_polyeval9, f_polyeval10};
33
34/// Modified Bessel of the first kind of order 1
35///
36/// Max ULP 0.5
37pub fn f_i1f(x: f32) -> f32 {
38    let ux = x.to_bits().wrapping_shl(1);
39    if ux >= 0xffu32 << 24 || ux == 0 {
40        // |x| == 0, |x| == inf, |x| == NaN
41        if ux == 0 {
42            // |x| == 0
43            return 0.;
44        }
45        if x.is_infinite() {
46            return if x.is_sign_positive() {
47                f32::INFINITY
48            } else {
49                f32::NEG_INFINITY
50            };
51        }
52        return x + f32::NAN; // x == NaN
53    }
54
55    let xb = x.to_bits() & 0x7fff_ffff;
56
57    if xb > 0x42b7d001 {
58        // x > 91.906261
59        return if x.is_sign_negative() {
60            f32::NEG_INFINITY
61        } else {
62            f32::INFINITY
63        };
64    }
65
66    static SIGN: [f64; 2] = [1., -1.];
67
68    let sign_scale = SIGN[x.is_sign_negative() as usize];
69
70    if xb <= 0x40f80000u32 {
71        // |x| <= 7.75
72        if xb <= 0x34000000u32 {
73            // |x| <= f32::EPSILON
74            // taylor series for I1(x) ~ x/2 + O(x^3)
75            return x * 0.5;
76        }
77        return i1f_small(f32::from_bits(xb), sign_scale) as f32;
78    }
79
80    i1f_asympt(f32::from_bits(xb), sign_scale)
81}
82
83/**
84Computes
85I1(x) = x/2 * (1 + 1 * (x/2)^2 + (x/2)^4 * P((x/2)^2))
86
87Generated by Woflram Mathematica:
88
89```text
90<<FunctionApproximations`
91ClearAll["Global`*"]
92f[x_]:=(BesselI[1,x]*2/x-1-1/2(x/2)^2)/(x/2)^4
93g[z_]:=f[2 Sqrt[z]]
94{err, approx}=MiniMaxApproximation[g[z],{z,{0.000000001,7.75},6,6},WorkingPrecision->60]
95poly=Numerator[approx][[1]];
96coeffs=CoefficientList[poly,z];
97TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
98poly=Denominator[approx][[1]];
99coeffs=CoefficientList[poly,z];
100TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
101```
102**/
103#[inline]
104fn i1f_small(x: f32, sign_scale: f64) -> f64 {
105    let dx = x as f64;
106    let x_over_two = dx * 0.5;
107    let x_over_two_sqr = x_over_two * x_over_two;
108    let x_over_two_p4 = x_over_two_sqr * x_over_two_sqr;
109
110    let p_num = f_estrin_polyeval7(
111        x_over_two_sqr,
112        f64::from_bits(0x3fb5555555555555),
113        f64::from_bits(0x3f706cdccca396c4),
114        f64::from_bits(0x3f23f9e12bdbba92),
115        f64::from_bits(0x3ec8e39208e926b2),
116        f64::from_bits(0x3e62e53b433c42ff),
117        f64::from_bits(0x3def7cb16d10fb46),
118        f64::from_bits(0x3d6747cd73d9d783),
119    );
120    let p_den = f_estrin_polyeval7(
121        x_over_two_sqr,
122        f64::from_bits(0x3ff0000000000000),
123        f64::from_bits(0xbfa2075f77b54885),
124        f64::from_bits(0x3f438c6d797c29f5),
125        f64::from_bits(0xbeda57e2a258c6da),
126        f64::from_bits(0x3e677e777c569432),
127        f64::from_bits(0xbdea9212a96babc1),
128        f64::from_bits(0x3d5e183186d5d782),
129    );
130    let p = p_num / p_den;
131
132    let p1 = f_fmla(0.5, x_over_two_sqr, 1.);
133    let p2 = f_fmla(x_over_two_p4, p, p1);
134    p2 * x_over_two * sign_scale
135}
136
137/**
138Asymptotic expansion for I1.
139
140Computes:
141sqrt(x) * exp(-x) * I1(x) = Pn(1/x)/Qn(1/x)
142hence:
143I1(x) = Pn(1/x)/Qm(1/x)*exp(x)/sqrt(x)
144
145Generated by Wolfram Mathematica:
146```text
147<<FunctionApproximations`
148ClearAll["Global`*"]
149f[x_]:=Sqrt[x] Exp[-x] BesselI[1,x]
150g[z_]:=f[1/z]
151{err, approx,err1}=MiniMaxApproximation[g[z],{z,{1/91.9,1/7.75},9,8},WorkingPrecision->60]
152poly=Numerator[approx];
153coeffs=CoefficientList[poly,z];
154TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
155poly=Denominator[approx];
156coeffs=CoefficientList[poly,z];
157TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
158```
159**/
160#[inline]
161fn i1f_asympt(x: f32, sign_scale: f64) -> f32 {
162    let dx = x as f64;
163    let recip = 1. / dx;
164    let p_num = f_polyeval10(
165        recip,
166        f64::from_bits(0x3fd9884533d43711),
167        f64::from_bits(0xc0309c047537243a),
168        f64::from_bits(0x4073bdb14a29bf68),
169        f64::from_bits(0xc0aaf9eca14d15af),
170        f64::from_bits(0x40d6c629318a9e42),
171        f64::from_bits(0xc0f7bee33088a4b0),
172        f64::from_bits(0x410d018cef093ee2),
173        f64::from_bits(0xc111f32b325d3fe4),
174        f64::from_bits(0x4100dddad80e0b42),
175        f64::from_bits(0xc0c96006c91a00e2),
176    );
177    let p_den = f_estrin_polyeval9(
178        recip,
179        f64::from_bits(0x3ff0000000000000),
180        f64::from_bits(0xc044a11d10bae889),
181        f64::from_bits(0x408843069497d993),
182        f64::from_bits(0xc0c058710de4b9b9),
183        f64::from_bits(0x40eb0d97f71420ae),
184        f64::from_bits(0xc10b55d181ef9ea1),
185        f64::from_bits(0x411f9413e1932a48),
186        f64::from_bits(0xc1213bff5bc7d2d6),
187        f64::from_bits(0x4105c53e92d9b9c0),
188    );
189    let z = p_num / p_den;
190    let e = core_expf(x);
191    let r_sqrt = j1f_rsqrt(dx);
192    (z * r_sqrt * e * sign_scale) as f32
193}
194
195#[cfg(test)]
196mod tests {
197    use super::*;
198
199    #[test]
200    fn test_i1f() {
201        assert!(f_i1f(f32::NAN).is_nan());
202        assert!(f_i1f(f32::INFINITY).is_infinite());
203        assert!(f_i1f(f32::NEG_INFINITY).is_infinite());
204        assert_eq!(f_i1f(0.), 0.);
205        assert_eq!(f_i1f(1.), 0.5651591);
206        assert_eq!(f_i1f(-1.), -0.5651591);
207        assert_eq!(f_i1f(9.), 1030.9147);
208        assert_eq!(f_i1f(-9.), -1030.9147);
209    }
210}
pxfm/bessel/i1f.rs

pxfm/bessel/
i1f.rs
