/* 
   ESI CS 2506. Correction de la feuille 3 
   Pour compiler ce fichier, il faut faire appel a la librairie mathematique avec l'option -lm
   Credit : certains des programmes presentes ont ete realises par Pierre Cohort et/ou Boris Velikson
*/

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>

/* La macro suivante simplie l'ecriture de irreductibleQ dans l'exerice V */
#define min(a,b) ((a)<(b)?(a):(b))



/* Exercice I */

/* on implemente d'abord la fonction f a partir de laquelle la
   suite est definie.*/

long f(int n)
{
  if (n%2==0)  /* n%2==0 equivaut a n pair */
    return n/2;
  else  return 3*n+1;
}   

long Syracuse(long u0)
     /* Cette fonction calcule les termes d'une suite de syracuse
	jusqu'a verifier la conjecture ou atteindre un indice maximal
	n_max fixe a l'avance. 
	ENTREE : u0 la condition initiale
	SORTIE : min(n(u0),n_max) ou n(u0) est le plus petit k pour lequel uk=1
     */
{
  const long n_max=1000000;
  long u=u0, n=0;
  while ((u!=1)&&(n<n_max)) {
    u=f(u);
    n++;
  }
  return n;
}  

void Collatz (long n)
     /* Test la conjecture pour u0 dans {1,...,n} 
	ENTREE : n 
     */
{
  long i;
  for (i=1;i<=n;i++){
    printf("u0 = %ld ; u%ld = 1\n",i,Syracuse(i));
  }
}

/* commentaire : il est interessant de visualiser la fonction Syracuse.
   Pour cela, on ecrit une fonction qui stocke les couples (u0,Syracuse(u0))
   dans un fichier. 

   Les fichiers seront vu plus tard dans le cours, on anticipe un peu 
*/

void writeFile(long n) 
     /* Test la conjecture pour u0 dans {1,...,n} 
	ENTREE : n 
	SORTIE : les couples (u0,Syracuse(u0)) dans le fichier syracuse.dat   

	Verifier que le programme a bien fonctionne en examinant
	le contenu de syracuse.dat avec la commande linux
	more syracuse.dat, ou en l'ouvrant avec emacs (dans ce cas, attention
	a ne pas modifier le fichier). Vous pouvez visualiser les resultats 
	graphiquement en utilisant gnuplot (lancer gnuplot puis, sous l'invite 
	de commandes gnuplot, tapez plot "syracuse.dat"), ou en utilisant
	matlab/scilab (en scilab : x=fscanfMat("syracuse.dat") puis 
	fonctions graphiques appliquees a x (voir aide du logiciel)
     */
{
  FILE* file=NULL;
  int i;

  file=fopen("syracuse.dat","w"); /* cree le fichier syracuse.dat */
  if (file==NULL){ return; } /* arret si probleme d'ouverture du fichier */ 
  
  for (i=1;i<=n;i++){
    /* ecriture dans le fichier syracuse.dat des couples (u0,Syracuse(u0))
       pour uo dans {0,...,n} */
    fprintf(file,"%d %ld\n",i,Syracuse(i));
  }
  fclose(file); /* fermeture du fichier */
}
    
int main_exercice1() {
  Collatz(1000);
  writeFile(1000);
  return 0;
}



/* Exercice II */

double g(double x)
{
  return x*x;
}

double point_gauche(double a, double b, double h)
     /* Integration par la mathode des rectangles point gauche
	ENTREE : les bornes de l'intervalle a et b (a<b)
	         pas, le pas d'integration
	SORTIE : une valeur approchee de l'integrale de f entre a et b
     */
{
  double x=a, I=0;

  while (x+h<=b) {
    I += g(x)*h;
    x += h;
  }
    
  return I;
}


double point_droit(double a, double b, double h)
     /* Integration par la mathode des rectangles point droit
	ENTREE : les bornes de l'intervalle a et b (a<b)
	         pas, le pas d'integration
	SORTIE : une valeur approchee de l'integrale de f entre a et b
     */
{
  double x=a, I=0;

  while (x+h<=b) {
    I += g(x+h)*h;
    x += h;
  }
    
  return I;
}


