6

我有一个程序打算使用 Chudnovsky 算法来近似 pi，但是我的方程中一个非常小的项一直被四舍五入为零。

这是算法：

import math
from decimal import *
getcontext().prec = 100

pi = Decimal(0.0)
C = Decimal(12/(math.sqrt(640320**3)))
k = 0
x = Decimal(0.0)
result = Decimal(0.0)
sign = 1
while k<10:
    r = Decimal(math.factorial(6*k)/((math.factorial(k)**3)*math.factorial(3*k)))
    s = Decimal((13591409+545140134*k)/((640320**3)**k))
    x += Decimal(sign*r*s)
    sign = sign*(-1)
    k += 1
result = Decimal(C*x)
pi = Decimal(1/result)


print Decimal(pi)

如果没有“小数”项，等式可能会更清晰。

import math

pi = 0.0
C = 12/(math.sqrt(640320**3))
k = 0
x = 0.0
result = 0.0
sign = 1
while k<10:
    r = math.factorial(6*k)/((math.factorial(k)**3)*math.factorial(3*k))
    s = (13591409+545140134*k)/((640320**3)**k)
    x += sign*r*s
    sign = sign*(-1)
    k += 1
result = C*x
pi = 1/result

print pi

问题在于“s”变量。对于 k>0，它总是为零。例如，在 k=1 时，s 应该等于大约 2.1e-9，但它只是零。因此，我在第一个 =0 之后的所有条款。如何让 python 计算 s 的确切值而不是将其四舍五入为 0？

4

4 回答 4

3

尝试：

s = Decimal((13591409+545140134*k)) / Decimal(((640320**3)**k))

你正在做的算术是原生 python - 通过允许 Decimal 对象执行你的除法，你应该消除你的错误。

那么，在计算r.

于 2013-04-11T01:01:20.237 回答
2

一些评论。

如果您使用的是 Python 2.x，则/返回整数结果。如果您想要 Decimal 结果，请先将至少一侧转换为 Decimal。

math.sqrt()仅返回约 16 位精度。由于您的 C 值只能精确到 ~16 位，因此您的最终结果只能精确到 16 位。

于 2013-04-11T01:02:35.450 回答
2

如果您在 Python 2.x 中进行数学运算，您可能应该将这一行放入每个模块中：

from __future__ import division

这改变了除法运算符的含义，以便在需要时返回一个浮点数以给出（更接近）精确的答案。历史行为是x / y返回一个intif both xand yare ints，这通常会迫使答案向下舍入。

如果需要，返回一个浮点数通常被认为是一种更好的方式来处理像 Python 这样鼓励鸭子类型的语言中的除法，因为你可以只担心你的数字的值，而不是为不同的类型获得不同的行为。

在 Python 3 中，这实际上是默认设置，但由于旧程序依赖于除法运算符的历史行为，因此人们认为这种更改过于向后不兼容，无法在 Python 2 中进行。这就是为什么您必须显式打开它的原因与__future__进口。我建议始终在可能进行任何数学运算的任何模块中添加该导入（或者只是任何模块，如果您可能会被打扰）。你几乎永远不会对它的存在感到不安，但是没有它是我不得不追逐的一些晦涩的错误的原因。

于 2013-04-11T02:13:46.603 回答
0

我觉得 's' 的问题在于所有项都是整数，因此您正在做整数数学。一个非常简单的解决方法是3.0在分母中使用。计算中只需要一个浮点数即可返回一个浮点数。

于 2013-04-11T01:03:03.320 回答