Algorithme de Babylone : une boucle sous toutes ses formes

une boucle ? une suite !

Qu'est-ce qu'une boucle?


Une héroine d'un conte pour enfants ?

Boucle d'or

Ou bien ceci ?

Pseudo-langage
r  r0
pour k de 1 jusque n faire
  r  f(r)
fpour
retourne r 

Mais non, bien sûr, C'est une suite:

$$r_{k+1}=f(r_k)=f^{k+1}(r_0) \ {\rm avec } \ r_0=r0$$

L'algorithme peut d'ailleurs être écrit récursivement:

Pseudo-langage
boucle(r : type 1, n : entier) : type 2
  si n <= 0 alors 
     retourne r
  sinon
     retourne boucle(f(r), n - 1)
  fsi
Orbite d'un système dynamique
L'orbite issue de $r_0$ pour la fonction $f$ est la suite $(r_n)$ définie par $r_n=f^n(r_0)$


En programmation fonctionnelle, on traduit une boucle par un pliage. Par exemple avec Haskell.

Haskell
Prelude> let boucle invariantBoucle ini ni nf pas = foldl (\ r k -> invariantBoucle r) ini [ni, ni + pas .. nf] 
Prelude> boucle (\ x -> 2*x) 1 1 10 1
1024

On peut avoir l'orbite avec la fonction iterate qui renvoie une liste infinie : c'est bien l'orbite au sens mathématique du terme.

Haskell
let orbite = iterate

Pour n'avoir que les n premiers éléments de l'orbite, on utilise la fonction take:

Haskell
Prelude> take 11 ( orbite (\ x -> 2*x) 1 )
[1,2,4,8,16,32,64,128,256,512,1024]

On peut avoir presque la même chose sur XCAS.
On utilise l'opérateur de composition itérée: f@@n est alors l'itérée n fois de la composée de f par elle-même.

Giac-XCAS
0>> orbite := (f,r0,n) -> (f@@n)(r0) 
1>> f := x -> 2*x
2>> orbite(f,1,10)
1024

Le programmeur voudrait savoir si sa boucle va effectivement renvoyer ce qui est
prévu dans sa spécification et surtout en combien de temps: s'il s'agit de
calculer la vitesse d'une voiture en fonction de la distance qui la sépare de la
voiture qui est en face, il faut faire vite.

exemple d'étude: l'algorithme de Babylone

Il y a presque 4000 ans, les petits Babyloniens calculaient la racine carrée de
2 avec 6 bonnes décimales en utilisant un algorithme élémentaire.

On peut lire $1 + \frac{24}{60} + \frac{51}{60^2}+ \frac{10}{60^3}\approx 1,414 213$

babylone

Pseudo-langage
r  1.0
pour k de 1 jusque n faire
  r  (r + 2/r) * 0.5
fpour
retourne r 

Vous pouvez le traduire instantanément dans n'importe quelle langage avec une boucle pour. Évitez autant que faire se peut les
tant que pour ne pas risquer de lancer des boucles infinies.


Avec XCAS, c'est immédiat: programmer en XCAS c'est écrire en pseudo-langage...mais pour de vrai !


Giac-XCAS
boucle(n):={
   r := 1;
   pour k de 1 jusque n faire
      r := (r + 2/r) * 0.5
   fpour;
   retourne r  
}:;

En récursif:

Pseudo-langage
boucle r n
  si n <= 0 alors 
     retourne r
  sinon
     retourne boucle((r + 2 / r) * 0.5, n - 1)
  fsi

La traduction sur XCAS est carément du mot à mot:

Giac-XCAS
boucle(r,n):={
  si n <= 0 alors
    retourne r
  sinon
    retourne boucle((r + 2 / r) * 0.5, n - 1);
  fsi
}:;

On peut aussi utiliser notre fonction boucle créée avec Haskell:

Haskell
Prelude> boucle (\x -> (x + 2 / x) * 0.5) 1 1 6 1
1.414213562373095