double point_milieu(double a, double b, double h)
     /* Integration par la mathode des rectangles point milieu
	ENTREE : les bornes de l'intervalle a et b (a<b)
	         pas, le pas d'integration
	SORTIE : une valeur approchee de l'integrale de f entre a et b
     */
{
  double x=a, I=0;

  while (x+h<=b) {
    I += g(x+.5*h)*h;
    x += h;
  }

  return I;
}


double trapeze(double a, double b, double h)
     /* Integration par la mathode des trapezes
	ENTREE : les bornes de l'intervalle a et b (a<b)
	         pas, le pas d'integration
	SORTIE : une valeur approchee de l'integrale de f entre a et b
     */
{
  double x=a, I=0;

  while (x+h<=b) {
    I += .5*(g(x)+g(x+h))*h;
    
    /* Une programmation plus astucieuse permettrait d'aller environ deux
       fois plus vite car g(x+h) a cette iteration est g(x) a l'iteration
       suivante */

    x += h;
  }

  return I;
}


double trapeze_optimise(double a, double b, double h)
     /* Integration par la mathode des trapezes
	ENTREE : les bornes de l'intervalle a et b (a<b)
	         pas, le pas d'integration
	SORTIE : une valeur approchee de l'integrale de f entre a et b
     */
{
  double x=a+h, aux1=g(a), aux2=g(a+h), I=0;
  
  /* Programmation plus astucieuse */

  do {
    I+=(aux1+aux2)*0.5*h;
    aux1=aux2;
    x+=h;
    aux2=g(x);
  } while (x<b);

  return I;
}


void integrales (double a, double b, double h)
{
  printf("Calculs pour h=%f\n",h);
  printf("Approximation par rectangle point gauche : \t%f\n",point_gauche(a,b,h));
  printf("Approximation par rectangle point droit :  \t%f\n",point_droit(a,b,h));
  printf("Approximation par rectangle point milieu : \t%f\n",point_milieu(a,b,h));
  printf("Approximation par trapeze :                \t%f\n",trapeze(a,b,h));
  printf("Approximation par trapeze optimise :       \t%f\n",trapeze_optimise(a,b,h));
  printf("\n");
}


int main_exercice2()
{
  double a=0, b=5;
  printf("Integrale de g entre %f et %f\n",a,b);
  printf("Resultat theorique :                       \t%f\n\n",(b*b*b-a*a*a)/3.);
  integrales(a,b,1);
  integrales(a,b,0.1);
  integrales(a,b,0.01);

  return 0;
}



/* Exercice III */

int NombreDeFaces(int n)
     /* Simule le lancer de n pieces 
	ENTREE : nombre de lancer n
	SORTIE : le nombre de "face" */
{
  int i, faces;
  faces=0;
  for (i=1;i<=n;i++)
    if (random()%2==1) /* random() fourni un nombre aleatoire */
      faces++; /* si ce nombre est pair con considere qu'on a pile sinon qu'on a face */
  return faces;
}

int experience_exercice3(int TypeExperience)
     /* Realise l'experience
	ENTREE : le type d'experience :
                 1 pour "3 faces exactement"
                 2 pour "3 faces au plus"
                 3 pour "3 faces au moins" 
	SORTIE : Elle retourne 1 si l'experience est reussie et 0 sinon */
{
  int res;

  assert ((TypeExperience>=1)&&(TypeExperience<=3));

  switch (TypeExperience) {
  case 1: res=NombreDeFaces(10)==3; break;
  case 2: res=NombreDeFaces(10)<=3; break;
  case 3: res=NombreDeFaces(10)>=3; break;
  }

  return(res);
}

int main_exercice3()
{
  int nombretotal=1e6;
  int i, j, reussi;

  for (i=1;i<=3;i++) {
    printf("%ie experience : ",i);
    reussi=0;
    for (j=1;j<=nombretotal;j++)
      reussi+=experience_exercice3(i); /* On compte les experiences reussies */
    printf("la probabilite est %0.0f\n",100.0*reussi/nombretotal);
  }

  return 0;
}

/* Resultas trouves en MF 100 : 
   1e experience : 15/128  (approx 0.12)
   2e experience : 11/64   (approx 0.17)
   3e experience : 121/128 (approx 0.95) */


