如何加速Leibniz级数收敛以高效计算π？

背景 #

Leibniz级数是用于近似π的无穷级数，虽能收敛到π，但收敛速度极慢。我曾实现多种π计算函数，其中Ramanujan公式收敛效率最高，仅3次迭代就能得到15位正确小数。

之后尝试通过添加修正项加速Leibniz级数：

• 1个修正项：2048次迭代时得到10位正确小数
• 2个修正项：638次迭代时得到15位正确小数
• 3个修正项：仅需131次迭代即可得到15位正确小数，但4个修正项会导致收敛变慢，154次迭代时出现下溢

我还实现了Euler变换、Aitken加速等收敛加速算法，但这些算法虽减少了所需项数，单步计算开销却大幅增加，导致总执行时间更长。

问题 #

寻求基于Leibniz级数维基页面的收敛加速方法及编程技巧，实现Leibniz级数计算π的加速，同时减少所需级数项数并降低执行时间。允许使用编译库，但禁止缓存等作弊手段。

相关代码及性能测试 #

import math
from itertools import cycle
from math import factorial


def Ramanujan(n):
    def term(k): return factorial(4 * k) * (1103 + 26390 * k) / \
        (factorial(k) ** 4 * 396 ** (4 * k))
    s = sum(term(i) for i in range(n))
    return 1 / (2 * 2 ** 0.5 * s / 9801)


sign = [1, -1]


def Leibniz_Series(n):
    return [1/i*s for i, s in zip(range(1, 2*n+1, 2), cycle(sign))]


def Leibniz_0(n):
    return 4 * sum(Leibniz_Series(n))


def Leibniz_1(n):
    correction = 1 / n*sign[n % 2]
    return Leibniz_0(n) + correction


def Leibniz_2(n):
    correction = 1 / (
        n + 1 / (
            4 * n
        )
    )*sign[n % 2]
    return Leibniz_0(n) + correction


def Leibniz_3(n):
    correction = 1 / (
        n + 1 / (
            4 * n + 4 / n
        )
    )*sign[n % 2]
    return Leibniz_0(n) + correction


def Leibniz_4(n):
    correction = 1 / (
        n + 1 / (
            4 * n + 4 / (
                n + 9 / n
            )
        )
    )*sign[n % 2]
    return Leibniz_0(n) + correction


def bincoeff(n, k):
    r = 1
    if (k > n):
        return 0
    for d in range(1, k + 1):
        r = r * n / d
        n -= 1
    return r


def abs(n):
    return n * (-1) ** (n < 0)


def Euler(arr):
    out = []

    for i in range(len(arr)):
        delta = 0
        for j in range(i + 1):
            coeff = bincoeff(i, j)
            delta += (-1) ** j * coeff * abs(arr[j])

        out.append(0.5 ** (i + 1) * delta)
    return out


def Leibniz_Euler(n):
    return 4 * sum(Euler(Leibniz_Series(n)))


def Leibniz_Aitken(n):
    series = Leibniz_Series(n)
    s0, s1, s2 = (sum(series[:n-i]) for i in (2, 1, 0))
    return 4 * (s2 * s0 - s1 * s1) / (s2 - 2*s1 + s0)


def LCP(s1, s2):
    i = 0
    for a, b in zip(s1, s2):
        if a != b:
            break
        i += 1
    return i


for func in (Ramanujan, Leibniz_Euler, Leibniz_Aitken, Leibniz_4, Leibniz_3, Leibniz_2, Leibniz_1, Leibniz_0):
    i = 12
    last = 0
    while (pi := func(i)) != math.pi:
        i += 1
        if i == 2048:
            if abs(pi - math.pi) <= 1e-6:
                print('计算结果接近π')
            else:
                print('方法收敛过慢')
            break
        if pi == last:
            print('浮点下溢')
            break
        last = pi
    print(
        f'函数: {func.__name__}, 迭代次数: {i}, 正确位数: {LCP(str(pi), str(math.pi))}')

性能测试输出 #

函数: Ramanujan, 迭代次数: 12, 正确位数: 17
浮点下溢
函数: Leibniz_Euler, 迭代次数: 52, 正确位数: 16
计算结果接近π
函数: Leibniz_Aitken, 迭代次数: 2048, 正确位数: 11
函数: Leibniz_4, 迭代次数: 154, 正确位数: 17
函数: Leibniz_3, 迭代次数: 131, 正确位数: 17
函数: Leibniz_2, 迭代次数: 638, 正确位数: 17
计算结果接近π
函数: Leibniz_1, 迭代次数: 2048, 正确位数: 12
方法收敛过慢
函数: Leibniz_0, 迭代次数: 2048, 正确位数: 4

执行时间测试 #

In [2]: %timeit Ramanujan(3)
每次循环3.7微秒 ± 229纳秒（7次运行，每次100000循环的均值±标准差）

In [3]: %timeit Leibniz_3(131)
每次循环17.4微秒 ± 123纳秒（7次运行，每次100000循环的均值±标准差）

In [4]: %timeit Leibniz_Aitken(2048)
每次循环295微秒 ± 10.8微秒（7次运行，每次1000循环的均值±标准差）

In [5]: %timeit Leibniz_Euler(52)
每次循环3.81毫秒 ± 430微秒（7次运行，每次100循环的均值±标准差）

内容的提问来源于stack exchange，提问作者Ξένη Γήινος