Les nouvelles technologies pour l’enseignement des mathématiques
Intégration des TICE dans l’enseignement des mathématiques

MathémaTICE, première revue en ligne destinée à promouvoir les TICE à travers l’enseignement des mathématiques.

Des décimales en cascade

Dans cet article, Julien Baldacci présente un algorithme simple et efficace pour obtenir un flot continu de décimales pour $\sqrt{2}$.
En annexe, Angelo Laplace a élargi la portée de l’algorithme et l’a appliqué au cas plus complexe de $\sqrt{14}$

Article mis en ligne le 2 janvier 2025
dernière modification le 13 janvier 2025

par Julien Baldacci, Angelo Laplace

N.D.L.R : L’article que voici a une histoire singulière. Proposé en novembre 2021 par Julien Baldacci (professeur au lycée du BTP Saint Lambert, Belleville, Paris), sa rédaction a été brutalement interrompue à la suite d’un problème personnel, alors même que l’article était en voie d’achèvement. Il a semblé regrettable au comité de rédaction de MathémaTICE qu’un article aussi prometteur, présentant un spectaculaire et puissant algorithme peu connu, reste à l’état de chantier dans nos archives.
Angelo Laplace, membre du comité de rédaction, a repris le fil de l’article, il en a comblé les parties inachevées et élargi le propos dans une annexe dont l’intérêt n’échappera à personne.
Il nous a paru naturel d’attribuer l’article à Julien Baldacci (qui en est l’initiateur et qui a rédigé l’essentiel du corps de l’article) et à Angelo Laplace (sans lequel cet article serait resté ignoré en dehors du comité de rédaction) qui y a ajouté une annexe consistante.


Des cascades (sans décimales...) dans le parc national des lacs de Plitvice (Croatie)


Par définition, un nombre irrationnel est un nombre qui ne peut pas s’écrire sous forme de fraction. Son développement décimal est donc infini non périodique, tout comme d’ailleurs est infini son développement en fraction continue. Ainsi la suite des décimales d’un nombre irrationnel est impossible à connaître dans sa totalité. La curiosité du mathématicien ne pourra donc être satisfaite que via un procédé théorique — un algorithme — pour déterminer un nombre arbitrairement grand de décimales. C’est le maximum qu’il puisse espérer.
Mais comme le souligne Bernard Ycart, le problème des algorithmes d’approximation a été moteur dans l’histoire des mathématiques. Une grande partie de l’analyse moderne en est issue. Plus précisément, la notion de fraction continue a bien été engendrée par l’approximation de racines carrées. Nous en verrons les liens dans notre propos et Bernard Ycart vous le conte à cette adresse.

Pour une constante mathématique usuelle donnée, plusieurs algorithmes permettent généralement d’atteindre l’objectif de calculer un grand nombre de décimales, plus ou moins rapidement. On se propose ici d’approcher $\sqrt{2}$ et de comparer deux méthodes.

La première illustre bien les problèmes que pose la détermination effective des décimales d’un irrationnel (ce n’est pas si simple). De plus, elle permet par comparaison de sentir toute l’originalité de la seconde méthode, appelée streaming, car les décimales coulent en flot continu comme un ruisseau. C’est la volonté de présenter cette méthode de streaming qui a motivé ce texte. On se contentera de ce seul exemple détaillé en profondeur : un algorithme simple et efficace pour obtenir un flot continu de décimales pour $\sqrt{2}$. Mais nous ouvrirons de larges perspectives en complément.


La cible : le nombre d’argent

Considérons la fonction $f$ définie par $f(x) = \dfrac{1}{2+x}$ sur $\mathbb{R}$ privé de $-2$.
L’équation $f(x) = x$ admet une unique solution positive $\alpha=\sqrt{2}-1\in[0,1]$.

C’est le nombre $\alpha$ qu’on va approcher, plutôt que $\sqrt{2}$ qui n’en diffère que du chiffre des unités.

Par définition, le nombre $\alpha$ vérifie la suite d’égalités

$$\alpha=\frac{1}{2+\alpha}=\frac{1}{2+\frac{1}{2+\alpha}}=\frac{1}{2+\frac{1}{2+\frac{1}{2+\alpha}}}=\cdots$$

Cette écriture de $\alpha$ en fraction continue ouvre la voie à la mise au point d’algorithmes pour l’approcher.

Notons au passage que cette écriture permet aussi de comprendre pourquoi on appelle parfois $\alpha$ le nombre d’argent, par analogie avec le nombre d’or $\varphi=\frac{1}{2}(1+\sqrt{5})$ dont l’écriture en fraction continue est encore plus simple :

$$\varphi=1+ \frac{1}{1+\frac{1}{1+\frac{1}{1+ \frac{1}{1+\cdots}}}}$$


Premier algorithme

Commençons par un premier algorithme, celui qui vient naturellement à l’esprit lorsqu’on contemple la suite d’égalités

$$\alpha=f(\alpha)=f(f(\alpha))=f(f(f(\alpha)))=\cdots$$

On peut se convaincre facilement que l’intervalle $[0 ;1]$ est stable par $f$. Considérons une suite $(u_n)$ définie par :

$$\begin{cases}u_{n+1} = f(u_n)\\u_0 \in [0 ;1]\end{cases}$$


Cette suite est à valeurs dans $[0 ;1]$.
De plus, pour tous $a,b \in [0 ;1]$, on a :

$$f(b)-f(a) = \frac{1}{2+b} - \frac{1}{2+a} = \frac{(2+a)-(2+b)}{(2+b)(2+a)} = \frac{a-b}{(2+b)(2+a)} $$

.
Et ainsi :

$$|f(b)-f(a)|\le\frac{|b-a|}{(2+a)(2+b)}\le\frac{|b-a|}{4}$$


