/*************************************************************************************************************************************/ /* */ /* 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); }