尝试定义 pi 的欧拉近似值之一，获取 'list 和 'int' 不支持的操作数类型

Question

我正在尝试定义一个函数，该函数将使用 Euler 的一种方法在 python 中逼近 pi。他的公式如下：

到目前为止我的代码是这样的：

def pi_euler1(n):
    numerator = list(range(2 , n))
    for i in numerator:
        j = 2
        while i * j <= numerator[-1]:
            if i * j in numerator:
                numerator.remove(i * j)
            j += 1
    for k in numerator:
        if (k + 1) % 4 == 0:
            denominator = k + 1
        else:
            denominator = k - 1
    #Because all primes are odd, both numbers inbetween them are divisible by 2,
    #and by extension 1 of the 2 numbers is divisible by 4
    term = numerator / denominator

我知道这是错误的，也不完整。我只是不太确定我之前提到的 TypeError 到底是什么意思。我只是坚持下去，我想创建一个术语列表，然后找到他们的产品。我在正确的路线上吗？

更新：我已经解决了这个问题，修复了由于 msconi 和 Johanc 而普遍存在的明显错误，现在使用以下代码：

import math
def pi_euler1(n):
    numerator = list(range(2 , 13 + math.ceil(n*(math.log(n)+math.log(math.log(n))))))
    denominator=[]
    for i in numerator:
        j = 2
        while i * j <= numerator[-1]:
            if (i * j) in numerator:
                numerator.remove(i * j)
            j += 1
    numerator.remove(2)
        for k in numerator:
            if (k + 1) % 4 == 0:
                denominator.append(k+1)
            else:
                denominator.append(k-1)
        a=1
        for i in range(n):
            a *= numerator[i] / denominator[i]
        return 4*a

这似乎可行，当我尝试绘制 pi 的半轴刻度误差图时，出现域错误，但我需要将范围的上限更改为 n+1因为 log(0) 未定义。谢谢大家

Answer 1

在 term = numerator / denominator 中，您将列表除以数字，这没有意义。在循环中将 k 除以分母，以便将 numerator 元素逐一用于方程的每个因数。然后你可以将它们重复乘以术语 term *= i / denominator，你在开始时将其初始化为 term = 1.

另一个问题是第一个循环，它不会给你第一个 n 个素数。例如，对于 n=3、list(range(2 , n)) = [2]。因此，您将获得的唯一素数是 2。

Answer 2

这里是经过一些小修改以使其正常工作的代码：

import math
def pi_euler1(n):
    lim = n * n + 4
    numerator = list(range(3, lim, 2))
    for i in numerator:
        j = 3
        while i * j <= numerator[-1]:
            if i * j in numerator:
                numerator.remove(i * j)
            j += 2
    euler_product = 1
    for k in numerator[:n]:
        if (k + 1) % 4 == 0:
            denominator = k + 1
        else:
            denominator = k - 1
        factor = k / denominator
        euler_product *= factor
    return euler_product * 4

print(pi_euler1(3))
print(pi_euler1(10000))
print(math.pi)

输出：

3.28125
3.148427801913721
3.141592653589793

备注：

您只需要奇数素数，因此您可以从奇数列表开始。
j可以从3开始，递增2步。实际上，j可以从i开始，因为[=17=的所有倍数] 小于 i*i 的已经被删除了。
一般来说，从要迭代的列表中删除元素是非常糟糕的做法。参见例如this post。在内部，Python 在它迭代的列表中使用索引。巧合的是，在这种特定情况下这不是问题，因为只有大于当前的数字才会被删除。
此外，从一个很长的列表中删除元素非常慢，因为每次都需要移动整个列表来填补空白。因此，最好使用两个单独的列表。
你没有计算结果，你也没有return它。
如您所见，这个公式收敛得非常慢。
评论中提到，之前的版本将n解释为最高素数的极限，而实际上n应该是素数的个数。我修改了代码来纠正这个问题。在上面的版本中有一个粗略的限制；下面的版本尝试对限制进行更严格的近似。

这是修改后的版本，没有从您迭代的列表中删除。它不是删除元素，而是标记它们。这样就快多了，所以可以在合理的时间内使用更大的n：

import math
def pi_euler_v3(n):
    if n < 3:
        lim = 6
    else:
        lim = n*n
        while lim / math.log(lim) / 2 > n:
            lim //= 2
    print(n, lim)

    numerator = list(range(3, lim, 2))
    odd_primes = []
    for i in numerator:
        if i is not None:
            odd_primes.append(i)
            if len(odd_primes) >= n:
                break
            j = i
            while i * j < lim:
                numerator[(i*j-3) // 2] = None
                j += 2
    if len(odd_primes) != n:
       print(f"Wrong limit calculation, only {len(odd_primes)} primes instead of {n}")
    euler_product = 1
    for k in odd_primes:
        denominator = k + 1 if k % 4 == 3 else k - 1
        euler_product *= k / denominator
    return euler_product * 4

print(pi_euler_v2(100000))
print(math.pi)

输出：

3.141752253548891
3.141592653589793

尝试定义 pi 的欧拉近似值之一，获取 'list 和 'int' 不支持的操作数类型

Trying to define one of Euler's approximations to pi, getting unsupported operand type(s) for 'list and 'int'

python

pi

approximation