Par suite, pour tout $n\in\mathbb{N}$, en prenant $b = u_{n}$ et $a = \alpha$ : $|u_{n+1}-\alpha|\le \frac{1}{4}|u_n-\alpha|$.
En utilisant cette inégalité en cascade, il vient facilement :

$$|u_n-\alpha|\le \frac{1}{4^n} \quad\underset{n\to+\infty}{\to}\quad 0.$$


Ainsi, pour tout terme initial $u_0 \in [0 ;1]$, la suite $(u_n)$ converge vers $\alpha$.

Une première façon d’approcher $\alpha$ consiste donc à calculer des termes de $(u_n)$. L’avant-dernière majoration montre qu’on gagne au moins 2 bits par étape, soit 3 décimales toutes les 5 étapes. La vitesse de convergence est tout à fait acceptable.

Cette méthode, avec $u_0 = 0$, revient à prendre pour valeurs approchées de $\alpha$ les réduites de son développement en fraction continue, c’est-à-dire à approcher $\alpha$ par $u_n=f^{(n)}(0)$ où l’on note

$$f^{(n)} = \underbrace{f\circ\cdots\circ f}_{n \text{ facteurs}} \text{ si } n > 0 \text{ et } f^{(0)} = \mathrm{id}.$$


En classe

Avec une classe de Terminale STI2D de niveau modeste mais curieuse, on s’est contenté de :

 Retravailler la résolution des équations du second degré, sur l’équation $x=\frac{1}{2+x}$. Il a fallu du temps et de l’aide pour comprendre que c’en est une et la ramener à la forme $ax^2+bx+c=0$ ;

 Admirer l’écriture de $\alpha$ en fraction continue, dont la beauté a frappé certains élèves ;

 Rechercher bien sûr les décimales de $\alpha$ (c’était l’objectif). En exploitant la monotonie de $f$, facile à dériver, on a écrit des encadrements de $\alpha$, obtenus à la main dans un premier temps afin de se familiariser avec la relation de récurrence, puis avec une calculatrice. Ci-dessous est reproduite la table des valeurs obtenue par les élèves, où le nombre $\alpha$ est encadré par $u_n=f^{(n)}(0)$ et $v_n=f^{(n)}(1)$. On vérifie que la convergence est rapide, mais que la calculatrice ne nous dévoile pas plus que quelques décimales.

On constatera qu’on a utilisé une calculatrice et non pas un langage de programmation. Dans Python par exemple, avec les réglages par défaut, il aurait été possible de dépasser la dizaine de décimales correctes, mais pas beaucoup plus. En effet, il y a deux difficultés en Python sur lesquelles nous allons revenir : ce qu’on approche, c’est $u_n$ et non pas $\alpha$. De plus, les nombres rationnels y sont mal calculés.

Avec les élèves, on ne s’est finalement pas écarté de l’exercice de seuil typique de BAC, où l’on approche en général une limite par une seule suite, par exemple $(u_n)$. Le contraire nous aurait forcé à aborder avec les élèves quantité de questions intéressantes mais très chronophages. Par exemple :

 Comment s’assurer que les décimales sont exactes ? C’est l’encadrement qui semble l’assurer. Mais ce ne sont que des encadrements et non des approximations.
 Et si par hasard on veut plus de décimales, comment s’y prendre ?

Un programme Python de ce type pourrait permettre d’aborder quelques uns de ces sujets qui fâchent :

  1. def Approx(seuil):
  2.     n=0, u=0
  3.     while 1/4^n>seuil:
  4.         u=1/(2+u)
  5.         n=n+1
  6.     return u

Télécharger

ou :

  1. (from math import sqrt)
  2. def approx(epsilon):
  3.     a=sqrt(2)-1
  4.     u=0
  5.     while abs(a-u)>epsilon:
  6.         u=1/(2+u)
  7.     return u

Télécharger

ou encore :

  1. def calcul_u(n):
  2.     u=0
  3.     for k in range(n):
  4.         u=1/(2+u)
  5.     return u

Télécharger

Programmes 1 et 2.
D’abord, le programme fournit certes un terme $u_n$ arbitrairement proche de $\alpha$, qui en valeur absolue n’en diffère pas de plus de $10^{-p}$. Mais les $p$ premières décimales de $u_n$ sont-elles bien celles de $\alpha$ ? Pas toujours. Il y a là une subtilité qu’il faut bien expliquer.
Ensuite, ce n’est pas l’irrationnel $\alpha$ qu’on développe mais le rationnel $u_n$ qui l’approche. Comme ils ne sont pas égaux, il est inutile de chercher à connaître plus de décimales que celles qui sont garanties.
Au delà des $p$ premières, les décimales sont a priori les mauvaises, avant, elles sont douteuses tant qu’on ne les a pas regardées. On ne le sait pas à l’avance.

Programme 2.
Enfin, pour un ordinateur représentant les nombres en base 2, du fait qu’une telle représentation est finie, les seuls rationnels qui sont représentés sans erreur sont ceux pour lesquels le dénominateur se réduit à une puissance de 2 dans leur écriture en fraction irréductible. La conséquence est qu’une machine est en général incapable de calculer sans arrondi $\frac{1}{2+x}$, et que les erreurs s’accumulant, la précision des calculs ne dépasse pas une ou deux dizaines de décimales.

Donc non seulement, on ne calcule pas le bon nombre, mais en plus on le calcule mal.

Bref, ce n’est pas si simple, les calculs sur machines sont si douteux, qu’il est plus rassurant de construire des encadrements (en faisant confiance aux premières décimales trouvées par la machine).


Trompeuses machines

Nous fiant à la calculatrice, on a ainsi pu déterminer les sept premières décimales après la virgule, puisqu’elles coïncident dans $u_{10}$ et $v_{10}$. Mais dans quelle mesure pouvons-nous faire confiance à la machine ? Et comment obtenir davantage de décimales ? Ces questions, qui ne sont pas si simples comme on va le voir, n’ont pas été abordées avec les élèves (car trop long et trop difficile).

