#include <math.h>
/* Ce programme cherche le maximum d'une fonction de deux variables  */
/* par methode de grimpe
/*     fonc(x,y) : fonction de deux variables, fournie par l'usager  */
/*     epsilon   : niveau de precision recherche                     */
/*     itermax   : nombre maximal d'iteration toleres                */
main() 
{
/* Declarations ---------------------------------------------------- */
  float fonc(float,float) ;      /* La fonction a maximiser */
  int k=0 ;                      /* Compteur d'iterations */
  int itermax=100 ;              /* Nombre maximal d'iterations */
  float epsilon=1.e-5  ;         /* Precision requise pour la racine */       
  float pi=3.1519926536 ;        /* la valeur de pi */       
  float f, fn, dx, dy, grad ;    /* Variables locales */
  float x, y, xn, yn, pas, del ; 
/* Entree/Sortie --------------------------------------------------- */
/* On fournit a l'execution le point de depart (x,y) de la grimpe    */
  printf ("SVP donner une position initiale en x :") ;
  scanf ("%f", &x) ;
  printf ("SVP donner une position initiale en y :") ;
  scanf ("%f", &y) ;
/* Executable ------------------------------------------------------ */
  pas = 0.1 ;                    /* Taille initiale du pas */
  f  = fonc(x,y) ; del=10.*epsilon ;
  printf ("Position initale du grimpeur: x= %f y= %f\n", x,y) ;
  while ( del > epsilon && k < itermax)  /* Boucle de grimpe */
    { 
      dx = fonc(x+pas,y)-fonc(x-pas,y) ; /* composante-x du gradient */
      dy = fonc(x,y+pas)-fonc(x,y-pas) ; /* composante-y du gradient */
      grad = sqrt(dx*dx+dy*dy) ;         /* grandeur du gradient */
      xn = x + pas*dx/grad ;             /* pas dans la direction-x */
      yn = y + pas*dy/grad ;             /* pas dans la direction-x */
      fn = fonc(xn,yn) ;
      if ( fn > f ) {                    /* On accepte le pas */
        x=xn ; y=yn ;                    /* Deplacement du grimpeur */
        del = fabs( fn-f ) ;             /* Pour le test de sortie */
        f=fn ;
        }
      else { pas=pas/2. ; }              /* On refuse le pas */
      k += 1 ;
    }
  printf ("Maximum a: x= %f y= %f\n fonction = %f\n apres %d iterations\n",
    x,y,fn,k) ;
}
/* Definition de la fonction dont on recherche le maximum */
float fonc (float x, float y)
{
  float f ;                              /* Variable locale */
  f = fabs( sin(x)*sin(y)/(x*y) );       /* La fonction qu'on calcule */
  return f ;
}