pxfm/bessel/
i1e.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::i0::bessel_rsqrt_hard;
30use crate::bessel::i0_exp;
31use crate::bessel::i1::i1_0_to_7p75;
32use crate::common::f_fmla;
33use crate::double_double::DoubleDouble;
34use crate::dyadic_float::{DyadicFloat128, DyadicSign};
35
36/// Modified exponentially scaled Bessel of the first kind of order 1
37///
38/// Computes exp(-|x|)*I1(x)
39pub fn f_i1e(x: f64) -> f64 {
40    let ux = x.to_bits().wrapping_shl(1);
41
42    if ux >= 0x7ffu64 << 53 || ux <= 0x7960000000000000u64 {
43        // |x| <= f64::EPSILON, |x| == inf, x == NaN
44        if ux <= 0x760af31dc4611874u64 {
45            // |x| <= 2.2204460492503131e-24
46            return x * 0.5;
47        }
48        if ux <= 0x7960000000000000u64 {
49            // |x| <= f64::EPSILON
50            // Power series of I1(x)*exp(-|x|) ~ x/2 - x^2/2 + O(x^3)
51            return f_fmla(x, -x * 0.5, x * 0.5);
52        }
53        if x.is_infinite() {
54            return 0.;
55        }
56        return x + f64::NAN; // x == NaN
57    }
58
59    let xb = x.to_bits() & 0x7fff_ffff_ffff_ffff;
60
61    static SIGN: [f64; 2] = [1., -1.];
62
63    let sign_scale = SIGN[x.is_sign_negative() as usize];
64
65    if xb < 0x401f000000000000u64 {
66        // |x| <= 7.75
67        let v_exp = i0_exp(-f64::from_bits(xb));
68        let vi1 = i1_0_to_7p75(f64::from_bits(xb));
69        let r = DoubleDouble::quick_mult(vi1, v_exp);
70        return f64::copysign(r.to_f64(), sign_scale);
71    }
72
73    i1e_asympt(f64::from_bits(xb), sign_scale)
74}
75
76/**
77Asymptotic expansion for I1.
78
79Computes:
80sqrt(x) * exp(-x) * I1(x) = Pn(1/x)/Qn(1/x)
81hence:
82I1(x)exp(-|x|) = Pn(1/x)/Qm(1/x)/sqrt(x)
83
84Generated by Wolfram Mathematica:
85```text
86<<FunctionApproximations`
87ClearAll["Global`*"]
88f[x_]:=Sqrt[x] Exp[-x] BesselI[1,x]
89g[z_]:=f[1/z]
90{err,approx}=MiniMaxApproximation[g[z],{z,{1/713.98,1/7.75},11,11},WorkingPrecision->120]
91poly=Numerator[approx][[1]];
92coeffs=CoefficientList[poly,z];
93TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
94poly=Denominator[approx][[1]];
95coeffs=CoefficientList[poly,z];
96TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
97```
98**/
99#[inline]
100fn i1e_asympt(x: f64, sign_scale: f64) -> f64 {
101    let dx = x;
102    let recip = DoubleDouble::from_quick_recip(x);
103    const P: [(u64, u64); 12] = [
104        (0xbc73a823f28a2f5e, 0x3fd9884533d43651),
105        (0x3cc0d5bb78e674b3, 0xc0354325c8029263),
106        (0x3d20e1155aaaa283, 0x4080c09b027c46a4),
107        (0xbd5b90dcf81b99c1, 0xc0bfc1311090c839),
108        (0xbd98f2fda9e8fa1b, 0x40f3bb81bb190ae2),
109        (0xbdcec960752b60da, 0xc1207c0bbbc31cd9),
110        (0x3dd3c9a299c9c41f, 0x414253e25c4584af),
111        (0xbde82e7b9a3e1acc, 0xc159a656aece377c),
112        (0x3e0d3d30d701a8ab, 0x416398df24c74ef2),
113        (0xbdf57b85ab7006e2, 0xc151fd119be1702b),
114        (0x3dd760928f4515fd, 0xc1508327e42639bc),
115        (0x3dc09e71bc648589, 0x4143e4933afa621c),
116    ];
117
118    let x2 = DoubleDouble::quick_mult(recip, recip);
119    let x4 = DoubleDouble::quick_mult(x2, x2);
120    let x8 = DoubleDouble::quick_mult(x4, x4);
121
122    let e0 = DoubleDouble::mul_add(
123        recip,
124        DoubleDouble::from_bit_pair(P[1]),
125        DoubleDouble::from_bit_pair(P[0]),
126    );
127    let e1 = DoubleDouble::mul_add(
128        recip,
129        DoubleDouble::from_bit_pair(P[3]),
130        DoubleDouble::from_bit_pair(P[2]),
131    );
132    let e2 = DoubleDouble::mul_add(
133        recip,
134        DoubleDouble::from_bit_pair(P[5]),
135        DoubleDouble::from_bit_pair(P[4]),
136    );
137    let e3 = DoubleDouble::mul_add(
138        recip,
139        DoubleDouble::from_bit_pair(P[7]),
140        DoubleDouble::from_bit_pair(P[6]),
141    );
142    let e4 = DoubleDouble::mul_add(
143        recip,
144        DoubleDouble::from_bit_pair(P[9]),
145        DoubleDouble::from_bit_pair(P[8]),
146    );
147    let e5 = DoubleDouble::mul_add(
148        recip,
149        DoubleDouble::from_bit_pair(P[11]),
150        DoubleDouble::from_bit_pair(P[10]),
151    );
152
153    let f0 = DoubleDouble::mul_add(x2, e1, e0);
154    let f1 = DoubleDouble::mul_add(x2, e3, e2);
155    let f2 = DoubleDouble::mul_add(x2, e5, e4);
156
157    let g0 = DoubleDouble::mul_add(x4, f1, f0);
158
159    let p_num = DoubleDouble::mul_add(x8, f2, g0);
160
161    const Q: [(u64, u64); 12] = [
162        (0x0000000000000000, 0x3ff0000000000000),
163        (0xbcb334d5a476d9ad, 0xc04a75f94c1a0c1a),
164        (0xbd324d58ed98bfae, 0x4094b00e60301c42),
165        (0x3d7c8725666c4360, 0xc0d36b9d28d45928),
166        (0x3d7f8457c2945822, 0x4107d6c398a174ed),
167        (0x3dbc655ea216594b, 0xc1339393e6776e38),
168        (0xbdebb5dffef78272, 0x415537198d23f6a1),
169        (0xbdb577f8abad883e, 0xc16c6c399dcd6949),
170        (0x3e14261c5362f109, 0x4173c02446576949),
171        (0x3dc382ededad42c5, 0xc1547dff5770f4ec),
172        (0xbe05c0f74d4c7956, 0xc165c88046952562),
173        (0xbdbf9145927aa2c7, 0x414395e46bc45d5b),
174    ];
175
176    let e0 = DoubleDouble::mul_add_f64(
177        recip,
178        DoubleDouble::from_bit_pair(Q[1]),
179        f64::from_bits(0x3ff0000000000000),
180    );
181    let e1 = DoubleDouble::mul_add(
182        recip,
183        DoubleDouble::from_bit_pair(Q[3]),
184        DoubleDouble::from_bit_pair(Q[2]),
185    );
186    let e2 = DoubleDouble::mul_add(
187        recip,
188        DoubleDouble::from_bit_pair(Q[5]),
189        DoubleDouble::from_bit_pair(Q[4]),
190    );
191    let e3 = DoubleDouble::mul_add(
192        recip,
193        DoubleDouble::from_bit_pair(Q[7]),
194        DoubleDouble::from_bit_pair(Q[6]),
195    );
196    let e4 = DoubleDouble::mul_add(
197        recip,
198        DoubleDouble::from_bit_pair(Q[9]),
199        DoubleDouble::from_bit_pair(Q[8]),
200    );
201    let e5 = DoubleDouble::mul_add(
202        recip,
203        DoubleDouble::from_bit_pair(Q[11]),
204        DoubleDouble::from_bit_pair(Q[10]),
205    );
206
207    let f0 = DoubleDouble::mul_add(x2, e1, e0);
208    let f1 = DoubleDouble::mul_add(x2, e3, e2);
209    let f2 = DoubleDouble::mul_add(x2, e5, e4);
210
211    let g0 = DoubleDouble::mul_add(x4, f1, f0);
212
213    let p_den = DoubleDouble::mul_add(x8, f2, g0);
214
215    let z = DoubleDouble::div(p_num, p_den);
216
217    let r_sqrt = DoubleDouble::from_rsqrt_fast(dx);
218
219    let r = z * r_sqrt;
220
221    let err = f_fmla(
222        r.hi,
223        f64::from_bits(0x3c40000000000000), // 2^-59
224        f64::from_bits(0x3ba0000000000000), // 2^-69
225    );
226    let up = r.hi + (r.lo + err);
227    let lb = r.hi + (r.lo - err);
228    if up == lb {
229        return f64::copysign(r.to_f64(), sign_scale);
230    }
231    i1e_asympt_hard(x, sign_scale)
232}
233
234/**
235Asymptotic expansion for I1.
236
237Computes:
238sqrt(x) * exp(-x) * I1(x) = Pn(1/x)/Qn(1/x)
239hence:
240I1(x)exp(-|x|) = Pn(1/x)/Qm(1/x)/sqrt(x)
241
242Generated by Wolfram Mathematica:
243```text
244<<FunctionApproximations`
245ClearAll["Global`*"]
246f[x_]:=Sqrt[x] Exp[-x] BesselI[1,x]
247g[z_]:=f[1/z]
248{err,approx}=MiniMaxApproximation[g[z],{z,{1/713.98,1/7.75},15,15},WorkingPrecision->120]
249poly=Numerator[approx][[1]];
250coeffs=CoefficientList[poly,z];
251TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
252poly=Denominator[approx][[1]];
253coeffs=CoefficientList[poly,z];
254TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
255```
256**/
257#[cold]
258#[inline(never)]
259fn i1e_asympt_hard(x: f64, sign_scale: f64) -> f64 {
260    static P: [DyadicFloat128; 16] = [
261        DyadicFloat128 {
262            sign: DyadicSign::Pos,
263            exponent: -129,
264            mantissa: 0xcc42299e_a1b28468_bea7da47_28f13acc_u128,
265        },
266        DyadicFloat128 {
267            sign: DyadicSign::Neg,
268            exponent: -124,
269            mantissa: 0xda979406_3df6e66f_cf31c3f5_f194b48c_u128,
270        },
271        DyadicFloat128 {
272            sign: DyadicSign::Neg,
273            exponent: -120,
274            mantissa: 0xd60b7b96_c958929b_cabe1d8c_3d874767_u128,
275        },
276        DyadicFloat128 {
277            sign: DyadicSign::Pos,
278            exponent: -113,
279            mantissa: 0xd27aad9a_8fb38d56_46ab4510_8479306e_u128,
280        },
281        DyadicFloat128 {
282            sign: DyadicSign::Neg,
283            exponent: -108,
284            mantissa: 0xe0167305_c451bd1f_d2f17b68_5c62e2ff_u128,
285        },
286        DyadicFloat128 {
287            sign: DyadicSign::Pos,
288            exponent: -103,
289            mantissa: 0x8f6d238f_c80d8e4a_08c130f6_24e1c925_u128,
290        },
291        DyadicFloat128 {
292            sign: DyadicSign::Neg,
293            exponent: -100,
294            mantissa: 0xfe32280f_2ea99024_d9924472_92d7ac8f_u128,
295        },
296        DyadicFloat128 {
297            sign: DyadicSign::Pos,
298            exponent: -96,
299            mantissa: 0xa48815ac_d265609f_da4ace94_811390b2_u128,
300        },
301        DyadicFloat128 {
302            sign: DyadicSign::Neg,
303            exponent: -93,
304            mantissa: 0x9ededfe5_833b4cc1_731efd5c_f8729c6c_u128,
305        },
306        DyadicFloat128 {
307            sign: DyadicSign::Pos,
308            exponent: -91,
309            mantissa: 0xe5b43203_2784ae6a_f7458556_0a8308ea_u128,
310        },
311        DyadicFloat128 {
312            sign: DyadicSign::Neg,
313            exponent: -89,
314            mantissa: 0xf5df279a_3fb4ef60_8d10adee_7dd2f47b_u128,
315        },
316        DyadicFloat128 {
317            sign: DyadicSign::Pos,
318            exponent: -87,
319            mantissa: 0xbdb59963_7a757ed1_87280e0e_7f93ca2b_u128,
320        },
321        DyadicFloat128 {
322            sign: DyadicSign::Neg,
323            exponent: -86,
324            mantissa: 0xc87fdea5_53177ca8_c91de5fb_3f8f78d3_u128,
325        },
326        DyadicFloat128 {
327            sign: DyadicSign::Pos,
328            exponent: -85,
329            mantissa: 0x846d16a7_4663ef5d_ad42d599_5bc726b8_u128,
330        },
331        DyadicFloat128 {
332            sign: DyadicSign::Neg,
333            exponent: -86,
334            mantissa: 0xb3ed2055_74262d95_389f33e4_2ac3774a_u128,
335        },
336        DyadicFloat128 {
337            sign: DyadicSign::Pos,
338            exponent: -88,
339            mantissa: 0xa511aa32_c18c34e4_3d029a90_a71b7a55_u128,
340        },
341    ];
342    static Q: [DyadicFloat128; 16] = [
343        DyadicFloat128 {
344            sign: DyadicSign::Pos,
345            exponent: -127,
346            mantissa: 0x80000000_00000000_00000000_00000000_u128,
347        },
348        DyadicFloat128 {
349            sign: DyadicSign::Neg,
350            exponent: -122,
351            mantissa: 0x877b771a_ad8f5fd3_5aacf5f9_f04ee9de_u128,
352        },
353        DyadicFloat128 {
354            sign: DyadicSign::Neg,
355            exponent: -118,
356            mantissa: 0x89475ecd_9c84361e_800c8a3a_c8af23bf_u128,
357        },
358        DyadicFloat128 {
359            sign: DyadicSign::Pos,
360            exponent: -111,
361            mantissa: 0x837d1196_cf2723f1_23b54da8_225efe05_u128,
362        },
363        DyadicFloat128 {
364            sign: DyadicSign::Neg,
365            exponent: -106,
366            mantissa: 0x8ae3aecb_15355751_a9ee12e5_a4dd9dde_u128,
367        },
368        DyadicFloat128 {
369            sign: DyadicSign::Pos,
370            exponent: -102,
371            mantissa: 0xb0886afa_bc13f996_ab45d252_75c8f587_u128,
372        },
373        DyadicFloat128 {
374            sign: DyadicSign::Neg,
375            exponent: -98,
376            mantissa: 0x9b37d7cd_b114b86b_7d14a389_26599aa1_u128,
377        },
378        DyadicFloat128 {
379            sign: DyadicSign::Pos,
380            exponent: -95,
381            mantissa: 0xc716bf54_09d5dd9f_bc16679b_93aaeca4_u128,
382        },
383        DyadicFloat128 {
384            sign: DyadicSign::Neg,
385            exponent: -92,
386            mantissa: 0xbe0cd82e_c8af8371_ab028ed9_c7902dd2_u128,
387        },
388        DyadicFloat128 {
389            sign: DyadicSign::Pos,
390            exponent: -89,
391            mantissa: 0x875f8d91_8ef5d434_a39d00f9_2aed3d2a_u128,
392        },
393        DyadicFloat128 {
394            sign: DyadicSign::Neg,
395            exponent: -87,
396            mantissa: 0x8e030781_5aa4ce7f_70156b82_8b216b7c_u128,
397        },
398        DyadicFloat128 {
399            sign: DyadicSign::Pos,
400            exponent: -86,
401            mantissa: 0xd4dd2687_92646fbd_5ea2d422_da64fc0b_u128,
402        },
403        DyadicFloat128 {
404            sign: DyadicSign::Neg,
405            exponent: -85,
406            mantissa: 0xd6d72ab3_64b4a827_0499af0f_13a51a80_u128,
407        },
408        DyadicFloat128 {
409            sign: DyadicSign::Pos,
410            exponent: -84,
411            mantissa: 0x828f4e8b_728747a9_2cebe54a_810e2681_u128,
412        },
413        DyadicFloat128 {
414            sign: DyadicSign::Neg,
415            exponent: -85,
416            mantissa: 0x91570096_36a3fcfb_6b936d44_68dda1be_u128,
417        },
418        DyadicFloat128 {
419            sign: DyadicSign::Pos,
420            exponent: -89,
421            mantissa: 0xf082ad00_86024ed4_dd31613b_ec41e3f8_u128,
422        },
423    ];
424
425    let recip = DyadicFloat128::accurate_reciprocal(x);
426
427    let mut p_num = P[15];
428    for i in (0..15).rev() {
429        p_num = recip * p_num + P[i];
430    }
431    let mut p_den = Q[15];
432    for i in (0..15).rev() {
433        p_den = recip * p_den + Q[i];
434    }
435    let z = p_num * p_den.reciprocal();
436    let r_sqrt = bessel_rsqrt_hard(x, recip);
437    (z * r_sqrt).fast_as_f64() * sign_scale
438}
439
440#[cfg(test)]
441mod tests {
442    use super::*;
443    #[test]
444    fn test_fi1e() {
445        assert_eq!(f_i1e(f64::EPSILON), 1.1102230246251563e-16);
446        assert_eq!(f_i1e(7.750000000757874), 0.13605110007443239);
447        assert_eq!(f_i1e(7.482812501363189), 0.13818116726273896);
448        assert_eq!(f_i1e(-7.750000000757874), -0.13605110007443239);
449        assert_eq!(f_i1e(-7.482812501363189), -0.13818116726273896);
450        assert!(f_i1e(f64::NAN).is_nan());
451        assert_eq!(f_i1e(f64::INFINITY), 0.);
452        assert_eq!(f_i1e(f64::NEG_INFINITY), 0.);
453        assert_eq!(f_i1e(0.01), 0.004950311047118276);
454        assert_eq!(f_i1e(-0.01), -0.004950311047118276);
455        assert_eq!(f_i1e(-9.01), -0.12716101566063667);
456        assert_eq!(f_i1e(9.01), 0.12716101566063667);
457        assert_eq!(f_i1e(763.), 0.014435579051182581);
458        assert_eq!(i1e_asympt_hard(9.01, 1.), 0.12716101566063667);
459    }
460}
pxfm/bessel/i1e.rs

pxfm/bessel/
i1e.rs