Pour déterminer plus de décimales, il faut changer d’environnement et utiliser un des programmes donnés ci-dessus, par exemple le second. Quitte à prendre $n$ assez grand, on peut en théorie approcher $\alpha$ d’aussi près qu’on veut par $u_n$. Mais les problèmes commencent dès qu’on passe à la pratique. Par exemple, la fonction Python approx suivante, assez typique d’une recherche de seuil au BAC, qui prétend trouver le terme $u_n$ d’indice minimum vérifiant $|u_n-\alpha|\le\varepsilon$, entre en boucle infinie dès qu’on l’appelle avec $\varepsilon=10^{-17}$ :

  1. from math import sqrt                                              
  2. def approx(epsilon):                                              
  3.     a=sqrt(2)-1                                                    
  4.     u=0                                                            
  5.     while abs(a-u)>epsilon:                                        
  6.         u=1/(2+u)                                                  
  7.     return u  

Télécharger

Que se passe-t-il ? Voici les valeurs de $(u_n)$ calculées par récurrence avec Python, avec les décimales qui diffèrent de celles de $\alpha$ en rouge :

u(5)=0.4142857142857143
u(10)=0.41421355164605467
u(15) = 0.41421356237468987
u(20) = 0.41421356237309487
u(25) = 0.4142135623730951
u(30) = 0.4142135623730951
u(35) = 0.4142135623730951

Les valeurs convergent, mais pas vers $\alpha$... Le nombre $\alpha$ est lui même mal approché par Python, qui affiche deux décimales fausses (0.41421356237309515). Avec cette valeur de $\alpha$ et les valeurs de $u_n$ qui stationnent selon la machine, l’écart $\alpha-u_n$ reste bloqué sur $5\times 10^{-17}$, et on ne sort jamais de la boucle.

De toute façon, tout est faux dans ces calculs, sur lesquels un doute indélébile pèse maintenant. Rappelons-le une machine, dans ce contexte, est en général incapable de calculer sans arrondi $\frac{1}{2+x}$, et puisque les erreurs s’accumulent, la précision des calculs ne dépasse pas une ou deux dizaines de décimales. Alors qu’un logiciel de calcul formel comme Xcas, lui, est capable de donner autant de décimales correctes que souhaité. Faire des calculs sur des flottants, c’est s’exposer à la malédiction des arrondis, et est particulièrement inadapté lorsqu’on cherche à calculer des décimales. Ajoutons que même en "trafiquant" Python pour avoir une précision supérieure, on ne fait que repousser le problème, et de toute façon, on ne calcule pas $\alpha$ mais un $u_n$, donc il est inutile de calculer plus de chiffres, qu’il n’en a en commun avec $\alpha$, ce qui n’est pas du tout facile à prévoir. Pour résumer donc, nous l’avons déjà dit, on ne calcule pas les bons nombres, et on les calcule mal. L’affaire semble mal engagée. Et histoire de s’en convaincre, voici un exemple encore plus pathologique !

Exemple :
Majorer l’erreur ne garantit aucune décimale. Voici une première difficulté liée à la base de numération. Supposons que nous voulions 12 décimales de $\alpha$. Certes, on sait trouver un rang $n$ à partir duquel $|u_n-\alpha|$ < $10^{-12}$ (en effet, on a vu que $|u_n-\alpha|$ < $\frac{1}{4^n}$, il suffit donc de résoudre une inéquation). Peut-on en conclure que $\alpha$ a des décimales communes avec $u_n$ ? Certainement pas avant d’avoir vu le nombre $u_n$. Imaginons par exemple que $\alpha=0,414 + r$ avec 0 < r < $10^{-20}$ et que l’algorithme donne $u_n=0,413 999 999 999 999=0,414-10^{-15}$. On a $|u_n-\alpha|=r+10^{-15}$ < $10^{-12}$. Pourtant, très peu de décimales de $\alpha$ se retrouvent dans $u_n$. C’est ce qui se passe si on a le malheur de tomber sur une série de 9 consécutifs. C’est gênant parce que pour un nombre de décimales requis donné, on ne peut pas savoir à l’avance combien de termes $u_n$ sont nécessaires pour les garantir.


Second algorithme (streaming)


Annonçons d’emblée que l’algorithme que nous allons présenter a été inventé pour résoudre explicitement ces problèmes. C’est un exemple de travail mathématique servant à résoudre un problème d’ordre informatique.
En 2005, l’informaticien théorique Jeremy Gibbons publie un article dans lequel il présente un algorithme qui résout les problèmes que nous avons rencontrés auparavant : il ne calcule pas une approximation de la cible, mais il calcule directement le vrai nombre. De plus, il n’utilise que des entiers, et ainsi il n’y a donc pas d’erreurs d’arrondis. Pour se débarrasser de ces difficultés, Gibbons n’a pas amélioré les machines, il a changé d’algorithme. Il y a de nombreux exemples de résultats mathématiques obtenus par l’intervention de l’informatique, nous avons ici un exemple de la situation contraire, où la réflexion mathématique vient résoudre un problème informatique.
Comme souvent en sciences, la découverte de Gibbons s’inscrit dans une généalogie de travaux précédents. Ici, l’histoire remonte à 1968 lorsque A. Sale découvre un algorithme dit compte-goutte, qui permet de calculer les décimales de la constante e, base des logarithmes népériens. Ce travail est repris par d’autres, jusqu’à ce que Gibbons réussisse à prendre les choses à l’envers et à rendre infini son algorithme.

On vient de voir que pour tout $x\in[0 ;1]$, la suite de terme général $f^{(n)}(x)$ converge vers $\alpha$. Autrement dit, la suite de fonctions $(f^{(n)})$ converge simplement sur $[0 ;1]$ vers la fonction constante égale à $\alpha$.

