Explanation of the use of 16384.0*floor(v/16384.0) in SciPy

Question

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.

Answer 1

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.

Answer 2

This is about maintaining precision of a double variable.

One step at a time:

i = v - 16384.0 * floor(v / 16384.0);

The expression v/16384.0 is calculated first.
The expression 16384 is a precision limit as we'll explore.

Next, the floor function is applied to the result of v/16384.0. The result is to truncate the division to whole numbers, like integer division but with floating point.

The whole number is multiplied by 16384.0 to convert the variable v to the nearest whole number divisible by 16384.

The final subtraction produces the remainder part when divided by 16384.

Any values smaller than 1/16384 will be truncated and not considered. So the precision of v is restricted to 1/16384.