Talk:Chudnovsky algorithm

Plagiarism

The materials which I saw in this footnote are first wrong and also stolen from another web site. Sunos 6 (talk | contribs) 05:15, 9 April 2008 (UTC)

Usage

How can the reader use this algorithm - from what point he can certainly know the n-th digit after the decimal dot is correct? 79.179.42.44 (talk) 21:16, 16 February 2012 (UTC)

The error will be approximately equal to the next term, so by estimating very roughly the size of the next term, you know up to where the approximation is correct. The factor (6k)!/(3k)!k!^3 grows by a factor 693, 982, 1147, 1252 for the first 5 terms, and ~ 1500 for the next 20 terms. This is to be divided by 262537412640768000, which yields a ratio of ~ 1.5e14 between subsequent terms. If this is not wrong, it should yield roughly 14 more digits at each step. — MFH:Talk 17:50, 14 March 2018 (UTC)

multiple of e ?

The link between 640320 and e^(pi sqrt 163) is given, but is there a simple explanation for 13591409 x 2e-7 = 2.7182818 ~ e? — MFH:Talk 17:32, 14 March 2018 (UTC)

Integer division in the Python code

In the Python section, on the line

   M = (K**3 - 16*K) * M // k**3

Why is this an integer division given that the whole calculation is in arbitrary-precision floats? — Preceding unsigned comment added by Eje211 (talk • contribs) 01:21, 27 May 2018 (UTC)

possible reason

