function [x1,f1] = recherche_section_doree(x0, f0, d, f, pas, ...
					   tol, tolbis, verb)
//     x0 : point de départ de la recherche
//     f0 : valeur de f en x0
//      d : direction de recherche
//      f : fonction à minimiser dans la direction d (g(t) = f(x0 + t*d))
//    pas : pas initial pour la phase 1
//    tol : tolérance sur l'intervalle [a,c]
// tolbis : tolérance sur les variations de g
//   verb : booléen (vrai pour afficher des informations)

  tau = 0.5*(1 + sqrt(5))
  tau2 = tau^2
  
  // 1- phase de recherche du triplet (a,b,c)
  if verb then printf("\n recherche lin phase 1"), end
  a = 0; ga = f0; c = pas; gc = f(x0 + c*d);
  if gc >= ga then
     while %t
	b = c/tau2; gb = f(x0 + b*d);
	if gb < ga then, break, end
	c = b; gc = gb
     end
  else // gc < ga
     b = c; gb = gc; 
     while %t
	c = a + (b-a)*tau2; gc = f(x0 + c*d);
	if gc > gb then, break, end
	a = b; ga = gb
	b = c; gb = gc
     end
  end
  
  // 2 - phase de réduction de l'intervalle par la section dorée
  if verb then printf("\n recherche lin phase 2"), end
  while  c-a > tol  |  max(ga,gc)-gb > tolbis*max(1,abs(gb))
     delta = c + a - b
     gdelta = f(x0 + delta*d);
     if gdelta < gb then
	if delta < b then, c = b; gc = gb; else, a = b; ga = gb; end
	b = delta; gb = gdelta
     else
	if delta < b then, a = delta; ga = gdelta; else, c = delta; gc = gdelta; end
     end
  end
  
  // mise à jour sortie
  x1 = x0 + b*d; f1 = gb
endfunction

function [xopt, fopt, gopt, iter] = gradient_pas_fixe(x0, f, grdf, tau, tol, itermax, verb)
   // gradient à pas fixe (sans enregistrement de l'historique)
   iter = 0
   xopt = x0; fopt = f(x0); gopt = grdf(x0); 
   
   while norm(gopt) > tol  & iter < itermax then
      iter = iter + 1
      xopt = xopt - tau*gopt;
      gopt = grdf(xopt)
      if verb then
	 printf("\n iteration = %d", iter) 
	 printf("\n ||g(x)|| = %g", norm(gopt))
	 printf("\n f(x) = %g\n", fopt)
      end
   end
endfunction

function [xopt, fopt, gopt, iter] = gradient_pas_optimal(x0, f, grdf, pas, tol, itermax, verb)
   iter = 0
   xopt = x0; fopt = f(x0); gopt = grdf(x0); 
   
   while norm(gopt) > tol  & iter < itermax then
      iter = iter + 1
      [xopt,fopt] = recherche_section_doree(xopt, fopt, -gopt, f, pas, 1e-7, 1e-10, %f)
      gopt = grdf(xopt)
      if verb then
	 printf("\n iteration = %d", iter) 
	 printf("\n ||g(x)|| = %g", norm(gopt))
	 printf("\n f(x) = %g\n", fopt)
      end
   end
endfunction


function [xopt, fopt, gopt, iter] = gradient_conjugue(x0, f, grdf, pas, tol, itermax, verb)
   iter = 0
   xopt = x0; fopt = f(x0); gopt = grdf(x0); d = -gopt; 
   norm2_gopt = d'*d;   
   
   while sqrt(norm2_gopt) > tol  & iter < itermax then
      iter = iter + 1
      [xopt,fopt] = recherche_section_doree(xopt, fopt, d, f, pas, 1e-7, 1e-10, %f)
      gopt_new = grdf(xopt)
      Beta = (gopt_new - gopt)'*gopt_new / norm2_gopt
      gopt = gopt_new; norm2_gopt = gopt'*gopt; 
      d = -gopt + Beta*d
      if verb then
	 printf("\n iteration = %d", iter) 
	 printf("\n ||g(x)|| = %g", norm(gopt))
	 printf("\n f(x) = %g\n", fopt)
      end
   end
endfunction

function [xopt, fopt, gopt, iter] = newton(x0, f, grdf, pas, tol, itermax, verb)
// x0   : point de départ
// f    : fonction à optimiser
// grdf : gradient de cette fonction
// pas  : le pas a utiliser au début de la recherche unidimensionnelle
//       pour le cas ou on fait une iteration de gradient
// tol  : la tolérance : on sort lorsque || grdf(x) || <= tol
// itermax : le nb max d’itérations
// verb : booléen (si vrai alors on affiche le numero d'itération, la valeur de f,
//        et de la norme du gradient)
//
   iter = 0;
   xopt = x0; fopt = f(x0); gopt = grdf(x0); 
   
   while norm(gopt) > tol & iter < itermax then
      iter = iter + 1;
      // approximation de la matrice hessienne par differences finies
      // en utilisant la fonction derivative
      H = derivative(grdf, xopt);
      d = - H\gopt;
      
      if d'*gopt < 0 then  // d est une direction de descente
	 meth = "newton";  // pour affichage si verb=%t
	 tau = 1;
	 while %t
	    xnew = xopt + tau*d;
	    fnew = f(xnew);
	    if ( fnew < fopt) then
	       break;
	    else
	       tau = tau/2;
	    end
	 end
	 xopt = xnew; fopt = fnew; gopt = grdf(xopt);
	 
      else  // d n'est pas une direction de descente => gradient à pas
            // optimal
	    meth = "grd_opt";   // pour affichage si verb=%t
	    [xopt,fopt] = recherche_section_doree(xopt, fopt, -gopt, f, pas, 1e-7, 1e-10, %f)
	    gopt = grdf(xopt)
      end

      if verb then
	 printf("\n iteration = %d", iter) 
	 printf("\n methode = %s", meth)
	 printf("\n ||g(x)|| = %g", norm(gopt))
	 printf("\n f(x) = %g\n", fopt)
      end

   end
endfunction