4

使用使用 Gauss-Legendre 算法计算 pi 的 python 脚本时出现此错误。在得到这个之前，您最多只能使用 1024 次迭代：

    C:\Users\myUsernameHere>python Desktop/piWriter.py
    End iteration: 1025
    Traceback (most recent call last):
      File "Desktop/piWriter.py", line 15, in <module>
        vars()['t' + str(sub)] = vars()['t' + str(i)] - vars()['p' + str(i)] * math.
    pow((vars()['a' + str(i)] - vars()['a' + str(sub)]), 2)
    OverflowError: long int too large to convert to float

这是我的代码：

import math

a0 = 1
b0 = 1/math.sqrt(2)
t0 = .25
p0 = 1

finalIter = input('End iteration: ')
finalIter = int(finalIter)

for i in range(0, finalIter):
        sub = i + 1
        vars()['a' + str(sub)] = (vars()['a' + str(i)] + vars()['b' + str(i)])/ 2
        vars()['b' + str(sub)] = math.sqrt((vars()['a' + str(i)] * vars()['b' + str(i)]))
        vars()['t' + str(sub)] = vars()['t' + str(i)] - vars()['p' + str(i)] * math.pow((vars()['a' + str(i)] - vars()['a' + str(sub)]), 2)
        vars()['p' + str(sub)] = 2 * vars()['p' + str(i)]
        n = i

pi = math.pow((vars()['a' + str(n)] + vars()['b' + str(n)]), 2) / (4 * vars()['t' + str(n)])
print(pi)

理想情况下，我希望能够插入一个非常大的数字作为迭代值，稍后再回来查看结果。

任何帮助表示赞赏！谢谢！

4

4 回答 4

4

浮点数只能表示最大为 sys.float_info.max 或 1.7976931348623157e+308 的数字。一旦你有一个超过 308 位（或更多）的 int，你就会被卡住。当 p1024 有 309 位时，您的迭代失败：

179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137216L

您必须为 pi 找到一种不同的算法，它不需要这么大的值。

实际上，您必须小心周围的浮点数，因为它们只是近似值。如果你修改你的程序来打印 pi 的逐次逼近，它看起来像这样：

2.914213562373094923430016933707520365715026855468750000000000
3.140579250522168575088244324433617293834686279296875000000000
3.141592646213542838751209274050779640674591064453125000000000
3.141592653589794004176383168669417500495910644531250000000000
3.141592653589794004176383168669417500495910644531250000000000
3.141592653589794004176383168669417500495910644531250000000000
3.141592653589794004176383168669417500495910644531250000000000

换句话说，仅在 4 次迭代之后，您的近似值就不再变得更好了。这是由于您使用的浮点数不准确，可能以1/math.sqrt(2). 计算 pi 的许多位需要非常仔细地理解数字表示。

于 2012-12-20T19:06:24.467 回答
1

如上一个答案所述，该float类型具有数字大小的上限。在典型实现中，sys.float_info.max是 1.7976931348623157e+308，它反映了 64 位浮点数中指数字段使用 10 位加号。（注意 1024*math.log(2)/math.log(10) 大约是 308.2547155599。）

您可以使用Decimal number 类型在指数大小上再增加半个十进制数。这是一个示例（从 ipython 解释器会话中截取）：

In [48]: import decimal, math    
In [49]: g=decimal.Decimal('1e12345')    
In [50]: g.sqrt()
Out[50]: Decimal('3.162277660168379331998893544E+6172')
In [51]: math.sqrt(g)
Out[51]: inf

这说明小数的sqrt()函数在比 更大的数字时正确执行math.sqrt()。

于 2012-12-20T19:26:07.493 回答
0

如上所述，获得大量数字会很棘手，但是看着所有这些数字vars会伤害我的眼睛。因此，在 (1) 将您的使用替换为vars字典和 (2) 使用**而不是数学函数之后，这是您的代码版本：

a, b, t, p = {}, {}, {}, {}
a[0] = 1
b[0] = 2**-0.5
t[0] = 0.25
p[0] = 1

finalIter = 4

for i in range(finalIter):
    sub = i + 1
    a[sub] = (a[i] + b[i]) / 2
    b[sub] = (a[i] * b[i])**0.5
    t[sub] = t[i] - p[i] * (a[i] - a[sub])**2
    p[sub] = 2 * p[i]
    n = i

pi_approx = (a[n] + b[n])**2 / (4 * t[n])

vars我没有用来玩游戏，而是使用字典来存储值（那里有官方 Python 教程的链接），这使您的代码更具可读性。您现在甚至可以看到一两个优化。

如评论中所述，您实际上不需要存储所有值，只需要存储最后一个值，但我认为更重要的是，您了解如何在不动态创建变量的情况下做事。除了 a dict，您也可以简单地将值附加到 a list，但列表始终是零索引的，您不能轻易“跳过”并在任意索引处设置值。在使用算法时，这有时会让人感到困惑，所以让我们从简单的开始。

无论如何，上面给了我

>>> print(pi_approx)
3.141592653589794
>>> print(pi_approx-math.pi)
8.881784197001252e-16
于 2012-12-20T21:25:21.250 回答
0

一个简单的解决方案是安装和使用现在支持 Python 3 的任意精度mpmath模块。但是，由于我完全同意 DSM 的观点，即您使用动态vars()创建变量是实现算法的一种不受欢迎的方式，因此我基于我对他重写您的代码的回答并[微不足道地]对其进行了修改以用于进行mpmath计算。

如果您坚持使用vars()，您可能会做类似的事情——尽管我怀疑它可能会更难，结果肯定更难阅读、理解和修改。

from mpmath import mpf  # arbitrary-precision float type

a, b, t, p = {}, {}, {}, {}
a[0] = mpf(1)
b[0] = mpf(2**-0.5)
t[0] = mpf(0.25)
p[0] = mpf(1)

finalIter = 10000

for i in range(finalIter):
    sub = i + 1
    a[sub] = (a[i] + b[i]) / 2
    b[sub] = (a[i] * b[i])**0.5
    t[sub] = t[i] - p[i] * (a[i] - a[sub])**2
    p[sub] = 2 * p[i]
    n = i

pi_approx = (a[n] + b[n])**2 / (4 * t[n])
print(pi_approx)  # 3.14159265358979
于 2013-09-25T18:11:19.310 回答