Winter: Pi avec un million de chiffres

Télécharger le projet
Commenter

Description

Voici le célèbre code minimal de 158 caractères écrit en C par Dik Winter, permettant de calculer et d'afficher les 2400 premiers chiffres décimaux de Pi:

int a=10000,b,c=8400,d,e,f[8401],g;
main(){for(;b-c;)f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}

Il est basé sur la formule d'Euler (voir: Formule de Euler) qui nécessite environ 3½ itérations par chiffre.
Pi/2 = 1 + 1/3 + 2/(3*5) + 2*3/(3*5*7) + 2*3*4/(3*5*7*9) + 2*3*4*5/(3*5*7*9*11) + ...

Le programme du fichier Winter0.cpp vérifie ce code: 
/* Code de Dik Winter (158 caractères):
int a=10000,b,c=8400,d,e,f[8401],g;main(){for(;b-c;)f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}
*/

#include <stdio.h>

int a=10000,b=0,c=8400,d,e=0,f[8401],g; // ajouté =0 pour a et e.
void main(){
  for(;b-c;)f[b++]=a/5;
  for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
    for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);

  getchar(); // ajouté
}

De manière plus académique, on peut écrire (fichier Winter1.cpp): 
#include <stdio.h>

int i, j, c, r, t, a[8401];

void main() {
  for (i=0; i<=8400; i++) a[i] = 2000;
  for (c=0, r=0, j=8400/14; j>0; j--) {
    for (i=j*14; i>0; i--) {
      t = (r*i + a[i]*10000);
      r = t / (2*i - 1);
      a[i] = t % (2*i - 1);
    }
    printf("%.4d", c + r/10000);
    c = r % 10000;
  }
  getchar();
}

Question: doit-on vraiment se limiter à 2400 décimales ?

Essayons d'abord de préciser la "signification" des constantes utilisées et changeons tous les int en unsigned int (fichier Winter2.cpp): 
#include <stdio.h>
#include <time.h>

typedef unsigned int uint;

const uint NbrChf = 32000;         // Nombre de chiffres
const uint ValBlk = 10000;         // Block de 4 chiffres
const uint NbrBlk = NbrChf/4;      // Nombre de blocks
const uint IteBlk = 14;            // Nombre d'itérations par block
const uint IteTot = NbrBlk*IteBlk; // Nombre d'itérations total

uint i, j, c, r, t, a[IteTot+1];

void main() {
  clock_t s=clock()-500;

  for (i=0; i<=IteTot; i++) a[i] = ValBlk/5;
  for (c=0, r=0, j=NbrBlk; j>0; j--) {
    for (i=j*IteBlk; i>0; i--) {
      t = (r*i + a[i]*ValBlk);
      r = t / (2*i - 1);
      a[i] = t % (2*i - 1);
    }
    printf("%.4d", c + r/ValBlk);
    c = r % ValBlk;
  }
  printf("nWinter2: 32 bits, %u chiffres, %u s", NbrChf, (clock()-s)/1000);
  getchar();
}
En posant NbrChf = 2400, on obtient le même résultat que ci-dessus.
Avec NbrChf = 32000, tout semble correct.
Par contre, ce n'est plus le cas pour NbrChf = 33000 !
Faites un contrôle partiel avec par exemple: Décimales de pi, e, phi, etc.

Jusqu'à maintenant, les blocks étaient de 4 chiffres pour que leur produit ne dépasse pas 32 bits.

Adaptons ce code aux 64 bits

Les ajustements principaux sont:
- Remplacer les unsigned int par unsigned long long.
- Block de 8 chiffres (ValBlk = 100000000).
- 28 itérations par block (IteBlk = 28).
- Adapter le format d'écriture. 
#include <stdio.h>
#include <time.h>

typedef unsigned long long ull;

const ull NbrChf = 1000000;       // Nombre de chiffres
const ull ValBlk = 100000000;     // Block de 8 chiffres
const ull NbrBlk = NbrChf/8;      // Nombre de blocks
const ull IteBlk = 28;            // Nombre d'itérations par block
const ull IteTot = NbrBlk*IteBlk; // Nombre d'itérations total

ull i, j, c, r, t, a[IteTot+1];

void main() {
  clock_t s=clock()-500;

  for (i=0; i<=IteTot; i++) a[i] = ValBlk/5;
  for (c=0, r=0, j=NbrBlk; j>0; j--) {
    for (i=j*IteBlk; i>0; i--) {
      t = (r*i + a[i]*ValBlk);
      r = t / (2*i - 1);
      a[i] = t % (2*i - 1);
    }
    printf("%.8d", c + r/ValBlk);
    c = r % ValBlk;
  }
  printf("nWinter3: 64 bits, %u chiffres, %u s", NbrChf, (clock()-s)/1000);
  getchar();
}

Le temps de calcul et d'affichage augmente avec le carré du nombre de chiffres NbrChf:
50'000: 8s, 100'000: 29 s, 200'000: 116 s, 400'000: 456 s, 1'000'000: 2819 s.

Revenons au code minimal de Winter.
Winter4: 32 bits, 32000 chiffres

unsigned int a=10000,b,c=112000,d,e,f[112001],g;
main(){for(;b-c;)f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}

Ce code a l'air de générer 32'000 chiffres corrects de Pi.

Winter5: 64 bits, 100000000 chiffres

unsigned long long a=100000000,b=0,c=3500000,d,e=0,f[3500001],g;
void main(){for(;b-c;)f[b++]=a/5;
for(;d=0,g=c*2;c-=28,printf("%.8d",e+d/a),e=d%a)
for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);

Les 1'000'000 chiffres de Pi calculés par ce programme semblent être corrects.

Attention: pour un million de chiffres, le temps d'exécution est d'environ une heure ! (55 minutes avec un processeur I5 3.2 GHz).

Bonne lecture.
Télécharger le projet

Codes Sources

S'inscrire maintenant

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.

Du même auteur (William VOIROL)