Premier problème : est-ce que la suite des valeurs calculées par la boucle

converge vers $\sqrt{2}$?###

Nous aurons besoin de deux théorèmes que nous admettrons et dont nous ne
donnerons que des versions simplifiées.


|Théorème de la limite monotone|
|---------------------------------------------|
|Toute suite croissante majorée (ou décroissante minorée) converge|




Un théorème fondamental est celui du point fixe. Il faudrait plutôt parler des
théorèmes du point fixe car il en existe de très nombreux avatars qui portent
les noms de mathématiciens renommés: Banach, Borel,
Brouwer, Kakutani, Kleene,…Celui qui nous
intéresserait le plus en informatique est celui de Knaster-Tarski. On
peut en trouver une présentation claire dans
Les démonstrations et les algorithmes.

Ils donnent tous des conditions d'existence de points fixes (des solutions de
$f(x)=x$ pour une certaine fonction).

Nous nous contenterons de la version light vue au lycée.



|Un théorème light du point fixe|
|------------------------------------------|
|Soit I un intervalle fermé de $\mathbb R$, soit $f$ une fonction de Ι vers I et soit
$(r_n)$ une suite d'éléments de I telle que $r_{n+1}=f(r_n)$ pour tout entier
naturel $n$. |
|SI $(r_n)$ est convergente ALORS sa limite est UΝ point fixe de $f$
appartenant à I.|


Cela va nous aider à étudier notre suite définie par $r_{n+1}= f(r_n)=\frac{r_n +
\frac{2}{r_n}}{2}$ et $r_0=1$.

On peut alors démontrer par récurrence que

  • pour tout entier naturel non nul $n$, on a $r_n \geqslant \sqrt{2}$
  • la suite est décroissante
  • $\sqrt{2}$ est l'unique point fixe positif de $f$

On peut alors conclure et être assuré que notre algorithme va nous permettre d'approcher $\sqrt{2}$.

Deuxième problème: quelle est la vitesse de convergence?

Ici, cela se
traduit par « combien de décimales obtient-t-on en plus à chaque
itération?».

Introduisons la notion d'ordre d'une suite:




Définition: Ordre d'une suite - Constante asymptotique d'erreur
Soit $(r_n)$ une suite convergeant vers $\ell$. S'il existe un entier $k>0$ tel que:

$$
\displaystyle\lim_{n\to +\infty}\frac{|r_{n+1}-\ell|}{|r_n-\ell|^k}=C
$$

avec $C\neq 0$ et $C\neq +\infty$ alors on dit que $(r_n)$ est d'ordre $k$ et
que C est la constante asymptotique d'erreur.|




Déterminons l'ordre et la constante asymptotique d'erreur de la suite
de Babylone donnant une approximation de $\sqrt{2}$.

On démontre que

$$
\left|r_{n+1} - \sqrt{2}\right| = \left| \frac{(r_n-\sqrt{2})^2}{2r_n}\right|
$$

Ainsi la suite est d'ordre 2 et sa constante asymptotique est $\frac{1}{2\sqrt{2}}$.

Ici, on peut même avoir une idée de la vitesse sans étude asymptotique. En effet, on a $r_n \geqslant 1$ donc

$$
\left|r_{n+1} - \sqrt{2}\right| \leqslant \left| (r_n-\sqrt{2})^2\right|
$$

Les élèves de terminale étant fort habitué(e)s à manipuler des pH et des concentrations de $H_3O^+$, on peut leur demander de démontrer
que le nombre de décimales d'un nombre positif plus petit que 1 est $-\log_{10}x$.

En posant $d_n=-\log_{10}|r_n-\sqrt{2}|$ on obtient donc que $d_{n+1} \geqslant 2d_n$: à chaque tour de boucle, on double au minimum le
nombre de bonnes décimales.

La précision de la machine étant de 16 décimales si les nombres sont codés en binaires sur 64 bits, il faut au maximum 6 itérations pour
obtenir la précision maximale.

