Les scripts d'Archimède

Tous les programmes du DM2, complétés, exécutés et vérifiés, suivis de la table promise par le document : les trente doublements en précision étendue. Les squelettes sont ceux du DM, les lignes en pointillés ont été remplies ; chaque bloc se copie tel quel dans votre éditeur. Convention d’affichage : virgule décimale dans le texte et les tableaux, point décimal dans les blocs de code et leurs sorties.

1. La boucle naïve (Q14) : le programme qui finit par dire que $\pi$ vaut $6$

Le squelette de Q14, complété avec la relation de Q2.d). Comme le document l’exige, la variable stockée est le côté $c$, jamais son carré : la version qui stocke $c^2$ raconte une autre histoire, et n’affiche jamais le $6$ du titre.

from math import sqrt

def archimede(kmax):
    c = 1.0           # cote de l'hexagone, c_6 = 1
    n = 6
    for k in range(kmax + 1):
        A = n * c / 2             # demi-perimetre inscrit
        print(k, n, A)
        c = sqrt(2 - sqrt(4 - c*c))   # relation de Q2.d), puis racine carree
        n = 2 * n

archimede(28)

Sortie :

0 6 3.0
1 12 3.1058285412302498
2 24 3.132628613281237
3 48 3.139350203046872
4 96 3.14103195089053
5 192 3.1414524722853443
6 384 3.141557607911622
7 768 3.141583892148936
8 1536 3.1415904632367617
9 3072 3.1415921060430483
10 6144 3.1415925165881546
11 12288 3.1415926186407894
12 24576 3.1415926453212157
13 49152 3.1415926453212157
14 98304 3.1415926453212157
15 196608 3.1415926453212157
16 393216 3.141593669849427
17 786432 3.1415923038117377
18 1572864 3.1416086962248038
19 3145728 3.1415868396550413
20 6291456 3.1416742650217575
21 12582912 3.1416742650217575
22 25165824 3.1430727401700396
23 50331648 3.1598061649411346
24 100663296 3.181980515339464
25 201326592 3.3541019662496847
26 402653184 4.242640687119286
27 805306368 6.0
28 1610612736 0.0

Tout y est : la belle convergence des premiers rangs, le gel des rangs $12$ à $15$, la dérive, puis l’effondrement. Au rang $26$, la machine affiche $4{,}2426\ldots$, c’est-à-dire $3\sqrt2$ : le côté calculé est devenu si faux que le polygone n’a plus rien à voir avec le cercle. Au rang $27$, elle annonce $6$. Au rang $28$, plus rien. Cette exécution reproduit le tableau de Q15 au dernier chiffre près, sur les dix lignes affichées par le DM et sur la colonne des écarts : vous pouvez la citer les yeux fermés.

2. La réparation par quantité conjuguée (Q18 et Q19)

Le squelette commun de l’annexe, complété dans ses deux branches. Une seule ligne change entre les deux versions, et cette ligne change tout.

from math import sqrt, pi

def table(kmax, version):
    c = 1.0                      # cote de l'hexagone
    n = 6
    for k in range(kmax + 1):
        A = n * c / 2
        print(k, n, A, abs(A - pi))
        if version == "naif":
            c = sqrt(2 - sqrt(4 - c*c))      # relation de Q2.d), voir Q14
        else:
            c = c / sqrt(2 + sqrt(4 - c*c))  # relation de Q18.a), voir Q19
        n = 2 * n

table(28, "stable")

Sortie :

