Python中如何避免向零舍入

Question

我有一个用Chudnovsky算法近似计算圆周率的程序，但方程中一个非常小的项一直被舍入为零。

这是算法：

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，但实际上它只是零。因此，我的第一个项之后的所有项都等于零。我如何让 python 计算 s 的精确值，而不是将其舍入为零？

- jb37

---

这是哪个版本的Python？ - Asad Saeeduddin

你尝试过s = float(13591409+545140134*k)/((640320**3)**k)吗？ - karthikr

您可能会对Nick Craig-Wood关于Chudnovsky算法的文章感兴趣。此外，如果您想要100位数字，则我不认为k < 10足够大。 - unutbu

@karthikr，我会把你的评论变成答案。这对我有用。 - AnthonyVO

你可以使用from decimal import Decimal代替from decimal import *。 - Solomon Ucko

Answer 1

尝试：

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

你正在进行的算术操作是使用原生Python实现的 - 通过允许Decimal对象执行你的除法，你应该能够消除错误。

当计算r时，你可以做同样的操作。

- Sajjan Singh

Answer 2

如果您正在使用Python 2.x进行数学运算，那么您应该在每个模块中添加以下行：

from __future__ import division

这改变了除法运算符的含义，使其在需要返回浮点数以给出（更接近）精确答案时返回浮点数。历史行为是当x和y都是int时，x / y返回一个int，这通常会强制向下舍入答案。

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

在Python 3中，这实际上是默认设置，但由于旧程序依赖于除法运算符的历史行为，人们认为这种更改过于不兼容Python 2。这就是为什么你必须明确地使用__future__导入它。我建议在任何可能进行任何数学运算的模块中始终添加该导入（或者如果你可以费心的话，在任何模块中都添加它）。你几乎永远不会因为它存在而感到沮丧，但没有它会导致我追踪一些晦涩的错误。

- Ben

Answer 3

几点评论。

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

math.sqrt()只返回约16位有效数字。由于C的值仅精确到约16位数字，因此最终结果仅精确到16位数字。

- casevh

---

有没有一种方法可以得到更精确的平方根？ - jb37

Decimal 实例有 sqrt() 方法，所以尝试 C = 12/(Decimal(640320**3).sqrt())。 - casevh

谢谢casevh，我之前不知道这个。 - jb37

Answer 4

我觉得“s”的问题在于所有的术语都是整数，因此你正在进行整数运算。一个非常简单的解决方法是在分母中使用3.0。只需要在计算中使用一个浮点数就可以返回一个浮点数。

- Nick Burns