/* Exercice IV - Programmation avec des structures */

struct perm {
	int p1,p2,p3,p4;
};

typedef struct perm perm;


int max4(int a,int b, int c, int d) {
	int m= 1;
	if  (a>m) m=a;
	if  (b>m) m=b;
	if  (c>m) m=c;
	if  (d>m) m=d;
	return m;
}


int min4(int a,int b, int c, int d) {
	int m= 4;
	if  (a<m) m=a;
	if  (b<m) m=b;
	if  (c<m) m=c;
	if  (d<m) m=d;
	return m;
}


perm creeP() {
	perm p;
	int a,b,c,d;
	do {
		printf(" tapez les images des 4 elements separes par des espaces "); 
		scanf("%d %d %d %d",&a, &b, &c, &d);
		if ((max4(a,b,c,d)>4) ||(min4(a,b,c,d)<1) ){
			printf("au moins un des nombres est hors limites [1..4]\n");	
			continue;
		}
		else if (a*b*c*d!=24) {
			printf("ce n'est pas bijectif!\n");
			continue;
		}
	} while (	(max4(a,b,c,d)>4) ||(min4(a,b,c,d)<1)||(a*b*c*d!=24));
	p.p1=a; p.p2=b; p.p3=c;p.p4=d;	
	return(p);
}


void afficheP(perm p) {
	printf("1\t2\t3\t4\n");
	printf("%d\t%d\t%d\t%d\n",p.p1,p.p2,p.p3,p.p4);
}


perm compP(perm q,perm p) {
	perm r;

	if (p.p1==1)  r.p1=q.p1;
	else if (p.p1==2)  r.p1=q.p2;
	else if (p.p1==3)  r.p1=q.p3;
	else  r.p1=q.p4;

	if (p.p2==1)  r.p2=q.p1;
	else if (p.p2==2)  r.p2=q.p2;
	else if (p.p2==3)  r.p2=q.p3;
	else  r.p2=q.p4;

	if (p.p3==1)  r.p3=q.p1;
	else if (p.p3==2)  r.p3=q.p2;
	else if (p.p3==3)  r.p3=q.p3;
	else  r.p3=q.p4;

	if (p.p4==1)  r.p4=q.p1;
	else if (p.p4==2)  r.p4=q.p2;
	else if (p.p4==3)  r.p4=q.p3;
	else  r.p4=q.p4;

	return (r);
}


int main_exercice4s() {
	perm pp, qq;
	pp=creeP();
	qq=creeP();
	afficheP(compP(qq,pp));	
	return 0;
}


/* Exercice IV - Programmation avec des tableaux */

typedef int* permutation;
/* Les elements consideres sont 1, 2, 3, 4 (classes d'equivalence) 
   L'indice i correspond a l'element i+1 */


permutation CreerPermutation()
/* Lit une permutation et la retourne */ 
{
  int i, j, x, entree_correcte;
  permutation P=(permutation)malloc(4*sizeof(int));

  for (i=0;i<4;i++) {

    do{

      /* Lecture de l'image de i */
      printf("Entrer l'image de %d \n",i+1);
      scanf("%d",&x);

      /* On verifie que cette entree ne conduit pas a un defaut d'injectivite */
      entree_correcte=1;

      if (i>0)
	for (j=0;j<i;j++)
	  if (P[j]==x-1) entree_correcte=0;

      if ((x<0)||(x>4)) entree_correcte=0;

    } while (!entree_correcte);

    P[i]=x-1;
  }

  return P;
}


void AfficherPermutation(permutation P)
/* Affiche la permutation donnee en parametre */
{
  int i;

  for (i=0;i<4;i++)
    printf("%d -> %d\n",i+1,P[i]+1);
}


permutation CompositionPermutation(permutation P1, permutation P2)
/* Composition des permutations P1 et P2, retourne P1 o P2 */ 
{
  int i;
  permutation P=(permutation)malloc(4*sizeof(int));

  for (i=0;i<4;i++)
    P[i]=P1[P2[i]];
  
  return(P);
}


