scipy.special.ellipkinc中的bug - 不完整的椭圆积分

Question

在使用SciPy的scipy.special.ellipeinc和ellipkinc时，似乎会出现一些数值不稳定的情况。例如，

>>> from scipy.special import ellipkinc
>>> ellipkinc(0.9272952180016123, 0.68359375000000011)
nan
>>> ellipkinc(0.9272952180016123, 0.6835937500000002)
2.0518660200390668
>>> ellipkinc(0.9272952180016123, 0.68359375)
1.0259330100195332
>>> ellipkinc(0.9272952180016123, 0.68359374)
1.0259330081166262

这种情况发生在k的平方乘以sin的平方（phi）接近0.3的时候，但这里的椭圆积分本身没有什么特别之处，所以推测是数值计算的问题。我对这个算法的内部工作原理了解不够，无法判断具体哪里出错了，那我该怎么办呢？我想到了： round(0.68359375000000011,8)，不过这样做肯定会让我的代码变慢吧？

Answer 1

（这更像是对问题的详细评论，而不是解决方案。）

这看起来像是ellipkinc中的一个错误。我发现函数返回nan时，只有一个浮点值，而当函数返回两倍的“正确”值时，周围有四个相邻的浮点值：

In [91]: phi = 0.9272952180016123

In [92]: mbad = 0.68359375000000011

In [93]: m = np.nextafter(mbad, 0)

In [94]: mvals = []

In [95]: for j in range(10):
   ....:     mvals.append(m)
   ....:     m = np.nextafter(m, 1)
   ....:     

In [96]: mvals
Out[96]: 
[0.68359375,
 0.68359375000000011,
 0.68359375000000022,
 0.68359375000000033,
 0.68359375000000044,
 0.68359375000000056,
 0.68359375000000067,
 0.68359375000000078,
 0.68359375000000089,
 0.683593750000001]

In [97]: ellipkinc(phi, mvals)
Out[97]: 
array([ 1.02593301,         nan,  2.05186602,  2.05186602,  2.05186602,
        2.05186602,  1.02593301,  1.02593301,  1.02593301,  1.02593301])