C’est pourquoi à la façon de Jeremy Gibbons, un informaticien théorique britannique, on se permettra d’écrire $\alpha$ comme composée infinie $\alpha= f\circ f \circ f \circ\cdots$ ou plus explicitement encore :

$$\alpha= \left(\dfrac{1}{2+x}\right)\circ \left(\dfrac{1}{2+x}\right) \circ \left(\dfrac{1}{2+x}\right) \circ\cdots$$


Ce genre d’écritures inspire Gibbons. II a imaginé une manière de composer un nombre infini de fonctions qui fournit au fur et à mesure autant de précision que le permet le nombre de facteurs utilisés, et qui commence par la gauche, ce qui permet de les utiliser toutes (du moins en théorie), au lieu de prescrire un nombre $n$ donné de facteurs et de négliger le reste. En 2004, il publie un algorithme qui appliqué ici, se résume à :

  1. Poser n=0
  2. Poser w(x)=x
  3. Répéter indéfiniment
  4.     Si w(0) et w (1) ont même partie entière d
  5.         Afficher d
  6.         Poser n=n+1
  7.         Poser w(x)=10(w(x)-d)
  8.     Sinon  
  9.         Poser w(x)=w(f(x))

Télécharger

C’est tout !

On constate que cet algorithme est court et simple, mais on devine mal ce qu’il fait à première vue. Il se trouve qu’il affiche les décimales de $\alpha$ les unes à la suite des autres. Voyons pourquoi :


Explication de texte

On note

$$\alpha = d_0,d_1 d_2d_3\dots \text{ le développement décimal de } \alpha$$


de sorte que $d_n$ est la $n$-ième décimale de $\alpha$ ($d_0 = 0$, $d_1 = 4$, $d_2 = 1$, $d_3 = 4$, etc.).

On remarque déjà que $w$ reste pendant toute la phase de calculs une fonction homographique (non constante) à coefficients entiers, même si $w$ va sans cesse changer d’expression. Ceci signifie qu’il existe des entiers $a, b,c,d$ avec $ad-bc\neq 0$ tels que $w(x) =\frac{ax + b }{cx+d}$ (pour tout réel $x\neq \frac{-d}{c}$). Essayons de nous en convaincre. Au départ, $w$ est initialisé à l’identité. C’est une fonction affine non constante, donc un cas particulier de fonction homographique. De plus, chaque boucle compose ou bien $w$ à gauche par une fonction affine non constante à coefficients entiers, ou bien $w$ à droite par la fonction homographique $f$. La composée de deux fonctions homographiques en étant encore une, on comprend que le caractère homographique de $w$ est conservé à tout moment de l’algorithme.

En tant que fonction homographique, $w$ est en particulier une fonction monotone. Un simple calcul de la dérivée permet de se convaincre en effet que son dénominateur est un carré et que son numérateur est une constante non nulle qui peut être positive ou négative.

Une autre propriété essentielle se conserve tout au long de l’algorithme et ceci malgré les changements successifs de la fonction $w$ :
pour tout entier naturel $n\in\mathbb{N}$, pour la fonction $w$ associée au rang $n$ et telle que les parties entières de $w(0)$ et $w(1)$ sont égales, on a

$$w(\alpha)=d_n,d_{n+1}d_{n+2}d_{n+3}\dots$$

Démonstration :
Notons $(P_n)$ cette propriété.
 $(P_0)$ est vraie à l’initialisation puisque $w(\alpha) = \alpha = d_0,d_1d_2d_3\dots$
 Supposons maintenant que $(P_n)$ est vraie pour un certain entier $n\in\mathbb{N}$.
Si $w(0)$ et $w(1)$ n’ont pas la même partie entière, $w$ prend d’après l’algorithme la valeur $w \circ f$, ce qui ne change pas la valeur de $w(\alpha)$ car $f(\alpha) = \alpha$. Le nombre $n$ lui n’évolue pas non plus (on y ajoute 1 seulement quand $w(0)$ et $w(1)$ ont même partie entière). La propriété $(P_{n})$ se conserve donc tant que $w(0)$ et $w(1)$ n’ont pas la même partie entière, bien que la fonction homographique $w$ évolue).
Si par contre $w(0)$ et $w(1)$ ont la même partie entière $d$, alors nécessairement $d$ est égal à $d_n$, car par monotonie de $w$, $w(0)$ et $w(1)$ encadrent $w(\alpha)=d_n,d_{n+1}d_{n+2}d_{n+3}\dots$ Dans ce cas, l’algorithme fait prendre à $n$ la valeur $n +1$ et à $w$ la valeur $10(w - d_n)$. Or lorsqu’on applique à $\alpha=d_n,d_{n+1}d_{n+2}d_{n+3}\dots$ la fonction $10(w - d_n)$, on obtient $d_{n+1},d_{n+2}d_{n+3}d_{n+4}\dots$, ce qui établit $(P_{n+1})$.

Retenons ce qui est apparu dans le raisonnement :
Si $w(0)$ et $w(1)$ ont même partie entière $d$, alors cette valeur commune est $d_n$, la $n$-ième décimale de $\alpha$. Ceci explique que l’algorithme teste l’égalité des parties entières, c’est le signal qu’une décimale de $\alpha$ est déterminée (qui est alors affichée). Il y a alors une instruction modifiant $w$ pour que $w(\alpha)$ ait $d_{n +1}$ pour partie entière, ce qui adapte $w$ à la recherche de la décimale suivante (bien sûr, on ne connaît pas $w(\alpha)$, mais on peut calculer $w(0)$ et $w(1)$ qui l’encadrent). L’entier $n$ s’interprète lui comme le nombre de décimales découvertes.

