Python 2.7 - 连分数展开 - 理解错误

Python 2.7 - Continued Fraction Expansion - Understanding the error

我编写了这段代码来使用欧几里德算法计算有理数 N 的连分式展开:

from __future__ import division

def contFract(N):
    while True:
        yield N//1
        f = N - (N//1)
        if f == 0:
            break
        N = 1/f

如果说 N 是 3.245,则函数永远不会结束,因为显然 f 永远不会等于 0。展开式的前 10 项是:

[3.0, 4.0, 12.0, 3.0, 1.0, 247777268231.0, 4.0, 1.0, 2.0, 1.0]

这显然是一个错误,因为实际的扩展只是:

[3;4,12,3,1] or [3;4,12,4]

这里的问题是什么引起的?这是某种舍入错误吗?

问题是您正在测试 f == 0(整数 0),这对于浮点数来说几乎从来都不是真的。所以循环永远持续下去。

相反,检查等同于 0 (which can be wrong sometimes) 的某些精度:

>>> from __future__ import division
>>>
>>> def contFract(N):
...     while True:
...         yield N//1
...         f = N - (N//1)
...         if f < 0.0001:  # or whatever precision you consider close enough to 0
...             break
...         N = 1/f
...
>>>
>>> list(contFract(3.245))
[3.0, 4.0, 12.0, 3.0, 1.0]
>>>

如果 f 可能为负数,请执行 -0.0001 < f < 0.0001abs(f) < 0.0001。这也被认为是不准确的,请参阅链接文章。

根据我的评论,使用 int(N) 而不是 N//1 因为它更清晰 - 稍微 效率较低:

>>> import timeit
>>> timeit.timeit('N = 2.245; N//1', number=10000000)
1.5497028078715456
>>> timeit.timeit('N = 2.245; int(N)', number=10000000)
1.7633858824068103

显然这个错误是由于 4.0 除以 1 的整数。4.0 在浮点表示中有一个错误(如果你对浮点表示有想法)。所以实际上 int(4.0) < 4。这导致 N//1 变为 3.0,分数 f 类似于 0.999999...。所以这永远不会结束。在每次迭代中打印 N 和 f 并尝试。你会明白的。

您正在使用 float 进行运算,不幸的是,有些数字无法用二进制表示。

有两种解决方法,首先 - 假设您的数字是 "close enough"(甚至新的 Python 3.5.2 引入了 math.isclose),或者您使用的是不同的浮动的实现,例如Decimal 你可以得到正确的结果。

N.b。这就是为什么对于所有金融系统,没有人使用浮点数,只使用 int/bigint 或小数。

 In [21] > N = decimal.Decimal('3.245')

 In [22] > while True:
    print 'N: %s' % (N//1,)
    f = N - N//1
    print 'f: %s' % f
    if f == 0:
        break
    N = 1/f

N: 3
f: 0.245
N: 4
f: 0.081632653061224489795918367
N: 12
f: 0.25000000000000000000000005
N: 3
f: 0.999999999999999999999999200
N: 1
f: 8.00E-25
N: 1250000000000000000000000
f: 0