Algorithms for computing digits of pi

The algorithm used in that question is not one of those that Wikipedia shows in its entry on approximations of pi.

Mathematically this algorithm can be expressed by the following recurrence relation:

which can be more easily visualized if we use the "expanded" version for the cases n=2 and n=3 for example:

The only reference I've found to this approximation of pi comes from an answer on Math.stackexchange , which itself was simply a copy of a Facebook post, which included the following image:

Therefore it looks like a geometric approximation, designed to be easy to visualize and not precisely efficient. There are much more efficient approaches of which I will highlight just a couple (because the topic of computing pi efficiently is certainly a world).

Gauss–Legendre approximation

This approximation converges much faster than the one used in this question. In the one used in this question, the number of correct decimals only increases by 1 (or sometimes 0) with each new iteration. On the other hand, in the Gauss-Legendre approximation, the number of correct digits doubles in each iteration, so that after a few iterations there is already a high number of decimal places.

The mathematical details are complex (they can be found on Wikipedia ), and here I will give only a Python implementation (adapted from this answer on SO ):

from decimal import Decimal, getcontext

num_iter = int(input("iteraciones? "))
getcontext().prec = 65
a = 1
b = 1/Decimal(2).sqrt()
t = 1/Decimal(4)
p = 1                
i = 0
while i < num_iter:
  an = (a + b) / 2
  b  = (a * b).sqrt()
  t -= p * (a - an) * (a - an)
  a = an
  pi = (a + b) * (a + b) / (4 * t)
  p  = 2*p
  i = i + 1

# Aproximacion
print(pi)  
# Valor exacto (hasta decimal 65, obtenido de WolframAlpha)
print("3.1415926535897932384626433832795028841971693993751058209749445923")

With only three iterations, many correct decimals are already obtained!

iteraciones? 3
3.141592653589793238 2795127748018639743812255048354469357873307022
3.141592653589793238 4626433832795028841971693993751058209749445923

And with one more iteration, the number of correct decimal places doubles.

iteraciones? 4
3.1415926535897932384626433832795028841971 146782836489215566171072
3.1415926535897932384626433832795028841971 693993751058209749445923

already with five or more we have the 64 correct digits (the last digit is never correct due to the rounding error caused by working only with 65 decimals)

In fact, a small variation of the above program can allow us to obtain any desired number of decimal places. It is enough to work with a precision of a couple more decimal places (for example, if we want 500 digits of pi, we work with 502) to eliminate this problem in the final rounding. And instead of iterating a preset number of times, we iterate until two successive iterations produce the same result, then stopping:

from decimal import Decimal, getcontext

n_digitos = int(input("Número de dígitos deseado? "))
getcontext().prec = n_digitos+2
a = 1
b = 1/Decimal(2).sqrt()
t = 1/Decimal(4)
p = 1                
i = 0
pi_anterior = None
while True:
  an = (a + b) / 2
  b  = (a * b).sqrt()
  t -= p * (a - an) * (a - an)
  a = an
  pi = (a + b) * (a + b) / (4 * t)
  if pi == pi_anterior:
    break
  p  = 2*p
  pi_anterior = pi
  i = i + 1

print("Obtenido resultado en {} iteraciones".format(i))
# Aproximacion
print(str(pi)[:-2])  # Quitar los dos últimos, que serían incorrectos  
# Valor exacto (hasta decimal 65, obtenido de WolframAlpha)
print("3.1415926535897932384626433832795028841971693993751058209749445923")

Número de dígitos deseado? 65
Obtenido resultado en 5 iteraciones
3.1415926535897932384626433832795028841971693993751058209749445923
3.1415926535897932384626433832795028841971693993751058209749445923

Let's go big!

