rv/misc/
bessel.rs

1const MAX_ITER: usize = 500;
2
3const BESSI0_COEFFS_A: [f64; 30] = [
4    -4.415_341_646_479_339_5E-18,
5    3.330_794_518_822_238_4E-17,
6    -2.431_279_846_547_955E-16,
7    1.715_391_285_555_133E-15,
8    -1.168_533_287_799_345_1E-14,
9    7.676_185_498_604_936E-14,
10    -4.856_446_783_111_929E-13,
11    2.955_052_663_129_64E-12,
12    -1.726_826_291_441_556E-11,
13    9.675_809_035_373_237E-11,
14    -5.189_795_601_635_263E-10,
15    2.659_823_724_682_386_6E-9,
16    -1.300_025_009_986_248E-8,
17    6.046_995_022_541_919E-8,
18    -2.670_793_853_940_612E-7,
19    1.117_387_539_120_103_7E-6,
20    -4.416_738_358_458_750_5E-6,
21    1.644_844_807_072_889_6E-5,
22    -5.754_195_010_082_104E-5,
23    1.885_028_850_958_416_5E-4,
24    -5.763_755_745_385_824E-4,
25    1.639_475_616_941_335_7E-3,
26    -4.324_309_995_050_576E-3,
27    1.054_646_039_459_499_8E-2,
28    -2.373_741_480_589_947E-2,
29    4.930_528_423_967_071E-2,
30    -9.490_109_704_804_764E-2,
31    1.716_209_015_222_087_7E-1,
32    -3.046_826_723_431_984E-1,
33    6.767_952_744_094_761E-1,
34];
35
36const BESSI0_COEFFS_B: [f64; 25] = [
37    -7.233_180_487_874_754E-18,
38    -4.830_504_485_944_182E-18,
39    4.465_621_420_296_76E-17,
40    3.461_222_867_697_461E-17,
41    -2.827_623_980_516_583_6E-16,
42    -3.425_485_619_677_219E-16,
43    1.772_560_133_056_526_3E-15,
44    3.811_680_669_352_622_4E-15,
45    -9.554_846_698_828_307E-15,
46    -4.150_569_347_287_222E-14,
47    1.540_086_217_521_41E-14,
48    3.852_778_382_742_142_6E-13,
49    7.180_124_451_383_666E-13,
50    -1.794_178_531_506_806_2E-12,
51    -1.321_581_184_044_771_3E-11,
52    -3.149_916_527_963_241_6E-11,
53    1.188_914_710_784_643_9E-11,
54    4.940_602_388_224_97E-10,
55    3.396_232_025_708_386_5E-9,
56    2.266_668_990_498_178E-8,
57    2.048_918_589_469_063_8E-7,
58    2.891_370_520_834_756_7E-6,
59    6.889_758_346_916_825E-5,
60    3.369_116_478_255_694_3E-3,
61    8.044_904_110_141_088E-1,
62];
63
64const BESSI1_COEFFS_A: [f64; 29] = [
65    2.777_914_112_761_046_4E-18,
66    -2.111_421_214_358_166E-17,
67    1.553_631_957_736_200_5E-16,
68    -1.105_596_947_735_386_2E-15,
69    7.600_684_294_735_408E-15,
70    -5.042_185_504_727_912E-14,
71    3.223_793_365_945_575E-13,
72    -1.983_974_397_764_943_6E-12,
73    1.173_618_629_889_090_1E-11,
74    -6.663_489_723_502_027E-11,
75    3.625_590_281_552_117E-10,
76    -1.887_249_751_722_829_4E-9,
77    9.381_537_386_495_773E-9,
78    -4.445_059_128_796_328E-8,
79    2.003_294_753_552_135_3E-7,
80    -8.568_720_264_695_455E-7,
81    3.470_251_308_137_678_5E-6,
82    -1.327_316_365_603_943_6E-5,
83    4.781_565_107_550_054E-5,
84    -1.617_608_158_258_967_4E-4,
85    5.122_859_561_685_758E-4,
86    -1.513_572_450_631_253_2E-3,
87    4.156_422_944_312_888E-3,
88    -1.056_408_489_462_619_7E-2,
89    2.472_644_903_062_651_6E-2,
90    -5.294_598_120_809_499E-2,
91    1.026_436_586_898_471E-1,
92    -1.764_165_183_578_340_6E-1,
93    2.525_871_864_436_336_5E-1,
94];
95
96#[allow(clippy::unreadable_literal)]
97#[allow(clippy::excessive_precision)]
98const BESSI1_COEFFS_B: [f64; 25] = [
99    7.51729631084210481353E-18,
100    4.41434832307170791151E-18,
101    -4.65030536848935832153E-17,
102    -3.20952592199342395980E-17,
103    2.96262899764595013876E-16,
104    3.30820231092092828324E-16,
105    -1.88035477551078244854E-15,
106    -3.81440307243700780478E-15,
107    1.04202769841288027642E-14,
108    4.27244001671195135429E-14,
109    -2.10154184277266431302E-14,
110    -4.08355111109219731823E-13,
111    -7.19855177624590851209E-13,
112    2.03562854414708950722E-12,
113    1.41258074366137813316E-11,
114    3.25260358301548823856E-11,
115    -1.89749581235054123450E-11,
116    -5.58974346219658380687E-10,
117    -3.83538038596423702205E-9,
118    -2.63146884688951950684E-8,
119    -2.51223623787020892529E-7,
120    -3.88256480887769039346E-6,
121    -1.10588938762623716291E-4,
122    -9.76109749136146840777E-3,
123    7.78576235018280120474E-1,
124];
125
126fn chbevl(x: f64, coeffs: &[f64]) -> f64 {
127    let mut b0 = coeffs[0];
128    let mut b1 = 0.0;
129    let mut b2 = 0.0;
130
131    coeffs.iter().skip(1).for_each(|c| {
132        b2 = b1;
133        b1 = b0;
134        b0 = x.mul_add(b1, *c) - b2;
135    });
136
137    0.5 * (b0 - b2)
138}
139
140/// Modified Bessel function, I<sub>0</sub>(x)
141#[must_use]
142pub fn i0(x: f64) -> f64 {
143    let ax = x.abs();
144
145    if ax <= 8.0 {
146        let y = ax.mul_add(0.5, -2.0);
147        ax.exp() * chbevl(y, &BESSI0_COEFFS_A)
148    } else {
149        ax.exp() * chbevl(32.0_f64.mul_add(ax.recip(), -2.0), &BESSI0_COEFFS_B)
150            / ax.sqrt()
151    }
152}
153
154/// Logarithm of Modified Bessel function, `log I<sub>0</sub>(x)`
155#[must_use]
156pub fn log_i0(x: f64) -> f64 {
157    let ax = x.abs();
158
159    if ax <= 8.0 {
160        let y = ax.mul_add(0.5, -2.0);
161        ax + chbevl(y, &BESSI0_COEFFS_A).ln()
162    } else {
163        0.5_f64.mul_add(
164            -ax.ln(),
165            ax + chbevl(32.0_f64.mul_add(ax.recip(), -2.0), &BESSI0_COEFFS_B)
166                .ln(),
167        )
168    }
169}
170
171/// Modified Bessel function, I<sub>1</sub>(x)
172#[must_use]
173pub fn i1(x: f64) -> f64 {
174    let z = x.abs();
175    let res = if z <= 8.0 {
176        let y = z.mul_add(0.5, -2.0);
177        chbevl(y, &BESSI1_COEFFS_A) * z * z.exp()
178    } else {
179        z.exp() * chbevl(32.0_f64.mul_add(x.recip(), -2.0), &BESSI1_COEFFS_B)
180            / z.sqrt()
181    };
182
183    res * x.signum()
184}
185
186/// An encounterable error when computing Bessel's I function
187#[derive(Clone, Debug, PartialEq, Eq)]
188pub enum BesselIvError {
189    /// The order, v, must be an integer if z is negative.
190    OrderNotIntegerForNegativeZ,
191    /// Arguments would lead to an overflow
192    Overflow,
193    /// Failed to converge
194    FailedToConverge,
195    /// Precision Error
196    PrecisionLoss,
197    /// Input parameters are outside the domain
198    Domain,
199}
200
201impl std::error::Error for BesselIvError {}
202
203#[cfg_attr(coverage_nightly, coverage(off))]
204impl std::fmt::Display for BesselIvError {
205    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
206        match self {
207            Self::OrderNotIntegerForNegativeZ => {
208                write!(f, "If z is negative, the order, v, must be an integer.")
209            }
210            Self::Overflow => {
211                write!(f, "Arguments would cause an overflow.")
212            }
213            Self::FailedToConverge => {
214                write!(f, "Failed to converge.")
215            }
216            Self::PrecisionLoss => {
217                write!(f, "Precision loss.")
218            }
219            Self::Domain => {
220                write!(f, "Parameters are outside domain.")
221            }
222        }
223    }
224}
225
226/// Modified Bessel function of the first kind of real order
227pub fn bessel_iv(v: f64, z: f64) -> Result<f64, BesselIvError> {
228    if v.is_nan() || z.is_nan() {
229        return Ok(f64::NAN);
230    }
231    let (v, t) = {
232        let t = v.floor();
233        if v < 0.0 && (t - v).abs() < f64::EPSILON {
234            (-v, -t)
235        } else {
236            (v, t)
237        }
238    };
239
240    let sign: f64 = if z < 0.0 {
241        // Return error if v is not an integer if x < 0
242        if (t - v).abs() > f64::EPSILON {
243            return Err(BesselIvError::OrderNotIntegerForNegativeZ);
244        }
245
246        if 2.0_f64.mul_add(-(v / 2.0).floor(), v) > f64::EPSILON {
247            -1.0
248        } else {
249            1.0
250        }
251    } else {
252        1.0
253    };
254
255    if z == 0.0 {
256        if v == 0.0 {
257            return Ok(1.0);
258        } else if v < 0.0 {
259            return Err(BesselIvError::Overflow);
260        }
261        return Ok(0.0);
262    }
263
264    let az = z.abs();
265    let res: f64 = if v.abs() > 50.0 {
266        // Use asymptotic expansion for large orders
267        bessel_ikv_asymptotic_uniform(v, az)?.0
268    } else {
269        // Use Temme's method for small orders
270        bessel_ikv_temme(v, az)?.0
271    };
272
273    Ok(res * sign)
274}
275
276const N_UFACTORS: usize = 11;
277const N_UFACTOR_TERMS: usize = 31;
278const ASYMPTOTIC_UFACTORS: [[f64; N_UFACTOR_TERMS]; N_UFACTORS] = [
279    [
280        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
281        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
282        0.0, 0.0, 1.0,
283    ],
284    [
285        0.0,
286        0.0,
287        0.0,
288        0.0,
289        0.0,
290        0.0,
291        0.0,
292        0.0,
293        0.0,
294        0.0,
295        0.0,
296        0.0,
297        0.0,
298        0.0,
299        0.0,
300        0.0,
301        0.0,
302        0.0,
303        0.0,
304        0.0,
305        0.0,
306        0.0,
307        0.0,
308        0.0,
309        0.0,
310        0.0,
311        0.0,
312        -0.208_333_333_333_333_34,
313        0.0,
314        0.125,
315        0.0,
316    ],
317    [
318        0.0,
319        0.0,
320        0.0,
321        0.0,
322        0.0,
323        0.0,
324        0.0,
325        0.0,
326        0.0,
327        0.0,
328        0.0,
329        0.0,
330        0.0,
331        0.0,
332        0.0,
333        0.0,
334        0.0,
335        0.0,
336        0.0,
337        0.0,
338        0.0,
339        0.0,
340        0.0,
341        0.0,
342        0.334_201_388_888_888_9,
343        0.0,
344        -0.401_041_666_666_666_7,
345        0.0,
346        0.070_312_5,
347        0.0,
348        0.0,
349    ],
350    [
351        0.0,
352        0.0,
353        0.0,
354        0.0,
355        0.0,
356        0.0,
357        0.0,
358        0.0,
359        0.0,
360        0.0,
361        0.0,
362        0.0,
363        0.0,
364        0.0,
365        0.0,
366        0.0,
367        0.0,
368        0.0,
369        0.0,
370        0.0,
371        0.0,
372        -1.025_812_596_450_617_3,
373        0.0,
374        1.846_462_673_611_111_2,
375        0.0,
376        -0.891_210_937_5,
377        0.0,
378        0.073_242_187_5,
379        0.0,
380        0.0,
381        0.0,
382    ],
383    [
384        0.0,
385        0.0,
386        0.0,
387        0.0,
388        0.0,
389        0.0,
390        0.0,
391        0.0,
392        0.0,
393        0.0,
394        0.0,
395        0.0,
396        0.0,
397        0.0,
398        0.0,
399        0.0,
400        0.0,
401        0.0,
402        4.669_584_423_426_247,
403        0.0,
404        -11.207_002_616_222_995,
405        0.0,
406        8.789_123_535_156_25,
407        0.0,
408        -2.364_086_914_062_5,
409        0.0,
410        0.112_152_099_609_375,
411        0.0,
412        0.0,
413        0.0,
414        0.0,
415    ],
416    [
417        0.0,
418        0.0,
419        0.0,
420        0.0,
421        0.0,
422        0.0,
423        0.0,
424        0.0,
425        0.0,
426        0.0,
427        0.0,
428        0.0,
429        0.0,
430        0.0,
431        0.0,
432        -28.212_072_558_200_244,
433        0.0,
434        84.636_217_674_600_74,
435        0.0,
436        -91.818_241_543_240_03,
437        0.0,
438        42.534_998_745_388_46,
439        0.0,
440        -7.368_794_359_479_631,
441        0.0,
442        0.227_108_001_708_984_38,
443        0.0,
444        0.0,
445        0.0,
446        0.0,
447        0.0,
448    ],
449    [
450        0.0,
451        0.0,
452        0.0,
453        0.0,
454        0.0,
455        0.0,
456        0.0,
457        0.0,
458        0.0,
459        0.0,
460        0.0,
461        0.0,
462        212.570_130_039_217_1,
463        0.0,
464        -765.252_468_141_181_6,
465        0.0,
466        1_059.990_452_527_999_9,
467        0.0,
468        -699.579_627_376_132_7,
469        0.0,
470        218.190_511_744_211_6,
471        0.0,
472        -26.491_430_486_951_554,
473        0.0,
474        0.572_501_420_974_731_4,
475        0.0,
476        0.0,
477        0.0,
478        0.0,
479        0.0,
480        0.0,
481    ],
482    [
483        0.0,
484        0.0,
485        0.0,
486        0.0,
487        0.0,
488        0.0,
489        0.0,
490        0.0,
491        0.0,
492        -1_919.457_662_318_406_8,
493        0.0,
494        8_061.722_181_737_308,
495        0.0,
496        -13_586.550_006_434_136,
497        0.0,
498        11_655.393_336_864_536,
499        0.0,
500        -5_305.646_978_613_405,
501        0.0,
502        1_200.902_913_216_352_5,
503        0.0,
504        -108.090_919_788_394_64,
505        0.0,
506        1.727_727_502_584_457_4,
507        0.0,
508        0.0,
509        0.0,
510        0.0,
511        0.0,
512        0.0,
513        0.0,
514    ],
515    [
516        0.0,
517        0.0,
518        0.0,
519        0.0,
520        0.0,
521        0.0,
522        20_204.291_330_966_15,
523        0.0,
524        -96_980.598_388_637_5,
525        0.0,
526        192_547.001_232_531_5,
527        0.0,
528        -203_400.177_280_415_55,
529        0.0,
530        122_200.464_983_017_47,
531        0.0,
532        -41_192.654_968_897_56,
533        0.0,
534        7_109.514_302_489_364,
535        0.0,
536        -493.915_304_773_088,
537        0.0,
538        6.074_042_001_273_483,
539        0.0,
540        0.0,
541        0.0,
542        0.0,
543        0.0,
544        0.0,
545        0.0,
546        0.0,
547    ],
548    [
549        0.0,
550        0.0,
551        0.0,
552        -242_919.187_900_551_33,
553        0.0,
554        1_311_763.614_662_977,
555        0.0,
556        -2_998_015.918_538_106,
557        0.0,
558        3_763_271.297_656_404,
559        0.0,
560        -2_813_563.226_586_534,
561        0.0,
562        1_268_365.273_321_624_8,
563        0.0,
564        -331_645.172_484_563_6,
565        0.0,
566        45_218.768_981_362_74,
567        0.0,
568        -2_499.830_481_811_209,
569        0.0,
570        24.380_529_699_556_064,
571        0.0,
572        0.0,
573        0.0,
574        0.0,
575        0.0,
576        0.0,
577        0.0,
578        0.0,
579        0.0,
580    ],
581    [
582        3_284_469.853_072_037_5,
583        0.0,
584        -19_706_819.118_432_22,
585        0.0,
586        50_952_602.492_664_63,
587        0.0,
588        -74_105_148.211_532_64,
589        0.0,
590        66_344_512.274_729_03,
591        0.0,
592        -37_567_176.660_763_35,
593        0.0,
594        13_288_767.166_421_82,
595        0.0,
596        -2_785_618.128_086_455,
597        0.0,
598        308_186.404_612_662_45,
599        0.0,
600        -13_886.089_753_717_039,
601        0.0,
602        110.017_140_269_246_74,
603        0.0,
604        0.0,
605        0.0,
606        0.0,
607        0.0,
608        0.0,
609        0.0,
610        0.0,
611        0.0,
612        0.0,
613    ],
614];
615
616///Compute Iv, Kv from (AMS5 9.7.7 + 9.7.8), asymptotic expansion for large v
617///
618/// Heavily inspired by
619/// <https://github.com/scipy/scipy/blob/1984f97749a355a6767cea55cad5d1dc6977ad5f/scipy/special/cephes/scipy_iv.c#L249>
620fn bessel_ikv_asymptotic_uniform(
621    v: f64,
622    x: f64,
623) -> Result<(f64, f64), BesselIvError> {
624    use std::f64::consts::PI;
625    let (v, sign) = (v.abs(), v.signum());
626
627    let z = x / v;
628    let t = z.mul_add(z, 1.0).sqrt().recip();
629    let t2 = t * t;
630    let eta = z.mul_add(z, 1.0).sqrt() + (z / (1.0 + t.recip())).ln();
631
632    let i_prefactor = (t / (2.0 * PI * v)).sqrt() * (v * eta).exp();
633    let mut i_sum = 1.0;
634
635    let k_prefactor = (PI * t / (2.0 * v)).sqrt() * (-v * eta).exp();
636    let mut k_sum = 1.0;
637
638    let mut divisor = v;
639    let mut term = 0.0;
640
641    for (n, item) in ASYMPTOTIC_UFACTORS
642        .iter()
643        .enumerate()
644        .take(N_UFACTORS)
645        .skip(1)
646    {
647        term = 0.0;
648        for k in
649            ((N_UFACTOR_TERMS - 1 - 3 * n)..(N_UFACTOR_TERMS - n)).step_by(2)
650        {
651            term *= t2;
652            term += item[k];
653        }
654        for _k in (1..n).step_by(2) {
655            term *= t2;
656        }
657        if n % 2 == 1 {
658            term *= t;
659        }
660
661        term /= divisor;
662        i_sum += term;
663        k_sum += if n % 2 == 0 { term } else { -term };
664
665        if term.abs() < f64::EPSILON {
666            break;
667        }
668        divisor *= v;
669    }
670
671    // check convergence
672    if term.abs() > 1E-3 * i_sum.abs() {
673        Err(BesselIvError::FailedToConverge)
674    } else if term.abs() > f64::EPSILON * i_sum.abs() {
675        Err(BesselIvError::PrecisionLoss)
676    } else {
677        let k_value = k_prefactor * k_sum;
678        let i_value = if sign > 0.0 {
679            i_prefactor * i_sum
680        } else {
681            i_prefactor.mul_add(
682                i_sum,
683                (2.0 / PI) * (PI * v).sin() * k_prefactor * k_sum,
684            )
685        };
686        Ok((i_value, k_value))
687    }
688}
689
690/// Compute I(v, x) and K(v, x) simultaneously by Temme's method, see
691/// Temme, Journal of Computational Physics, vol 19, 324 (1975)
692/// Heavily inspired by
693/// <https://github.com/scipy/scipy/blob/1984f97749a355a6767cea55cad5d1dc6977ad5f/scipy/special/cephes/scipy_iv.c#L532>
694#[allow(clippy::many_single_char_names)]
695pub(crate) fn bessel_ikv_temme(
696    v: f64,
697    x: f64,
698) -> Result<(f64, f64), BesselIvError> {
699    use std::f64::consts::PI;
700    let (v, reflect) = if v < 0.0 { (-v, true) } else { (v, false) };
701
702    let n = v.round();
703    let u = v - n;
704    let n = n as isize;
705
706    if x < 0.0 {
707        return Err(BesselIvError::Domain);
708    } else if x == 0.0 {
709        return Err(BesselIvError::Overflow);
710    }
711
712    let w = x.recip();
713    let (ku, ku_1) = if x <= 2.0 {
714        temme_ik_series(u, x)?
715    } else {
716        cf2_ik(u, x)?
717    };
718
719    let mut prev = ku;
720    let mut current = ku_1;
721    for k in 1..=n {
722        let kf = k as f64;
723        let next = 2.0 * (u + kf) * current / x + prev;
724        prev = current;
725        current = next;
726    }
727
728    let kv = prev;
729    let kv1 = current;
730
731    let lim = (4.0_f64.mul_add(v * v, 10.0) / (8.0 * x)).powi(3) / 24.0;
732
733    let iv = if lim < 10.0 * f64::EPSILON && x > 100.0 {
734        bessel_iv_asymptotic(v, x)?
735    } else {
736        let fv = cf1_ik(v, x)?;
737        w / kv.mul_add(fv, kv1)
738    };
739
740    if reflect {
741        let z = u + ((n % 2) as f64);
742        Ok(((2.0 / PI).mul_add((PI * z).sin() * kv, iv), kv))
743    } else {
744        Ok((iv, kv))
745    }
746}
747
748/// Modified Bessel functions of the first and second kind of fractional order
749///
750/// Calculate K(v, x) and K(v+1, x) by method analogous to
751/// Temme, Journal of Computational Physics, vol 21, 343 (1976)
752#[allow(clippy::many_single_char_names)]
753fn temme_ik_series(v: f64, x: f64) -> Result<(f64, f64), BesselIvError> {
754    use crate::consts::EULER_MASCERONI;
755    use crate::misc::gammafn;
756    use std::f64::consts::PI;
757    /*
758     * |x| <= 2, Temme series converge rapidly
759     * |x| > 2, the larger the |x|, the slower the convergence
760     */
761    debug_assert!(x.abs() <= 2.0);
762    debug_assert!(v.abs() <= 0.5);
763
764    let gp = gammafn(v + 1.0) - 1.0;
765    let gm = gammafn(1.0 - v) - 1.0;
766
767    let a = (x / 2.0).ln();
768    let b = (v * a).exp();
769    let sigma = -a * v;
770    let c = if v.abs() < 2.0 * f64::EPSILON {
771        1.0
772    } else {
773        (PI * v).sin() / (PI * v)
774    };
775    let d = if sigma.abs() < f64::EPSILON {
776        1.0
777    } else {
778        sigma.sinh() / sigma
779    };
780    let gamma1 = if v.abs() < f64::EPSILON {
781        -EULER_MASCERONI
782    } else {
783        (0.5 / v) * (gp - gm) * c
784    };
785    let gamma2 = (2.0 + gp + gm) * c / 2.0;
786
787    let mut p = (gp + 1.0) / (2.0 * b);
788    let mut q = (gm + 1.0) * b / 2.0;
789    let mut f = d.mul_add(-a * gamma2, sigma.cosh() * gamma1) / c;
790    let mut h = p;
791    let mut coef = 1.0;
792    let mut sum = coef * f;
793    let mut sum1 = coef * h;
794
795    for k in 1..MAX_ITER {
796        let kf = k as f64;
797        f = kf.mul_add(f, p + q) / kf.mul_add(kf, -v * v);
798        p /= kf - v;
799        q /= kf + v;
800        h = kf.mul_add(-f, p);
801        coef *= x * x / (4.0 * kf);
802        sum += coef * f;
803        sum1 += coef * h;
804
805        if (coef * f).abs() < sum.abs() * f64::EPSILON {
806            return Ok((sum, 2.0 * sum1 / x));
807        }
808    }
809
810    Err(BesselIvError::FailedToConverge)
811}
812
813/// Calculate K(v, x) and K(v+1, x) by evaluating continued fraction
814/// z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x), see
815/// Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987)
816#[allow(clippy::many_single_char_names)]
817fn cf2_ik(v: f64, x: f64) -> Result<(f64, f64), BesselIvError> {
818    use std::f64::consts::PI;
819    /*
820     * Steed's algorithm, see Thompson and Barnett,
821     * Journal of Computational Physics, vol 64, 490 (1986)
822     */
823    debug_assert!(x.abs() > 1.0);
824
825    let mut a = v.mul_add(v, -0.25);
826    let mut b = 2.0 * (x + 1.0);
827    let mut d = b.recip();
828
829    let mut delta = d;
830    let mut f = d;
831    let mut prev = 0.0;
832    let mut cur = 1.0;
833    let mut q = -a;
834    let mut c = -a;
835    let mut s = q.mul_add(delta, 1.0);
836
837    for k in 2..MAX_ITER {
838        let kf = k as f64;
839        a -= 2.0 * (kf - 1.0);
840        b += 2.0;
841        d = a.mul_add(d, b).recip();
842        delta *= b.mul_add(d, -1.0);
843        f += delta;
844
845        let t = (b - 2.0).mul_add(-cur, prev) / a;
846        prev = cur;
847        cur = t;
848        c *= -a / kf;
849        q += c * t;
850        s += q * delta;
851
852        if (q * delta).abs() < s.abs() * f64::EPSILON / 2.0 {
853            let kv = (PI / (2.0 * x)).sqrt() * (-x).exp() / s;
854            let kv1 = kv * v.mul_add(v, -0.25).mul_add(f, 0.5 + v + x) / x;
855            return Ok((kv, kv1));
856        }
857    }
858    Err(BesselIvError::FailedToConverge)
859}
860
861/// Evaluate continued fraction fv = I_(v+1) / `I_v`, derived from
862/// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73 */
863#[allow(clippy::many_single_char_names)]
864fn cf1_ik(v: f64, x: f64) -> Result<f64, BesselIvError> {
865    /*
866     * |x| <= |v|, CF1_ik converges rapidly
867     * |x| > |v|, CF1_ik needs O(|x|) iterations to converge
868     */
869
870    /*
871     * modified Lentz's method, see
872     * Lentz, Applied Optics, vol 15, 668 (1976)
873     */
874
875    const TOL: f64 = f64::EPSILON;
876    let tiny: f64 = f64::MAX.sqrt().recip();
877    let mut c = tiny;
878    let mut f = tiny;
879    let mut d = 0.0;
880
881    for k in 1..MAX_ITER {
882        let kf = k as f64;
883        let a = 1.0;
884        let b = 2.0 * (v + kf) / x;
885        c = b + a / c;
886        d = a.mul_add(d, b);
887        if c == 0.0 {
888            c = tiny;
889        }
890        if d == 0.0 {
891            d = tiny;
892        }
893        d = d.recip();
894        let delta = c * d;
895        f *= delta;
896
897        if (delta - 1.0).abs() <= TOL {
898            return Ok(f);
899        }
900    }
901
902    Err(BesselIvError::FailedToConverge)
903}
904
905/// Compute `I_v` from (AMS5 9.7.1), asymptotic expansion for large |z|
906///  `I_v` ~ exp(x)/sqrt(2 pi x) ( 1 + (4*v*v-1)/8x + (4*v*v-1)(4*v*v-9)/8x/2! + ...)
907///  Heavily inspired by
908///  <https://github.com/scipy/scipy/blob/1984f97749a355a6767cea55cad5d1dc6977ad5f/scipy/special/cephes/scipy_iv.c#L145>
909fn bessel_iv_asymptotic(v: f64, x: f64) -> Result<f64, BesselIvError> {
910    let prefactor = x.exp() / (2.0 * std::f64::consts::PI * x).sqrt();
911
912    if prefactor.is_infinite() {
913        Ok(x)
914    } else {
915        let mu = 4.0 * v * v;
916        let mut sum: f64 = 1.0;
917        let mut term: f64 = 1.0;
918        let mut k: usize = 1;
919
920        while term.abs() > f64::EPSILON * sum.abs() {
921            let kf = k as f64;
922            let factor = 2.0_f64
923                .mul_add(kf, -1.0)
924                .mul_add(-(2.0_f64.mul_add(kf, -1.0)), mu)
925                / (8.0 * x)
926                / kf;
927            if k > 100 {
928                return Err(BesselIvError::FailedToConverge);
929            }
930            term *= -factor;
931            sum += term;
932            k += 1;
933        }
934        Ok(sum * prefactor)
935    }
936}
937
938#[cfg(test)]
939mod tests {
940    use super::*;
941
942    const TOL: f64 = 1E-12;
943
944    #[test]
945    fn chbevl_val() {
946        assert::close(
947            chbevl(2.23525, &BESSI0_COEFFS_A),
948            0.139_251_866_282_289,
949            TOL,
950        );
951    }
952
953    #[test]
954    fn bessi0_small() {
955        assert::close(i0(3.74), 9.041_496_849_012_773, TOL);
956        assert::close(i0(-3.74), 9.041_496_849_012_773, TOL);
957        assert::close(i0(8.0), 427.564_115_721_804_74, TOL);
958    }
959
960    #[test]
961    fn bessi0_large() {
962        assert::close(i0(8.1), 469.500_606_710_121_4, TOL);
963        assert::close(i0(10.0), 2_815.716_628_466_254, TOL);
964    }
965
966    #[test]
967    fn bessi1_small() {
968        assert::close(i1(3.74), 7.709_894_215_253_694, TOL);
969        assert::close(i1(-3.74), -7.709_894_215_253_694, TOL);
970        assert::close(i1(0.0024), 0.001_200_000_864_000_207_2, TOL);
971        assert::close(i1(8.0), 399.873_136_782_559_9, TOL);
972    }
973
974    #[test]
975    fn bessi1_large() {
976        assert::close(i1(8.1), 439.484_308_910_358_44, TOL);
977        assert::close(i1(10.0), 2_670.988_303_701_255, TOL);
978    }
979
980    #[test]
981    fn bessel_iv_basic_limits() {
982        assert::close(bessel_iv(0.0, 0.0).unwrap(), 1.0, TOL);
983        assert::close(bessel_iv(1.0, 0.0).unwrap(), 0.0, TOL);
984    }
985
986    #[test]
987    fn bessel_iv_high_order() {
988        assert::close(
989            bessel_iv(60.0, 40.0).unwrap(),
990            0.071_856_419_684_526_32,
991            TOL,
992        );
993    }
994
995    #[test]
996    fn bessel_iv_low_order() {
997        assert::close(
998            bessel_iv(0.0, 1.0).unwrap(),
999            1.266_065_877_752_008_4,
1000            TOL,
1001        );
1002        assert::close(
1003            bessel_iv(0.0, 10.0).unwrap(),
1004            2_815.716_628_466_254_4,
1005            TOL,
1006        );
1007
1008        assert::close(
1009            bessel_iv(1.0, 10.0).unwrap(),
1010            2_670.988_303_701_254,
1011            TOL,
1012        );
1013        assert::close(
1014            bessel_iv(20.0, 10.0).unwrap(),
1015            0.000_125_079_973_564_494_78,
1016            TOL,
1017        );
1018    }
1019
1020    #[test]
1021    fn cf1_ik_checks() {
1022        assert::close(cf1_ik(0.0, 10.0).unwrap(), 0.948_599_825_954_845_8, TOL);
1023        assert::close(
1024            cf1_ik(10.0, 10.0).unwrap(),
1025            0.389_913_883_928_382_94,
1026            TOL,
1027        );
1028        assert::close(
1029            cf1_ik(60.0, 5.0).unwrap(),
1030            0.040_916_097_908_833_04,
1031            TOL,
1032        );
1033    }
1034
1035    #[test]
1036    fn cf2_ik_checks() {
1037        let (k1, k2) = cf2_ik(0.0, 2.0).unwrap();
1038        assert::close(k1, 0.113_893_872_749_533_53, TOL);
1039        assert::close(k2, 0.139_865_881_816_522_54, TOL);
1040
1041        let (k1, k2) = cf2_ik(5.0, 5.0).unwrap();
1042        assert::close(k1, 3.270_627_371_203_162e-2, TOL);
1043        assert::close(k2, 8.067_161_323_456_37e-2, TOL);
1044    }
1045
1046    #[test]
1047    fn temme_ik_series_checks() {
1048        let res = temme_ik_series(0.0, 0.0);
1049        assert!(res.is_err());
1050
1051        let (k1, k2) = temme_ik_series(0.0, 1.0).unwrap();
1052        assert::close(k1, 4.210_244_382_407_083_4e-1, TOL);
1053        assert::close(k2, 6.019_072_301_972_346e-1, TOL);
1054
1055        let (k1, k2) = temme_ik_series(0.5, 2.0).unwrap();
1056        assert::close(k1, 1.199_377_719_680_612_3e-1, TOL);
1057        assert::close(k2, 1.799_066_579_520_924_3e-1, TOL);
1058    }
1059
1060    #[test]
1061    fn bessel_ikv_temme_checks() {
1062        let (i, k) = bessel_ikv_temme(0.0, 1.0).unwrap();
1063        assert::close(i, 1.266_065_877_752_008_4, TOL);
1064        assert::close(k, 0.421_024_438_240_708_34, TOL);
1065
1066        let (i, k) = bessel_ikv_temme(5.0, 2.0).unwrap();
1067        assert::close(i, 0.009_825_679_323_131_702, TOL);
1068        assert::close(k, 9.431_049_100_596_468, TOL);
1069
1070        let (i, k) = bessel_ikv_temme(20.0, 2.0).unwrap();
1071        assert::close(i, 4.310_560_576_109_548E-19, TOL);
1072        assert::close(k, 5.770_856_852_700_242_4E16, TOL);
1073
1074        let (i, k) = bessel_ikv_temme(20.0, 2.0).unwrap();
1075        assert::close(i, 4.310_560_576_109_548E-19, TOL);
1076        assert::close(k, 5.770_856_852_700_242_4E16, TOL);
1077
1078        let (i, k) = bessel_ikv_temme(1.0, 10.0).unwrap();
1079        assert::close(i, 2_670.988_303_701_254, TOL);
1080        assert::close(k, 1.864_877_345_382_558_5E-5, TOL);
1081    }
1082
1083    #[test]
1084    fn bessel_ikv_asymptotic_uniform_checks() {
1085        let (i, k) = bessel_ikv_asymptotic_uniform(60.0, 40.0).unwrap();
1086        assert::close(i, 7.185_641_968_452_632e-2, TOL);
1087        assert::close(k, 9.649_278_749_222_319e-2, TOL);
1088
1089        let (i, k) = bessel_ikv_asymptotic_uniform(100.0, 60.0).unwrap();
1090        assert::close(i, 2.883_277_090_649_164e-7, TOL);
1091        assert::close(k, 1.487_001_275_494_647_4e4, TOL);
1092    }
1093
1094    #[test]
1095    fn proptest_i0_vs_log_i0() {
1096        use proptest::prelude::*;
1097
1098        proptest!(|(x: f64)| {
1099            // Skip NaN/infinite inputs
1100            prop_assume!(x.is_finite());
1101
1102            let i0_val = i0(x);
1103            let log_i0_val = log_i0(x);
1104
1105            // i0 should equal exp(log_i0) within tolerance
1106            // Only check when values are in reasonable range to avoid overflow
1107            if i0_val.is_finite() && log_i0_val.is_finite() && log_i0_val < 700.0 {
1108                assert::close(i0_val.ln(), log_i0_val, 1e-10);
1109            }
1110        });
1111    }
1112}
rv/misc/bessel.rs

rv/misc/
bessel.rs