0 6 3.0 0.14159265358979312
1 12 3.105828541230249 0.03576411235954424
2 24 3.1326286132812378 0.008964040308555354
3 48 3.1393502030468667 0.0022424505429263775
4 96 3.1410319508905093 0.000560702699283766
5 192 3.1414524722854615 0.00014018130433157694
6 384 3.1415576079118575 3.504567793566338e-05
7 768 3.1415838921483177 8.76144147543556e-06
8 1536 3.1415904632280496 2.1903617435370393e-06
9 3072 3.141592105999271 5.475905222596111e-07
10 6144 3.1415925166921563 1.3689763678215172e-07
11 12288 3.141592619365383 3.422440997269405e-08
12 24576 3.1415926450336897 8.556103381351932e-09
13 49152 3.1415926514507664 2.1390267335164026e-09
14 98304 3.141592653055036 5.347571274683105e-10
15 196608 3.141592653456103 1.336899480008924e-10
16 393216 3.14159265355637 3.342304211173541e-11
17 786432 3.141592653581437 8.355982572538778e-12
18 1572864 3.1415926535877032 2.0898838215543947e-12
19 3145728 3.14159265358927 5.231370892033738e-13
20 6291456 3.1415926535896617 1.3145040611561853e-13
21 12582912 3.1415926535897594 3.375077994860476e-14
22 25165824 3.1415926535897842 8.881784197001252e-15
23 50331648 3.1415926535897905 2.6645352591003757e-15
24 100663296 3.1415926535897922 8.881784197001252e-16
25 201326592 3.141592653589793 0.0
26 402653184 3.141592653589793 0.0
27 805306368 3.141592653589793 0.0
28 1610612736 3.141592653589793 0.0

Même suite mathématique, même machine, mêmes seize chiffres de précision : mais plus aucune soustraction de deux nombres voisins, donc plus d’annulation catastrophique. À partir du rang $25$, l’écart affiché est exactement nul : la machine est arrivée au flottant le plus proche de $\pi$, et elle n’en bougera plus. La suite exacte, elle, continue de croître, comme le montre la table de la section 5.

3. Les moyennes de Q9, et la moyenne pondérée de Q27

Le second programme de l’annexe : les deux demi-périmètres se calculent l’un à partir de l’autre par moyenne harmonique et moyenne géométrique, sans aucune soustraction dangereuse. C’est lui qui fournit les valeurs en pleine précision demandées par Q13, Q26 et Q27.

from math import sqrt, pi

def moyennes(kmax):
    A, B = 3.0, 2*sqrt(3)
    for k in range(kmax + 1):
        S = (2*A + B)/3          # la moyenne ponderee de Q27
        print(k, A, B, B - A, S, abs(S - pi))
        B = 2*A*B/(A + B)        # moyenne harmonique, Q9.a)
        A = sqrt(A*B)            # moyenne geometrique, Q9.b)

moyennes(8)

Sortie :

0 3.0 3.4641016151377544 0.4641016151377544 3.154700538379251 0.013107884789457902
1 3.1058285412302493 3.215390309173473 0.10956176794322348 3.142349130544657 0.0007564769548640271
2 3.132628613281238 3.1596599420975005 0.027031328816262246 3.1416390562199923 4.640263019917157e-05
3 3.139350203046867 3.146086215131435 0.006736012084567644 3.14159554040839 2.886818597058749e-06
4 3.1410319508905093 3.1427145996453683 0.0016826487548589064 3.1415928338087955 1.802190023880712e-07
5 3.141452472285462 3.141873049979824 0.00042057769436221193 3.1415926648502492 1.1260456123096674e-08
6 3.1415576079118575 3.1416627470568486 0.00010513914499110655 3.1415926542935213 7.03728186834951e-10
7 3.141583892148318 3.1416101766046896 2.6284456371428178e-05 3.1415926536337753 4.398215125434035e-11
8 3.14159046322805 3.141597034321526 6.571093476015477e-06 3.1415926535925416 2.7484681197620375e-12

Notez l’avantage signalé par Q20.c) : à chaque rang, on dispose des deux bornes, donc d’un encadrement certifié de $\pi$, ce que la version de la section 2 ne donne pas. Et la colonne de $S$ illustre l’accélération de Snell : comparez sa dernière colonne à la colonne $B - A$ du même rang.

4. La suite de Muller (partie IV et Q31)

Le programme de l’annexe, tel quel : trois lignes de calcul, et deux arithmétiques au choix.

from fractions import Fraction

def muller(nmax, exact):
    if exact:
        u, v = Fraction(2), Fraction(-4)
    else:
        u, v = 2.0, -4.0
    for n in range(1, nmax):
        u, v = v, 111 - 1130/v + 3000/(u*v)
        print(n + 1, float(v))

Sortie de muller(30, False), en flottants :

