/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        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...                                                                      */
                    }
          }



Copyright © Jean-François Colonna, 2021-2021.
Copyright © CMAP (Centre de Mathématiques APpliquées) UMR CNRS 7641 / Ecole Polytechnique, 2021-2021.