Correction de l’algorithme. Il reste à prouver que l’algorithme finit toujours par produire une décimale. C’est bien ce qui arrive, quel que soit l’état des variables car quitte à composer $w$ à droite par assez de facteurs $f$, $w(0)$ et $w(1)$ finissent bien par avoir même partie entière. Plus précisément, il existe un entier $k\in\mathbb{N}$ tel que $w(f^{(k)}(0))$ et $w(f^{(k)}(1))$ ont même partie entière. En effet, si ce n’était pas le cas, il existerait pour tout $k\in\mathbb{N}$ un entier $n_k$ coincé entre $w(f^{(k)}(0))$ et $w(f^{(k)}(1))$. Or, par continuité de $w$, ces deux nombres sont les termes de deux suites convergeant vers $w(\alpha)$. Donc la suite $(n_k)_{k\in\mathbb{N}}$ convergerait aussi vers $w(\alpha)$ d’après le théorème des gendarmes. Mais une suite d’entiers qui converge est stationnaire, donc $w(\alpha)$ serait un entier $m$ et $\alpha = w^{-1}(m)$ serait rationnel, ce qui constitue une contradiction avec la définition même de $\alpha$.

Ceci explique le choix fait par l’algorithme de composer $w$ à droite par $f$ lorsque $w(0)$ et $w(1)$ ont des parties entières distinctes, c’est-à-dire lorsqu’il reste une ambiguïté sur la décimale cherchée. Cela revient à pousser l’indice $k$ jusqu’à ce que l’écart entre $f^{(k)}(0)$ et $f^{(k)}(1)$ soit juste assez réduit pour que leurs images par $w$ ne diffèrent plus par les parties entières.

Par comparaison avec les suites du premier algorithme, du type $f^{(k)}(x_0)$ avec $x_0\in[0 ;1]$, qui approchent $\alpha$ par itération de $f$, l’algorithme de Gibbons ne consomme pour chaque décimale que le strict nécessaire parmi les facteurs $f$ pour le produire. Cette opposition entre les phases de consommation de l’information encodée dans une écriture (instruction $w$ prend la valeur $w\circ f$) et de production d’une information dans une autre (instructions Afficher $d$ et $w$ prend la valeur $10(w-d)$) est caractéristique de sa méthode de streaming.

L’algorithme de Gibbons est :
Simple : Il tient en quelques lignes et s’exécute presque à la main.
Rapide : Au moins 3 décimales toutes les 5 étapes, comme pour les suites sur lesquelles il est basé.
Sans fin : Alors qu’avec des approximations par des suites, on doit se fixer une précision décimale au préalable, cet algorithme produit indéfiniment des digits sans jamais s’arrêter. C’est unique.
Exact : Ici, lorsqu’un digit est produit, c’est le bon, ce n’est pas celui d’une approximation de $\alpha$. De plus, on ne travaille qu’avec des fractions, donc des entiers. C’est important en programmation, car des erreurs d’approximation apparaissent quasi-systématiquement lorsqu’on travaille avec des flottants.


Programmation en Python

1. $w(x)$ est de la forme $\frac{ax+b}{cx+d}$
On représentera $w$ dans le code par quatre variables entières $a, b, c , d,$ respectivement initialisées à 0, 1, 1, 2. Les changements d’état de $w$ se font en modifiant les coefficients $a, b, c, d$ sans jamais définir $w$.
2. Pour le test, inutile d’utiliser la fonction partie entière. La partie entière de $w(0)$ est le quotient de $b$ par $d$, et celle de $w(1)$ est le quotient de $a + b$ par $c + d$. Ceci permet d’écrire la condition du test en n’utilisant que des entiers. Ce qui est un énorme avantage.
3. Pour y voir quelque chose, il faut renoncer au streaming infini (ça défile trop vite), se contenter de disons 10 000 digits, et en demander l’affichage par lignes de 40.

Version minimaliste

  1. def streamrac2():
  2.     a,b,c,d=0,1,1,2
  3.     while True:
  4.         y=b//d
  5.         if y==(a+b)//(c+d):
  6.             print(y)
  7.             a=10*(a-c*y)
  8.             b=10*(b-d*y)
  9.         else:
  10.             a,b=b,a+2*b
  11.             c,d=d,c+2*d

Télécharger

Version avec affichage plus lisible

  1. def streamrac2(N):
  2.   # Affiche la suite des N premiers digits de sqrt(2)-1 par streaming
  3.   L=""                   # Chaîne de digits à afficher.
  4.   k=0                    # Nombre de digits trouvés.
  5.   a,b,c,d=0,1,1,2        # Initialisation de w(x)=(ax+b)/(cx+d).
  6.   while k<N:            # While True pour streaming infini
  7.     y=b//d               # Candidat digit = part. entière de w(0).
  8.     if y==(a+b)//(c+d):  # Cond. pour produire le digit y.
  9.       a=10*(a-c*y)       # M.à.J : Compos. par 10(x-y) à gauche.
  10.       b=10*(b-d*y)
  11.       L=L+str(y)
  12.       k=k+1
  13.       if len(L)>=20:     # Affichage et vidange tous les 20 digits.
  14.         print(L)
  15.         L=""
  16.     else:                # Consommation d'un facteur 1/(2+x) :
  17.       a,b=b,a+2*b        # M.à.J : Compos. par 1/(2+x) à droite.
  18.       c,d=d,c+2*d

Télécharger

Il reste alors à tester l’algorithme dans sa version la plus lisible avec l’aide de EduPython par exemple. Voici l’affichage (instantané) obtenu pour N =40.

