/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        C A L C U L   I N D I V I D U E L   D E S   C H I F F R E S   D E   P I   E N   B A S E   1 6  :                           */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        ATTENTION :                                                                                                                */
/*                                                                                                                                   */
/*                    Cette methode est decrite a la                                                                                 */
/*                  page 54 du numero 1, volume 19 de                                                                                */
/*                  The Mathematical Intelligencer.                                                                                  */
/*                  Ce programme est a compiler avec                                                                                 */
/*                  'cc -64 -O3' sur 'SYSTEME_SGO200A1_IRIX_CC'                                                                      */
/*                  par exemple.                                                                                                     */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        Author of '$xtc/pi_chif16.01$c' :                                                                                          */
/*                                                                                                                                   */
/*                    Jean-Francois COLONNA (LACTAMME, AAAAMMJJhhmmss).                                                              */
/*                                                                                                                                   */
/*************************************************************************************************************************************/

#include  <stdio.h>

extern    double    fmod();
extern    double    log();
extern    double    pow();

#define   UN        1.0
#define   DEUX      2.0
#define   QUATRE    4.0
#define   SEPT      7.0
#define   HUIT      8.0

#define   SEIZE     16

#define   Entier    long int

#define   INTE(x)   ((Entier)(x))
#define   FLOT(x)   ((double)(x))
#define   MODU(x,m) ((x) % (m))

#define   EXPONENTIELLE(exponentielle,base,exposant,modulo)                                                                             \
                    {                                                                                                                   \
                    Entier    Cbase=base;                                                                                               \
                    Entier    Cexposant=exposant;                                                                                       \
                    Entier    Cmodulo=modulo;                                                                                           \
                                                                                                                                        \
                    exponentielle=1;                                                                                                    \
                                                                                                                                        \
                    if        (Cexposant > 0)                                                                                           \
                              {                                                                                                         \
                              Entier    puissance_de_2;                                                                                 \
                              Entier    iterer=0;                                                                                       \
                                                                                                                                        \
                              puissance_de_2=INTE(pow(DEUX,FLOT(INTE(log(FLOT(Cexposant))/log(DEUX)))));                                \
                                                                                                                                        \
                              while     (iterer == 0)                                                                                   \
                                        {                                                                                               \
                                        if        (Cexposant >= puissance_de_2)                                                         \
                                                  {                                                                                     \
                                                  exponentielle = MODU(Cbase*exponentielle,Cmodulo);                                    \
                                                  Cexposant = Cexposant - puissance_de_2;                                               \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  }                                                                                     \
                                                                                                                                        \
                                        puissance_de_2 = puissance_de_2 / INTE(DEUX);                                                   \
                                                                                                                                        \
                                        if        (puissance_de_2 >= 1)                                                                 \
                                                  {                                                                                     \
                                                  exponentielle = MODU(exponentielle*exponentielle,Cmodulo);                            \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  iterer++;                                                                             \
                                                  }                                                                                     \
                                        }                                                                                               \
                              }                                                                                                         \
                    else                                                                                                                \
                              {                                                                                                         \
                              }                                                                                                         \
                                                                                                                                        \
                    if        (exponentielle < 0)                                                                                       \
                              {                                                                                                         \
                              printf("\n exponentielle negative : %ld puissance %ld modulo %ld = %ld"                                   \
                                    ,base                                                                                               \
                                    ,exposant                                                                                           \
                                    ,modulo                                                                                             \
                                    ,exponentielle                                                                                      \
                                     );                                                                                                 \
                                        /* Ceci peut se produire lorsque les operations entieres donnent des resulats ne tenant      */ \
                                        /* pas dans 'Entier'...                                                                      */ \
                              }                                                                                                         \
                    else                                                                                                                \
                              {                                                                                                         \
                              }                                                                                                         \
                    }                                                                                                                   \
                                        /* Calcul de :                                                                               */ \
                                        /*                                                                                           */ \
                                        /*                                      exposant                                             */ \
                                        /*                  exponentielle = base        modulo                                       */ \
                                        /*                                                                                           */ \
                                        /* quelle que soit (ou presque) la "taille"...                                               */