2 18.5
3 9.378378378378379
4 7.801152737752169
5 7.154414480975333
6 6.806784736924811
7 6.592632768721792
8 6.449465934053933
9 6.348452060746624
10 6.274438662728116
11 6.218696768582163
12 6.17585385581539
13 6.142627170481006
14 6.120248704570159
15 6.166086559598099
16 7.235021165534931
17 22.062078463525793
18 78.57557488787224
19 98.34950312216536
20 99.8985692661829
21 99.99387098890278
22 99.99963038728635
23 99.99997773067949
24 99.99999865921669
25 99.99999991932181
26 99.99999999514776
27 99.99999999970828
28 99.99999999998246
29 99.99999999999893
30 99.99999999999993

Sortie de muller(30, True), en fractions exactes (extrait : début et fin) :

2 18.5
3 9.378378378378379
4 7.801152737752162
5 7.154414480975249
6 6.806784736923633
7 6.592632768704439
8 6.449465933790288
9 6.348452056654357
10 6.274438598216328
...
25 6.014174914550819
26 6.011784587871333
27 6.0098012392984845
28 6.008154378912229
29 6.006786093031206
30 6.00564868877142

Jusqu’au rang $13$, les deux colonnes se ressemblent. Au rang $14$, la version flottante commence à remonter ; au rang $17$, elle s’enfuit ; à partir du rang $22$, elle s’installe en $100$ avec une régularité parfaite, sans le moindre signe d’alarme. La version exacte descend tranquillement vers $6$. Les deux sorties reproduisent chiffre pour chiffre le tableau de Q24. Quant au prix payé par l’arithmétique exacte, c’est l’objet de la question Q31 : exécutez le mode exact en affichant la fraction $u_{30}$ elle-même plutôt que sa valeur décimale, la réponse vaut le détour.

5. La table des trente doublements

Le DM l’a promise, la voici : la suite $(A_k)$ des demi-périmètres inscrits, du rang $0$ (l’hexagone) au rang $29$, calculée en arithmétique décimale à $60$ chiffres, hors d’atteinte de l’annulation catastrophique. Le script est le suivant : $\pi$ y est obtenu par la formule de Machin (vérifié sur les $59$ premiers chiffres publiés par ailleurs, le soixantième portant l’arrondi du calcul), la récurrence est la forme stable de Q18.a), et la table a été recalculée une seconde fois, indépendamment, par la récurrence des moyennes de Q9, avec un accord à $10^{-59}$ près, la pleine précision du calcul.

from decimal import Decimal, getcontext, ROUND_HALF_EVEN

getcontext().prec = 60

def pi_machin():
    """pi = 16 arctan(1/5) - 4 arctan(1/239), serie de Taylor en Decimal."""
    def arctan_inv(x):
        x = Decimal(x)
        terme, total, k, signe = 1/x, Decimal(0), 0, 1
        while terme != 0:
            total += signe * terme / (2*k + 1)
            terme, k, signe = terme/(x*x), k + 1, -signe
        return total
    return 16*arctan_inv(5) - 4*arctan_inv(239)

PI = pi_machin()

c = Decimal(1)                  # cote de l'hexagone, exact
n = 6
q = Decimal("1.000000000000000")     # gabarit : 15 decimales
for k in range(30):
    A = n * c / 2
    print(k, n, A.quantize(q, rounding=ROUND_HALF_EVEN), "{:.1e}".format(float(PI - A)))
    c = c / (2 + (4 - c*c).sqrt()).sqrt()    # relation stable de Q18.a)
    n = 2 * n

La sortie, reformatée, donne la table suivante. Les valeurs sont arrondies à quinze décimales ; l’écart $\pi - A_k$ est donné avec deux chiffres significatifs.

