c - SciPy中16384.0*floor(v/16384.0)的使用说明

Question

当我发现这个函数时，我正在研究 scipy 源代码 ( 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 = v - 16384.0 * floor(v/16384.0) 这行感到困惑。我知道 16384=16*1024，但我想知道为什么在当前代码中使用它。

该代码用于贝塞尔函数的计算，特别是Y_v(z)。当 v 为负数且接近整数时，标准反射公式(包括 cos(pi*nu) 和 sin(pi*nu))会产生无效结果。我敢说它是用来平息这种情况的，但我不知道它是如何以及为什么起作用的。

Answer 1

代码想要检查已知为整数的 float 的奇偶性。

在实践中，重要的是将整数阶的情况与任意阶的情况分开处理，因为 sin/cos 反射公式最终可能会接近 0*inf 的情况，这会破坏数值精度。

直接转换 (int)v 可能会溢出，因此在 i = v 之前取除以某个偶数的余数；如果 (i & 1).

为什么是 16384.0 而不是 2.0？我想，重点是减少浮点运算中的精度损失，从而弄乱结果—— float 的奇偶校验可能来自尾数的最后一位。这个技巧取自其他有点旧的代码(从 80 年代开始)，并且可能是在考虑非 IEEE fp 算法的情况下编写的；假设 IEEE fp，我想有更好的方法来做到这一点。