Número de dígitos deseado? 1000
Obtenido resultado en 9 iteraciones
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899
8628034825342117067982148086513282306647093844609550582231725359408128481117450284
1027019385211055596446229489549303819644288109756659334461284756482337867831652712
0190914564856692346034861045432664821339360726024914127372458700660631558817488152
0920962829254091715364367892590360011330530548820466521384146951941511609433057270
3657595919530921861173819326117931051185480744623799627495673518857527248912279381
8301194912983367336244065664308602139494639522473719070217986094370277053921717629
3176752384674818467669405132000568127145263560827785771342757789609173637178721468
4409012249534301465495853710507922796892589235420199561121290219608640344181598136
2977477130996051870721134999999837297804995105973173281609631859502445945534690830
2642522308253344685035261931188171010003137838752886587533208381420617177669147303
5982534904287554687311595628638823537875937519577818577805321712268066130019278766
11195909216420198
3.1415926535897932384626433832795028841971693993751058209749445923

Algorithms based on getting figures

If the previous algorithm is based on approximations (each iteration produces a number closer to pi), there are other algorithms that are based on obtaining numbers of pi. The main advantage of this type of algorithm would be that they do not require as much storage space, because once a digit has been computed, it is printed and it is not necessary to keep all the previous digits in memory when computing the next one. This is an advantage over methods such as Gauss-Legendre seen before, in which if we want to obtain a precision of one million digits, we must work with a million digits from the beginning.

One very notable algorithm is the one based on the Bailey-Borwein-Plouffe formula , since it allows computing the digit that appears in position N without the need to compute the previous N-1. Unfortunately that algorithm produces the digits in hexadecimal base, and not in the base 10 that we are used to. Knowing the Nth digit in hexadecimal does not help us to know a specific digit of the decimal expansion. What is possible is to write pi with N digits (hexadecimal), which would be an approximation of pi, and convert that number to decimal, to obtain another approximation in base 10, which would be correct with a precision of 2 ^-4*N

However, this saving the first N hexadecimal digits to pass them to decimal spoils that advantage of "not having to compute or store the previous digits".

Other methods of obtaining pi "digit by digit" are known globally as Spigot algorithms . One of them (which I've adapted from this answer on SO ), claims to be quite fast, although it's not at all obvious to understand what it does (not even reading the article it's based on), so it could even be wrong, although it seems to give correctly the first 65 decimals:

def pi_digits():
  k,a,b,a1,b1 = 2,4,1,12,4
  while True:
    p,q,k = k * k, 2 * k + 1, k + 1
    a,b,a1,b1 = a1, b1, p*a + q*a1, p*b + q*b1
    d,d1 = a/b, a1/b1
    while d == d1:
        yield int(d)
        a,a1 = 10*(a % b), 10*(a1 % b1)
        d,d1 = a//b, a1//b1

num_decimales = 65
for i, digit in enumerate(pi_digits()):
  print(digit, end="")
  if i==0: print(".", end="")  # El punto decimal
  if i==num_decimales-1: break

print()
# Valor exacto para comparación
print("3.1415926535897932384626433832795028841971693993751058209749445923")

Execution result:

3.1415926535897932384626433832795028841971693993751058209749445923
3.1415926535897932384626433832795028841971693993751058209749445923

This method does not converge as fast as the one seen before, since to obtain 65 decimal places it has iterated 65 times. The function pi_digits()is an infinite generator, each time it is called it returns the next digit of pi. The one I foruse it from has a breakto stop after 65 calls (otherwise it would be infinite).

The advantage of this method is that it does not depend on precision nor does it require decimal, since it works with integers, and that it is more efficient in memory since it does not need to retain all the digits of pi that it has computed so far, but only a few variables ( p,q,k,a,b,a1,1,d,d1) . . In return it is slower from a certain decimal.

If we put it to compute 1000 digits we see that it obtains exactly the same result as the Gauss-Legendre algorithm, and it takes twice as long. If we try to get 100,000 digits of pi, Gauss-Legendre takes 21s, while the Spigot algorithm takes on the order of 3 minutes. Each additional digit takes a little longer.