/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        T E S T   D U   P R O B L E M E   ' v $xiii/tri_image$FON 20040206173024 '                                                 */
/*        S U R   ' $LACT15 ' ,   ' $LACT16 '   E T   ' $CMAP28 '  :                                                                 */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        Nota :                                                                                                                     */
/*                                                                                                                                   */
/*                    Sur les MACHINEs ou les options suivantes de                                                                   */
/*                  '$Cc' sont disponibles, le probleme disparait :                                                                  */
/*                                                                                                                                   */
/*                                      -ffloat-store                                                                                */
/*                                                                                                                                   */
/*                  (disponible partout, elle reduit de 80 a 64 bits                                                                 */
/*                  la precision) ou :                                                                                               */
/*                                                                                                                                   */
/*                                      -mfpmath=sse,387 -msse2                                                                      */
/*                                                                                                                                   */
/*                  (utilise simultanement les instructions flottantes                                                               */
/*                  des jeux 'standard 387' et 'SSE', elle est disponible                                                            */
/*                  sur '$LACT16', mais pas sur '$LACT15').                                                                          */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        Author of '$Dbugs/APC$D/LinuxDebian$D/GCC$D/flottant.02$c' :                                                               */
/*                                                                                                                                   */
/*                    John F. Colonna (LACTAMME, 20040219104426).                                                                    */
/*                                                                                                                                   */
/*************************************************************************************************************************************/

#include  < stdio.h >

extern    double    sin();

typedef   union     flottant  {
                              double    d;
                              struct    i         {
                                                  int       i1;
                                                  int       i2;
                                                  }         i;
                              }         flottant;
#define   pd(x)     (x.d)
#define   pi1(x)    (x.i.i1)
#define   pi2(x)    (x.i.i2)
                                        /* Outils destines a l'edition de nombre flottants 64 bits sous forme de 8+8 chiffres        */
                                        /* hexa-decimaux.                                                                            */

#define   dimY      3
#define   Ymin      0
#define   Ymax      (Ymin+dimY-1)

#define   dimX      7
#define   Xmin      0
#define   Xmax      (Xmin+dimX-1)

#define   DM(t)               t[dimY*dimX]
#define   IM(t,x,y)           (*((t)+((y-Ymin)*(dimX)+(x-Xmin))))
#define   PM(c,min,max)       (c=(min) ; c <= (max) ; c++)
                                        /* Definition at acces a des matrices sous forme de vecteurs.                                */

#define   ABSO(x)   (((x) > 0) ? (x) : -(x))
#define   MIN2(a,b) (((a) < (b)) ? (a) : (b))
#define   MAX2(a,b) (((a) > (b)) ? (a) : (b))
#define   SOUA(a,b) ABSO((a) - (b))

#define   Amplification(x)    (x)
                                        /* Voir le commentaire du 20040220100843...                                                  */