On en déduit que :
$\alpha \simeq$ 0,414213562373095048801688724209698078569 avec l’assurance d’avoir 39 décimales correctes.
Ou encore : $\sqrt{2} \simeq$ 1,414213562373095048801688724209698078569 toujours avec la même assurance. On en trouvera une confirmation sur cette page.
Par curiosité, on peut pousser un peu l’ordinateur pour voir le temps de calcul pour 400 décimales, puis 4000 etc. Il suffit alors d’importer "time" et d’ajouter quelques lignes à l’algorithme qui devient ainsi :

  1. import time
  2. def stream2rac2(N):
  3.   start_time = time.time()
  4.   # Affiche la suite des N premiers digits de sqrt(2)-1 par streaming
  5.   L=""                   # Chaîne de digits à afficher.
  6.   k=0                    # Nombre de digits trouvés.
  7.   a,b,c,d=0,1,1,2        # Initialisation de w(x)=(ax+b)/(cx+d).
  8.   while k<N:            # While True pour streaming infini
  9.     y=b//d               # Candidat digit = part. entière de w(0).
  10.     if y==(a+b)//(c+d):  # Cond. pour produire le digit y.
  11.       a=10*(a-c*y)       # M.à.J : Compos. par 10(x-y) à gauche.
  12.       b=10*(b-d*y)
  13.       L=L+str(y)
  14.       k=k+1
  15.       if len(L)>=20:     # Affichage et vidange tous les 20 digits.
  16.         print(L)
  17.         L=""
  18.     else:                # Consommation d'un facteur 1/(2+x) :
  19.       a,b=b,a+2*b        # M.à.J : Compos. par 1/(2+x) à droite.
  20.       c,d=d,c+2*d
  21.   end_time = time.time()
  22.   execution_time = end_time - start_time
  23.   print(f"Temps d'exécution : {execution_time} secondes")

Télécharger

Sur ma machine plutôt vieillotte, une milliseconde suffit pour afficher 400 décimales.
Le tableau suivant montre l’étendue des talents de l’algorithme :

nombre de décimales 100 1 000 10 000 100 000 1 000 000
Temps de calcul en s <$1 \times 10^{-3}$ $5 \times 10^{-3}$ 0,2 15,1 1 440

On constate qu’il faut 1 440 secondes soit 24 minutes pour calculer un million de décimales de $\sqrt{2}$ que l’on sait toutes justes, sur une machine assez ancienne. C’est époustouflant ! Les décimales pleuvent réellement en cascade !


Conclusion
Nous l’avons vu, le traditionnel exercice de bac qui permet d’estimer une valeur approchée d’un nombre comme $\sqrt{2}$ ou le nombre d’argent a des côtés insatisfaisants. On s’en contente pour l’examen. Mais pourtant, les décimales issues des procédures d’utilisation de suites récurrentes ne produisent pas à coup sûr des décimales correctes du nombre cherché, même avec Python qui éprouve des difficultés avec les flottants. Toutes ces questions sont capitales pour le calcul numérique, mais malheureusement faute de temps dans les classes, il est assez illusoire de les aborder sereinement. Pourtant l’algorithme de streaming de Gibbons est d’une remarquable efficacité. Mais en expliquer les méandres et les motivations serait chronophage au lycée et en dehors des espaces balisés par le programme.
On peut se demander d’ailleurs si cet algorithme est applicable à d’autres nombres. Évidemment oui, au moins aux nombres irrationnels dits quadratiques, c’est-à- dire ceux qui sont solutions d’une équation du second degré, et avec un développement en fraction continue 1-périodique. Cela permet alors de définir une fonction $f$ unique pour la composition. On pense naturellement aux radicaux et au nombre d’or $\varphi$, dont nous parlions en début d’article. Il est en effet lui aussi solution d’une équation du second degré (à savoir $x²-x-1=0$) et un développement en fraction continue est donné par :

$$\varphi=1+ \frac{1}{1+\frac{1}{1+\frac{1}{1+ \frac{1}{1+\cdots}}}}$$


Bonne nouvelle, on peut calculer sans grande difficulté les décimales de $\varphi - 1$.
Pour ce faire, il suffit d’adapter l’algorithme ainsi. On remplace la fonction $f(x) = \dfrac{1}{2+x}$ par $f(x) = \dfrac{1}{1+x}$.
Dans Python, on initialise alors les quatre variables entières a, b, c , d à respectivement 0, 1, 1, 1.
La partie "Else" de la boucle est également à modifier. Elle devient :

else:                # Consommation d'un facteur 1/(1+x) :
     a,b=b,a+1*b        # M.à.J : Compos. par 1/(1+x) à droite.
     c,d=d,c+1*d

Et on laisse agir !
On obtient environ 0,6180339887498948482. Ce qui donne le nombre d’or $\varphi \simeq 1,6180339887498948482$.
L’algorithme de streaming n’est donc pas l’algorithme d’un seul nombre. On peut facilement adapter celui-ci à l’usage que l’on veut faire.
Et si l’on s’attaque à des nombres irrationnels quadratiques mais 2-périodiques ? C’est un peu plus technique, mais cela fonctionne tout aussi bien, on peut le voir avec le nombre $\sqrt{3}$ dont voici le développement en fraction continue :

$$\sqrt{3}=1+ \frac{1}{1+\frac{1}{2+\frac{1}{1+ \frac{1}{2+\cdots}}}}$$


Autrement dit, la composition à droite alterne entre les deux fonctions $f(x) = \dfrac{1}{1+x}$ et $g(x) = \dfrac{1}{2+x}$. Il convient donc d’introduire un compteur de compositions à droite. On le note $p$. Lorsque celui-ci aura une valeur paire, on compose par la fonction $g$ et sinon par la fonction $f$. Cela permet d’obtenir les décimales du nombre $\sqrt{3} - 1$ avec le programme Python suivant :

  1. def streamrac3(N):
  2.   # Affiche la suite des N premiers digits de sqrt(3)-1 par streaming
  3.   L=""                   # Chaîne de digits à afficher.
  4.   k=0                    # Nombre de digits trouvés.
  5.   p=1                    # Compteur de compositions à droite.
  6.   a,b,c,d=0,1,1,1      # Initialisation de w(x)=(ax+b)/(cx+d).
  7.   while k<N:            # While True pour streaming infini
  8.     y=b//d               # Candidat digit = part. entière de w(0).
  9.     if y==(a+b)//(c+d):  # Cond. pour produire le digit y.
  10.       a=10*(a-c*y)       # M.à.J : Compos. par 10(x-y) à gauche.
  11.       b=10*(b-d*y)
  12.       L=L+str(y)
  13.       k += 1
  14.       if len(L)>=20:     # Affichage et vidange tous les 20 digits.
  15.         print(L)
  16.         L=""
  17.     else:
  18.       p=p+1             # M.à.J du compteur de compositions à droite.
  19.       if p%2 != 0 :            #test de parité pour la composition à droite par l'une des 2 fonctions :
  20.          a,b=b,a+1*b        # M.à.J : Compos. par 1/(1+x) à droite.
  21.          c,d=d,c+1*d
  22.       else:
  23.          a,b=b,a+2*b        # M.à.J : Compos. par 1/(2+x) à droite.
  24.          c,d=d,c+2*d

