Python decimal.InvalidOperation Ошибка с использованием большого числа

На питоне я написал программу для вычисления значения числа пи с использованием алгоритма Чудновского. Он работает с числами до 800. Однако, если я использую число выше 800, он возвращает ошибку decimal.InvalidOperation. Это мой код:

from math import factorial
from decimal import Decimal, getcontext

getcontext().prec = 1000

pi_input = input("How Many Digits Of Pi Are To Be Represented: ")
num = int(pi_input)

def cal(n):
    t = Decimal(0)
    pi = Decimal(0)
    deno = Decimal(0)

    for k in range(n):
        t = ((-1) ** k) * (factorial(6 * k)) * (13591409 + 545140134 * k)
        deno = factorial(3 * k) * (factorial(k) ** 3) * (640320 ** (3 * k))
        pi += Decimal(t) / Decimal(deno)
        
    pi = pi * Decimal(12) / Decimal(640320 ** Decimal(1.5))
    pi = 1 / pi

    return round(pi, n)

print(cal(num))

Кто-нибудь может мне помочь?

Какая строка вызывает исключение?

President James K. Polk 12.12.2020 22:55

Я не вижу недопустимых операций, и у меня все работает нормально для num=800. Хотя вы могли бы значительно ускорить его.

President James K. Polk 12.12.2020 23:11

Вы можете сравнить истинное значение Pi с from mpmath import mp и mp.dps = num; print(mp.pi). Вы также можете выполнить полный расчет с помощью mpmath, который уже реализует функцию factorial. mpmath работает намного быстрее, чем Decimal

JohanC 12.12.2020 23:12

@PresidentJamesK.Polk Извините, вы правы. math.factorial отлично работает с целыми числами. round(pi, 1000) выдает ошибку, когда getcontext().prec = 1000.

JohanC 12.12.2020 23:28
Почему в Python есть оператор "pass"?
Почему в Python есть оператор "pass"?
Оператор pass в Python - это простая концепция, которую могут быстро освоить даже новички без опыта программирования.
Некоторые методы, о которых вы не знали, что они существуют в Python
Некоторые методы, о которых вы не знали, что они существуют в Python
Python - самый известный и самый простой в изучении язык в наши дни. Имея широкий спектр применения в области машинного обучения, Data Science,...
Основы Python Часть I
Основы Python Часть I
Вы когда-нибудь задумывались, почему в программах на Python вы видите приведенный ниже код?
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
Алиса и Боб имеют неориентированный граф из n узлов и трех типов ребер:
Оптимизация кода с помощью тернарного оператора Python
Оптимизация кода с помощью тернарного оператора Python
И последнее, что мы хотели бы показать вам, прежде чем двигаться дальше, это
Советы по эффективной веб-разработке с помощью Python
Советы по эффективной веб-разработке с помощью Python
Как веб-разработчик, Python может стать мощным инструментом для создания эффективных и масштабируемых веб-приложений.
1
4
231
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Я бы не рекомендовал ваш метод. Попробуйте использовать модуль gmpy следующим образом:

import sys
import math
import os.path
from gmpy2 import mpz, sqrt
from time import time, sleep

print_pi = input("Should The Value Of Pi Be Printed: ").lower()
print_pi_bool = 0

if "y" in print_pi:
    print_pi_bool = True
elif "n" in print_pi:
    print_pi = False
else:
    print("Incorrect Input. Please Try Again.")
    sys.exit()

def sqrt(n, one):
    """
    Return the square root of n as a fixed point number with the one
    passed in.  It uses a second order Newton-Raphson convgence.  This
    doubles the number of significant figures on each iteration.
    """
    # Use floating point arithmetic to make an initial guess
    floating_point_precision = 10**16
    n_float = float((n * floating_point_precision) // one) / floating_point_precision
    x = (int(floating_point_precision * math.sqrt(n_float)) * one) // floating_point_precision
    n_one = n * one
    count = 0
    while 1:
        count += 1
        print(count)
        x_old = x
        x = (x + n_one // x) // 2
        if x == x_old:
            break
    return x

def pi_chudnovsky_bs(digits):
    """
    Compute int(pi * 10**digits)

    This is done using Chudnovsky's series with binary splitting
    """
    C = 640320
    C3_OVER_24 = C**3 // 24
    def bs(a, b):
        """
        Computes the terms for binary splitting the Chudnovsky infinite series

        a(a) = +/- (13591409 + 545140134*a)
        p(a) = (6*a-5)*(2*a-1)*(6*a-1)
        b(a) = 1
        q(a) = a*a*a*C3_OVER_24

        returns P(a,b), Q(a,b) and T(a,b)
        """
        if b - a == 1:
            # Directly compute P(a,a+1), Q(a,a+1) and T(a,a+1)
            if a == 0:
                Pab = Qab = mpz(1)
            else:
                Pab = mpz((6*a-5)*(2*a-1)*(6*a-1))
                Qab = mpz(a*a*a*C3_OVER_24)
            Tab = Pab * (13591409 + 545140134*a) # a(a) * p(a)
            if a & 1:
                Tab = -Tab
        else:
            # Recursively compute P(a,b), Q(a,b) and T(a,b)
            # m is the midpoint of a and b
            m = (a + b) // 2
            # Recursively calculate P(a,m), Q(a,m) and T(a,m)
            Pam, Qam, Tam = bs(a, m)
            # Recursively calculate P(m,b), Q(m,b) and T(m,b)
            Pmb, Qmb, Tmb = bs(m, b)
            # Now combine
            Pab = Pam * Pmb
            Qab = Qam * Qmb
            Tab = Qmb * Tam + Pam * Tmb
        return Pab, Qab, Tab
    # how many terms to compute
    DIGITS_PER_TERM = math.log10(C3_OVER_24/6/2/6)
    N = int(digits/DIGITS_PER_TERM + 1)
    # Calclate P(0,N) and Q(0,N)
    P, Q, T = bs(0, N)
    one_squared = mpz(10)**(2*digits)
    sqrtC = sqrt(10005*one_squared, one_squared)
    return (Q*426880*sqrtC) // T

# The last 5 digits or pi for various numbers of digits
check_digits = {
        100 : 70679,
       1000 :  1989,
      10000 : 75678,
     100000 : 24646,
    1000000 : 58151,
   10000000 : 55897,
}

if __name__ == "__main__":
    digits = 100
    pi = pi_chudnovsky_bs(digits)
    #raise SystemExit
    for log10_digits in range(1,11):
        digits = 10**log10_digits
        start =time()
        pi = pi_chudnovsky_bs(digits)
        pi = str(pi)
        pi = pi[:(len(str(pi)) // 2) + 1]
        length = int(len(str(pi))) - 1
        print("Chudnovsky Binary Splitting Using GMPY: Digits",f"{digits:,}","\n--------------",time()-start,"--------------")
        print("Length Of " + f"{length:,}" + " Digits")

Это мой первый ответ, поэтому я надеюсь, что он поможет!

Ааа спасибо! Теперь это работает. Думаю, я выброшу другой

The Pilot Dude 12.01.2021 16:13

Другие вопросы по теме