// Conjecture de Syracuse
function u= syracuse(a,n)
    u = a
    for i=1:n
        if modulo(u,2)==0 then
            u = u/2
        else  
            u = 3*u +1
        end
    end 

endfunction

function [maxsuite,nite]= syracusebis(a)
    u = a
    maxsuite = a
    nite = 0

    while u<>1 
        if modulo(u,2)==0 then
            u = u/2
        else  
            u = 3*u +1
        end
        maxsuite = max(maxsuite,u)
        nite = nite +1
    end 

endfunction

// calcul de la fonction dont on evalue l'aire
   function y = f1(x)
       y = x^(2)
    
endfunction

//  fonction dont on evalue l'aire pour calculer approx de pi
 function y = f2(x)
       y = sqrt(1- x^2)
    
endfunction


// Monte Carlo par esperance
   // calcul d'une integrale//aire
   // n tirages aleatoires 
   function E = mtesperance(n,f)
    E = 0
    Points = []
    Images = []
      for i=1:n
            x = rand(1,1,'uniform')    // tirage aleatoire loi uniforme
            E= E + f(x)        // somme avec les valeurs precedentes
                Points = [Points,x]
                Images = [Images,f(x)]
                
      end
      E = E/n
      plot2d(Points,Images,style=-1,  frameflag=6,axesflag=5)
      X = [0:0.001:1]
      Y = f(X)
      plot2d(X,Y,2)
      legends(['nuage de points'    'fonction'],[-1,2],2) 
endfunction


// Monte Carlo par rejet
   // calcul d'une integrale//aire
   // n tirages aleatoires d'un point (x,y) compris dans [0,1]*[0,1]
   function E = mtrejet(n,f)
    Ndessous = 0
    Xdessus = []
    Ydessus = []
     Xdessous = []
    Ydessous = []
      for i=1:n
            x = rand(1,1,'uniform')    // tirage aleatoire loi uniforme
            y = rand(1,1,'uniform')    // tirage aleatoire loi uniforme    
            if f(x)>y then  // si le point se trouve en dessous de la courbe de la fonction considérée
                Ndessous = Ndessous+1
                Xdessous = [Xdessous,x]  
            Ydessous = [Ydessous,y]    
        else
            Xdessus = [Xdessus,x]  
            Ydessus = [Ydessus,y]    
            end
              
      end
      
      E = Ndessous/n
      plot2d(Xdessus,Ydessus,style=-4,  frameflag=6,axesflag=5)
      plot2d(Xdessous,Ydessous,style=-5,  frameflag=6,axesflag=5)
      X = [0:0.001:1]
      Y = f(X)
      plot2d(X,Y,2)
    
  legends(['nuage de points au dessus'  'nuage de points en dessous'  'fonction'],[-4,-5,2],2) 
   
endfunction




