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.0001
或 abs(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
我编写了这段代码来使用欧几里德算法计算有理数 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.0001
或 abs(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