void      fonction(argumentA,argumentB)
double    argumentA;
double    argumentB;
                                        /* Ces deux arguments sont destines a forcer 'cX=1' ci-dessous sans perturber l'optmisation. */
          {
          int       X,Y,YY;
          double    DM(A1),DM(TM),DM(A2),DM(A1N);
          double    minA1=+1e308,maxA1=-1e308;

          for       PM(Y,Ymin,Ymax)
                    {
                    for       PM(X,Xmin,Xmax)
                              {
                              double    valeurX=((1.0+sin(10.0*((double)(X-Xmin))/((double)(Xmax-Xmin))))/2.0);
                              double    valeurY=((1.0+sin(10.0*((double)(Y-Ymin))/((double)(Ymax-Ymin))))/2.0);
                                        /* Le 20040219091351, je note sur '$LACT15' que le probleme disparait si l'on remplace       */
                                        /* la constante multiplicative 10.0 ci-dessus par 3.0. Or la difference entre les deux       */
                                        /* fichiers assembleur correspondant est :                                                   */
                                        /*                                                                                           */
                                        /*                  .long     0x0,0x40240000                (10.0)                           */
                                        /*                                                                                           */
                                        /* different de :                                                                            */
                                        /*                                                                                           */
                                        /*                  .long     0x0,0x40080000                (3.0)                            */
                                        /*                                                                                           */
                                        /* c'est-a-dire la constante elle-meme. Ainsi, l'optimisation ne depend pas de cette         */
                                        /* constante, alors que le comportement en depend. Peut-etre le processeur flottant est-il   */
                                        /* en cause, cette defaillance ne se manifestant que lorsque l'optimisation est activee et   */
                                        /* que donc certaines instructions sont utilisees alors qu'elles ne le sont pas lorsqu'il    */
                                        /* n'y a pas optimisation...                                                                 */

                              IM(A1,X,Y) = valeurY;
                              minA1=MIN2(minA1,valeurY);
                              maxA1=MAX2(maxA1,valeurY);

                              IM(A2,X,Y) = valeurX;
                              IM(TM,X,Y) = valeurX;
                                        /* On notera que :                                                                           */
                                        /*                                                                                           */
                                        /*                  IM(TM,X,Ymin)=...=IM(TM,X,Y)=...=IM(TM,X,Ymax)              \-/ Y        */
                                        /*                                                                                           */
                              }
                    }

          for       PM(Y,Ymin,Ymax)
                    {
                    for       PM(X,Xmin,Xmax)
                              {
                              IM(A1N,X,Y) = (IM(A1,X,Y)-minA1)/(maxA1-minA1);
                                        /* A priori, etant donnees les definitions de {valeurX,valeurY}, les valeurs de 'IM(A1,X,Y)' */
                                        /* sont deja dans [0,1] sans que les bornes soient attteintes. Cette renormalisation fait    */
                                        /* que les bornes [0,1] sont atteintes ; la supprimer fait disparaitre le probleme. C'est    */
                                        /* peut-etre plus les valeurs flottantes alors generees dans 'IM(A1N,X,Y)' que le code de    */
                                        /* renormalisation qui engendre le probleme...                                               */
                              }
                    }

          for       PM(Y,Ymin,Ymax)
                    {
                    for       PM(X,Xmin,Xmax)
                              {
                              int       cX=(int)((argumentA*(Xmin+((Xmax-Xmin)*IM(A1N,X,Y))))+argumentB);
                                        /* On ne peut forcer "betement" :                                                            */
                                        /*                                                                                           */
                                        /*                  int       cX=(int)(Xmin+1)                                               */
                                        /*                                                                                           */
                                        /* sans faire disparaitre le probleme, d'ou {argumentA,argumentB}...                         */
                              int       cY=Ymin;
                              double    ncA2=IM(A2,X,Y);
                              double    dm=+1e308;

                              if        ((cX < Xmin) || (cX > Xmax) || (cY < Ymin) || (cY > Ymax))
                                        /* Ce test ne peut a priori jamais etre vrai etant donne que :                               */
                                        /*                                                                                           */
                                        /*                  cX=Xmin+1           (<= Xmax)                                            */
                                        /*                  cY=Ymin                                                                  */
                                        /*                                                                                           */
                                        /* par construction...                                                                       */
                                        {
                                        printf("les coordonnees du produit generalise inverse a droite sont incorrectes\n");
                                        printf("X=%d : doit etre dans [%d,%d]\n",cX,Xmin,Xmax);
                                        fflush(stdout);
                                        printf("Y=%d : doit etre dans [%d,%d]\n",cY,Ymin,Ymax);
                                        fflush(stdout);
                                        /* On notera que supprimer, par exemple, l'un des 'fflush(...)'s fait disparaitre le         */
                                        /* probleme...                                                                               */
                                        }
                              else
                                        {
                                        }

                              for       PM(YY,Ymin,Ymax)
                                        /* Boucle dite "interne" destinee a recherche le minimum 'dm' de l'ensemble des 'dc' (qui    */
                                        /* rappelons-le sont tous egaux entre-eux) decrit par 'YY'.                                  */
                                        {
                                        double    nctm=IM(TM,cX,YY);
                                        double    dc;

                                        dc = SOUA(nctm,ncA2);
                                        /* On se souviendra que :                                                                    */
                                        /*                                                                                           */
                                        /*                  IM(TM,X,Ymin)=...=IM(TM,X,YY)=...=IM(TM,X,Ymax)             \-/ YY       */
                                        /*                                                                                           */
                                        /* ainsi, a l'interieur de cette boucle interne 'PM(YY,Ymin,Ymax)', on a :                   */
                                        /*                                                                                           */
                                        /*                  nctm=constante                                              \-/ YY       */
                                        /*                                                                                           */
                                        /* enfin, 'ncA2' etant exterieur a cette boucle, on a donc :                                 */
                                        /*                                                                                           */
                                        /*                  dc=constante                                                \-/ YY       */
                                        /*                                                                                           */
                                        /* et le test suivant qui recherche le minimum ('dm') de l'ensemble des 'dc' decrit par      */
                                        /* la boucle interne 'PM(YY,Ymin,Ymax)' doit etre toujours faux, sauf la premiere fois       */
                                        /* pour :                                                                                    */
                                        /*                                                                                           */
                                        /*                  YY=Ymin                                                                  */
                                        /*                                                                                           */
                                        /* 'dm' etant initialise avec un tres grand nombre flottant (+1e308)...                      */

                                        if        (dm > dc)
                                        /* Le 20040303134627, je note que si l'on remplace ce test par :                             */
                                        /*                                                                                           */
                                        /*                  if        (FfIFGT(dm,dc))                                                */
                                        /*                                                                                           */
                                        /* en compilant grace a la commande :                                                        */
                                        /*                                                                                           */
                                        /*                  cc        -O3 flottant.02$c                                           \  */
                                        /*                            -lm                                                         \  */
                                        /*                            $xbg/*$SO $xbii/*$SO $xbipf/*$SO $xbmf/*$SO $xbin/*$SO         */
                                        /*                                                                                           */
                                        /* cela se met a fonctionner correctement ('v $xig/fonct$vv$FON 20040303111309'), comme      */
                                        /* je m'y attendais...                                                                       */
                                                  {
                                                  if        (YY > Ymin)
                                        /* En toute logique, ce test ne peut jamais etre vrai...                                     */
                                                            {
                                                            flottant  fdm,fdc,fdiff;

                                                            pd(fdm) = Amplification(dm);
                                                            pd(fdc) = Amplification(dc);
                                        /* La procedure 'Amplification(...)' a ete introduite le 20040220100843 afin de montrer      */
                                        /* qu'en fait 'dm' et 'dc' n'etaient pas egaux : ils l'etaient sur 64 bits (voir la sortie   */
                                        /* hexa-decimale), mais pas sur 80 bits (d'ou l'option '-ffloat-store'...).                  */
                                        /*                                                                                           */
                                        /* En definissant :                                                                          */
                                        /*                                                                                           */
                                        /*                  #define   Amplification(x)    (1e20*(x))                                 */
                                        /*                                                                                           */
                                        /* on obtient (pour la premiere anomalie rencontree) :                                       */
                                        /*                                                                                           */
                                        /*                  X=0002 Y=0000 YY=0001 :                                                  */
                                        /*                  dm=59298796031362523136.000000000000000000 (267f0da7,4409b77d)           */
                                        /*                                                                     #                     */
                                        /*                  dc=59298796031362514944.000000000000000000 (267f0da6,4409b77d)           */
                                        /*                                                                     #                     */
                                        /*                                                                                           */
                                        /*                  dm-dc=0.000000000000000056 (00000000,3c900000)                           */
                                        /*                                                                                           */
                                        /* ou l'on voit bien (voir les "#"s) que 'dm' et 'dc' etaient en fait bien differents, en    */
                                        /* notant bien que 'diff' est egal a 'dm-dc' et non pas a 'pd(fdm)-pd(fdc)', c'est-a-dire    */
                                        /* la difference n'est pas amplifiee par 'Amplification(...)' a cause de la remarque faite   */
                                        /* dans le commentaire suivant...                                                            */
                                                            pd(fdiff) = dm - dc;
                                        /* ATTENTION, ici a cause de '$LACT16', il est impossible d'ecrire :                         */
                                        /*                                                                                           */
                                        /*                  pd(fdiff) =  pd(fdm) -  pd(fdc);                                         */
                                        /*                                                                                           */
                                        /* sous peine de voir disparaitre le probleme...                                             */

                                                            printf("X=%04d Y=%04d YY=%04d :\n"
                                                                  ,X,Y,YY
                                                                   );

                                        /* ATTENTION : on notera que couper en deux le 'printf(...)' qui suit fait disparaitre       */
                                        /* le probleme !                                                                             */
                                                            printf("dm   =%.18f (0x%08x,0x%08x)\ndc   =%.18f (0x%08x,0x%08x)\ndm-dc=%.18f (0x%08x,0x%08x)\n"
                                                                  ,pd(fdm),pi1(fdm),pi2(fdm)
                                                                  ,pd(fdc),pi1(fdc),pi2(fdc)
                                                                  ,pd(fdiff),pi1(fdiff),pi2(fdiff)
                                                                   );
                                        /* ATTENTION : on notera que couper en deux le 'printf(...)' qui precede fait disparaitre    */
                                        /* le probleme !                                                                             */
                                                            }
                                                  else
                                                            {
                                                            }

                                                  dm=dc;
                                                  cY=YY;
                                                  }
                                        else
                                                  {
                                                  }
                                        }
                              }
                    }
          }

main()
          {
          double    parametreA=0.0;
          double    parametreB=(double)MIN2(Xmin+1,Xmax);
                                        /* Ces deux arguments sont destines a forcer 'cX=1' ci-dessus sans perturber l'optmisation.  */

          fonction(parametreA,parametreB);
          }



Copyright (c) Jean-François Colonna, 2004-2012.
Copyright (c) France Telecom R&D and CMAP (Centre de Mathématiques APpliquées) UMR CNRS 7641 / Ecole Polytechnique, 2004-2012.