int main_exercice4t()
{
  permutation P1, P2, P3;

  printf("Entrer une permutation P1 :\n");
  P1=CreerPermutation();

  printf("Entrer une permutation P2 :\n");
  P2=CreerPermutation();

  printf("P1 o P2 =\n");
  P3=CompositionPermutation(P1,P2);
  AfficherPermutation(P3);

  return 0;
}


/* Exercice V */

typedef struct rationnel { int N, D; } Q;

Q creerQ (int num, int den)
     /* ENTREE : le numerateur num et le denominateur den!=0
        SORTIE : le nombre rationnel num/den */
{
  Q resultat;
  assert(den!=0);
  resultat.N=num;
  resultat.D=den;
  return(resultat);
}

void afficherQ (Q r)
     /* ENTREE : le nombre rationnel r
        Affichage de r sur le standard output */
{
  printf("%d/%d",r.N,r.D);
  /* Le buffer n'est pas flushe pour permettre l'ecriture de plusieur */
  /* nombres rationnels sur une meme ligne, si necessaire. */
}

double approximerQ (Q r)
     /* ENTREE : le nombre rationnel r
        SORTIE : une approximation reele de r */
{
  return((double)r.N/r.D); 
  /* Attention a ne pas oublie le type casting  */
  /* pour ne pas faire une division entiere */
}

int egalQ (Q r1, Q r2)
     /* ENTREE : deux nombres rationnels r1 et r2
        SORTIE : 1 si r1=r2 et 0 sinon */
{
  return(r1.N*r2.D==r1.D*r2.N);
  /* Attention a bien faire les produits croises et ne pas simplement */
  /* comparer numerateurs et denomniateurs (2/3=4/6) */
}

Q produitQ (Q r1, Q r2)
     /* ENTREE : deux nombres rationnels r1 et r2
        SORTIE : r1*r2 */
{
  Q resultat;
  resultat.N=r1.N*r2.N;
  resultat.D=r1.D*r2.D;
  return(resultat);
}

Q divisionQ (Q r1, Q r2)
     /* ENTREE : deux nombres rationnels r1 et r2!=0
        SORTIE : r1/r2 */
{
  Q resultat;
  assert(r2.N!=0);
  resultat.N=r1.N*r2.D;
  resultat.D=r1.D*r2.N;
  return(resultat);
}

Q sommeQ (Q r1, Q r2)
     /* ENTREE : deux nombres rationnels r1 et r2
        SORTIE : r1+r2 */
{
  Q resultat;
  resultat.D=r1.D*r2.D; /* on met au meme denominateur */
  resultat.N=r2.N*r1.D+r1.N*r2.D;
  return(resultat);
}

Q irreductibleQ (Q r)
     /* ENTREE : un nombre rationnel r
        SORTIE : r sous forme irreductible */
{
  int n;

  /* On commence par regarder si 2 est un facteur commun 
     au numerateur et au denominateur */
  n=2; 

  while (n<=min(r.N,r.D)) {
    if ((r.N%n==0)&&(r.D%n==0))
      {
	/* si c'est le cas on simplifie */
	r.N/=n;
	r.D/=n;
      }
    /* sinon on passe a l'entier suivant  */
    /* (le nombre premier suivant serait suffisant) */
    else n++;
  }

  return r;
}

int main_exercice5()
{
  Q r, r1, r2, r3, r4;

  r1=creerQ(1,2);
  r2=creerQ(1,4);
  r3=creerQ(2,9);
  r4=creerQ(1,7);
  r=divisionQ(sommeQ(r1,r2),sommeQ(r3,r4));
  r=irreductibleQ(r);
  printf("r=");
  afficherQ(r);
  printf(" et sa valeur approchee est %f\n",approximerQ(r));
  return 0;
}



/* main */

int main()
{
  int n;

  printf("Feuille d'exercices 3\n");
  printf("Choix de l'exercice : ");
  scanf("%d",&n);

  switch (n) {
  case 1: main_exercice1(); break;
  case 2: main_exercice2(); break;
  case 3: main_exercice3(); break;
  case 4: main_exercice4s(); main_exercice4t(); break;
  case 5: main_exercice5(); break;
  }

  return 0;
}