/*
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 */
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_exercice3()
{
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;
}
/* Exercice IV */
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()
/* Realise l'experience "3 faces exactement"
SORTIE : Elle retourne 1 si l'experience est reussie et 0 sinon */
{
int res;
res=NombreDeFaces(10)==3;
return(res);
}
int main_exercice4()
{
int nombretotal=1e6;
int i, j, reussi;
reussi=0;
for (j=1;j<=nombretotal;j++)
reussi+=experience(i); /* On compte les experiences reussies */
printf("la probabilite est %.2f\n",(float)reussi/nombretotal);
return 0;
}
/*
Resultat trouve en MF 1302 :
15/128 (approx 0.12)
*/
/* Exercice V */
/* Type fonctions homographiques */
typedef struct H_struct { double a, b, c, d; } H;
H CreerH (double a, double b, double c, double d)
/* Creation d'une fonction homographique
ENTREE : les coefficients a, b, c, d
SORTIE : la fonction homographique correspondante
*/
{
H f;
f.a = a;
f.b = b;
f.c = c;
f.d = d;
return f;
}
int EstAffineH (H f)
/* Predicat indiquant si f est affine
ENTREE : la fonction homographique f
SORTIE : 1 si f est affine, 0 sinon
*/
{
return((f.b==0)&&(f.d==0));
}
double PointSingulierH (H f)
/* Determination du point singulier d'une fonction homographique non constante
ENTREE : une fonction homographique f non constante
SORTIE : le reel ou elle n'est pas defini
*/
{
assert(f.c!=0);
return(-f.d/f.c);
}
double CalculerH (H f, double x)
/* Evaluation d'une fonction homographique en un point
ENTREE : la fonction f et le reel x
SORTIE : f(x)
*/
{
double numerateur = f.a*x+f.b;
double denominateur = f.c*x+f.d;
assert(denominateur!=0);
return(numerateur/denominateur);
}
H InverseH (H f)
/* Calcul de la reciproque d'une fonction homographique
ENTREE : la fonction homographique f
SORTIE : son inverse
*/
{
H finv;
finv.a = -f.d;
finv.b = f.b;
finv.c = f.c;
finv.d = -f.a;
return finv;
}
H ComposeH (H f1, H f2)
/* Calcul de la composee de deux fonctions homographiques
ENTREE : f1 et f2
SORTIE : f1 o f2
*/
{
H fcomp;
double a1 = f1.a;
double b1 = f1.b;
double c1 = f1.c;
double d1 = f1.d;
double a2 = f2.a;
double b2 = f2.b;
double c2 = f2.c;
double d2 = f2.d;
fcomp.a = a1*a2+b1*c2;
fcomp.b = a1*b2+b1*d2;
fcomp.c = c1*a2+d1*c2;
fcomp.d = c1*b2+d1*d2;
return fcomp;
}
int EstEgalH (H f1, H f2)
/* Predicat d'egalite de deux fonctions homographiques
ENTREE : les deux fonctions f1 et f2
SORTIE : 1 si les fonctions sont egales, 0 sinon
*/
{
double a1 = f1.a;
double b1 = f1.b;
double c1 = f1.c;
double d1 = f1.d;
double a2 = f2.a;
double b2 = f2.b;
double c2 = f2.c;
double d2 = f2.d;
return ((a1*c2==a2*c1)&&(b1*c2+a1*d2==a2*d1+c1*b2)&&(b1*d2==b2*d1));
}
int main_exercice5()
{
H f1, f2, g;
double a=1, b=2, c=3, d=4;
double x=2.1;
f1 = CreerH(a,b,c,d);
f2 = CreerH(a,0,c,0);
if (EstAffineH(f1)) printf("f1 est affine Affine\n");
if (EstAffineH(f2)) printf("f2 est affine Affine\n");
if ((f1.b!=0)&&(f1.d!=0)) printf("f1 n'est pas definie en %f\n",PointSingulierH(f1));
printf("f1(%f)=%f\n",x,CalculerH(f1,x));
f2 = InverseH(f1);
g = ComposeH(f1,f2);
if (EstEgalH(g,CreerH(1,0,0,1))) printf("Tout va bien\n");
else printf("Il y a un probleme !\n");
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_exercice4(); break;
case 5: main_exercice5(); break;
}
return 0;
}