Explanation of the use of 16384.0*floor(v/16384.0) in SciPy
c, scientific-computing, scipy
Solution
The code wants to check the parity of a floating point number that is known to be an integer.
In practice, it is important to handle the integer order case separately from the arbitrary order case, since the sin/cos reflection formula can end up close to a 0*inf situation which destroys numerical precision.
Direct cast `(int)v` can overflow, so a remainder with division by some even number is taken before `i = v; if (i & 1)`.
Why 16384.0 rather than 2.0? The point, I suppose, was to reduce loss of precision in the floating point operations messing up the result --- the parity of the floating point number may come from the last digits in the mantissa. This trick is taken from other somewhat old code (from the 80s), and probably has been written with non-IEEE fp arithmetic in mind; assuming IEEE fp, I suppose there are better ways to do this.
Problem
I was studying the scipy source code when I found this function (1) ``` static int reflect_jy(npy_cdouble *jy, double v) { /* NB: Y_v may be huge near negative integers -- so handle exact * integers carefully */ int i; if (v != floor(v)) return 0; i = v - 16384.0 * floor(v / 16384.0); if (i & 1) { jy->real = -jy->real; jy->imag = -jy->imag; } return 1; ``` } I am puzzled by the line `i = v - 16384.0 * floor(v / 16384.0)`. I am aware that 16384=16*1024, but I would want to know why this is used in the present code. The code is used in the computation of Bessel functions, specifically `Y_v(z)`. When `v` is negative and near an integer, the standard reflection formulae (which include cos(pi*nu) and sin(pi*nu)) produce invalid results. I venture that this it is used to quell this, but I don't know how and why this works.