试图定义欧拉对 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

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

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 开始，因为所有小于 i*i 的 i 的倍数都已被删除。
• 一般来说，从您正在迭代的列表中删除元素是非常糟糕的做法。参见例如this post。在内部，Python 在它迭代的列表中使用索引。巧合的是，在这种特定情况下这不是问题，因为只会删除大于当前的数字。
• 此外，从很长的列表中删除元素非常慢，因为每次都需要移动完整的列表以填补空白。因此，最好使用两个单独的列表。
• 您没有计算结果产品，也没有退货。
• 如您所见，此公式收敛速度非常慢。
• 如 cmets 中所述，以前的版本将 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

Answer 2

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

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