Approximer pi avec une suite de Ramanujan

Srinisiva Ramanujan (22 décembre 1887 – 26 avril 1920) est un mathématicien autodidacte indien, auteur génial de nombreuses formules très originales.

Godfroy Hardy – mathématicien britannique qui a fait connaître les travaux de Ramanujan – rapporte cette anecdote :

« Je me souviens que j'allais le voir une fois alors qu'il était malade, à Putney. J'avais pris un taxi portant le numéro 1729 et je remarquais que ce nombre me semblait peu intéressant, ajoutant que j'espérais que ce ne fût pas mauvais signe.
— Non, me répondit-il, c'est un nombre très intéressant : c'est le plus petit nombre décomposable en somme de deux cubes de deux manières différentes. »
En effet, $1729 = 9^3 + 10^3 = 1^3 + 12^3$.

Il a produit des formules permettant de calculer $\pi$. En voici une :

$$\frac{1}{\pi} = \frac{2\sqrt{2}}{9801} \sum_{k=0}^{\infty} \frac{(4k)! (1\,103 + 26\,390k)}{(k!)^4 \cdot 396^{4k}}$$

Il faut bien dire que cette formule fait un peu peur et parait magique.

On ne prend pas spécialement de précautions, on fait les calculs normalement en Python. Puisque les résultats sont des nombres à virgule, on utilisera des float.

from math import sqrt # racine carrée pour le calcul de racine(2)

La formule de Ramanujan es utile parce qu'elle converge vite. On ne peut bien sûr pas faire la somme à l'infini et donc on ne peut pas obtenir une valeur exacte de $\pi$, mais on atteint des approximations excellentes même en ne sommant que jusqu'à un N faible.

• Avec N == 0 on a déjà 6 bons chiffres après la virgule,
• avec N == 1 on en a 14 !

Une telle formule n'a d'intérêt que pour calculer $\pi$ avec une très grande précision. Alors nous allons devoir utiliser autre chose que les simples float. En effet, on sait que les float ont une précision limitée. On va devoir utiliser un autre mode de représentation des nombres.

Pour cela nous allons utiliser le module decimal qui permet de stocker des nombres décimaux avec une précision arbitraire et sans le problème habituel de l'encodage binaire – rappelez vous, que se passe-t-il en Python quand on écrit 0.1 + 0.2 ?

from decimal import Decimal, getcontext

C = getcontext() # permettra de régler des paramètres
C.prec = 1024    # par exemple, une précision à 1024 chiffres

# un exemple d'utilisation
x = Decimal(3)   # x est égal au nombre 3 mais stocké au format Decimal
x.sqrt()         # renvoie la racine carrée de x, donc 3, avec une précision à 1024 chiffres !
3*x              # renvoie le Decimal égal à 3*x, donc 9 ici
0.1*x            # à ne surtout pas faire : si on utilise Decimal, c'est justement pour éviter
                 # le problème des float qui ne codent qu'imparfaitement les nombres décimaux

dixieme = D(0.1) # surtout pas ! comme vous devez le savoir, quand on écrit 0.1,
                 # la machine produit un nombre codé en binaire qui ne fait qu'approximer 0.1
                 # et ce nombre est 0.100000000000000005551...
dixieme = D('0.1') # oui ! on laisse à Decimal le soin de comprendre correctement '0.1'

• Ajouter les lignes 1 à 4 (import et réglage de contexte) de ce qui précède,
• supprimer l'import de sqrt de la bibliothèque math qui ne servira plus,
• vous pouvez laisser votre méthode facto inchangée,
• dans la fonction terme, convertissez 1103 en Decimal(1103) et de même pour 26390 et 396
• Dans le calcul final, au lieu de sqrt(2) remplacez par Decimal(2).sqrt().

Vous pouvez recommencer et apprécier le résultat.