I tried changing '//' to '/' and got an error because '/' is not supported by Decimal objects. The good news is that I checked the output of the longest of the runs, and compared them with the output of a completely different algorithm, and they were identical. The code (based on https://gist.github.com/markhamilton1/9716714) is much simpler, but also pretty fast (they both run in about a second on my hardware).

   q, r, t, k, n, l = 1, 0, 1, 1, 3, 3                                                                                               
   while True:
       if 4*q+r-t < n*t:
           yield n
           nr = 10*(r-n*t)
           n  = ((10*(3*q+r))//t)-10*n
           q  *= 10
           r  = nr
       else:
           nr = (2*q+r)*l
           nn = (q*(7*k)+2+(r*l))//(t*l)
           q  *= k
           t  *= l
           l  += 2
           k += 1
           n  = nn
           r  = nr

kogorman (talk) 00:26, 10 October 2018 (UTC)

How to choose parameters

It's not clear to me what parameters to choose if I want to run this beyond the most precise example. What's the best relationship among maxk, prec and disp?

kogorman (talk) 00:34, 10 October 2018 (UTC)

Incorrect Algorithm

Please correct me if I am wrong, but I believe the formulae for M_q+1, K₀, and K_q+1 are incorrect as written:

${\begin{alignedat}{4}M_{q+1}&=M_{q}\cdot \left({\frac {(12q+2)(12q+6)(12q+10)}{(q+1)^{3}}}\right)\,\,&&{\textrm {where}}\,\,M_{0}&&=1\\[4pt]\end{alignedat}}$

The computation of M_q can be further optimized by introducing an additional term K_q as follows:

${\begin{alignedat}{4}K_{q+1}&=K_{q}+12\,\,&&{\textrm {where}}\,\,K_{0}&&=6\\[4pt]M_{q+1}&=M_{q}\cdot \left({\frac {K_{q}^{3}-16K_{q}}{(q+1)^{3}}}\right)\,\,&&{\textrm {where}}\,\,M_{0}&&=1\\[12pt]\end{alignedat}}$

---

I believe the formulae should instead be:

${\begin{alignedat}{4}M_{q+1}&=M_{q}\cdot \left({\frac {(12q-2)(12q-6)(12q-10)}{q^{3}}}\right)\,\,&&{\textrm {where}}\,\,M_{0}&&=1\\[4pt]\end{alignedat}}$

The computation of M_q can be further optimized by introducing an additional term K_q as follows:

${\begin{alignedat}{4}K_{q+1}&=K_{q}+12\,\,&&{\textrm {where}}\,\,K_{0}&&=-6\\[4pt]M_{q+1}&=M_{q}\cdot \left({\frac {K_{q}^{3}-16K_{q}}{q^{3}}}\right)\,\,&&{\textrm {where}}\,\,M_{0}&&=1\\[12pt]\end{alignedat}}$

By making the above changes, I am able to correctly calculate π in my python program.

ScowlingFlavoredJava (talk) 18:43, 4 June 2021 (UTC)

I think the original formula is correct. How could you even calculate M1 from M0 with your formula? It means that you need to divide by 0. 31.208.187.87 (talk) 20:57, 15 June 2022 (UTC)

No. For M1, the value of q=1, so you're dividing by 1 (from 1^3).

The original formula is incorrect. The "q + 1" term is incorrect and should be changed to "q". I discovered this issue independently when trying to implement this algorithm in a bc script.

By making the above changes, I am able to correctly calculate π in my bc script. 24.217.141.228 (talk) 22:49, 20 December 2022 (UTC)

Are you sure? If I wanted to evaluate M₁ using ScowlingFlavoredJava's formula, I would need to begin by replacing each occurrence of q with 0, in order to obtain M₁ on the left hand side and a function of M₀ on the right hand side. This results in a division by zero.

If instead I substitute q=1 as you suggest, I obtain M₂ = M₁ . (...), which isn't useful because the value of M₂ isn't known yet.

ScowlingFlavoredJava's formula would however be completely valid with a slight change of indices:

M_q = M_q-1 . (12q-2)(12q-6)(12q-10)/q³, defined for q >= 1

This is exactly what the Python code below is doing - q is incremented before the next value of M is calculated. E.g. on the first iteration of the loop, q gets assigned the value 1 and then M gets assigned the calculated value M₁.

This way of expressing recurrence is mathmatically sound, but it's just syntactically inconsistent with the L_q+1 and X_q+1 recurrences. You can substitute q -> q+1 to fix this, which after some simplification arrives you back at the (12q+2)(12q+6)(12q+10)/(q+1)³ formula.

I believe the original equation with denominator (q+1)³ is correct, the Python code is also correct, and ScowlingFlavoredJava and 24.217.141.228 have both made the same off-by-one error in their reasoning about the mathematical expression. kierdavis (talk) 20:16, 10 February 2023 (UTC)

Proposed Python code

I have implemented the above changes in a python module should you wish to try it out. Please see below:

#!/usr/bin/env python3
import decimal
from decimal import Decimal as D  # For internal use only


def pi(precision: int = decimal.getcontext().prec,
       rounding: decimal.Context.rounding = decimal.ROUND_FLOOR
       ) -> decimal.Decimal:
    """Return pi as a decimal.Decimal to the specified precision.

    This function uses the Chudnovsky algorithm
    (https://wikipedia.org/wiki/Chudnovsky_algorithm) and implements
    acceleration by using the Mq (and Kq), Lq, and Xq recurrence
    patterns.

    Note that the first argument is the precision NOT the number of
    decimal places. The relationship is:
        decimal places = precision + 1

    Positional/keyword arguments:
        precision -- the precision to which pi is calculated
                (default decimal.getcontext().prec)
        rounding -- how to round the final value of pi
                (default decimal.ROUND_FLOOR)
    """
    # Add 2 to the precision to ensure the accuracy of the answer
    precision += 2
    decimal.getcontext().prec = precision

    # Use incremental calculations.
    # Set initial values for q, C, L, X, K, and M
    q = D(0)
    C = D(426880)*D(10005).sqrt()
    L = D(13591409)
    X = D(1)
    K = D(-6)
    M = D(1)
    running_sum = D(0)
    pi = D('NaN')

    # Iteration 0
    lastpi = pi
    running_sum += M*L/X
    pi = C * running_sum**D(-1)

    # Loop until pi and lastpi converge to the current precision
    while pi - lastpi:
        lastpi = pi
        q += D(1)
        L += D(545140134)
        X *= D(-262537412640768000)
        K += D(12)
        M *= (K**D(3) - D(16)*K) / q**D(3)
        running_sum += M*L/X
        pi = C * running_sum**D(-1)

    # Reduce precision back to the specified precision
    precision -= 2
    decimal.getcontext().prec = precision

    # quantize pi to sepcified precision with specified rounding
    return pi.quantize(D('1e-'+str(precision-1)), rounding)


if __name__ == '__main__':
    print(pi(50))
    # Expected output:
    #
    # >>> print(pi(50))
    # 3.1415926535897932384626433832795028841971693993751

ScowlingFlavoredJava (talk) 20:00, 4 June 2021 (UTC)

Thanks. Works. The one in the main page is not useful. 2001:559:B:BE96:6188:BBA:F027:9892 (talk) 16:51, 17 November 2023 (UTC)

Circular references

The source for the complexity of the algorithm directly cites this article as a reference. --Tbjep (talk) 17:26, 12 December 2021 (UTC)

vandalism protection

due to the title of the algorithm containing the word "chud" in it, people have decided to vandalize this article more frequently than usual. please put relevant protections in place to prevent this page from being vandalized in the future. 66.103.200.96 (talk) 22:32, 12 November 2023 (UTC)

There have been maybe two incidents of vandalism since the start of 2020. That is not frequent vandalism. —David Eppstein (talk) 23:11, 12 November 2023 (UTC)