Haskell
Prelude> take 7 ( orbite (\ x -> (x + 2/x) * 0.5) 1 )
[1.0,1.5,1.4166666666666665,1.4142156862745097,1.4142135623746899,1.414213562373095,1.414213562373095]
Prelude> sqrt 2
1.4142135623730951

Ou, si vous préférez avec XCAS:

Giac-XCAS
racine(a,xk,n):={
   si n == 0 alors 
      xk 
   sinon
      racine(a,(xk + a / xk) * 0.5, n - 1);
   fsi
}:;

ou toujours avec XCAS/giac mais en ligne de commande:

Giac-XCAS
0>> DIGITS := 16:;
"Done"
// Time 0
1>> racine(a,xk,n) := si n <= 0 alors xk; sinon racine(a,(xk + a / xk) * 0.5, n - 1); fsi;
// Interprète racine
// Succès lors de la compilation racine
racine: définition récursive
"Done"
// Time 0
2>> racine(2,1,5)
1.414213562373095
// Time 0

ou encore plus simplement avec l'opérateur de composition itérée @@.
Cela correspond donc au pliage de la programmation fonctionnelle : on peut tout faire avec XCAS...
Ainsi:

Giac-XCAS
0>> DIGITS := 16
1>> f := x -> (x + 2/x) * 0.5
2>> racine := n -> (f@@n)(1)
3>> racine(6)
1.414213562373095

Avec Python:

Python
def racine(a,xk,n):
    """
    Calcule une approximation de a à partir d'une première valeur par défaut xk
    On effectue n itérations de l'algorithme de Babylone (ou de Héron)
    """
    if n == 0:
        return xk
    else:
        return racine(a, (xk + a / xk) * 0.5 , n - 1)

Une application : on peut ainsi calculer rapidement sur un petit système
embarqué (voire un smart phone) la vitesse de la
voiture en fonction de la distance la séparant de la voiture qui la précède.

Voici par exemple une application disponible sur iPhone:






safety



En multiprécision

Il se peut que vous ayez besoin d'aller au-delà de la précision de la machine (environs 16 décimales pour 64 bits).
Vous pouvez vous tourner vers des logiciels travaillant en précision infinie naturellement comme par exemple XCAS/giac. Il suffit
alors de fixer DIGITS à une valeur voulue:

Giac-XCAS
0>> DIGITS := 100:;
1>> racine(a,xk,n):= si n <= 0 alors xk; sinon racine(a,(xk + a / xk) * 0.5, n - 1); fsi;
2>> racine(2,1,8)
// avec le doublement de la précision, 8 devrait suffire :
1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573
// on compare avec la valeur donnée par l'approximation native de XCAS:
3>> evalf(sqrt(2),100)
1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573

En fait, XCAS/giac utilise la bibliothèque C MPFR de calcul multi-précision mise au point par une équipe française et qui est devenue la référence mondiale.

On peut l'utiliser via Sage dont la syntaxe de programmation est celle de Python. On précise la taille de la mantisse
utilisée dans la représentation des nombres: 53 en 64 bits ou 113 en 128 bits.

Sage
On travaille sur 128 bits (on donne la précision = taille de la mantisse)
RN = RealField(113)  # utilise MPFR de C
 
def racine(a,xk,n):
    """
    Calcule une approximation de a à partir d'une première valeur par défaut xk
    On effectue n itérations de l'algorithme de Babylone (ou de Héron)
    On travaille sur une mantisse de prec bits fixée globalement via la bibliothèque MPFR
    """
    if n == 0:
        return RN(xk)
    else:
        return racine(RN(a), RN(xk + a / xk) * RN(0.5) , n - 1)

7 itérations suffisent:

Sage
sage: racine(2,1,7)
1.41421356237309504880168872420970
sage: sqrt(2).n(113) 
1.41421356237309504880168872420970

On pourra se référer à ce poly de cours pour plus de précisions sur le traitement des nombres à virgule flottante.

courtesy of webmatter.de