Comment calculer le 10’000’000’000’000’000’000 ème terme de la suite de Fibonacci

Tombé l’autre jour sur un problème idiot [1] mais intéressant : calculer le 10¹⁹ ème terme de la suite de Fibonacci. Idiot parce que ça ne sert à rien. Intéressant parce que ça sous-entend qu’il existe une manière de calculer le n-ième terme de cette suite définie par récurrence sans calculer tous les termes précédents. En effet, calculer les termes les uns après les autres prendrait dans les 300’000 ans à raison d’une microseconde par terme…

Un nombre d’or, mais flottant

Comme même les ésotéristes le savent, la suite de Fibonacci est liée au nombre d’or. A partir de ce fait, Moivre, Euler et Binet ont indépendamment obtenu ce qu’on appelle aujourd’hui la formule de Binet, et qui donne directement le n-ième terme de la suite:

\(\mathcal F_n=\frac1{\sqrt5}(\varphi^n-\varphi’^n)\), avec \( \varphi=\frac{1+\sqrt5}2\), et \( \varphi’=-\frac1\varphi\) .

En pratique, le terme en \( \varphi’\) devient rapidement négligeable et il suffit de chercher l’entier le plus proche de \(\frac{\varphi^n}{\sqrt5}\). Par exemple pour n=50, on peut calculer très vite \(\mathcal F_{50}\approx\frac{\varphi^{50}}{\sqrt5}\) = 12586269025

Mais on se heurte rapidement au problème de la précision de calcul en nombres flottants sur nos ordinateurs : \(\varphi\) étant irrationnel (mais pas aussi transcendant que \(\pi\) …), il faut calculer \(\varphi^n\) avec au moins autant de décimales que \(\mathcal F_n\) comporte de chiffres. Dès n=71 on dépasse les 53 bits de la mantisse des nombres « extended precision » et les dernières décimales des termes calculés sont faux.

Dit autrement, à partir de n=70 la division \(\mathcal F_n/\mathcal F_{n-1}\) donne au moins autant de décimales de \( \varphi\) que son calcul en nombres flottants

La Matrice salvatrice

Pour une méthode de calcul en nombres entiers atteignant (presque) la rapidité de la formule de Binet, il faut utiliser l’idée d’un certain Brenner : écrire la récurrence de Fibonacci sous forme matricielle [2]. La matrice toute simple \(Q=\begin{pmatrix}1&1\\1&0\end{pmatrix}\) permet d’obtenir un nouveau terme de la série de Fibonacci en la multipliant par un vecteur formé des deux termes précédents:

\(\begin{pmatrix}1&1\\1&0\end{pmatrix}\begin{pmatrix}\mathcal F_{n-1}\\\mathcal F_{n-2}\end{pmatrix}=\begin{pmatrix}\mathcal F_{n}\\\mathcal F_{n-1}\end{pmatrix}\)

En enchaînant les multiplications matricielles, on obtient le n-ième terme à partir des deux premiers (0,1) ainsi :

\(\begin{pmatrix}1&1\\1&0\end{pmatrix}^{n-1}\begin{pmatrix}1\\0\end{pmatrix}=\begin{pmatrix}\mathcal F_{n}\\\mathcal F_{n-1}\end{pmatrix}\)

En fait on retrouve les termes de la suite directement dans la matrice  \(Q^n = \begin{pmatrix}\mathcal F_{n+1}&\mathcal F_{n}\\\mathcal F_{n}&\mathcal F_{n-1}\end{pmatrix}\).

L’algorithme de l'exponentiation rapide permet d’élever la matrice Q à la puissance n en effectuant log₂(n) multiplications de matrices 2×2, soit environ 63 pour n=10¹⁹. Ultra rapide, et facilement généralisable à d’autres formules de récurrence !

Malheureusement, la fonction matrix_power de la librairie Python numpy souffre d’une limitation, voire d’un bug qui empêche de l’utiliser dès n=71 car elle utilise les nombres flottants. Il m’a donc fallu la réécrire, en résolvant une petite contrainte technologique au passage.

Vers les pétaoctets et au delà …

Il se trouve que le 10¹⁹-ième terme de la suite de Fibonacci est trop grand pour tenir dans le plus super des ordinateurs. En effet, la longueur du n-ième terme de la suite s’obtient facilement en prenant le log (base 10)  de la formule de Binet : \(\mathcal L_n=\lceil n\log{\varphi}\rceil\) , soit grosso-modo n*0.208987640249978733769272089237555416822459239918210953539…

Vérifions. Pour n=1000, on obtient:

\(\mathcal F_{1000}\)=43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875 qui comporte bien pile 209 chiffres.

Donc \(\mathcal F_{10^{19}}\) comporte 2089876402499787337 chiffres … Il faudrait dans les 867 pétaoctets de RAM (de préférence…) pour stocker ce nombre …

Et modulo 1000000007 ?

Fort heureusement, le problème idiot avait un petit détail en prime : il fallait calculer le résultat modulo 1000000007, et ça, ça change tout.

Ce nombre est premier (pour que ça soit plus intéressant), mais surtout il est inférieur à 2³⁰, ce qui fait que la multiplication modulo 1000000007 ne nécessite que des opérations sur des entiers 32 bits (signés), ultra rapide.

Voici donc quelques fonctions Python tirées du module math2 de ma librairie Goulib qui effectuent:

• le produit matriciel, optionnellement modulaire
• l’exponentiation rapide d’une matrice, optionnellement modulaire, et qui corrige le bug de numpy.matrix_power
• le calcul quasi instantané du n-ième terme de la série de Fibonacci, avec une option modulaire vivement recommandée pour de grands n

Avec ça on obtient fibonacci(int(1E19),1000000007) = 647754067 ce qui nous fait une belle jambe,  mais nous a appris pas mal de choses intéressantes, non ?

Prochaine étape : trouver la période de Pisano correspondante …

Références

1. Nth Fibonacci number for n as big as 10^19? sur StackOverflow (et réponse de will qui m’a inspiré cet article)
2. Weisstein, Eric W. « Fibonacci Q-Matrix. » From MathWorld–A Wolfram Web Resource.
3. Arthur Charpentier « Fibonacci, les lapins, le nombre d’or et les calculs actuariels », 2016 sur Freakonometrics (mention spéciale pour l’image de la spirale…)