// version générique des deux méthodes MC
// conçue pour une dimension arbotraire
// calcule les intégrales I=Int(P,f)
// où P est un pavé  de dimension d
// et f(x) est une fonction sur R^d

funcprot(0);
// fonction constante pour tester
function y=f0(x) // attention! x est un vecteur de dimension d
  y=1;
endfunction

//une autre fonction pour tester
// ici x est de dimension 2
function y=f1(x)
  y=x(1)+x(2);//f(x,y)=x+y
endfunction

// fonction du problème 4
// ici x est de dimension 2
function z=fprob4(x)
  z=(1-x(2))*(2+sin(5*%pi*x(1)));//x(1)=x, x(2)=y
endfunction

// Fonction qui calcule l'intégrale par méthode de moyenne
// dim: dimension de l'espace de x
// a,b: vecteurs de dimension dim
// les bornes du pavé: a: point en bas à gauche, b: en heut à droite
// f: le nom de la fonction à intégrer (son argument doit etre un vecteur de dimension dim)
// N: nombre de points aléatoires à générer

function I=MC_moyenne(dim,a,b,f,N)
  U=rand(dim,N); // N nombres aléatoires distribués uniformément entre 0 et 1
  S=0;

  for i=1:N
    S=S+f(a+(b-a).*U(:,i)); // on cumule la somme pour la moyenne
  end
 
  volumePave=1;
  for i=1:dim
    volumePave=volumePave*(b(i)-a(i));
  end

  I=volumePave*S/N;// intégrale = moyenne * longueur de l'intervalle
endfunction


// Fonction qui calcule l'intégrale par méthode de fléchettes
// a,b: les bornes de l'intervalle
// f: le nom de la fonction à intégrer
// N: nombre de points aléatoires  à générer
// M: un majorant de la fonction  sur l'intervalle 
function I=MC_flechettes(dim,a,b,f,M,N)

// on génère  deux séries de nombres aléatoires uniformes
// U pour les coordonnées x(i) et V pour les coordonnées y(i)

  U=rand(dim,N); 
  V=rand(1,N);
  
  cpt=0;
  for i=1:N
    X=a+(b-a).*U(:,i);
    Y=M*V(i);
    if(Y<=f(X))  // on teste si le point de coordonnées (x(i), y(i)) se trouve sous le graphe de la fonction intégrée
      cpt=cpt+1;
    end
  end
  
// cpt/N  donne une approximation de la probabilité de tomber sous le graphe de la fonction

volumePave=1;
for i=1:dim
  volumePave=volumePave*(b(i)-a(i));
end
  I=M*volumePave*cpt/N;
endfunction

dim=2;
N=1000000;
 
// test avec la fonction constante : volume d'un pavé [0,1]x[2,5]

a=[0 ;2];
b=[1 ;5];
M=1.3;

I1=MC_moyenne(dim,a,b,f0,N);

I2=MC_flechettes(dim,a,b,f0,M,N);

I=3.0;

e1=abs(I-I1);
e2=abs(I-I2);

// pour avoir un affichage console formaté on peut utiliser en scilab
// la fonction printf de C

printf("\n ******test1 - fonction constante********\n valeur exacte : %f\n méthode de moyenne : %f\n méthode de flechetes  : %f\n erreur  de méthode de moyenne %f\n erreur de méthode de fléchettes : %f\n", I,I1,I2,e1,e2);

// test avec une fonction affine simple f(x,y)=x+y
a=[0 ;0];
b=[1 ;1];
M=2.2;

I1=MC_moyenne(dim,a,b,f1,N);

I2=MC_flechettes(dim,a,b,f1,M,N);

I=1.0;

e1=abs(I-I1);
e2=abs(I-I2);

// pour avoir un affichage console formaté on peut utiliser en scilab
// la fonction printf de C

printf("\n ******test2 - fn affine (x+y)********\n valeur exacte : %f\n méthode de moyenne : %f\n méthode de fléchettes  : %f\n erreur de méthode de moyenne %f\n erreur de méthode de fléchettes : %f\n", I,I1,I2,e1,e2);

// test avec la fonc du prob 4
a=[0 ;0];
b=[1 ;1];
M=3.0;

I1=MC_moyenne(dim,a,b,fprob4,N);

I2=MC_flechettes(dim,a,b,fprob4,M,N);

I=1+1/(5*%pi);

e1=abs(I-I1);
e2=abs(I-I2);

// pour avoir un affichage console formaté on peut utiliser en scilab
//  la fonction printf de C

printf("\n ******test3 - fn du probleme 4********\n valeur exacte : %f\n méthode de moyenne : %f\n méthode de fléchettes  : %f\n erreur de méthode de moyenne %f\n erreur de méthode de fléchettes : %f\n", I,I1,I2,e1,e2);

//*****************************************************************
//*****************************************************************
//*****************************************************************
//**************************Problème 5*****************************
//*****************************************************************
//*****************************************************************
//*****************************************************************

// test pour une sphère de rayon 1
function t=test1(x)
  if(x(1)^2+x(2)^2+x(3)^2<1)
    t=1;
  else
    t=0;
  end
endfunction

//test pour une sphère et un cylindre
function t=test2(x)
  if((x(1)^2+x(2)^2+x(3)^2<9)&((x(1)-2)^2+x(2)^2<4)&(x(3)<4)&(x(3)>-2))
    t=1;
  else
    t=0;
  end
endfunction

// Fonction qui calcule l'intégrale  par méthode de fléchettes
// a,b: les bornes de l'intervalle
// f: le nom de la fonction à intégrer
// N: nombre de points aléatoires à générer
// M: un majorant de la fonction sur l'intervalle 
function I=MC_flechettesBis(dim,a,b,test,N)
  
  // on génère "dim" séries de nombres aléatoires uniformes
  // U pour les coordonnées des x(i), 
  // U est une matrice de "dim" lignes et de "N" colonnes
  
  U=rand(dim,N); 

  cpt=0;
  for i=1:N
    X=a+(b-a).*U(:,i);

    if(test(X)==1)  // on teste  si le point de coordonnées (x(i), ) se trouve dans le volume de la fonction intégrée
      cpt=cpt+1;
    end
  end

// cpt/N donne une approximation de la probabilité de tomber sous le graphe de la fonction

  volumePave=1;
  for i=1:dim
    volumePave=volumePave*(b(i)-a(i));
  end
  I=volumePave*cpt/N;
endfunction

a=[-1;-1;-1];
b=[1;1;1];
//N=100000;
dim=3;

I2=MC_flechettesBis(dim,a,b,test1,N);

I=4*%pi/3;

e2=abs(I-I2);

// pour avoir un affichage console formaté on peut utiliser en scilab
//  la fonction printf de C

printf("\n ******test4 - sphère rayon 1********\n valeur exacte : %f\n méthode de fléchettes  : %f\n erreur de méthode de fléchettes : %f\n", I,I2,e2);

a=[-3,-3,-3]';
b=[4;3;4];
//N=1000000;
dim=3;

I2=MC_flechettesBis(dim,a,b,test2,N);

// pour avoir un affichage console formaté on peut utiliser en scilab
// la fonction printf de C

printf("\n ******test5 - sphère rayon 3 et cylindre ********\n méthode de flechetes  : %f \n", I2 );
disp(N)