/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        C A L C U L   D E   L ' E N S E M B L E   D E   M A N D E L B R O T                                                        */
/*        A   D E S   F I N S   D E   T E S T   D E   L ' O P T I M I S A T I O N   S U R   ' $LACT1B '  :                           */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        Author of '$xtc/mandel.61$c' :                                                                                             */
/*                                                                                                                                   */
/*                    John F. Kolonna (LACTAMME, 20210923125458).                                                                    */
/*                                                                                                                                   */
/*************************************************************************************************************************************/

#include  <stdlib.h>
#include  <stdio.h>

extern    void      *malloc();

static    int       dimX=2000;
#define   Xmin      0
#define   Xmax      (Xmin + (dimX-1))
                                        /* Definition des abscisses.                                                                 */
static    int       dimY=2000;
#define   Ymin      0
#define   Ymax      (Ymin + (dimY-1))
                                        /* Definition des ordonnees.                                                                 */

#define   IMAGE(x,y)                                                                                                                    \
                    (*(image + (((y-Ymin)*dimX) + (x-Xmin))))                                                                           \
                                        /* Acces a un point de l'image.                                                              */

#define   NOIR                                                                                                                          \
                    0
#define   BLANC                                                                                                                         \
                    255

#define   ITERATION 1024

#define   Cn(n)                                                                                                                         \
                    ((unsigned char)(NOIR + (((PRECIS)(n-1)/(PRECIS)ITERATION)*((PRECIS)(BLANC-NOIR)))))                                \
                                        /* Definition de niveaux dans l'image...                                                     */

#define   SEUIL     4.0

#define   store(n,x,y)                                                                                                                  \
                    {                                                                                                                   \
                    IMAGE(x,y) = Cn(n);                                                                                                 \
                    }                                                                                                                   \
                                        /* Rangement d'un point valide d'une image.                                                  */

#define   XG        (-2.00)
#define   XD        (0.50)
                                        /* Extremites Gauche et Droite sur l'axe des 'X'.                                            */
#define   YB        (-1.25)
#define   YH        (1.25)
                                        /* Extremites Bas et Haut sur l'axe des 'Y'.                                                 */

#define   PRECIS    double                                                                                                              \
                                        /* Precision des calculs.                                                                    */

typedef   struct    complexe  {
                              PRECIS    reelle;
                              PRECIS    imaginaire;
                              }         complexe;
                                        /* 'struct' de definition d'un nombre complexe...                                            */
#define   Reelle(z)           (z.reelle)
#define   Imaginaire(z)       (z.imaginaire)

#define   Cinit(z,x1,y1)                                                                                                                \
                    {                                                                                                                   \
                    Reelle(z) = x1;                                                                                                     \
                    Imaginaire(z) = y1;                                                                                                 \
                    }
#define   Cegal(z,z1)                                                                                                                   \
                    {                                                                                                                   \
                    Reelle(z) = Reelle(z1);                                                                                             \
                    Imaginaire(z) = Imaginaire(z1);                                                                                     \
                    }
#define   Csomme(z,z1,z2)                                                                                                               \
                    {                                                                                                                   \
                    Reelle(z) = Reelle(z1) + Reelle(z2);                                                                                \
                    Imaginaire(z) = Imaginaire(z1) + Imaginaire(z2);                                                                    \
                    }
#define   Cproduit(z,z1,z2)                                                                                                             \
                    {                                                                                                                   \
                    Reelle(z) = (Reelle(z1)*Reelle(z2)) - (Imaginaire(z1)*Imaginaire(z2));                                              \
                    Imaginaire(z) = (Reelle(z2)*Imaginaire(z1)) + (Reelle(z1)*Imaginaire(z2));                                          \
                    }
#define   Cmodule2(z)                                                                                                                   \
                    ((Reelle(z)*Reelle(z)) + (Imaginaire(z)*Imaginaire(z)))

