/*************************************************************************************************************************************/ /* */ /* P A R A D O X E D E B E R T R A N D : */ /* */ /* */ /* */ /* Author of '$xtc/Paradoxe.01$c' : */ /* */ /* Jean-Francois COLONNA (LACTAMME, 20041211205842). */ /* */ /*************************************************************************************************************************************/ #include "INCLUDES.01.I" /* Introduit le 20051116094850... */ extern double drand48(); extern double cos(); extern double sin(); extern double tan(); extern double sqrt(); #define ITERATIONS \ 10000000 #define RAYON_DU_CERCLE \ 1.0 /* Definition du cercle. */ #define X_DOMAINE \ 8.0 #define Y_DOMAINE \ 8.0 #define RAYON_DU_DOMAINE \ 8.0 /* Definition du domaine a priori circulaire (si ce parametre 'RAYON_DU_DOMAINE' est */ /* inferieur a 'min(abs(X_DOMAINE),abs(Y_DOMAINE))'. Plus ce domaine est grand par rapport */ /* au disque 'RAYON_DU_CERCLE' puis la probabilite s'approche de 0.5... */ #define RANDOM(inf,sup) \ ((((sup)-(inf))*drand48())+(inf)) main() { double cote_du_triangle_equilateral_inscrit_dans_le_cercle=2*RAYON_DU_CERCLE*sin(PI/3); /* Seuil de 'longueur' qui est la longueur du cote du triangle equilateral inscrit dans */ /* le cercle... */ int n; /* Index des iterations... */ int nombre_total=0; int nombre_de_cas=0; /* Pour l'analyse frequentielle... */ if (MIN2(ABSO(X_DOMAINE),ABSO(Y_DOMAINE)) < RAYON_DU_DOMAINE) { printf("ATTENTION : le domaine n'est pas circulaire...\n"); } else { } if (RAYON_DU_CERCLE > RAYON_DU_DOMAINE) { printf("ATTENTION : le domaine est inclus dans le cercle...\n"); } else { } for (n=1 ; n<=ITERATIONS ; n++) { double x=RANDOM(-X_DOMAINE,+X_DOMAINE); double y=RANDOM(-Y_DOMAINE,+Y_DOMAINE); /* Un point arbitraire de la droite... */ if (sqrt(EXP2(x) + EXP2(y)) < RAYON_DU_DOMAINE) /* Le domaine est a priori considere circulaire... */ { double theta=RANDOM(-PI/2,+PI/2); /* Orientation de la droite... */ double a; double b; /* Equation de la droite... */ double A; double B; double C; double delta; /* Coefficients et discriminant de l'equation du second degre d'intersection entre la */ /* droite et le cercle... */ a = tan(theta); b = y - a*x; /* Equation de la droite... */ A = 1 + EXP2(a); B = 2*a*b; C = EXP2(b) - EXP2(RAYON_DU_CERCLE); delta = EXP2(B) - 4*A*C; /* Etude de l'intersection entre la droite et le cercle... */ if (delta > 0) /* Seul le cas strictement positif est interessant car il faut evidemment deux points... */ { double x1,y1; double x2,y2; /* Points d'intersection... */ double longueur_de_la_corde; /* Longueur du segment {{x1,y1},{x2,y2}}. */ nombre_total++; /* Comptage des cas ou la droite coupe le cercle. */ x1 = (-B - sqrt(delta)) / (2*A); y1 = a*x1 + b; x2 = (-B + sqrt(delta)) / (2*A); y2 = a*x2 + b; longueur_de_la_corde = sqrt(EXP2(x2-x1) + EXP2(y2-y1)); if (longueur_de_la_corde > cote_du_triangle_equilateral_inscrit_dans_le_cercle) { nombre_de_cas++; /* Comptage des cas ou la droite coupe le cercle de facon telle que la corde soit plus */ /* grande que le cote du triangle equilateral inscrit dans le cercle... */ } else { } } else { } } else { } } printf("cote....... = %f\n",cote_du_triangle_equilateral_inscrit_dans_le_cercle); printf("cercle..... = %f\n",RAYON_DU_CERCLE); printf("domaine.... = [%+f,%+f]x[%+f,%+f] (rayon=%f)\n",-X_DOMAINE,+X_DOMAINE,-Y_DOMAINE,+Y_DOMAINE,RAYON_DU_DOMAINE); printf("probabilite = %f (=%d/%d)\n",((double)nombre_de_cas)/((double)nombre_total),nombre_de_cas,nombre_total); /* Evidemment, si le domaine est inscrit dans le cercle et deux fois plus petit (rayon=0.5), */ /* alors la probabilite est egale a 1... */ }