$k$	$n = 6 \times 2^k$	demi-périmètre inscrit $A_k$	écart $\pi - A_k$
$0$	$6$	$3{,}000000000000000$	$1{,}4 \times 10^{-1}$
$1$	$12$	$3{,}105828541230249$	$3{,}6 \times 10^{-2}$
$2$	$24$	$3{,}132628613281238$	$9{,}0 \times 10^{-3}$
$3$	$48$	$3{,}139350203046867$	$2{,}2 \times 10^{-3}$
$4$	$96$	$3{,}141031950890510$	$5{,}6 \times 10^{-4}$
$5$	$192$	$3{,}141452472285462$	$1{,}4 \times 10^{-4}$
$6$	$384$	$3{,}141557607911858$	$3{,}5 \times 10^{-5}$
$7$	$768$	$3{,}141583892148318$	$8{,}8 \times 10^{-6}$
$8$	$1\,536$	$3{,}141590463228050$	$2{,}2 \times 10^{-6}$
$9$	$3\,072$	$3{,}141592105999272$	$5{,}5 \times 10^{-7}$
$10$	$6\,144$	$3{,}141592516692157$	$1{,}4 \times 10^{-7}$
$11$	$12\,288$	$3{,}141592619365384$	$3{,}4 \times 10^{-8}$
$12$	$24\,576$	$3{,}141592645033691$	$8{,}6 \times 10^{-9}$
$13$	$49\,152$	$3{,}141592651450768$	$2{,}1 \times 10^{-9}$
$14$	$98\,304$	$3{,}141592653055037$	$5{,}3 \times 10^{-10}$
$15$	$196\,608$	$3{,}141592653456104$	$1{,}3 \times 10^{-10}$
$16$	$393\,216$	$3{,}141592653556371$	$3{,}3 \times 10^{-11}$
$17$	$786\,432$	$3{,}141592653581438$	$8{,}4 \times 10^{-12}$
$18$	$1\,572\,864$	$3{,}141592653587704$	$2{,}1 \times 10^{-12}$
$19$	$3\,145\,728$	$3{,}141592653589271$	$5{,}2 \times 10^{-13}$
$20$	$6\,291\,456$	$3{,}141592653589663$	$1{,}3 \times 10^{-13}$
$21$	$12\,582\,912$	$3{,}141592653589761$	$3{,}3 \times 10^{-14}$
$22$	$25\,165\,824$	$3{,}141592653589785$	$8{,}2 \times 10^{-15}$
$23$	$50\,331\,648$	$3{,}141592653589791$	$2{,}0 \times 10^{-15}$
$24$	$100\,663\,296$	$3{,}141592653589793$	$5{,}1 \times 10^{-16}$
$25$	$201\,326\,592$	$3{,}141592653589793$	$1{,}3 \times 10^{-16}$
$26$	$402\,653\,184$	$3{,}141592653589793$	$3{,}2 \times 10^{-17}$
$27$	$805\,306\,368$	$3{,}141592653589793$	$8{,}0 \times 10^{-18}$
$28$	$1\,610\,612\,736$	$3{,}141592653589793$	$2{,}0 \times 10^{-18}$
$29$	$3\,221\,225\,472$	$3{,}141592653589793$	$5{,}0 \times 10^{-19}$

Trois lectures de cette table, pour finir.

L’écart est divisé par $4$ à chaque ligne. C’est la vitesse mesurée en Q13, visible ici sur vingt-neuf lignes d’affilée. À partir du rang $24$, les quinze décimales affichées sont celles de $\pi$ et ne bougent plus ; l’écart, lui, continue sa descente géométrique, ligne après ligne.

Cette table n’est pas celle de Q15, et c’est tout le sujet. Le tableau de Q15 montre ce que la machine affiche en double précision ; cette table montre ce que la suite vaut. Confrontez-les : au rang $10$, la machine affichait $3{,}141592516588155$ quand la vraie valeur est $3{,}141592516692157$, divergence dès la dixième décimale, déjà causée par l’annulation catastrophique. Aux rangs $12$ à $15$, la machine restait figée sur $3{,}141592645321216$ pendant que la suite exacte, strictement croissante comme l’a démontré Q10, continuait d’avancer. Au-delà, plus rien ne coïncide. En revanche, le programme naïf de la section 1, réexécuté pour cette page, reproduit le tableau de Q15 chiffre pour chiffre : le DM disait vrai sur ce que dit la machine, et la machine avait tort sur ce que vaut la suite.

La borne de van Ceulen est au bout de cette logique. La dernière ligne enferme $\pi$ à $5 \times 10^{-19}$ près avec trois milliards de côtés. Trente-cinq décimales demandent beaucoup plus, et c’est l’objet de Q13.c), de Q29, et de l’histoire racontée par la pierre de Leyde.