使用itertools和Leibniz级数逼近Pi的值



编写一个程序,通过对以下级数的项求和来近似ππ的值:

(4/1) - (4/3) + (4/5) - (4/7) + (4/9) - (4/11) + ……

这个代码出了什么问题

def main():
n=eval(input("Enter N: "))
x,y=0,0
for i,j in itertools.product ((1,1+4*round((n//2)),4),(-3,-3-4*int((n/2)),-4)):
x=x+(4/i)
y=y+(4/j)
print(x+y)

此代码有什么问题?

很多:

  • 可读性:它没有解释其目的
  • 可用性:用户必须输入N,而不知道它有什么好处
  • 可维护性:for循环有六个常量,称为";幻数";没有任何解释
  • 可读性:莱布尼茨系列真的很简单。为什么代码如此复杂
  • 正确性:您实现了main(),但从不调用它
  • 安全:eval()很危险

要解决这些问题:

  • 添加注释、代码的作用、公式的描述等
  • 改进问题以供用户输入
  • 去掉复杂的代码,实现一个简单的求和算法
  • 使用尽可能接近的公式以避免混淆。4*已经是一个拉伸
  • 使用int()而不是eval()从用户输入中获取数字

建议结果:

# Approximation of Pi using the Leibniz series
# See https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80
# k is used as described on Wikipedia
# n is the number of approximation steps. n → ∞ to get close approximations
n = int(input("Enter the number of approximation steps:"))
pi = 0
for k in range(0, n):
pi +=  4* (-1)**k / (2*k+1)
print(pi)

对于10.000.000步,它给出3.14159255,所以你真的需要大量的N.


如果你想追求输出的精度,那么使用@Patrick Artner建议的fractions是个好主意。但这确实降低了计算速度。这里有一个比较:

N      | FP                   | Fractions            | Diff
-------+----------------------+----------------------+---------
10000  | 3.1414926535900 345  | 3.1414926535900 434  | + 89E-16
20000  | 3.141542653589824 8  | 3.141542653589824 3  | -  5E-16
30000  | 3.14155932025646 20  | 3.14155932025646 93  | + 73E-16
40000  | 3.14156765358979 85  | 3.14156765358979 70  | - 15E-16
50000  | 3.1415726535897 814  | 3.1415726535897 950  | +136E-16
60000  | 3.1415759869231 020  | 3.1415759869231 277  | +257E-16
70000  | 3.141578367875 4820  | 3.141578367875 5083  | +263E-16
80000  | 3.1415801535897 496  | 3.1415801535897 936  | +440E-16
90000  | 3.1415815424786 238  | 3.1415815424786 824  | +586E-16
100000 | 3.1415826535897 198  | 3.1415826535897 935  | +740E-16
110000 | 3.141583562680 6436  | 3.141583562680 7027  | +591E-16

浮动数学和四舍五入将打断您的脖子降低近似的质量。

您可以使用分数模块:

import fractions

while True:
try:
n = int(input("Amount of terms used for estimation of pi by Leibniz series: "))
if n > 0:
break
except ValueError:
print("Try again and be sensible about your inputs.")
sign = 1
pi = fractions.Fraction(0,1)
for nn in range(n):
term = fractions.Fraction(sign * 4, 2*nn + 1)
print(f"Approx nr: {nn+1})   pi before: {str(pi):>20}         approx refinement: {('+' if sign > 0 else '') + str(term):>20}", end="")
pi += term
print(f"   new pi:{pi} ({float(pi):2.8})")
sign *= -1
print(float(pi))

获取(输入9(:

Approx nr: 1)   pi before:                    0         approx refinement:                   +4   new pi:4 (4.0)
Approx nr: 2)   pi before:                    4         approx refinement:                 -4/3   new pi:8/3 (2.6666667)
Approx nr: 3)   pi before:                  8/3         approx refinement:                 +4/5   new pi:52/15 (3.4666667)
Approx nr: 4)   pi before:                52/15         approx refinement:                 -4/7   new pi:304/105 (2.8952381)
Approx nr: 5)   pi before:              304/105         approx refinement:                 +4/9   new pi:1052/315 (3.3396825)
Approx nr: 6)   pi before:             1052/315         approx refinement:                -4/11   new pi:10312/3465 (2.9760462)
Approx nr: 7)   pi before:           10312/3465         approx refinement:                +4/13   new pi:147916/45045 (3.2837385)
Approx nr: 8)   pi before:         147916/45045         approx refinement:                -4/15   new pi:135904/45045 (3.0170718)
Approx nr: 9)   pi before:         135904/45045         approx refinement:                +4/17   new pi:2490548/765765 (3.2523659)
3.252365934718876

我认为您打算使用ziprange,而不是itertools.product,它将为您提供的两个序列的每个组合创建一个成对的序列。

例如,

  • 产物((1,2(,(3,4((的产率为(1,3(
  • 而zip(((1,2(,(3,4((产生(1,3(和(2,4(,像zip一样一个接一个地匹配位置

更正后的代码如下,一些注释也很有用:

def main():
n=eval(input("Enter N: "))
x,y=0,0
for i,j in zip(range(1,1+4*round((n//2)),4), range(-3,-3-4*int((n/2)),-4)):
x=x+(4/i) # sum of positive parts
y=y+(4/j) # sum of negative parts
print(x+y)

最新更新