# What algorithm does (or did?) Excel use for Bessel functions that is discontinuous at x=8?

Writing this comment reminded me of something I noticed years ago about evaluating Bessel functions of the first kind $$J_n(x)$$ in Excel. (BESSELJ)

I don't use Excel now but at the time I'd checked MacOS and Windows computers with several different released versions of Excel and they all had the same problem.

This was circa 2015.

When the argument passed through $$x=8$$, there was a step discontinuity of order $$10^{-5}$$ for $$n<9$$.

Plotting the first difference will show that the Excel values are step-wise discontinuous at $$8$$.

Question: Is it possible to know which algorithm had this behavior?

Example from Excel:

   x         J1             J3            J5
7.9997,  0.234593634,  -0.291131046,  0.185841208,
7.9998,  0.234607873,  -0.291131434,  0.185819063,
7.9999,  0.234622108,  -0.291131819,  0.185796918,
8.0000,  0.234628387,  -0.291125242,  0.185775021,
8.0001,  0.234642617,  -0.291125622,  0.185752874,
8.0002,  0.234656846,  -0.291126,     0.185730726,
8.0003,  0.234671072,  -0.291126376,  0.185708577


edit: per request:

n    x=7.99999952316284            x=8             x=8.00000095367431
--   ------------------    ------------------      ------------------
0     0.17165091838403      0.171650807317694       0.171650583550973
1     0.2346362739686       0.234628386507074       0.234628522235498
2    -0.112991846395527    -0.112993710690925      -0.112993459984572
3    -0.291132200533783    -0.291125241852537      -0.112993459984572

• Without the actual implementation, you can only guess. My guess is, that the integer order Bessel functions are computed from $J_0, J_1$ with the recurrence relation. And (and least) in the $J_1$ implementation there is a switch between two approximations. (E.g. netlib.org/fdlibm has a switch at $x=8$). Do the even-order functions are 'discontinuous' too? Can you check with real order Bessel functions, e.g. $J_{3\pm0.00001}(8)?$ Commented Sep 28, 2018 at 10:18
• The starting point for me for computing special functions has always been the DLMF. According to that (dlmf.nist.gov/10.74#i), it could be a switch between a power series and an asymptotic expansion, but obviously that doesn't prove anything. Commented Sep 28, 2018 at 10:39
• Out of curiosity, can you try 7.999999523162841796875 and 8.00000095367431640625 (printing the output to full precision), the two Float32 numbers next to 8? Commented Sep 28, 2018 at 18:02
• @Kirill the jump is definitely where you thought: BESSELJ(7.99999952316284,3)=-0.291132200533783000 but BESSELJ(8.00000095367431,3)=-0.291125245492850000 forcing digits until they're all zero in excel 2013 for windows. The latter matches BESSELJ(8,3) to 8 dp Commented Oct 2, 2018 at 14:53
• @uhoh, luckily I can truthfully say that my knowledge is so outdated as to be of little use to anyone - which is good as I provide similar assistance for a lot of other software Commented Oct 7, 2018 at 14:05

It is clear that Microsoft Excel's implementation of the BesselJ function for $$N=1$$ is identical to this algorithm for arguments $$X < 8$$:

double ax,z;
double xx,y,ans,ans1,ans2;

if ((ax=fabs(x)) < 8.0) {
y=x*x;
ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
+y*(99447.43394+y*(376.9991397+y*1.0))));
ans=ans1/ans2;


You can implement that equation in Excel and compare it to Excel's result for BESSELJ(X,1). They are identical -- and I would note, identically erroneous for all values $$0 < X < 8$$, which proves Excel used this equation. For instance, compared to the exact result, the errors of BESSELJ(7,1) and $$X=7$$ plugged into the above are both approximately $$2.24457 \times 10^{-9}$$.

However, Excel encoded something else for values of $$X\geq8$$, and the error is considerably worse (e.g. 7.96012E-6 for $$X=8$$). The following code gives errors that are four orders of magnitude more accurate:

If ax=fabs(x) >= 8:

z=8.0/ax;
y=z*z;
xx=ax-2.356194491;
ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+y*(0.2457520174e-5+y*(-0.240337019e-6))));
ans2=0.04687499995+y*(-0.2002690873e-3
+y*(0.8449199096e-5+y*(-0.88228987e-6
+y*0.105787412e-6)));
ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
if (x < 0.0) ans = -ans;


I thought perhaps Excel attempted to implement the above and made a typo in one of the values, but I tried varying each of the terms in turn and was unable to replicate the values BESSELJ(X,1)` computes.