#define   CUMUL(serie,terme)                                                                                                            \
                    serie = CORRIGE(serie + terme);

#define   SERIE(serie,rang,increment)                                                                                                   \
                    {                                                                                                                   \
                    Entier    base=SEIZE;                                                                                               \
                    Entier    exposant;                                                                                                 \
                    Entier    modulo;                                                                                                   \
                    Entier    k;                                                                                                        \
                    Entier    precision=10;                                                                                             \
                                        /* Parametre arbitraire permettant de prendre en compte quelques (=precision) termes         */ \
                                        /* de la serie apres le rang 'rang'...                                                       */ \
                                                                                                                                        \
                    serie=0;                                                                                                            \
                                                                                                                                        \
                    for       (k=0 ; k<=rang ; k++)                                                                                     \
                              {                                                                                                         \
                              Entier    exponentielle;                                                                                  \
                                                                                                                                        \
                              exposant = rang-k;                                                                                        \
                                                                                                                                        \
                              modulo = (8*k)+(increment);                                                                               \
                              EXPONENTIELLE(exponentielle,base,exposant,modulo);                                                        \
                                                                                                                                        \
                              CUMUL(serie,FLOT(exponentielle)/FLOT(modulo));                                                            \
                              }                                                                                                         \
                                                                                                                                        \
                    for       (k=(rang+1) ; k<=(rang+precision) ; k++)                                                                  \
                              {                                                                                                         \
                              exposant = rang-k;                                                                                        \
                                                                                                                                        \
                              modulo = (8*k)+(increment);                                                                               \
                                                                                                                                        \
                              CUMUL(serie,pow(FLOT(base),FLOT(exposant))/FLOT(modulo));                                                 \
                              }                                                                                                         \
                    }                                                                                                                   \
                                        /* Calcul de la serie :                                                                      */ \
                                        /*                                                                                           */ \
                                        /*                                   k=infini                                                */ \
                                        /*                                     _____                                                 */ \
                                        /*                                     \               1                                     */ \
                                        /*                serie(increment) =   /     ----------------------                          */ \
                                        /*                                     -----    k                                            */ \
                                        /*                                      k=0   16 (8.k + increment)                           */ \
                                        /*                                                                                           */

#define   CORRIGE(x)                                                                                                                    \
                    fmod((x),UN)

#define   DECIMALES(x)                                                                                                                  \
                    INTE(pow(FLOT(SEIZE),HUIT)*(x))

#define   EDITE(valeur,titre)                                                                                                           \
                    {                                                                                                                   \
                    char      chaine[]="12345678";                                                                                      \
                    sprintf(chaine,"%lx",INTE(DECIMALES(valeur)));                                                                      \
                    printf("%s = %.4s\n",titre,chaine);                                                                                 \
                    }                                                                                                                   \
                                        /* Edition en hexa-decimal d'une valeur...                                                   */

main()
     {
     Entier    rangD=999999;
     Entier    rangA=999999;
     Entier    rang;
     Entier    facteur1=+4,increment1=1;
     Entier    facteur2=-2,increment2=4;
     Entier    facteur3=-1,increment3=5;
     Entier    facteur4=-1,increment4=6;
     double    serie1;
     double    serie2;
     double    serie3;
     double    serie4;
     double    serie;
rangD=-1;
rangA=-1;

     for       (rang=rangD ; rang<=rangA ; rang++)
               {
               printf("position=%ld\n",rang+1);

               SERIE(serie1,rang,increment1);
               EDITE(serie1,"      serie 1");

               SERIE(serie2,rang,increment2);
               EDITE(serie2,"      serie 2");

               SERIE(serie3,rang,increment3);
               EDITE(serie3,"      serie 3");

               SERIE(serie4,rang,increment4);
               EDITE(serie4,"      serie 4");

               serie = CORRIGE((facteur1*serie1) + (facteur2*serie2) + (facteur3*serie3) + (facteur4*serie4));
                                        /* Calcul de pi :                                                                            */
                                        /*                                                                                           */
                                        /*                  pi = 4.serie(1) - 2.serie(4) - 1.serie(5) - 1.serie(6)                   */
                                        /*                                                                                           */

               EDITE(serie,"           pi");
               }
     }