/*************************************************************************************************************************************/
/* */
/* C O M P R E H E N S I O N D U P L U S P E T I T P R O G R A M M E D E C A L C U L D E P I : */
/* ( ' P O U R L A S C I E N C E ' , M A I 1 9 9 4 ) */
/* */
/* */
/* Author of '$xtc/pi_mini.04$c' : */
/* */
/* Unknown and Jean-Francois COLONNA (LACTAMME, AAAAMMJJhhmmss). */
/* */
/*************************************************************************************************************************************/
extern double pow();
#define PRECISION \
(2) \
/* Le nombre de chiffre par paquet (4) peut etre reduit, mais pas augmente... */
#define DECIMALES \
(1*PRECISION) \
/* Nombre de decimales calculees (ce parametre peut etre change sans probleme, mais en */ \
/* restant un multiple de 'PRECISION'...). ATTENTION, la liste des decimales inclut la */ \
/* partie entiere (3). */
#define RAPPORT_TERMES_DECIMALES \
(14) \
/* Constante representant le "rapport" entre le nombre de termes calcules et le nombre */ \
/* de decimales alors obtenues ; elle peut donc etre augmentee, mais pas reduite. Elle */ \
/* donne le nombre de termes d'une certaine serie qui semblent etre calcules "en parallele". */ \
/* Chacun d'eux est calcule pour une precision donnee (la boucle 'while' externe) dans un */ \
/* element du tableau 'termes' ; par exemple on aura : */ \
/* */ \
/* termes[14] <-O-> 1/27 = 1/(2x14-1) */ \
/* termes[13] <-O-> 1/25 = 1/(2x13-1) */ \
/* (...) */ \
/* termes[02] <-O-> 1/03 = 1/(2x02-1) */ \
/* termes[01] <-O-> 1/01 = 1/(2x01-1) */ \
/* */ \
/* dans la boucle 'while' interne, auquel cas, cela ressemble a la formule de Wallis, sauf */ \
/* qu'il faudrait que cela soit des carres, et d'autre part cette methode converge tres */ \
/* lentement alors que ce programme semble relativement rapide... */
#define INDICE_DE_GAUCHE \
0 \
/* Indice de depart d'un tableau. */
#define TAILLE \
(((DECIMALES * RAPPORT_TERMES_DECIMALES) / PRECISION) + 1) \
/* Taille du tableau 'termes'. */
#define P \
(1*precision)
#define DP \
(2*precision)
/* Precisions d'edition des entiers... */
#define EDIT(fonction) \
{ \
if (editer != 0) \
{ \
fonction; \
} \
else \
{ \
} \
} \
/* Fonction d'edition conditionnelle. */
main()
{
int editer=1;
/* Controleur d'edition (0=NON, 1=OUI). */
int base=(int)pow((double)10,(double)PRECISION);
int precision = PRECISION;
/* Nombre de chiffre par paquet... */
int paquet=0;
/* Paquet courant de 'precision' decimales. */
int termes[TAILLE];
/* Tableau contenant les differents termes de la serie (et non pas la suite des decimales) */
/* qui sont calcules d'une part dans l'ordre inverse (N, N-1, ..., 2, 1) dans la boucle */
/* 'while' interne, et d'autre part a un certain niveau de precision grace a la boucle */
/* 'while' externe. */
int indice_de_gauche=INDICE_DE_GAUCHE;
/* Indice d'indexation du tableau 'termes'. */
int indice_de_droite=TAILLE-1;
int retenue;
/* Pseudo-retenue propagee de la droite vers la gauche... */
int inverse_du_coefficient;
/* Il s'agit de l'inverse du coefficient '1/(2n+1)', mais malheureusement, comme il n'y */
/* nulle part de soustractions il ne doit pas s'agir du developpement en serie des */
/* fonctions arc-tangentes puisque cette serie est "alternee"... */
EDIT(printf("\n initialisation du processus\n"));
for (indice_de_gauche=INDICE_DE_GAUCHE ; indice_de_gauche < TAILLE ; indice_de_gauche++)
{
termes[indice_de_gauche] = (int)base*(2.0/10.0);
EDIT(printf(" %0*d",P,termes[indice_de_gauche]));
/* L'initialisation est la meme partout (2 a 'base' pres) car en fait le tableau 'termes' */
/* ne correspond pas aux decimales mais aux differents termes de la serie pour une precision */
/* donnee (et fixee par la boucle 'while' externe) ; voir les commentaires relatifs a la */
/* constante 'RAPPORT_TERMES_DECIMALES'. La division par 10 est destinee a mettre la partie */
/* entiere ('3') comme premiere decimale, c'est-a-dire a calculer pi/10. */
}
while ((inverse_du_coefficient=indice_de_droite*2) != 0)
{
int repeter=1;
retenue = 0;
indice_de_gauche = indice_de_droite;
EDIT(printf("\n boucle exterieure, indice de droite=%0*d",P,indice_de_droite));
while (repeter == 1)
/* Exemple d'imbrication en demandant 2*4 decimales : */
/* */
/* 14 28 */
/* ________________ */
/* \ _______ | */
/* \\ 5926 | 3141 | */
/* \\ | | */
/* \\ | | */
/* \\ | | G = indice_de_gauche */
/* \\ | | D = indice_de_droite */
/* \\ | | */
/* \\| | */
/* \ | */
/* \ | */
/* \ | */
/* G --> \ | <-- D */
/* \ | */
/* \ | */
/* \ | */
/* \| */
/* */
/* le grand triangle correspondant a la premiere boucle 'while' externe, alors que le */
/* petit triangle correspond a la seconde (et derniere). */
{
int Sretenue;
inverse_du_coefficient = inverse_du_coefficient - 2;
EDIT(printf("\n boucle interieure, indice de gauche=%0*d",P,indice_de_gauche));
EDIT(printf(" diviseur=%0*d",P,inverse_du_coefficient+1));
EDIT(printf(" retenue=%0*d",DP,retenue));
retenue = retenue + (termes[indice_de_gauche]*base);
EDIT(printf("+(%0*d*%0*d) --> %0*d"
,DP,termes[indice_de_gauche]
,DP,base
,DP,retenue
)
);
/* Cette multiplication par 'base' permet de cumuler les differents termes successifs de */
/* la serie utilisee, ainsi a la fin de chaque boucle 'while' interne, le dernier element */
/* calcule a gauche est le paquet courant de decimales. En ce qui concerne la nature de */
/* la serie utilisee, voir le programme 'v $xtc/pi_mini.02$c'. */
EDIT(Sretenue = retenue);
termes[indice_de_gauche] = retenue%(inverse_du_coefficient+1);
retenue = retenue/(inverse_du_coefficient+1);
EDIT(printf("/%0*d --> %0*d"
,DP,inverse_du_coefficient+1
,DP,retenue
)
);
indice_de_gauche = indice_de_gauche - 1;
if (indice_de_gauche == INDICE_DE_GAUCHE)
{
repeter = 0;
EDIT(printf(" %*s %*s"
,DP,""
,DP,""
)
);
}
else
{
retenue = retenue * indice_de_gauche;
EDIT(printf("*%0*d --> %0*d"
,DP,indice_de_gauche
,DP,retenue
)
);
}
EDIT(printf(" termes[%0*d]=%0*d%%%0*d=%0*d"
,P,indice_de_gauche+1
,DP,Sretenue
,DP,inverse_du_coefficient+1
,DP,termes[indice_de_gauche+1]
)
);
}
indice_de_droite = indice_de_droite - RAPPORT_TERMES_DECIMALES;
EDIT(printf("\n paquet courant=%0*d+(%0*d/%0*d)=",DP,paquet,DP,retenue,DP,base));
printf("%.*d",precision,paquet+(retenue/base));
paquet = retenue % base;
EDIT(printf("\n paquet=%0*d%%%0*d=%0*d",DP,retenue,DP,base,DP,paquet));
}
EDIT(printf("\n"));
}