pxfm/bessel/
i1.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::common::f_fmla;
32use crate::double_double::DoubleDouble;
33use crate::dyadic_float::{DyadicFloat128, DyadicSign};
34use crate::exponents::rational128_exp;
35use crate::polyeval::{f_estrin_polyeval5, f_polyeval6};
36
37/// Modified Bessel of the first kind of order 1
38///
39/// Max found ULP 0.5
40pub fn f_i1(x: f64) -> f64 {
41    let ux = x.to_bits().wrapping_shl(1);
42
43    if ux >= 0x7ffu64 << 53 || ux <= 0x7960000000000000u64 {
44        // |x| <= f64::EPSILON, |x| == inf, x == NaN
45        if ux <= 0x760af31dc4611874u64 {
46            // Power series of I1(x) ~ x/2 + O(x^3)
47            // |x| <= 2.2204460492503131e-24
48            return x * 0.5;
49        }
50        if ux <= 0x7960000000000000u64 {
51            // |x| <= f64::EPSILON
52            // Power series of I1(x) ~ x/2 + x^3/16 + O(x^4)
53            const A0: f64 = 1. / 2.;
54            const A1: f64 = 1. / 16.;
55            let r0 = f_fmla(x, x * A1, A0);
56            return r0 * x;
57        }
58        if x.is_infinite() {
59            return if x.is_sign_positive() {
60                f64::INFINITY
61            } else {
62                f64::NEG_INFINITY
63            };
64        }
65        return x + f64::NAN; // x == NaN
66    }
67
68    let xb = x.to_bits() & 0x7fff_ffff_ffff_ffff;
69
70    if xb >= 0x40864fe69ff9fec8u64 {
71        // |x| >= 713.9876098185423
72        return if x.is_sign_negative() {
73            f64::NEG_INFINITY
74        } else {
75            f64::INFINITY
76        };
77    }
78
79    static SIGN: [f64; 2] = [1., -1.];
80
81    let sign_scale = SIGN[x.is_sign_negative() as usize];
82
83    if xb < 0x401f000000000000u64 {
84        // |x| <= 7.75
85        return f64::copysign(i1_0_to_7p75(f64::from_bits(xb)).to_f64(), sign_scale);
86    }
87
88    i1_asympt(f64::from_bits(xb), sign_scale)
89}
90
91/**
92Computes
93I1(x) = x/2 * (1 + 1 * (x/2)^2 + (x/2)^4 * P((x/2)^2))
94
95Polynomial generated by Wolfram Mathematica:
96```text
97<<FunctionApproximations`
98ClearAll["Global`*"]
99f[x_]:=(BesselI[1,x]*2/x-1-1/2(x/2)^2)/(x/2)^4
100g[z_]:=f[2 Sqrt[z]]
101{err,approx}=MiniMaxApproximation[g[z],{z,{0.000000001,7.5},9,9},WorkingPrecision->60]
102poly=Numerator[approx][[1]];
103coeffs=CoefficientList[poly,z];
104TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
105poly=Denominator[approx][[1]];
106coeffs=CoefficientList[poly,z];
107TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
108```
109**/
110#[inline]
111pub(crate) fn i1_0_to_7p75(x: f64) -> DoubleDouble {
112    let half_x = x * 0.5; // this is exact
113    let eval_x = DoubleDouble::from_exact_mult(half_x, half_x);
114
115    const P: [(u64, u64); 5] = [
116        (0x3c55555555555555, 0x3fb5555555555555),
117        (0x3c1253e1df138479, 0x3f7304597c4fbd4c),
118        (0x3bcec398b7059ee9, 0x3f287b5b01f6b9c0),
119        (0xbb7354e7c92c4f77, 0x3ed21de117470d10),
120        (0xbb1d35ac0d7923cc, 0x3e717f3714dddc04),
121    ];
122
123    let ps_num = f_estrin_polyeval5(
124        eval_x.hi,
125        f64::from_bits(0x3e063684ca1944a4),
126        f64::from_bits(0x3d92ac4a0e48a9bb),
127        f64::from_bits(0x3d1425988b0b0aea),
128        f64::from_bits(0x3c899839e74ddefc),
129        f64::from_bits(0x3bed8747bcdd1e61),
130    );
131
132    let mut p_num = DoubleDouble::mul_f64_add(eval_x, ps_num, DoubleDouble::from_bit_pair(P[4]));
133    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[3]));
134    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[2]));
135    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[1]));
136    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[0]));
137
138    const Q: [(u64, u64); 4] = [
139        (0x0000000000000000, 0x3ff0000000000000),
140        (0xbc3e59afb81ac7ea, 0xbf9c4848e0661d70),
141        (0x3bd62fa3dbc1a19c, 0x3f38a9eafcd7e674),
142        (0x3b6f4688b9eab7d0, 0xbecbfdec51454533),
143    ];
144
145    let ps_den = f_polyeval6(
146        eval_x.hi,
147        f64::from_bits(0x3e56e7cde9266a32),
148        f64::from_bits(0xbddc316dff4a672f),
149        f64::from_bits(0x3d5a43aaee30ebb5),
150        f64::from_bits(0xbcd1fb023f4f1fa0),
151        f64::from_bits(0x3c4089ede324209f),
152        f64::from_bits(0xbb9f64f47ba69604),
153    );
154
155    let mut p_den = DoubleDouble::mul_f64_add(eval_x, ps_den, DoubleDouble::from_bit_pair(Q[3]));
156    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[2]));
157    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[1]));
158    p_den = DoubleDouble::mul_add_f64(eval_x, p_den, f64::from_bits(0x3ff0000000000000));
159
160    let p = DoubleDouble::div(p_num, p_den);
161
162    let eval_sqr = DoubleDouble::quick_mult(eval_x, eval_x);
163
164    let mut z = DoubleDouble::mul_f64_add_f64(eval_x, 0.5, 1.);
165    z = DoubleDouble::mul_add(p, eval_sqr, z);
166
167    let x_over_05 = DoubleDouble::from_exact_mult(x, 0.5);
168
169    let r = DoubleDouble::quick_mult(z, x_over_05);
170
171    let err = f_fmla(
172        r.hi,
173        f64::from_bits(0x3c40000000000000), // 2^-59
174        f64::from_bits(0x3be0000000000000), // 2^-65
175    );
176    let ub = r.hi + (r.lo + err);
177    let lb = r.hi + (r.lo - err);
178    if ub == lb {
179        return r;
180    }
181    i1_0_to_7p5_hard(x)
182}
183
184// Polynomial generated by Wolfram Mathematica:
185// I1(x) = x/2 * (1 + 1 * (x/2)^2 + (x/2)^4 * P((x/2)^2))
186// <<FunctionApproximations`
187// ClearAll["Global`*"]
188// f[x_]:=(BesselI[1,x]*2/x-1-1/2(x/2)^2)/(x/2)^4
189// g[z_]:=f[2 Sqrt[z]]
190// {err,approx}=MiniMaxApproximation[g[z],{z,{0.000000001,7.5},9,9},WorkingPrecision->60]
191// poly=Numerator[approx][[1]];
192// coeffs=CoefficientList[poly,z];
193// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
194// poly=Denominator[approx][[1]];
195// coeffs=CoefficientList[poly,z];
196// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
197#[cold]
198#[inline(never)]
199pub(crate) fn i1_0_to_7p5_hard(x: f64) -> DoubleDouble {
200    const ONE_OVER_4: f64 = 1. / 4.;
201    let eval_x = DoubleDouble::quick_mult_f64(DoubleDouble::from_exact_mult(x, x), ONE_OVER_4);
202
203    const P: [(u64, u64); 10] = [
204        (0x3c55555555555555, 0x3fb5555555555555),
205        (0x3c1253e1df138479, 0x3f7304597c4fbd4c),
206        (0x3bcec398b7059ee9, 0x3f287b5b01f6b9c0),
207        (0xbb7354e7c92c4f77, 0x3ed21de117470d10),
208        (0xbb1d35ac0d7923cc, 0x3e717f3714dddc04),
209        (0xba928dee24678e32, 0x3e063684ca1944a4),
210        (0xba36aa59912fcbed, 0x3d92ac4a0e48a9bb),
211        (0x39bad76f18b5ad37, 0x3d1425988b0b0aea),
212        (0xb923a6bab6928df4, 0x3c899839e74ddefc),
213        (0x3864356cdfa7b321, 0x3bed8747bcdd1e61),
214    ];
215
216    let mut p_num = DoubleDouble::mul_add(
217        eval_x,
218        DoubleDouble::from_bit_pair(P[9]),
219        DoubleDouble::from_bit_pair(P[8]),
220    );
221    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[7]));
222    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[6]));
223    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[5]));
224    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[4]));
225    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[3]));
226    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[2]));
227    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[1]));
228    p_num = DoubleDouble::mul_add(eval_x, p_num, DoubleDouble::from_bit_pair(P[0]));
229
230    const Q: [(u64, u64); 10] = [
231        (0x0000000000000000, 0x3ff0000000000000),
232        (0xbc3e59afb81ac7ea, 0xbf9c4848e0661d70),
233        (0x3bd62fa3dbc1a19c, 0x3f38a9eafcd7e674),
234        (0x3b6f4688b9eab7d0, 0xbecbfdec51454533),
235        (0x3af0fb4a17103ef8, 0x3e56e7cde9266a32),
236        (0xba71755779c6d4bd, 0xbddc316dff4a672f),
237        (0x39cf8ed8d449e2c6, 0x3d5a43aaee30ebb5),
238        (0x39704e900a373874, 0xbcd1fb023f4f1fa0),
239        (0xb8e33e87e4bffbb1, 0x3c4089ede324209f),
240        (0x380fb09b3fd49d5c, 0xbb9f64f47ba69604),
241    ];
242
243    let mut p_den = DoubleDouble::mul_add(
244        eval_x,
245        DoubleDouble::from_bit_pair(Q[9]),
246        DoubleDouble::from_bit_pair(Q[8]),
247    );
248    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[7]));
249    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[6]));
250    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[5]));
251    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[4]));
252    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[3]));
253    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[2]));
254    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[1]));
255    p_den = DoubleDouble::mul_add(eval_x, p_den, DoubleDouble::from_bit_pair(Q[0]));
256
257    let p = DoubleDouble::div(p_num, p_den);
258
259    let eval_sqr = DoubleDouble::quick_mult(eval_x, eval_x);
260
261    let mut z = DoubleDouble::mul_f64_add_f64(eval_x, 0.5, 1.);
262    z = DoubleDouble::mul_add(p, eval_sqr, z);
263
264    let x_over_05 = DoubleDouble::from_exact_mult(x, 0.5);
265
266    DoubleDouble::quick_mult(z, x_over_05)
267}
268
269/**
270Asymptotic expansion for I1.
271
272Computes:
273sqrt(x) * exp(-x) * I1(x) = Pn(1/x)/Qn(1/x)
274hence:
275I1(x) = Pn(1/x)/Qm(1/x)*exp(x)/sqrt(x)
276
277Generated by Wolfram Mathematica:
278```text
279<<FunctionApproximations`
280ClearAll["Global`*"]
281f[x_]:=Sqrt[x] Exp[-x] BesselI[1,x]
282g[z_]:=f[1/z]
283{err,approx}=MiniMaxApproximation[g[z],{z,{1/713.98,1/7.75},11,11},WorkingPrecision->120]
284poly=Numerator[approx][[1]];
285coeffs=CoefficientList[poly,z];
286TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
287poly=Denominator[approx][[1]];
288coeffs=CoefficientList[poly,z];
289TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
290```
291**/
292#[inline]
293fn i1_asympt(x: f64, sign_scale: f64) -> f64 {
294    let dx = x;
295    let recip = DoubleDouble::from_quick_recip(x);
296    const P: [(u64, u64); 12] = [
297        (0xbc73a823f28a2f5e, 0x3fd9884533d43651),
298        (0x3cc0d5bb78e674b3, 0xc0354325c8029263),
299        (0x3d20e1155aaaa283, 0x4080c09b027c46a4),
300        (0xbd5b90dcf81b99c1, 0xc0bfc1311090c839),
301        (0xbd98f2fda9e8fa1b, 0x40f3bb81bb190ae2),
302        (0xbdcec960752b60da, 0xc1207c0bbbc31cd9),
303        (0x3dd3c9a299c9c41f, 0x414253e25c4584af),
304        (0xbde82e7b9a3e1acc, 0xc159a656aece377c),
305        (0x3e0d3d30d701a8ab, 0x416398df24c74ef2),
306        (0xbdf57b85ab7006e2, 0xc151fd119be1702b),
307        (0x3dd760928f4515fd, 0xc1508327e42639bc),
308        (0x3dc09e71bc648589, 0x4143e4933afa621c),
309    ];
310
311    let x2 = DoubleDouble::quick_mult(recip, recip);
312    let x4 = DoubleDouble::quick_mult(x2, x2);
313    let x8 = DoubleDouble::quick_mult(x4, x4);
314
315    let e0 = DoubleDouble::mul_add(
316        recip,
317        DoubleDouble::from_bit_pair(P[1]),
318        DoubleDouble::from_bit_pair(P[0]),
319    );
320    let e1 = DoubleDouble::mul_add(
321        recip,
322        DoubleDouble::from_bit_pair(P[3]),
323        DoubleDouble::from_bit_pair(P[2]),
324    );
325    let e2 = DoubleDouble::mul_add(
326        recip,
327        DoubleDouble::from_bit_pair(P[5]),
328        DoubleDouble::from_bit_pair(P[4]),
329    );
330    let e3 = DoubleDouble::mul_add(
331        recip,
332        DoubleDouble::from_bit_pair(P[7]),
333        DoubleDouble::from_bit_pair(P[6]),
334    );
335    let e4 = DoubleDouble::mul_add(
336        recip,
337        DoubleDouble::from_bit_pair(P[9]),
338        DoubleDouble::from_bit_pair(P[8]),
339    );
340    let e5 = DoubleDouble::mul_add(
341        recip,
342        DoubleDouble::from_bit_pair(P[11]),
343        DoubleDouble::from_bit_pair(P[10]),
344    );
345
346    let f0 = DoubleDouble::mul_add(x2, e1, e0);
347    let f1 = DoubleDouble::mul_add(x2, e3, e2);
348    let f2 = DoubleDouble::mul_add(x2, e5, e4);
349
350    let g0 = DoubleDouble::mul_add(x4, f1, f0);
351
352    let p_num = DoubleDouble::mul_add(x8, f2, g0);
353
354    const Q: [(u64, u64); 12] = [
355        (0x0000000000000000, 0x3ff0000000000000),
356        (0xbcb334d5a476d9ad, 0xc04a75f94c1a0c1a),
357        (0xbd324d58ed98bfae, 0x4094b00e60301c42),
358        (0x3d7c8725666c4360, 0xc0d36b9d28d45928),
359        (0x3d7f8457c2945822, 0x4107d6c398a174ed),
360        (0x3dbc655ea216594b, 0xc1339393e6776e38),
361        (0xbdebb5dffef78272, 0x415537198d23f6a1),
362        (0xbdb577f8abad883e, 0xc16c6c399dcd6949),
363        (0x3e14261c5362f109, 0x4173c02446576949),
364        (0x3dc382ededad42c5, 0xc1547dff5770f4ec),
365        (0xbe05c0f74d4c7956, 0xc165c88046952562),
366        (0xbdbf9145927aa2c7, 0x414395e46bc45d5b),
367    ];
368
369    let e0 = DoubleDouble::mul_add_f64(
370        recip,
371        DoubleDouble::from_bit_pair(Q[1]),
372        f64::from_bits(0x3ff0000000000000),
373    );
374    let e1 = DoubleDouble::mul_add(
375        recip,
376        DoubleDouble::from_bit_pair(Q[3]),
377        DoubleDouble::from_bit_pair(Q[2]),
378    );
379    let e2 = DoubleDouble::mul_add(
380        recip,
381        DoubleDouble::from_bit_pair(Q[5]),
382        DoubleDouble::from_bit_pair(Q[4]),
383    );
384    let e3 = DoubleDouble::mul_add(
385        recip,
386        DoubleDouble::from_bit_pair(Q[7]),
387        DoubleDouble::from_bit_pair(Q[6]),
388    );
389    let e4 = DoubleDouble::mul_add(
390        recip,
391        DoubleDouble::from_bit_pair(Q[9]),
392        DoubleDouble::from_bit_pair(Q[8]),
393    );
394    let e5 = DoubleDouble::mul_add(
395        recip,
396        DoubleDouble::from_bit_pair(Q[11]),
397        DoubleDouble::from_bit_pair(Q[10]),
398    );
399
400    let f0 = DoubleDouble::mul_add(x2, e1, e0);
401    let f1 = DoubleDouble::mul_add(x2, e3, e2);
402    let f2 = DoubleDouble::mul_add(x2, e5, e4);
403
404    let g0 = DoubleDouble::mul_add(x4, f1, f0);
405
406    let p_den = DoubleDouble::mul_add(x8, f2, g0);
407
408    let z = DoubleDouble::div(p_num, p_den);
409
410    let e = i0_exp(dx * 0.5);
411    let r_sqrt = DoubleDouble::from_rsqrt_fast(dx);
412
413    let r = DoubleDouble::quick_mult(z * r_sqrt * e, e);
414
415    let err = f_fmla(
416        r.hi,
417        f64::from_bits(0x3c40000000000000), // 2^-59
418        f64::from_bits(0x3ba0000000000000), // 2^-69
419    );
420    let up = r.hi + (r.lo + err);
421    let lb = r.hi + (r.lo - err);
422    if up == lb {
423        return f64::copysign(r.to_f64(), sign_scale);
424    }
425    i1_asympt_hard(x, sign_scale)
426}
427
428/**
429Asymptotic expansion for I1.
430
431Computes:
432sqrt(x) * exp(-x) * I1(x) = Pn(1/x)/Qn(1/x)
433hence:
434I1(x) = Pn(1/x)/Qm(1/x)*exp(x)/sqrt(x)
435
436Generated by Wolfram Mathematica:
437```text
438<<FunctionApproximations`
439ClearAll["Global`*"]
440f[x_]:=Sqrt[x] Exp[-x] BesselI[1,x]
441g[z_]:=f[1/z]
442{err,approx}=MiniMaxApproximation[g[z],{z,{1/713.98,1/7.75},15,15},WorkingPrecision->120]
443poly=Numerator[approx][[1]];
444coeffs=CoefficientList[poly,z];
445TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
446poly=Denominator[approx][[1]];
447coeffs=CoefficientList[poly,z];
448TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
449```
450**/
451#[cold]
452#[inline(never)]
453fn i1_asympt_hard(x: f64, sign_scale: f64) -> f64 {
454    static P: [DyadicFloat128; 16] = [
455        DyadicFloat128 {
456            sign: DyadicSign::Pos,
457            exponent: -129,
458            mantissa: 0xcc42299e_a1b28468_bea7da47_28f13acc_u128,
459        },
460        DyadicFloat128 {
461            sign: DyadicSign::Neg,
462            exponent: -124,
463            mantissa: 0xda979406_3df6e66f_cf31c3f5_f194b48c_u128,
464        },
465        DyadicFloat128 {
466            sign: DyadicSign::Neg,
467            exponent: -120,
468            mantissa: 0xd60b7b96_c958929b_cabe1d8c_3d874767_u128,
469        },
470        DyadicFloat128 {
471            sign: DyadicSign::Pos,
472            exponent: -113,
473            mantissa: 0xd27aad9a_8fb38d56_46ab4510_8479306e_u128,
474        },
475        DyadicFloat128 {
476            sign: DyadicSign::Neg,
477            exponent: -108,
478            mantissa: 0xe0167305_c451bd1f_d2f17b68_5c62e2ff_u128,
479        },
480        DyadicFloat128 {
481            sign: DyadicSign::Pos,
482            exponent: -103,
483            mantissa: 0x8f6d238f_c80d8e4a_08c130f6_24e1c925_u128,
484        },
485        DyadicFloat128 {
486            sign: DyadicSign::Neg,
487            exponent: -100,
488            mantissa: 0xfe32280f_2ea99024_d9924472_92d7ac8f_u128,
489        },
490        DyadicFloat128 {
491            sign: DyadicSign::Pos,
492            exponent: -96,
493            mantissa: 0xa48815ac_d265609f_da4ace94_811390b2_u128,
494        },
495        DyadicFloat128 {
496            sign: DyadicSign::Neg,
497            exponent: -93,
498            mantissa: 0x9ededfe5_833b4cc1_731efd5c_f8729c6c_u128,
499        },
500        DyadicFloat128 {
501            sign: DyadicSign::Pos,
502            exponent: -91,
503            mantissa: 0xe5b43203_2784ae6a_f7458556_0a8308ea_u128,
504        },
505        DyadicFloat128 {
506            sign: DyadicSign::Neg,
507            exponent: -89,
508            mantissa: 0xf5df279a_3fb4ef60_8d10adee_7dd2f47b_u128,
509        },
510        DyadicFloat128 {
511            sign: DyadicSign::Pos,
512            exponent: -87,
513            mantissa: 0xbdb59963_7a757ed1_87280e0e_7f93ca2b_u128,
514        },
515        DyadicFloat128 {
516            sign: DyadicSign::Neg,
517            exponent: -86,
518            mantissa: 0xc87fdea5_53177ca8_c91de5fb_3f8f78d3_u128,
519        },
520        DyadicFloat128 {
521            sign: DyadicSign::Pos,
522            exponent: -85,
523            mantissa: 0x846d16a7_4663ef5d_ad42d599_5bc726b8_u128,
524        },
525        DyadicFloat128 {
526            sign: DyadicSign::Neg,
527            exponent: -86,
528            mantissa: 0xb3ed2055_74262d95_389f33e4_2ac3774a_u128,
529        },
530        DyadicFloat128 {
531            sign: DyadicSign::Pos,
532            exponent: -88,
533            mantissa: 0xa511aa32_c18c34e4_3d029a90_a71b7a55_u128,
534        },
535    ];
536    static Q: [DyadicFloat128; 16] = [
537        DyadicFloat128 {
538            sign: DyadicSign::Pos,
539            exponent: -127,
540            mantissa: 0x80000000_00000000_00000000_00000000_u128,
541        },
542        DyadicFloat128 {
543            sign: DyadicSign::Neg,
544            exponent: -122,
545            mantissa: 0x877b771a_ad8f5fd3_5aacf5f9_f04ee9de_u128,
546        },
547        DyadicFloat128 {
548            sign: DyadicSign::Neg,
549            exponent: -118,
550            mantissa: 0x89475ecd_9c84361e_800c8a3a_c8af23bf_u128,
551        },
552        DyadicFloat128 {
553            sign: DyadicSign::Pos,
554            exponent: -111,
555            mantissa: 0x837d1196_cf2723f1_23b54da8_225efe05_u128,
556        },
557        DyadicFloat128 {
558            sign: DyadicSign::Neg,
559            exponent: -106,
560            mantissa: 0x8ae3aecb_15355751_a9ee12e5_a4dd9dde_u128,
561        },
562        DyadicFloat128 {
563            sign: DyadicSign::Pos,
564            exponent: -102,
565            mantissa: 0xb0886afa_bc13f996_ab45d252_75c8f587_u128,
566        },
567        DyadicFloat128 {
568            sign: DyadicSign::Neg,
569            exponent: -98,
570            mantissa: 0x9b37d7cd_b114b86b_7d14a389_26599aa1_u128,
571        },
572        DyadicFloat128 {
573            sign: DyadicSign::Pos,
574            exponent: -95,
575            mantissa: 0xc716bf54_09d5dd9f_bc16679b_93aaeca4_u128,
576        },
577        DyadicFloat128 {
578            sign: DyadicSign::Neg,
579            exponent: -92,
580            mantissa: 0xbe0cd82e_c8af8371_ab028ed9_c7902dd2_u128,
581        },
582        DyadicFloat128 {
583            sign: DyadicSign::Pos,
584            exponent: -89,
585            mantissa: 0x875f8d91_8ef5d434_a39d00f9_2aed3d2a_u128,
586        },
587        DyadicFloat128 {
588            sign: DyadicSign::Neg,
589            exponent: -87,
590            mantissa: 0x8e030781_5aa4ce7f_70156b82_8b216b7c_u128,
591        },
592        DyadicFloat128 {
593            sign: DyadicSign::Pos,
594            exponent: -86,
595            mantissa: 0xd4dd2687_92646fbd_5ea2d422_da64fc0b_u128,
596        },
597        DyadicFloat128 {
598            sign: DyadicSign::Neg,
599            exponent: -85,
600            mantissa: 0xd6d72ab3_64b4a827_0499af0f_13a51a80_u128,
601        },
602        DyadicFloat128 {
603            sign: DyadicSign::Pos,
604            exponent: -84,
605            mantissa: 0x828f4e8b_728747a9_2cebe54a_810e2681_u128,
606        },
607        DyadicFloat128 {
608            sign: DyadicSign::Neg,
609            exponent: -85,
610            mantissa: 0x91570096_36a3fcfb_6b936d44_68dda1be_u128,
611        },
612        DyadicFloat128 {
613            sign: DyadicSign::Pos,
614            exponent: -89,
615            mantissa: 0xf082ad00_86024ed4_dd31613b_ec41e3f8_u128,
616        },
617    ];
618
619    let recip = DyadicFloat128::accurate_reciprocal(x);
620
621    let mut p_num = P[15];
622    for i in (0..15).rev() {
623        p_num = recip * p_num + P[i];
624    }
625    let mut p_den = Q[15];
626    for i in (0..15).rev() {
627        p_den = recip * p_den + Q[i];
628    }
629    let z = p_num * p_den.reciprocal();
630    let r_sqrt = bessel_rsqrt_hard(x, recip);
631    let f_exp = rational128_exp(x);
632    (z * r_sqrt * f_exp).fast_as_f64() * sign_scale
633}
634
635#[cfg(test)]
636mod tests {
637    use super::*;
638    #[test]
639    fn test_fi1() {
640        assert_eq!(
641            f_i1(0.0000000000000000000000000000000006423424234121),
642            3.2117121170605e-34
643        );
644        assert_eq!(f_i1(7.750000000757874), 315.8524811496668);
645        assert_eq!(f_i1(7.482812501363189), 245.58002285881892);
646        assert_eq!(f_i1(-7.750000000757874), -315.8524811496668);
647        assert_eq!(f_i1(-7.482812501363189), -245.58002285881892);
648        assert!(f_i1(f64::NAN).is_nan());
649        assert_eq!(f_i1(f64::INFINITY), f64::INFINITY);
650        assert_eq!(f_i1(f64::NEG_INFINITY), f64::NEG_INFINITY);
651        assert_eq!(f_i1(0.01), 0.005000062500260418);
652        assert_eq!(f_i1(-0.01), -0.005000062500260418);
653        assert_eq!(f_i1(-9.01), -1040.752038018038);
654        assert_eq!(f_i1(9.01), 1040.752038018038);
655    }
656}
pxfm/bessel/i1.rs

pxfm/bessel/
i1.rs
