/*************************************************************************************************************************************/
/* */
/* I N T E R P O L A T I O N D U N - I E M E D E G R E : */
/* */
/* */
/* Principe : */
/* */
/* On se donne un polynome de */
/* degre N en t, soit : */
/* */
/* k=N */
/* _____ */
/* \ */
/* \ k */
/* P(t) = / a .t */
/* /____ k */
/* */
/* k=0 */
/* */
/* de derivee : */
/* */
/* k=N */
/* _____ */
/* \ */
/* \ k-1 */
/* P'(t) = / k.a .t */
/* /____ k */
/* */
/* k=0 */
/* */
/* Ce polynome va etre defini sur [O,E] */
/* (pour 'Origine' et 'Extremite'), et */
/* l'on se donne les valeurs suivantes : */
/* */
/* P(O) */
/* P'(O) */
/* P(E) */
/* P'(E) */
/* */
/* Quatre indices {I,J,K,L} sont tires */
/* au sort dans [0,N] (en entier...). */
/* Les valeurs des coefficients 'a' pour */
/* des indices autres que {I,J,K,L} sont */
/* tirees au sort : */
/* */
/* a = RANDOM */
/* k */
/* */
/* avec : */
/* */
/* k # {I,J,K,L} */
/* */
/* En ce qui concerne les valeurs des quatre */
/* coefficients 'a' pour les indices {I,J,K,L}, */
/* elles sont calculees de la facon suivante ; */
/* soit : */
/* */
/* P (O) */
/* p */
/* */
/* P'(O) */
/* p */
/* */
/* P (E) */
/* p */
/* */
/* P'(E) */
/* p */
/* */
/* les valeurs respectives de P(O), P'(O) */
/* P(E) et P'(E) en omettant les monomes */
/* de rang {I,J,K,L}. On a alors le systeme */
/* lineaire suivant de 4 equations a 4 */
/* inconnues : */
/* */
/* I J K L */
/* a .O + a .O + a .O + a .O = P(O) - P (O) */
/* I J K L p */
/* */
/* I J K L */
/* a .E + a .E + a .E + a .E = P(E) - P (E) */
/* I J K L p */
/* */
/* I-1 J-1 K-1 L-1 */
/* I.a .O + J.a .O + K.a .O + L.a .O = P'(O) - P'(O) */
/* I J K L p */
/* */
/* I-1 J-1 K-1 L-1 */
/* I.a .E + J.a .E + K.a .E + L.a .E = P'(E) - P'(E) */
/* I J K L p */
/* */
/* qui, resolu avec une simple methode de */
/* Cramer, donnera la valeur des coefficients */
/* 'a' pour les indices {I,J,K,L}... */
/* */
/* */
/* Author of '$xtc/interpolN.03$c' : */
/* */
/* Jean-Francois Colonna (LACTAMME, 1997MMJJhhmmss). */
/* */
/*************************************************************************************************************************************/
#include "INCLUDES.01.I"
/* Introduit le 20051116101250... */
extern double drand48();
#define RANDOM \
(1.0*(drand48()-0.5)) \
/* Les coefficients 'a' sont aleatoires dans [-0.5,+0.5] (sauf ceux d'indice {I,J,K,L}). */
#define FONCTION_ORIGINE \
0
#define DERIVEE_ORIGINE \
0
/* Definition de la fonction Origine. */
#define FONCTION_EXTREMITE \
1
#define DERIVEE_EXTREMITE \
0
/* Definition de la fonction Extremite. */
#define EXTENSION \
1.0 \
/* Extension du segment [FONCTION_ORIGINE,FONCTION_EXTREMITE] pour la visualisation qui */ \
/* aura donc lieu dans [FONCTION_ORIGINE-EXTENSION,FONCTION_EXTREMITE+EXTENSION]. */
#define T0 \
(0.0)
#define T1 \
(+1.0)
#define PAS \
0.01
/* Definition du parametre d'interpolation. */
#define DEGRE_N \
5 \
/* Degre du polynome. */
#define DEGRE_0 \
0
#define DEGRE_1 \
(DEGRE_0+1)
#define DEGRE_2 \
(DEGRE_1+1)
#define DEGRE_3 \
(DEGRE_2+1)
/* Premiers coefficients du polynome. */
int max2(a,b)
int a,b;
{
return(MAX2(a,b));
}
/* Cette fonction est destinee a eviter les effets de bord a prevoir si 'MAX2(...)' est */
/* appele avec comme argument 'drand48()' qui risque alors d'etre appele plusieurs fois */
/* de suite sans redonner (evidemment...) le meme resultat puisqu'il est aletaoire... */
int min2(a,b)
int a,b;
{
return(MIN2(a,b));
}
/* Cette fonction est destinee a eviter les effets de bord a prevoir si 'MIN2(...)' est */
/* appele avec comme argument 'drand48()' qui risque alors d'etre appele plusieurs fois */
/* de suite sans redonner (evidemment...) le meme resultat puisqu'il est aletaoire... */
#define CHOIX(exposant_IJKL,forcer_exposant_IJKL) \
{ \
if (tirer_au_sort_les_exposants == 1) \
{ \
while (selection_IJKL[exposant_IJKL=min2((int)(DEGRE_0+(drand48()*(DEGRE_N+1))),DEGRE_N)] != 0) \
/* On notera l'ecriture 'DEGRE_N+1' et non pas 'DEGRE_N' que la logique reclame, et ce afin */ \
/* de garantir l'accessibilite de 'DEGRE_N'... */ \
{ \
} \
} \
else \
{ \
exposant_IJKL = forcer_exposant_IJKL; \
} \
\
selection_IJKL[exposant_IJKL]++; \
}
extern void *malloc();
static int dimX=0;
#define Xmin 0
#define Xmax (Xmin + (dimX-1))
/* Definition des abscisses. */
static int dimY=0;
#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 store(n,x,y) \
{ \
IMAGE(MAX2(MIN2(x,Xmax),Xmin),MAX2(MIN2(y,Ymax),Ymin)) = n; \
} \
/* Rangement d'un point valide d'une image. */
double puis(x,n)
double x;
int n;
{
int k;
double monome=((n>=0) ? 1 : 0);
for (k=1 ; k<=n ; k++)
{
monome=monome*x;
}
return(monome);
}
/* Calcul de 'x' a la puissance 'n'. */
static double coefficients[DEGRE_N+1];
/* Definition des N+1 coefficients du polynome 'P(t)'. */
double Cpolynome(x,degre)
double x;
int degre;
{
int exposant;
double polynome=0;
for (exposant=degre ; exposant>=DEGRE_0 ; exposant--)
{
polynome = (polynome*x)+coefficients[exposant];
}
return(polynome);
}
/* Calcul de la valeur du polynome en 'x' par la methode de Horner. */
double Cderivee(x,degre)
double x;
int degre;
{
int exposant;
double derivee=0;
for (exposant=degre ; exposant>DEGRE_0 ; exposant--)
{
derivee = (derivee*x)+(exposant*coefficients[exposant]);
}
return(derivee);
}
/* Calcul de la valeur de la derivee du polynome en 'x' par la methode de Horner. */
main()
{
int visualiser=1;
/* #0 : visualiser, =0 : editer les coefficients du polynome. */
int sauter;
int graine=1;
/* En fait indique un nombre d'iterations initiales sur 'drand48()' permettant de faire */
/* generer ainsi des suites de nombres aleatoires differentes. */
int tirer_au_sort_les_exposants=1;
int forcer_exposant_I=DEGRE_N;
int forcer_exposant_J=DEGRE_2;
int forcer_exposant_K=DEGRE_1;
int forcer_exposant_L=DEGRE_0;
/* #0 : tirer au sort les exposants {I,J,K,L}, =0 : les forcer avec les 4 valeurs ci-dessus. */
int tirer_au_sort_les_coefficients=1;
double valeur_commune_des_coefficients=1;
/* #0 : tirer au sort les coefficients des indices differents de {I,J,K,L}, =0 : les forcer */
/* avec la valeur commune 'valeur_commune_des_coefficients'. */
double fO=FONCTION_ORIGINE,dO=DERIVEE_ORIGINE;
double fE=FONCTION_EXTREMITE,dE=DERIVEE_EXTREMITE;
/* Definition des conditions aux limites de la fonction sur [O,E]. */
double t;
/* Variable 't' courante. */
int exposant;
/* Exposant courant. */
int selection_IJKL[DEGRE_N+1];
int exposant_I,exposant_J,exposant_K,exposant_L;
/* Exposants {I,J,K,L} aleatoires dans [0,N] et liste destinee a eviter que deux de ces */
/* exposants ne soient identiques... */
double polynome,derivee;
/* Valeur courante du polynome et de sa derivee. */
double polynome_partiel_O,derivee_partielle_O,polynome_partiel_E,derivee_partielle_E;
/* Valeur des polynomes et derivees "partiels" aux points {O,E}. */
double determinant=0;
/* Valeur du determinant du systeme cramerien. La valeur nulle indiquant qu'il n'a pas */
/* encore ete calcule. */
double membreD1,membreD2,membreD3,membreD4;
/* Valeurs des quatre membres de droite des quatre equations lineaires a resoudre... */
int x,y;
/* Definition des coordonnees de l'image. */
unsigned char *image;
/* Definition de l'image a generer... */
if (DEGRE_N < DEGRE_3)
{
printf("\n cette methode exige au moins le degre %d",DEGRE_3);
}
else
{
}
for (sauter=1 ; sauter<=graine ; sauter++)
{
double a_sauter=drand48();
}
if (visualiser == 0)
{
}
else
{
Get(dimX,"dimX");
Get(dimY,"dimY");
/* Recuperation des dimensions en 'X' et en 'Y' 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++)
{
store(NOIR,x,y);
/* Initialisation de l'image a generer... */
}
}
}
while (determinant==0)
/* Tant que le determinant est nul (et donc en particulier la premiere fois), on va */
/* choisir aleatoirement quatre indices {I,J,K,L}. */
{
for (exposant=DEGRE_0 ; exposant<=DEGRE_N ; exposant++)
{
selection_IJKL[exposant] = 0;
/* Aucun des exposants aleatoires (I,J,K,L) n'a encore ete selectionne... */
}
CHOIX(exposant_I,forcer_exposant_I);
CHOIX(exposant_J,forcer_exposant_J);
CHOIX(exposant_K,forcer_exposant_K);
CHOIX(exposant_L,forcer_exposant_L);
/* Choix aleatoire des indices {I,J,K,L} des quatre coefficients qui sont calcules a partir */
/* des autres... */
for (exposant=DEGRE_0 ; exposant<=DEGRE_N ; exposant++)
{
if ( (exposant!=exposant_I)
&& (exposant!=exposant_J)
&& (exposant!=exposant_K)
&& (exposant!=exposant_L)
)
{
if (tirer_au_sort_les_coefficients == 0)
{
coefficients[exposant]=valeur_commune_des_coefficients;
}
else
{
coefficients[exposant]=RANDOM;
}
/* Initialisation arbitraire de tous les coefficients 'a' du polynome sauf les quatre */
/* d'indices {I,J,K,L} que l'on va calculer ci-apres a l'aide d'un systeme lineaire. */
}
else
{
coefficients[exposant]=0;
/* Ainsi les coefficients d'indices {I,J,K,L} ne participeront pas ci-apres au calcul de */
/* 'polynome_partiel_?' et de 'derivee_partielle_?'... */
}
}
polynome_partiel_O = Cpolynome(T0,DEGRE_N);
polynome_partiel_E = Cpolynome(T1,DEGRE_N);
/* Calcul de la valeur du polynome "partiel" en {O,E} par la methode de Horner. */
derivee_partielle_O = Cderivee(T0,DEGRE_N);
derivee_partielle_E = Cderivee(T1,DEGRE_N);
/* Et de la valeur de la derivee du polynome "partiel" en {O,E} par la methode de Horner. */
determinant = DET4(puis(T0,exposant_I)
,puis(T0,exposant_J)
,puis(T0,exposant_K)
,puis(T0,exposant_L)
,puis(T1,exposant_I)
,puis(T1,exposant_J)
,puis(T1,exposant_K)
,puis(T1,exposant_L)
,exposant_I*puis(T0,exposant_I-1)
,exposant_J*puis(T0,exposant_J-1)
,exposant_K*puis(T0,exposant_K-1)
,exposant_L*puis(T0,exposant_L-1)
,exposant_I*puis(T1,exposant_I-1)
,exposant_J*puis(T1,exposant_J-1)
,exposant_K*puis(T1,exposant_K-1)
,exposant_L*puis(T1,exposant_L-1)
);
/* Calcul du determinant du systeme lineaire a resoudre. */
}
membreD1 = fO-polynome_partiel_O;
membreD2 = fE-polynome_partiel_E;
membreD3 = dO-derivee_partielle_O;
membreD4 = dE-derivee_partielle_E;
/* Calcul des quatre membres de droite. */
coefficients[exposant_I] = DET4(membreD1
,puis(T0,exposant_J)
,puis(T0,exposant_K)
,puis(T0,exposant_L)
,membreD2
,puis(T1,exposant_J)
,puis(T1,exposant_K)
,puis(T1,exposant_L)
,membreD3
,exposant_J*puis(T0,exposant_J-1)
,exposant_K*puis(T0,exposant_K-1)
,exposant_L*puis(T0,exposant_L-1)
,membreD4
,exposant_J*puis(T1,exposant_J-1)
,exposant_K*puis(T1,exposant_K-1)
,exposant_L*puis(T1,exposant_L-1)
)/determinant;
coefficients[exposant_J] = DET4(puis(T0,exposant_I)
,membreD1
,puis(T0,exposant_K)
,puis(T0,exposant_L)
,puis(T1,exposant_I)
,membreD2
,puis(T1,exposant_K)
,puis(T1,exposant_L)
,exposant_I*puis(T0,exposant_I-1)
,membreD3
,exposant_K*puis(T0,exposant_K-1)
,exposant_L*puis(T0,exposant_L-1)
,exposant_I*puis(T1,exposant_I-1)
,membreD4
,exposant_K*puis(T1,exposant_K-1)
,exposant_L*puis(T1,exposant_L-1)
)/determinant;
coefficients[exposant_K] = DET4(puis(T0,exposant_I)
,puis(T0,exposant_J)
,membreD1
,puis(T0,exposant_L)
,puis(T1,exposant_I)
,puis(T1,exposant_J)
,membreD2
,puis(T1,exposant_L)
,exposant_I*puis(T0,exposant_I-1)
,exposant_J*puis(T0,exposant_J-1)
,membreD3
,exposant_L*puis(T0,exposant_L-1)
,exposant_I*puis(T1,exposant_I-1)
,exposant_J*puis(T1,exposant_J-1)
,membreD4
,exposant_L*puis(T1,exposant_L-1)
)/determinant;
coefficients[exposant_L] = DET4(puis(T0,exposant_I)
,puis(T0,exposant_J)
,puis(T0,exposant_K)
,membreD1
,puis(T1,exposant_I)
,puis(T1,exposant_J)
,puis(T1,exposant_K)
,membreD2
,exposant_I*puis(T0,exposant_I-1)
,exposant_J*puis(T0,exposant_J-1)
,exposant_K*puis(T0,exposant_K-1)
,membreD3
,exposant_I*puis(T1,exposant_I-1)
,exposant_J*puis(T1,exposant_J-1)
,exposant_K*puis(T1,exposant_K-1)
,membreD4
)/determinant;
/* Resolution du systeme des qu'il n'est plus indetermnine... */
if (visualiser == 0)
{
if ( (Cpolynome(T0,DEGRE_N) != fO)
|| (Cderivee(T0,DEGRE_N) != dO)
|| (Cpolynome(T1,DEGRE_N) != fE)
|| (Cderivee(T1,DEGRE_N) != dE)
)
{
printf("\n erreur d'interpolation :");
printf("\n fO : attendue=%f calculee=%f",fO,Cpolynome(T0,DEGRE_N));
printf("\n dO : attendue=%f calculee=%f",dO,Cderivee(T0,DEGRE_N));
printf("\n fE : attendue=%f calculee=%f",fE,Cpolynome(T1,DEGRE_N));
printf("\n dE : attendue=%f calculee=%f",dE,Cderivee(T1,DEGRE_N));
}
else
{
}
printf("\n I=%d",exposant_I);
printf("\n J=%d",exposant_J);
printf("\n K=%d",exposant_K);
printf("\n L=%d",exposant_L);
printf("\n\n");
for (exposant=DEGRE_N ; exposant>=DEGRE_0 ; exposant--)
{
printf("+(%g*t^%d)",coefficients[exposant],exposant);
}
/* Edition du polynome selectionne lorsque cela est demande. */
printf("\n");
}
else
{
for (t=T0 ; t<=T1 ; t=t+PAS)
{
int x,y;
polynome = Cpolynome(t,DEGRE_N);
/* Calcul du polynome par la methode de Horner. */
x = Xmin+((Xmax-Xmin)*((t-T0)/(T1-T0)));
y = Ymin+((Ymax-Ymin)*((polynome-(fO-EXTENSION))/((fE+EXTENSION)-(fO-EXTENSION))));
store(BLANC,x,y);
/* Visualisation du polynome selectionne lorsque cela est demande. */
}
write(1,image,dimX*dimY);
/* Sortie de l'image... */
}
}