在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

代码想要检查已知为整数的浮点数的奇偶校验.

在实践中,重要的是与任意顺序情况分开处理整数阶情况,因为sin/cos反射公式可以最终接近破坏数值精度的0*inf情况.

直接投射(int)v可以溢出,因此之前采用除以某个偶数的余数i = v; if (i & 1).

为什么16384.0而不是2.0？我想,重点是减少浮点运算中精度损失的结果 - 浮点数的奇偶校验可能来自尾数中的最后一位数.这个技巧来自其他有点旧代码(来自80年代),并且可能是用非IEEE fp算法编写的; 假设IEEE fp,我想有更好的方法来做到这一点.