Télécharger

Voici alors un affichage pour 100 digits :
07320508075688772935
27446341505872366942
80525381038062805580
69794519330169088000
37081146186757248575
A rapprocher bien évidemment de $\sqrt{3} \simeq$ 1,732050807568877....

Et quand on connaît le théorème de Lagrange (Un irrationnel est quadratique si et seulement si son développement en fraction continue est périodique à partir d’un certain rang), cela ouvre de belles perspectives de calcul !


Annexe d’Angelo Laplace pour les calculateurs acharnés :

Dans ce qui suit, on utilise les notations traditionnelles pour le développement en fraction continue périodique des nombres irrationnels quadratiques.
Ainsi, pour le nombre d’or, on écrit : $\varphi = [\overline{1}]$ pour indiquer que le nombre 1 se répète à l’infini dans le développement qui s’écrit : $\varphi=1+ \frac{1}{1+\frac{1}{1+\frac{1}{1+ \frac{1}{1+\cdots}}}}$.

Le tableau suivant donne des informations sur le développement en fraction continue des premiers radicaux non entiers, afin de réaliser facilement les modifications dans le programme Python de streaming. A chaque fois, le programme calculera le nombre : $\sqrt{x}$ - partie entière de $\sqrt{x}$, autrement dit, uniquement la "partie décimale" du nombre ciblé. A noter que par exemple $\sqrt{8}$ n’est pas traité car ce nombre vaut également $2\sqrt{2}$.

Cible en fraction continue période initialisation composition à droite
$\sqrt{2}$ [1,$\overline{2}$] 1-périodique $a,b,c,d=0,1,1,2$ $a,b=b,a+2*b$ et $c,d=d,c+2*d$
$\sqrt{3}$ [1,$\overline{1,2}$] 2-périodique $a,b,c,d=0,1,1,1$
p=p+1
     if p%2 != 0 :            
        a,b=b,a+1*b      
        c,d=d,c+1*d
     else:
        a,b=b,a+2*b        
        c,d=d,c+2*d
$\sqrt{5}$ [2,$\overline{4}$] 1-périodique $a,b,c,d=0,1,1,4$ $a,b=b,a+4*b$ et $c,d=d,c+4*d$
$\sqrt{6}$ [2,$\overline{2,4}$] 2-périodique $a,b,c,d=0,1,1,2$
p=p+1
     if p%2 != 0 :            
        a,b=b,a+2*b      
        c,d=d,c+2*d
     else:
        a,b=b,a+4*b        
        c,d=d,c+4*d
$\sqrt{7}$ [2,$\overline{1,1,1,4}$] 4-périodique $a,b,c,d=0,1,1,1$
p=p+1
     if p%4 != 0 :            
        a,b=b,a+1*b      
        c,d=d,c+1*d
     else:
        a,b=b,a+4*b        
        c,d=d,c+4*d
$\sqrt{10}$ [3,$\bar{6}$] 1-périodique $a,b,c,d=0,1,1,6$ $a,b=b,a+6*b$ et $c,d=d,c+6*d$
$\sqrt{11}$ [3,$\overline{3,6}$] 2-périodique $a,b,c,d=0,1,1,3$
p=p+1
     if p%2 != 0 :            
        a,b=b,a+3*b      
        c,d=d,c+3*d
     else:
        a,b=b,a+6*b        
        c,d=d,c+6*d
$\sqrt{13}$ [3,$\overline{1,1,1,1,6}$] 5-périodique $a,b,c,d=0,1,1,1$
p=p+1
     if p%5 != 0 :            
        a,b=b,a+1*b      
        c,d=d,c+1*d
     else:
        a,b=b,a+6*b        
        c,d=d,c+6*d

Le caractère répétitif des modifications du programme nous a fait envisager de changer quelque peu l’aspect du programme de base et d’inclure davantage de paramètres. Ainsi, vous trouverez ci-dessous deux programmes tout prêts pour calculer en cascade les décimales des nombres irrationnels 1-périodiques et 2-périodiques.
De bonnes surprises ne sont pas à exclure, par exemple, $\sqrt{7}$ est certes 4-périodique mais son développement en fraction continue ne contient dans sa période que 2 chiffres différents, trois 1 successifs et un 4 qui se répètent tous les 4 chiffres. Aussi, il était aisé de modifier le programme prévu pour les nombres 2-périodiques pour l’adapter à $\sqrt{7}$. Le nombre $\sqrt{13}$ est lui 5-périodique mais sa période est elle aussi tout à fait coopérative !
Les premières grosses complications arrivent avec $\sqrt{14}$ qui est 4-périodique, de développement [3,$\overline{1,2,1,6}$].

