SheetJS bessel besseli (formulajs BESSELI) wrong for x >> n, and for x < 0
bessel 1.0.2 (npm, ~1.2M/week) besseli(100, 3) returns 8.21e41; the exact value is 1.026e42 (scipy.special.iv). besseli(574.7, 5) is 65% low. besselj/bessely/besselk are fine to ~4e-7. It is Numerical Recipes bessi: the Miller backward recurrence starts at m = n + sqrt(40n), which is too low once x is several times m (it starts failing near x=20 for n=2 and x=80 for n=10). Fix: use max(n, |x|) in place of n when choosing m. Separately, bessel0/bessel1 test x <= 3.75 instead of |x| <= 3.75, so every negative x below -3.75 gets the small-argument polynomial (this is the cause behind SheetJS/bessel#5).
Reproduced on npm bessel 1.0.2 (node 22) against SciPy 1.18.1 iv. besseli(100, 3) = 8.206e41 vs 1.0263e42 (20% low), and besseli(574.7, 5) = 2.224e247 vs 6.322e247 (65% low). At the onset points errors are still small: (20, 2) is 6e-5 relative, (15, 2) 7e-6, (80, 10) 1.7e-5, (60, 10) 3e-7. Positive small-argument values are fine, e.g. besseli(5, 0) to 2.6e-9. The negative-x bug is source-confirmed: bessel0 tests
x <= 3.75(line 185), and bessel1 testsx < 3.75. It grows fast with |x|. For besseli(x, 0): x = -5 gives 27.238 vs 27.240 (7e-5), -10 gives 2545 vs 2816 (10% low), -20 gives 3.30e6 vs 4.36e7 (13x low), and -50 gives 1.5e11 vs 2.9e20. besseli(x, 1) behaves the same. Spot checks: besselj(30, 2) and besselk(5, 1) agree to 9e-10 and 1.5e-8. I didn't test the max(n, |x|) fix.Follow-up: I tested the fix after all. Replacing n with max(n, |x|) when choosing m (line 206, in a copy of 1.0.2) brings every x >> n case to the bessel0 polynomial's own accuracy. Relative errors: (100, 3) 4.7e-7, (574.7, 5) 1.7e-7, (20, 2) 7e-8, (80, 10) 4.6e-7, (60, 10) 3.8e-7. So the patched besseli is good to about 5e-7, not to double precision.
Same wrong values from =BESSELI(100,3) in HyperFormula 3.5.0, fast-formula-parser and hot-formula-parser 4.0.0, since all three use the bessel package: 8.2059e41 vs the exact 1.0263e42.