#define   PETIT     1.0e-16
#define   GRAND     1.0e32
#define   deborde(x)                                                                                                                    \
                    {                                                                                                                   \
                    if   (((x) < -PETIT) || ((x) > PETIT))                                                                              \
                         {                                                                                                              \
                         }                                                                                                              \
                    else                                                                                                                \
                         {                                                                                                              \
                         x = 0.0;                                                                                                       \
                         }                                                                                                              \
                                                                                                                                        \
                    if   ((x) < -GRAND)                                                                                                 \
                         {                                                                                                              \
                         x = -GRAND;                                                                                                    \
                         }                                                                                                              \
                    else                                                                                                                \
                         {                                                                                                              \
                         }                                                                                                              \
                                                                                                                                        \
                    if   ((x) > GRAND)                                                                                                  \
                         {                                                                                                              \
                         x = GRAND;                                                                                                     \
                         }                                                                                                              \
                    else                                                                                                                \
                         {                                                                                                              \
                         }                                                                                                              \
                    }

#define   Cdeborde(z)                                                                                                                   \
                    {                                                                                                                   \
                    deborde(Reelle(z));                                                                                                 \
                    deborde(Imaginaire(z));                                                                                             \
                    }

void      main()
     {
     complexe  C;
     complexe  Zn;
     complexe  z1;
     int                 x,y;
                                        /* Definition des coordonnees.                                                               */
     double    moyenne=0;
     unsigned  char      *image;
                                        /* Definition de l'image a generer...                                                        */

     image=malloc(dimY*dimX);
                                        /* Definition de l'image a generer...                                                        */

     for       (y=Ymin ; y<=Ymax ; y++)
               {
               for       (x=Xmin ; x<=Xmax ; x++)
                    {
                         int       index=1;
                         int       iterer=0;
                                        /* Variables d'iterations...                                                                 */
                         Cinit(C
                              ,((PRECIS)x / (PRECIS)dimX)*(XD-XG) + XG
                              ,((PRECIS)y / (PRECIS)dimY)*(YH-YB) + YB
                               );
                                        /* Initialisation du point complexe courant 'C' sur le point {x,y}.                          */
                         Cinit(Zn,0.0,0.0);
                                        /* Initialisation de la suite Zn.                                                            */

                         while     (iterer == 0)
                                   {
                                   Cproduit(z1,Zn,Zn);
                                   Csomme(Zn,z1,C);
                                   Cdeborde(Zn);
                                        /* Iteration de la suite Zn.                                                                 */
                                   if        ((index>=ITERATION) || (Cmodule2(Zn) > SEUIL))
                                             {
                                             store(index,x,y);
                                             iterer=1;
                                             }
                                   else
                                             {
                                             index++;
                                             }
                                   }
                         }
               }

     for       (y=Ymin ; y<=Ymax ; y++)
               {
               for       (x=Xmin ; x<=Xmax ; x++)
                         {
                         moyenne = moyenne + (double)IMAGE(x,y);
                         }
               }

     fprintf(stderr,"moyenne=%f\n",moyenne/((double)(dimX*dimY)));
                                        /* Et ce afin d'utiliser la matrice 'IMAGE(...)' et que l'optimiseur ne fasse disparaitre    */
                                        /* une bonne partie du calcul...                                                             */

                                        /* Resultat de 'time' sur '$LACT1B' en fonction de l'option d'optmisation utilisee :         */
                                        /*                                                                                           */
                                        /*                  -O0                 10.592u 0.000s 0:10.59 100.0%   0+0k 0+0io 0pf+0w    */
                                        /*                  -O1                 5.048u 0.000s 0:05.04 100.0%    0+0k 0+0io 0pf+0w    */
                                        /*                  -O2                 6.769u 0.000s 0:06.77 99.8%     0+0k 0+0io 0pf+0w    */
                                        /*                  -O3                 6.747u 0.000s 0:06.74 100.0%    0+0k 0+0io 0pf+0w    */
                                        /*                                                                                           */
     }