Irrationnels 1-périodiques :
Voici un programme, utilisant davantage de variables, adapté aux nombres irrationnels 1-périodiques, comme le nombre d’or et $\sqrt{2}$ :

  1. def streamrac(N,a,b,c,d,p):
  2.   # Affiche la suite des N premiers digits de la cible par streaming, la période est p.
  3.   L=""                   # Chaîne de digits à afficher.
  4.   k=0                    # Nombre de digits trouvés
  5.   while k<N:            # While True pour streaming infini
  6.     y=b//d               # Candidat digit = part. entière de w(0).
  7.     if y==(a+b)//(c+d):  # Cond. pour produire le digit y.
  8.       a=10*(a-c*y)       # M.à.J : Compos. par 10(x-y) à gauche.
  9.       b=10*(b-d*y)
  10.       L=L+str(y)
  11.       k=k+1
  12.       if len(L)>=20:     # Affichage et vidange tous les 20 digits.
  13.         print(L)
  14.         L=""
  15.     else:                # Consommation d'un facteur 1/(p+x) :
  16.       a,b=b,a+p*b        # M.à.J : Compos. par 1/(p+x) à droite.
  17.       c,d=d,c+p*d

Télécharger

Dans celui-ci, le paramètre $p$ représente le chiffre qui se répète infiniment dans le développement en fraction continue du nombre ciblé.
Dans le cas du nombre d’or, les paramètres sont (N,0,1,1,1,1) et dans le cas du nombre d’argent, ce sont (N,0,1,1,2,2). Voici l’affichage obtenu :

Il est conforme à ce que nous espérions. L’affichage du 0 est à remplacer par la partie entière du nombre souhaité.
L’avantage de cette façon de programmer est qu’il n’est plus nécessaire de modifier les fonctions de composition dans le corps du programme, il suffit d’appeler la fonction "streamrac" avec des paramètres différents quand on change de nombre cible.

Irrationnels 2-périodiques :
Ici les paramètres $p$ et $q$ représentant les deux chiffres dans l’ordre de la période.

  1. def streamrac(N,a,b,c,d,p,q):   #cas 2-périodique, de période p,q
  2.   # Affiche la suite des N premiers digits de la cible par streaming
  3.   L=""                   # Chaîne de digits à afficher.
  4.   k=0                    # Nombre de digits trouvés
  5.   i=1                   #initialisation du compteur
  6.   while k<N:             # While True pour streaming infini
  7.     y=b//d               # Candidat digit = part. entière de w(0).
  8.     if y==(a+b)//(c+d):  # Cond. pour produire le digit y.
  9.       a=10*(a-c*y)       # M.à.J : Compos. par 10(x-y) à gauche.
  10.       b=10*(b-d*y)
  11.       L=L+str(y)
  12.       k=k+1
  13.       if len(L)>=20:     # Affichage et vidange tous les 20 digits.
  14.         print(L)
  15.         L=""
  16.     else:
  17.         i=i+1               # M.à.J du compteur
  18.         if i%2 != 0 :       # si le compteur est impair
  19.            a,b=b,a+p*b      # Compos. à droite par 1/(x+p)
  20.            c,d=d,c+p*d
  21.         else:
  22.            a,b=b,a+q*b      # si le compteur est pair
  23.            c,d=d,c+q*d      # Compos. à droite par 1/(x+q)

Télécharger

D’après ce qui précède, l’appel de la fonction "streamrac" est : streamrac(20,0,1,1,2,2,4) pour calculer 20 décimales de $\sqrt{6}$ ou streamrac(20,0,1,1,1,1,2) pour 20 décimales de $\sqrt{3}$.

Le cas de $\sqrt{14}$

$\sqrt{14}$ = [3,$\overline{1,2,1,6}$].
Un nouveau paramètre $r$ est nécessaire. Une nouvelle boucle If est nécessaire également. Lorsque le compteur de composition est impair, il faut composer par $\frac{1}{1+x}$. Ce sera le premier test préalable à la composition à droite. Si le compteur est pair, un second test sera nécessaire pour trancher si le compteur de composition est un multiple de 4, auquel cas on compose par $\frac{1}{6+x}$ et sinon par $\frac{1}{2+x}$.
En réalité, pour des nombres de plus en plus grands sous le radical, le problème revient souvent à optimiser les tests de composition à droite dans la boucle réservée.

Voici le programme utilisé pour $\sqrt{14}$ :

  1. def streamrac(N,a,b,c,d,p,q,r):   #sqrt(14)
  2.   # Affiche la suite des N premiers digits de la cible par streaming
  3.   L=""                   # Chaîne de digits à afficher.
  4.   k=0                    # Nombre de digits trouvés
  5.   i=1                    #initialisation du compteur
  6.   while k<N:             # While True pour streaming infini
  7.     y=b//d               # Candidat digit = part. entière de w(0).
  8.     if y==(a+b)//(c+d):  # Cond. pour produire le digit y.
  9.       a=10*(a-c*y)       # M.à.J : Compos. par 10(x-y) à gauche.
  10.       b=10*(b-d*y)
  11.       L=L+str(y)
  12.       k=k+1
  13.       if len(L)>=20:     # Affichage et vidange tous les 20 digits.
  14.         print(L)
  15.         L=""
  16.     else:
  17.         i=i+1               # M.à.J du compteur
  18.         if i%2 != 0 :       # si le compteur est impair
  19.            a,b=b,a+p*b      # Compos. à droite par 1/(x+p)
  20.            c,d=d,c+p*d
  21.         else:
  22.             if i%4 != 0 :         # si le compteur n'est pas multiple de 4
  23.                  a,b=b,a+q*b      # Compos. à droite par 1/(x+q)
  24.                  c,d=d,c+q*d
  25.             else :                 #si le compteur est un multiple de 4
  26.                     a,b=b,a+r*b    # Compos. à droite par 1/(x+r)
  27.                     c,d=d,c+r*d

Télécharger

Il convient alors d’appeler : streamrac(40,0,1,1,1,1,2,6).

Les décimales de chacune de nos cibles pleuvent en cascade depuis le début de l’article. Théoriquement en continu et infiniment...Les possibilités de nombres à cibler sont toutes aussi infinies.
Le lecteur est invité à poursuivre la cascade au-delà de $\sqrt{14}$ si le cœur lui en dit.