Как делать быстрее вычисление числа пи по Лейбницу

Рейтинг: 3Ответов: 3Опубликовано: 14.02.2023
def pi(*imp):
    if not imp:
        imp = 4
    else:
        imp = imp[0]
    rel_imp = '1'
    for n in range(imp + 1):
        rel_imp = rel_imp + '0'
    imp = int(rel_imp)
    del rel_imp
    return_pi = 0
    int_set_pi = 1
    for n in range(imp):
        return_pi += 1 / int_set_pi
        if int_set_pi > 0:
            int_set_pi = (int_set_pi + 2) * -1
        else:
            int_set_pi = (int_set_pi - 2) * -1
    return round(4 * return_pi, len(str(imp)) - 2)

Нужно как-то ускорить. Оно до 6-ого точного знака считает секунд 10, до 7-ого знака уже 120 секунд. Пытался добавить параллельное вычисление по кусочкам (asyncio), рационализм в памяти (tracemallok), сохранение уже подсчитанного (open). Но оно не помогает.

Ответы

▲ 0Принят

Странноватый у вас код, я некоторые вещи поменял, чтобы Numba смогла заработать, ну и лишнее выкинул. Вот вам до 8 знака за пару секунд:

from numba import njit

@njit
def pi(imp):
    if not imp:
        imp = 4
    imp = 10**(imp+1)
    return_pi = 0
    int_set_pi = 1
    for n in range(imp):
        return_pi += 1 / int_set_pi
        if int_set_pi > 0:
            int_set_pi = (int_set_pi + 2) * -1
        else:
            int_set_pi = (int_set_pi - 2) * -1
    return round(4 * return_pi, len(str(imp)) - 2)

pi(8)

Вывод с таймингом:

3.14159265
Wall time: 1.79 s

Numba очень хорошо ускоряет простые математические вычисления в длинных циклах.

▲ 2

Можно использовать формулу Бэйли-Борвина-Плаффа, о которой можно почитать здесь

Реализация:

import decimal

def pi_bbp(precision):
    decimal.getcontext().prec = precision + 1
    pi = decimal.Decimal(0)
    for k in range(precision):
        pi += (decimal.Decimal(1)/(16**k))*((decimal.Decimal(4)/(8*k+1))-(decimal.Decimal(2)/(8*k+4))-(decimal.Decimal(1)/(8*k+5))-(decimal.Decimal(1)/(8*k+6)))
    return pi

Запуск вычисления 1000 знаков после запятой:

print(pi_bbp(1000))

Время выполнения для 1000 знаков: 0.04690, 10000 знаков: 28.88430

▲ 1

Ваш код, по сути, вычисляет арктангенс единицы. Это можно значительно ускорить, используя тот факт что arctg(1/2) + arctg(1/3) = arctg(1)

Разложение в ряд для арктангенса выглядит так:

его простая реализация на python

def arctg(x, eps=0.0):
    res = 0
    k = 0

    while True:
        new = res + (-1)**k * x**(2*k+1)/(2*k+1)
        if abs(res - new) <= eps:
            return res
        res = new
        k += 1

вычисления pi:

print(4*(arctg(1/2) + arctg(1/3)))