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.
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.