/*************************************************************************************************************************************/
/* */
/* C A L C U L D E P I : */
/* */
/* */
/* Nota Important : */
/* */
/* J'ai un peu honte de la qualite et du style */
/* de ce programme ; il convient donc que je m'en */
/* explique ici. Il s'agit la d'un programme ecrit */
/* en 1970 sur un ordinateur CII 10070 (et oui, */
/* je BULLais deja...) en FORTRAN ETENDU (d'ou */
/* les affectations imbriquees entre autres...) ; */
/* ces extensions par rapport a F77 faisaient que */
/* je ne pouvais pas l'executer directement sur */
/* le CII SPS 9 (a pardon, je voulais dire le */
/* TELEMECANIQUE SPS 9 (a pardon, je voulais dire le */
/* SEMS SPS 9 (a pardon, je voulais dire le BULL SPS 9 */
/* (a pardon... (((c'est recursif afin d'etre */
/* prevoyant quant au futur))) ...)))). */
/* Donc, pour m'eviter une reecriture fastidieuse */
/* (et oui, la BULLe...), j'ai decide de le transformer */
/* simplement en un programme C : pour cela, j'ai */
/* mis des ";" en bout de ligne, remplace les boucles */
/* "DO" par des boucles "FOR",... Et ca marche, mais */
/* bien sur au prix d'abominations, telles des "GOTO"s... */
/* J'espere que vous voudrez bien me pardonner ces */
/* offenses, et me donner mon SPS9 quotidien... */
/* */
/* */
/* Author of '$xrp/pi.F1$K' : */
/* */
/* Jean-Francois Colonna (LACTAMME, 19????????????). */
/* */
/*************************************************************************************************************************************/
#undef goto
/*************************************************************************************************************************************/
/* */
/* C O N S T A N T E S D E B A S E : */
/* */
/*************************************************************************************************************************************/
#define zero ((DOOBLE)ZERO)
#define memoire 200000 \
/* Taille des tableaux DOOBLE 'terme' et 'serie'. */
/*= #include "DEFINITIONS.I" */
/* */
/*= #ifdef inclure_DEFINIT_DEF */
/*= # include gener_DEFINIT_1_DEF */
/*= #Aifdef inclure_DEFINIT_DEF */
/*= #Eifdef inclure_DEFINIT_DEF */
/* */
/*= #include gener_DEFINIT_2_DEF */
/*= #define decimales 10 */
/*= double taille = W; */
/*= int index = UNDEF; */
/*= for (index = BEGIN_AT_0 ; index <= decimales ; index++) */
/*= { */
/*= taille = taille * BASE10; */
/*= } */
/*= DEFINED("decimales",decimales,"nombre de chiffres significatifs par mot."); */
/*= DEFINEF("taille",taille,"taille du paquet de decimales contenu dans un mot memoire."); */
#define mots 10 \
/* Nombre de mots (paquet de decimales) par ligne. */
#define nombred 10 \
/* Nombre de decimales a imprimer par paquet, */
#define nombrem 10 \
/* Nombre de paquets de decimales a imprimer par ligne. */
#define e10 (1.0 * taille)
#define c32 (3.2 * taille)
#define c2 2.0
#define c3 3.0
#define c6 (c2 * c3)
#define c10 decimales
#define c19 (c10 + c10 - BEGIN_AT_0)
#define c410 (4.0 * taille)
#define c25 (5.0 * 5.0) \
/* Pour calculer ARCTG(1/5). */
#define c15625 (c25 * c25 * c25)
#define c239 239.0 \
/* Pour calculer ARCTG(1/239). */
#define c57121 (c239 * c239)
#define DOOBLE double
/*===================================================================================================================================*/
main()
/*-----------------------------------------------------------------------------------------------------------------------------------*/
{
int mindex={UNDEF};
/* Nombre de paquets de decimales a calculer. */
int mindexm1={UNDEF};
int index0={UNDEF};
int ibande={UNDEF};
int index={UNDEF};
/* Index d'acces aux tableaux 'a' et 's'. */
DOOBLE termec={UNDEF};
DOOBLE reste1={UNDEF};
DOOBLE reste2={UNDEF};
DOOBLE reste3={UNDEF};
DOOBLE angle2={UNDEF};
DOOBLE mangle2={UNDEF};
DOOBLE angle6={UNDEF};
DOOBLE mangle6={UNDEF};
DOOBLE quot1={UNDEF};
DOOBLE quot2={UNDEF};
DOOBLE quot3={UNDEF};
DOOBLE cumul1={UNDEF};
DOOBLE cumul2={UNDEF};
DOOBLE cumul3={UNDEF};
DOOBLE cumul4={UNDEF};
DOOBLE reste4={UNDEF};
DOOBLE divis1={UNDEF};
DOOBLE divis2={UNDEF};
DOOBLE mdivis1={UNDEF};
DOOBLE mdivis2={UNDEF};
DOOBLE divis3={UNDEF};
DOOBLE mdivis3={UNDEF};
DOOBLE termem={UNDEF};
DOOBLE seriem={UNDEF};
DOOBLE reste={UNDEF};
DOOBLE capacite={e10};
DOOBLE mcapacite={-e10};
int calcul={UNDEF};
DOOBLE terme[memoire+BEGIN_AT_0];
DOOBLE serie[memoire+BEGIN_AT_0];
DOOBLE q={UNDEF};
int chiffre={UNDEF};
int index1={UNDEF};
int index2={UNDEF};
char ligne[mots][decimales];
static char chiffres[BASE10]={'0','1','2','3','4','5','6','7','8','9'};
/*..............................................................................................................................*/
/*************************************************************************************************************************************/
/* */
/* I N I T I A L I S A T I O N S : */
/* */
/*************************************************************************************************************************************/
cumul1 = (double)W;
for (index = BEGIN_AT_0 ; index <= decimales ; index++)
{
cumul1 = cumul1 * (DOOBLE)BASE10;
}
if (cumul1 != (DOOBLE)taille)
{
CAL1(Prer0(" ATTENTION : les parametres 'decimales' et 'taille' sont incoherents !!! \n"));
exit();
}
CAL2(Prin1("\n Veuillez entrer le nombre minimal de paquets de %d decimales desirees : ",decimales));
scanf("%d",&mindex);
if (mindex >= memoire)
{
CAL1(Prer0(" ATTENTION : trop de decimales sont demandees !!! \n"));
exit();
}
if (mindex < zero)
{
CAL1(Prer0("\n ATTENTION, la valeur absolue de votre reponse va etre utilisee..."));
mindex = -mindex;
}
mindexm1 = (mindex = mots + (mots * (mindex / mots))) - I;
/* Afin que le nombre de decimales soit un multiple du nombre de chiffres... */
calcul = mindex * decimales;
CAL2(Prin1("\n\n Calcul de PI avec %d decimales.",calcul));
Bclock("",VRAI);
serie[INDEX0] = terme[INDEX0] = c32;
for (index = INDEX0 + I ; index <= mindexm1 ; index++)
{
serie[index] = terme[index] = zero;
/* Initialisation des tableaux de cumul. */
}
seriem = termem = zero;
index0 = INDEX0 + I;
mdivis1 = c3;
mangle6 = -(angle6 = c15625);
mangle2 = -(angle2 = c25);
/*************************************************************************************************************************************/
/* */
/* C A L C U L D E L A S E R I E A R C T G ( 1 / 5 ) : */
/* */
/*************************************************************************************************************************************/
eti23:
if (terme[index0 - I] == zero) index0 = index0 + I;
if ((ibande = ((divis3 = - (mdivis3 = (mdivis2 = - (divis2 = (divis1 = - (mdivis1 = mdivis1 - c6)) + c2)) - c2)) + c19) / c10)
>= mindex) goto eti200;
reste1 = capacite * (mangle6 * (termec = (terme[index0 - I] = AINT((cumul4 = terme[index0 - I]) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT(termec / divis1)) + termec);
reste3 = capacite * (mdivis2 * (quot2 = AINT(termec / divis2)) + termec);
reste4 = capacite * (mdivis3 * (quot3 = AINT(termec / divis3)) + termec);
serie[index0 - I] = serie[index0 - I] + (mangle2 * quot1 + quot2) * angle2 - quot3;
for (index = index0 ; index <= ibande ; index++)
{
reste1 = capacite * (mangle6 * (terme[index] = (termec = AINT((cumul4 = terme[index] + reste1) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT((cumul1 = termec + reste2) / divis1)) + cumul1);
reste3 = capacite * (mdivis2 * (quot2 = AINT((cumul2 = termec + reste3) / divis2)) + cumul2);
reste4 = capacite * (mdivis3 * (quot3 = AINT((cumul3 = termec + reste4) / divis3)) + cumul3);
serie[index] = serie[index] + (mangle2 * quot1 + quot2) * angle2 - quot3;
}
for (index = ibande + I ; index <= mindexm1 ; index++)
{
reste2 = capacite * (mdivis1 * (quot1 = AINT(reste2 / divis1)) + reste2);
reste3 = capacite * (mdivis2 * (quot2 = AINT(reste3 / divis2)) + reste3);
reste4 = capacite * (mdivis3 * (quot3 = AINT(reste4 / divis3)) + reste4);
serie[index] = serie[index] + (mangle2 * quot1 + quot2) * angle2 - quot3;
}
quot2 = AINT(reste3 / divis2);
quot3 = AINT(reste4 / divis3);
seriem = seriem + (mangle2 * AINT(reste2 / divis1) + quot2) * angle2 - quot3;
if (terme[index0 - I] == zero) index0 = index0 + I;
if ((ibande = ((divis3 = - (mdivis3 = (mdivis2 = - (divis2 = (divis1 = - (mdivis1 = mdivis1 - c6)) + c2)) - c2)) + c19) / c10)
>= mindex) goto eti201;
reste1 = capacite * (mangle6 * (termec = (terme[index0 - I] = AINT((cumul4 = terme[index0 - I]) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT(termec / divis1)) + termec);
reste3 = capacite * (mdivis2 * (quot2 = AINT(termec / divis2)) + termec);
reste4 = capacite * (mdivis3 * (quot3 = AINT(termec / divis3)) + termec);
serie[index0 - I] = serie[index0 - I] + (mangle2 * quot1 + quot2) * mangle2 + quot3;
for (index = index0 ; index <= ibande ; index++)
{
reste1 = capacite * (mangle6 * (terme[index] = (termec = AINT((cumul4 = terme[index] + reste1) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT((cumul1 = termec + reste2) / divis1)) + cumul1);
reste3 = capacite * (mdivis2 * (quot2 = AINT((cumul2 = termec + reste3) / divis2)) + cumul2);
reste4 = capacite * (mdivis3 * (quot3 = AINT((cumul3 = termec + reste4) / divis3)) + cumul3);
serie[index] = serie[index] + (mangle2 * quot1 + quot2) * mangle2 + quot3;
}
for (index = ibande + I ; index <= mindexm1 ; index++)
{
reste2 = capacite * (mdivis1 * (quot1 = AINT(reste2 / divis1)) + reste2);
reste3 = capacite * (mdivis2 * (quot2 = AINT(reste3 / divis2)) + reste3);
reste4 = capacite * (mdivis3 * (quot3 = AINT(reste4 / divis3)) + reste4);
serie[index] = serie[index] + (mangle2 * quot1 + quot2) * mangle2 + quot3;
}
quot2 = AINT(reste3 / divis2);
quot3 = AINT(reste4 / divis3);
seriem = seriem + (mangle2 * AINT(reste2 / divis1) + quot2) * mangle2 + quot3;
goto eti23;
eti200:
reste1 = capacite * (mangle6 * (termec = (terme[index0 - I] = AINT((cumul4 = terme[index0 - I]) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT(termec / divis1)) + termec);
reste3 = capacite * (mdivis2 * (quot2 = AINT(termec / divis2)) + termec);
reste4 = capacite * (mdivis3 * (quot3 = AINT(termec / divis3)) + termec);
serie[index0 - I] = serie[index0 - I] + (mangle2 * quot1 + quot2) * angle2 - quot3;
for (index = index0 ; index <= mindexm1 ; index++)
{
reste1 = capacite * (mangle6 * (terme[index] = (termec = AINT((cumul4 = terme[index] + reste1) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT((cumul1 = termec + reste2) / divis1)) + cumul1);
reste3 = capacite * (mdivis2 * (quot2 = AINT((cumul2 = termec + reste3) / divis2)) + cumul2);
reste4 = capacite * (mdivis3 * (quot3 = AINT((cumul3 = termec + reste4) / divis3)) + cumul3);
serie[index] = serie[index] + (mangle2 * quot1 + quot2) * angle2 - quot3;
}
quot2 = AINT(((termem = AINT((termem + reste1) / angle6)) + reste3) / divis2);
quot3 = AINT((termem + reste4) / divis3);
seriem = seriem + (mangle2 * AINT((termem + reste2) / divis1) + quot2) * angle2 - quot3;
if (terme[index0 - I] != zero) goto eti890;
if (index0 == mindex) goto eti2500;
index0 = index0 + I;
eti890:
mdivis3 = - (divis3 = (divis2 = - (mdivis2 = (mdivis1 = - (divis1 = divis1 + c6)) - c2)) + c2);
eti201:
reste1 = capacite * (mangle6 * (termec = (terme[index0 - I] = AINT((cumul4 = terme[index0 - I]) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT(termec / divis1)) + termec);
reste3 = capacite * (mdivis2 * (quot2 = AINT(termec / divis2)) + termec);
reste4 = capacite * (mdivis3 * (quot3 = AINT(termec / divis3)) + termec);
serie[index0 - I] = serie[index0 - I] + (mangle2 * quot1 + quot2) * mangle2 + quot3;
for (index = index0 ; index <= mindexm1 ; index++)
{
reste1 = capacite * (mangle6 * (terme[index] = (termec = AINT((cumul4 = terme[index] + reste1) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT((cumul1 = termec + reste2) / divis1)) + cumul1);
reste3 = capacite * (mdivis2 * (quot2 = AINT((cumul2 = termec + reste3) / divis2)) + cumul2);
reste4 = capacite * (mdivis3 * (quot3 = AINT((cumul3 = termec + reste4) / divis3)) + cumul3);
serie[index] = serie[index] + (mangle2 * quot1 + quot2) * mangle2 + quot3;
}
quot2 = AINT(((termem = AINT((termem + reste1) / angle6)) + reste3) / divis2);
quot3 = AINT((termem + reste4) / divis3);
seriem = seriem + (mangle2 * AINT((termem + reste2) / divis1) + quot2) * mangle2 + quot3;
if (terme[index0 - I] != zero) goto eti8906;
if (index0 == mindex) goto eti2510;
index0 = index0 + I;
eti8906:
mdivis3 = - (divis3 = (divis2 = - (mdivis2 = (mdivis1 = - (divis1 = divis1 + c6)) - c2)) + c2);
goto eti200;
eti2500:
divis3 = (divis2 = (divis1 = divis1 + c6) + c2) + c2;
quot2 = AINT((termem = AINT(termem / angle6)) / divis2);
quot3 = AINT(termem / divis3);
seriem = seriem + (mangle2 * AINT(termem / divis1) + quot2) * mangle2 + quot3;
eti2501:
if (termem == zero) goto eti25;
divis3 = (divis2 = (divis1 = divis1 + c6) + c2) + c2;
quot2 = AINT((termem = AINT(termem / angle6)) / divis2);
quot3 = AINT(termem / divis3);
seriem = seriem + (mangle2 * AINT(termem / divis1) + quot2) * angle2 - quot3;
if (termem == zero) goto eti25;
divis3 = (divis2 = (divis1 = divis1 + c6) + c2) + c2;
quot2 = AINT((termem = AINT(termem / angle6)) / divis2);
quot3 = AINT(termem / divis3);
seriem = seriem + (mangle2 * AINT(termem / divis1) + quot2) * mangle2 + quot3;
goto eti2501;
eti2510:
divis3 = (divis2 = (divis1 = divis1 + c6) + c2) + c2;
quot2 = AINT((termem = AINT(termem / angle6)) / divis2);
quot3 = AINT(termem / divis3);
seriem = seriem + (mangle2 * AINT(termem / divis1) + quot2) * angle2 - quot3;
eti2511:
if (termem == zero) goto eti25;
divis3 = (divis2 = (divis1 = divis1 + c6) + c2) + c2;
quot2 = AINT((termem = AINT(termem / angle6)) / divis2);
quot3 = AINT(termem / divis3);
seriem = seriem + (mangle2 * AINT(termem / divis1) + quot2) * mangle2 + quot3;
if (termem == zero) goto eti25;
divis3 = (divis2 = (divis1 = divis1 + c6) + c2) + c2;
quot2 = AINT((termem = AINT(termem / angle6)) / divis2);
quot3 = AINT(termem / divis3);
seriem = seriem + (mangle2 * AINT(termem / divis1) + quot2) * angle2 - quot3;
goto eti2511;
eti9200:
seriem = seriem + AINT((termem = AINT(termem / angle6)) / divis1);
eti9201:
if (termem == zero) goto eti92;
divis1 = divis1 + c2;
seriem = seriem - AINT((termem = AINT(termem / angle6)) / divis1);
if (termem == zero) goto eti92;
divis1 = divis1 + c2;
seriem = seriem + AINT((termem = AINT(termem / angle6)) / divis1);
goto eti9201;
eti9210:
seriem = seriem - AINT((termem = AINT(termem / angle6)) / divis1);
eti9211:
if (termem == zero) goto eti92;
divis1 = divis1 + c2;
seriem = seriem + AINT((termem = AINT(termem / angle6)) / divis1);
if (termem == zero) goto eti92;
divis1 = divis1 + c2;
seriem = seriem - AINT((termem = AINT(termem / angle6)) / divis1);
goto eti9211;
/*************************************************************************************************************************************/
/* */
/* C A L C U L D E L A S E R I E A R C T G ( 1 / 2 3 9 ) : */
/* */
/*************************************************************************************************************************************/
eti25:
reste1 = c410;
mangle6 = -(angle6 = c239);
for (index = INDEX0 ; index <= mindexm1 ; index++)
{
reste1 = capacite * (mangle6 * (terme[index] = AINT(reste1 / angle6)) + reste1);
serie[index] = serie[index] - terme[index];
}
seriem = seriem - (termem = AINT(reste1 / angle6));
divis1 = (DOOBLE)W;
mangle6 = -(angle6 = c57121);
index0 = INDEX0 + I;
eti10:
mdivis1 = - (divis1 = divis1 + c2);
if (terme[index0 - I] != zero) goto eti91;
if (index0 == mindex) goto eti9200;
index0 = index0 + I;
eti91:
reste1 = capacite * (mangle6 * (terme[index0 - I] = (termec = AINT((cumul4 = terme[index0 - I]) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT(termec / divis1)) + termec);
serie[index0 - I] = serie[index0 - I] + quot1;
for (index = index0 ; index <= mindexm1 ; index++)
{
reste1 = capacite * (mangle6 * (terme[index] = (termec = AINT((cumul4 = terme[index] + reste1) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT((cumul1 = termec + reste2) / divis1)) + cumul1);
serie[index] = serie[index] + quot1;
}
seriem = seriem + AINT(((termem = AINT((termem + reste1) / angle6)) + reste2) / divis1);
mdivis1 = - (divis1 = divis1 + c2);
if (terme[index0 - I] != zero) goto eti2;
if (index0 == mindex) goto eti9210;
index0 = index0 + I;
eti2:
reste1 = capacite * (mangle6 * (terme[index0 - I] = (termec = AINT((cumul4 = terme[index0 - I]) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT(termec / divis1)) + termec);
serie[index0 - I] = serie[index0 - I] - quot1;
for (index = index0 ; index <= mindexm1 ; index++)
{
reste1 = capacite * (mangle6 * (terme[index] = (termec = AINT((cumul4 = terme[index] + reste1) / angle6))) + cumul4);
reste2 = capacite * (mdivis1 * (quot1 = AINT((cumul1 = termec + reste2) / divis1)) + cumul1);
serie[index] = serie[index] - quot1;
}
seriem = seriem - AINT(((termem = AINT((termem + reste1) / angle6)) + reste2) / divis1);
goto eti10;
eti92:
index = mindex;
serie[index] = seriem;
mindex = mindexm1;
eti35:
if (serie[index] < zero) goto eti37;
serie[index] = (reste = AINT(serie[index] / capacite)) * mcapacite + serie[index];
if ((index = index - I) < zero) goto eti42;
serie[index] = serie[index] + reste;
goto eti35;
eti37:
serie[index] = capacite * (reste = AINT(((DOOBLE)I) + (serie[index] / mcapacite))) + serie[index];
if ((index = index - I) < zero) goto eti42;
serie[index] = serie[index] - reste;
goto eti35;
/*************************************************************************************************************************************/
/* */
/* I M P R E S S I O N D U R E S U L T A T : */
/* */
/*************************************************************************************************************************************/
eti42:
Eclock("Duree du calcul de PI",VRAI);
CAL2(Prin1("\n\n\n %1d,",INTE(reste)));
for (index = INDEX0 ; index < mindex ; index = index + mots)
{
for (index1 = INDEX0 ; index1 <= mots - BEGIN_AT_0 ; index1++)
{
for (index2 = INDEX0 ; index2 <= decimales - BEGIN_AT_0 ; index2++)
{
q = AINT(serie[index + index1] / BASE10);
chiffre = INTE(serie[index + index1] - (q * BASE10));
serie[index + index1] = q;
ligne[index1][decimales - BEGIN_AT_0 - index2] = chiffres[chiffre];
}
}
for (index1 = INDEX0 ; index1 <= mots - BEGIN_AT_0 ; index1++)
{
for (index2 = INDEX0 ; index2 <= decimales - BEGIN_AT_0 ; index2++)
{
CAL2(Prin1("%1c",ligne[index1][index2]));
cumul1 = ((index + index1) * decimales) + index2;
reste2 = cumul1 - ((nombrem * nombred) * AINT(cumul1 / (nombrem * nombred)));
reste3 = cumul1 - (nombred * AINT(cumul1 / nombred));
if (reste3 == (DOOBLE)(nombred - BEGIN_AT_0))
{
CAL2(Prin0(" "));
}
if (reste2 == (DOOBLE)((nombrem * nombred) - BEGIN_AT_0))
{
CAL2(Prin0("\n "));
}
}
}
}
CAL2(Prin0("\n"));
}
Copyright © Jean-François Colonna, 2019-2021.
Copyright © CMAP (Centre de Mathématiques APpliquées) UMR CNRS 7641 / Ecole Polytechnique, 2019-2021.