LE CALCUL DE PI

Tous d'abord, nous allons travailler à une précision de 1000 décimales.

dec := 1000

I. Les algorithmes géométriques

C'est tout d'abord l'algorithme d'Archimède, célèbre grec. Il a postulé que :

6*(2^n)*(Subscript[U, n])->π,  (Subscript[U, n])=\(TraditionalForm\`\@\(\(1\/2\) \((1 - \@\(1 - x²\))\)\)\)

archi1[x_] := N[Sqrt[1/2*(1 - Sqrt[1 - x^2])], dec]archi[n_] := 6*2^n*Last[NestList[archi1, 1/2, n]]

On va prendre le temps mis pour calculer à 10 itérations, 20, 50, 100, 200, 500, ainsi que la précision. On construit deux listes architps et archires pour pouvoir afficher des graphes à la fin.

r1 = Timing[archi[10]]a1 = %[[1]]b1 = %%[[2]]r2 = Timing[archi[20]]a2 = %[[1]]b2 = %%[[2]]r3 = Timing[archi[50]]a3 = %[[1]]b3 = %%[[2]]r4 = Timing[archi[100]]a4 = %[[1]]b4 = %%[[2]]r5 = Timing[archi[200]]a5 = %[[1]]b5 = %%[[2]]r6 = Timing[archi[500]]a6 = %[[1]]b6 = %%[[2]]

architps
       = {{10, a1 /. Second -> 1}, {20, a2 /. Second -> 1}, {50, 
          a3 /. Second -> 1}, {100, a4 /. Second -> 1}, {200, 
          a5 /. Second -> 1}, {500, a6 /. Second -> 1}}

archires
       = {{Log[10], Log[Abs[N[Pi, 1000] - b1]]}, {Log[20], 
          Log[Abs[N[Pi, 1000] - b2]]}, {Log[50], 
          Log[Abs[N[Pi, 1000] - b3]]}, {Log[100], 
          Log[Abs[N[Pi, 1000] - b4]]}, {Log[200], 
          Log[Abs[N[Pi, 1000] - b5]]}, {Log[500], 
          Log[Abs[N[Pi, 1000] - b6]]}};

ListPlot[archires, PlotRange -> All, PlotJoined -> True]

ListPlot[architps, PlotRange -> All, PlotJoined -> True]

Ensuite, c'est une application des suites de Fibonnaci, et on est déjà passé au Moyen-Age, voire à la Renaissance.

π=\(TraditionalForm
      \`4*\(∑\+\(k = 0\)\%∞\((\(-1\))\)\^k\/\(2  k + 1\)*
            \((\((\((1 + \@5)\)\/2)\)\^\(\(-2\) k - 1\) + 
                \((\((1 + \@5)\)\/2)\)\^\(\(-6\) k - 3\))\)\)\)

fibo1[k_] := 
      (-1)^k/(2  k + 1)*
        (((1 + Sqrt[5])/2)^(-2 k - 1) + 
            ((1 + Sqrt[5])/2)^(-6 k - 3))\(	fibo[n_] := N[4*Sum[fibo1[k], {k, 0, n}], dec]\)

r1 = Timing[fibo[10]]a1 = %[[1]]b1 = %%[[2]]r2 = Timing[fibo[20]]a2 = %[[1]]b2 = %%[[2]]r3 = Timing[fibo[50]]a3 = %[[1]]b3 = %%[[2]]r4 = Timing[fibo[100]]a4 = %[[1]]b4 = %%[[2]]r5 = Timing[fibo[200]]a5 = %[[1]]b5 = %%[[2]]r6 = Timing[fibo[500]]a6 = %[[1]]b6 = %%[[2]]

fibotps
       = {{10, a1 /. Second -> 1}, {20, a2 /. Second -> 1}, {50, 
          a3 /. Second -> 1}, {100, a4 /. Second -> 1}, {200, 
          a5 /. Second -> 1}, {500, a6 /. Second -> 1}}

fibores
       = {{Log[10], Log[Abs[N[Pi, 1000] - b1]]}, {Log[20], 
          Log[Abs[N[Pi, 1000] - b2]]}, {Log[50], 
          Log[Abs[N[Pi, 1000] - b3]]}, {Log[100], 
          Log[Abs[N[Pi, 1000] - b4]]}, {Log[200], 
          Log[Abs[N[Pi, 1000] - b5]]}, {Log[500], 
          Log[Abs[N[Pi, 1000] - b6]]}};

ListPlot[fibores, PlotRange -> All, PlotJoined -> True]

ListPlot[fibotps, PlotRange -> All, PlotJoined -> True]

II. Les algorithmes analytiques

On va s'intéresser à la formule de Gosper, qui indique que :

π=\(TraditionalForm
      \`\(9\/\@3\) 
        \(∑\+\(k = 1\)\%∞\(\((\(n!\))\) ²\)\/\(n*
              \(\((2  n)\)!\)\)\)\)

gosper2[n_] := (n!)^2/(n*(2  n)!)gosper[n_] := N[9/Sqrt[3]*Sum[gosper2[k], {k, 1, n}], dec]

r1 = Timing[gosper[10]]a1 = %[[1]]b1 = %%[[2]]r2 = Timing[gosper[20]]a2 = %[[1]]b2 = %%[[2]]r3 = Timing[gosper[50]]a3 = %[[1]]b3 = %%[[2]]r4 = Timing[gosper[100]]a4 = %[[1]]b4 = %%[[2]]r5 = Timing[gosper[200]]a5 = %[[1]]b5 = %%[[2]]r6 = Timing[gosper[500]]a6 = %[[1]]b6 = %%[[2]]

gospertps
       = {{10, a1 /. Second -> 1}, {20, a2 /. Second -> 1}, {50, 
          a3 /. Second -> 1}, {100, a4 /. Second -> 1}, {200, 
          a5 /. Second -> 1}, {500, a6 /. Second -> 1}}

gosperres
       = {{Log[10], Log[Abs[N[Pi, 1000] - b1]]}, {Log[20], 
          Log[Abs[N[Pi, 1000] - b2]]}, {Log[50], 
          Log[Abs[N[Pi, 1000] - b3]]}, {Log[100], 
          Log[Abs[N[Pi, 1000] - b4]]}, {Log[200], 
          Log[Abs[N[Pi, 1000] - b5]]}, {Log[500], 
          Log[Abs[N[Pi, 1000] - b6]]}};

ListPlot[gosperres, PlotRange -> All, PlotJoined -> True]

ListPlot[gospertps, PlotRange -> All, PlotJoined -> True]

Ensuite, c'est la fameuse formule de Leibniz, démontrée à partir du développement d'arctan.

π=8*\(TraditionalForm
      \`∑\+\(k = 0\)\%∞ 1
          \/\(\((4  k + 1)\) \((4  k + 3)\)\)\)

leibniz1[n_] := 1/((4*n + 1)*(4*n + 3))leibniz[n_] := N[8*Sum[leibniz1[k], {k, 0, n}], dec]

r1 = Timing[leibniz[10]]a1 = %[[1]]b1 = %%[[2]]r2 = Timing[leibniz[20]]a2 = %[[1]]b2 = %%[[2]]r3 = Timing[leibniz[50]]a3 = %[[1]]b3 = %%[[2]]r4 = Timing[leibniz[100]]a4 = %[[1]]b4 = %%[[2]]r5 = Timing[leibniz[200]]a5 = %[[1]]b5 = %%[[2]]r6 = Timing[leibniz[500]]a6 = %[[1]]b6 = %%[[2]]

leibniztps
       = {{10, a1 /. Second -> 1}, {20, a2 /. Second -> 1}, {50, 
          a3 /. Second -> 1}, {100, a4 /. Second -> 1}, {200, 
          a5 /. Second -> 1}, {500, a6 /. Second -> 1}}

leibnizres
       = {{Log[10], Log[Abs[N[Pi, 1000] - b1]]}, {Log[20], 
          Log[Abs[N[Pi, 1000] - b2]]}, {Log[50], 
          Log[Abs[N[Pi, 1000] - b3]]}, {Log[100], 
          Log[Abs[N[Pi, 1000] - b4]]}, {Log[200], 
          Log[Abs[N[Pi, 1000] - b5]]}, {Log[500], 
          Log[Abs[N[Pi, 1000] - b6]]}};

ListPlot[leibnizres, PlotRange -> All, PlotJoined -> True]

ListPlot[leibniztps, PlotRange -> All, PlotJoined -> True]

III. Les Algorithmes Modernes

C'est tout d'abord l'algorithe de Ramanjuan :

\(TraditionalForm
    \`1\/π = 
      ∑\+\(k = 0\)\%∞\(\(\((\(\((2  n)\)!\))\)\^3\) 
            \((42  n + 5)\)\)\/\(\(2\^\(12  n + 4\)\) \((\(n!\))\)\^6\)\)

raman1[n_] := 
      ((2  n)!)^3*
        (42*n + 5)/(2^(12  n + 4)*(n!)^6)raman2[n_] := N[Sum[raman1[k], {k, 0, n}], dec]raman[n_] := N[1/raman2[n], dec]

r1 = Timing[raman[10]]a1 = %[[1]]b1 = %%[[2]]r2 = Timing[raman[20]]a2 = %[[1]]b2 = %%[[2]]r3 = Timing[raman[50]]a3 = %[[1]]b3 = %%[[2]]r4 = Timing[raman[100]]a4 = %[[1]]b4 = %%[[2]]r5 = Timing[raman[200]]a5 = %[[1]]b5 = %%[[2]]r6 = Timing[raman[500]]a6 = %[[1]]b6 = %%[[2]]

ramantps
       = {{10, a1 /. Second -> 1}, {20, a2 /. Second -> 1}, {50, 
          a3 /. Second -> 1}, {100, a4 /. Second -> 1}, {200, 
          a5 /. Second -> 1}, {500, a6 /. Second -> 1}}

ramanres
       = {{Log[10], Log[Abs[N[Pi, 1000] - b1]]}, {Log[20], 
          Log[Abs[N[Pi, 1000] - b2]]}, {Log[50], 
          Log[Abs[N[Pi, 1000] - b3]]}, {Log[100], 
          Log[Abs[N[Pi, 1000] - b4]]}, {Log[200], 
          Log[Abs[N[Pi, 1000] - b5]]}, {Log[500], 
          Log[Abs[N[Pi, 1000] - b6]]}};

ListPlot[ramanres, PlotRange -> All, PlotJoined -> True]

ListPlot[ramantps, PlotRange -> All, PlotJoined -> True]

Et enfin, un algorithme à convergence ultra rapide, celui de Borwein. Dès la dixième itération, il donne Pi à 1000 décimales près, mais le calcul est très long.

\(TraditionalForm
      \`a\_n = \@a\_\(n - 1\)/2 + 1\/\(2 \@ a\_\(n - 1\)\)\),(Subscript[a, 0]) = (Sqrt[2])
\(TraditionalForm
      \`b\_n = b\_\(n - 1\)*
          \(1 + a\_\(n - 1\)\)\/\(a\_\(n - 1\) + b\_\(n - 1\)\)\), (Subscript[b, 0]) = 0
(Subscript[p, n]) = (Subscript[p, n - 1])*\(TraditionalForm
      \`\(b\_n\) \(1 + a\_n\)\/\(1 + b\_n\) -> π\), (Subscript[p, 0]) = 2 + (Sqrt[2])

a[0] := Sqrt[2]b[0] := 0p[0] := 2 + Sqrt[2]aa1[x_] := N[Sqrt[x]/2 + 1/(2*Sqrt[x]), dec]bb1[x_, y_] := N[Sqrt[y]*(1 + x)/(x + y), dec]aa[n_] := Last[NestList[aa1, Sqrt[2], n]]bb[n_] := Last[FoldList[bb1, 0, NestList[aa1, Sqrt[2], n - 1]]]c[n_] := Table[N[bb[i]*(1 + aa[i]), dec]/(1 + bb[i]), {i, n}]p1[x_, y_] := N[x*y, dec]borwein[n_] := Last[FoldList[p1, 2 + Sqrt[2], c[n]]]

r1 = Timing[borwein[10]]a1 = %[[1]]b1 = %%[[2]]

Abs[b1 - N[Pi, 1000]]

On va maintenant comparer la précision des différents algorithmes après 500 itérations.

prec
       = {{1, archires[[6]][[2]]}, {2, 
          fibores[[6]][[2]]}, {3, 
          gosperres[[6]][[2]]}, {4, 
          leibnizres[[6]][[2]]}, {5, 
          ramanres[[6]][[2]]}}

ListPlot[prec,  PlotRange -> All,  PlotJoined -> True]

<< Graphics`MultipleListPlot`

MultipleListPlot[architps, fibotps, gospertps, leibniztps, ramantps, 
      PlotRange -> All, PlotJoined -> True, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0], RGBColor[0, 0, 1], 
          RGBColor[0, 0, 0], RGBColor[1, 0, 1]}, 
      PlotLegend -> {Archimede, Fibonnaci, Gosper, Leibniz, Ramanjuan}]