- MonteCarlo.gif (223.55 Kio) Consulté 5561 fois
http://forum.system-cfg.com/viewtopic.p ... 27#p116811
je me suis décidé à tester la méthode de Monte-Carlo sur le MO5 et en Forth bien entendu.
Analysant le jet d’une pièce sur un carrelage, Buffon disait dans son « Essai d’Arithmétique morale » en 1778 « le sort du premier joueur fera à celui du second, comme la superficie de la couronne circonscrite est à la superficie de la figure inscrite; cela peut se démontrer aisément… »
C’est clair et limpide !
Soit un carré de coté 200, un quart de cercle de rayon 200 ayant pour origine l’un des sommets, un point de coordonnées X Y – X Y compris entre 0 et 200 -, la probabilité de ce point de se trouver à l’intérieur du ¼ de cercle « est comme la superficie de l’un est à la superficie de l’autre »
Même si elle est moins glamour, l’avantage de la méthode de Monte-Carlo est qu’elle ne peut pas être suspectée de triche : PI apparait miraculeusement à partir de la formule de Pythagore qui vérifie si le point est à l’intérieur ou à l’extérieur du cercle. Alors que la simulation informatique de l’aiguille de Buffon fait généralement intervenir des radians et de la trigonométrie que l’on peut accuser de cacher PI.
Malheureusement le langage Forth n’est pas le plus adapté pour ce genre de calculs (avec les décimales en plus). Mais on peut quand même s’amuser un peu.
Wikipedia nous donne les clefs pour la mettre en œuvre.
Comme le Forth travaille avec des entiers l’approximation de PI sera donnée de la façon suivante
PI(approximé)= 10000 * 4 * Nombre de points dans le ¼ de cercle/ nombre d’essais.
Ce programme n’est qu’une ébauche pour 32000 essais mais on pourrait aller plus loin :
- DO … LOOP nous limite à des entiers signés < 32767 on pourrait prévoir un test sur des entiers non signés ramenant à 65536
- La fonction RANDOM est elle vraiment aléatoire ?
- Il y a quelque part une accumulation de données (peut être un reste de division …) mais comme je n’ai pas le temps de faire un debug, j’ai opté pour la solution des paresseux : la remise à plat de la pile SP à chaque LOOP par SP! si quelqu’un le trouve merci de le signaler.
- Les couples X Y dont X2+Y2=R2 sont considérés rouges, donc inclus … Doivent-ils l’être ?
Code : Tout sélectionner
( Début programme)
: TASK ;
( Mode décimal)
DECIMAL
( a b c d . . . c d a b )
: 2SWAP ROT >R ROT R> ;
(a b . . . a b a b)
: DDUP OVER OVER ;
(dun . . . ) (imprime double unsigned : 32 bits non signés)
: DU. SWAP OVER <# #S #> TYPE DROP ;
( d1 n . . . t=d1*n ) (double*simple = triple)
: T* DUP ROT U* >R >R U* 0 R> R> D+ ;
(t n . . . d1= t/n) (Triple/simple=double)
: T/ >R I U/ SWAP ROT 0 I U/ SWAP ROT R> U/ SWAP DROP 0 2SWAP SWAP D+ ;
(ici RND initial = valeur de la variable système HERE – On peut changer par n’importe quelle valeur)
0 VARIABLE RND HERE RND !
(constante double)
: DCONSTANT <BUILDS SWAP , , DOES> DUP @ SWAP 2+ @ ;
( a b . . . b)
: NIP SWAP DROP ;
(Formule aléatoire)
: RANDOM RND @ 31421 * 6927 + DUP RND ! ;
(N . . . 0<Aléatoire<N)
: CHOOSE RANDOM U* NIP ;
(Coordonnées du point)
0 VARIABLE XA
0 VARIABLE YA
(Nombre de points rouges dans le ¼ de cercle)
0 VARIABLE ROUGE
(création de constantes pour gain de temps durant l’exécution)
(Constante simple)
20000 CONSTANT 20K
(Constantes doubles)
40000. DCONSTANT 40K
400000000. DCONSTANT 400M
(100X 100Y . . . ) (Impression du point X/100 Y/100)
: PLOT INKG >R 100 / R> 100 / PSET ;
(dX dY . . . dX2+dY2) (Pythagore)
: DIAGONALE DUP U* ROT DUP U* D+ ;
(. . . . .) (Génère un couple, calcule le carré de sa distance par DIAGONALE compare avec 400M et met à jour ROUGE si
le point est à l’intérieur)
: ESSAI 20K CHOOSE 20K CHOOSE DDUP DIAGONALE DMINUS 400M D+ 0< SWAP DROP 0= IF 1 ROUGE +! 5 PLOT ENDIF ;
(Nbre essais . . .) ( 40000 * Nbre de rouges / nombre d’essais)
: PRINTPI >R 40K ROUGE @ DUP 9 27 LOCATE . T* R> T/ 17 28 LOCATE DU. ;
(Génération de Random)
: GENALEAT CLS 4 9 LOCATE ." LA METHODE DE MONTE-CARLO" 6 5 LOCATE ." DETERMINATION DE LA VALEUR DE PI" 20 15 LOCATE
." DOMINIQUE CONTANT 2016" 11 5 LOCATE ." GENERATION DES VALEURS ALEATOIRES" 100 0 DO LOOP 14 9 LOCATE
." LAISSER UN CERTAIN TEMPS" 15 5 LOCATE ." PUIS APPUYER SUR LA BARRE ESPACE" BEGIN 1 RND +! KEYF 25 = UNTIL ;
(. . . . ) (écran principal)
: PRESENTATION GENALEAT CLS 5 INKG 201 0 319 199 BOX 1 26 LOCATE ." PI:METHODE DE" 2 27 LOCATE ." MONTE CARLO"
4 27 LOCATE ." ===========" 8 26 LOCATE ." POINTS RGES:" 11 26 LOCATE ." ESSAIS:" 14 27 LOCATE ." ==========="
16 26 LOCATE ." VALEUR DE PI:" 21 28 LOCATE ."30/04/2016" ;
( . . . ) (Barre d’espace pour arrêt)
: TOUCHE KEYF 25 = IF BEGIN KEYF 25 = 0= UNTIL KEY DROP ENDIF ;
(Lancement du programme)
: MONTE-CARLO PRESENTATION 10000 1 DO ESSAI I DUP 12 27 LOCATE . PRINTPI TOUCHE SP! LOOP ;- Chargement.jpg (13.56 Kio) Consulté 5561 fois
Edit : Correction fichier vide
Je pense que c'est un générateur "modulo-congruentiel". Cela se voit par l'apparition d'alignements diagonaux quand il y a pas mal de points affichés. Les générateurs "modulo-congruentiels" sont simples et rapides. C'est pourquoi ils étaient pas mal employés dans les langages des années 80. Hélas, en pratique, ils sont assez mauvais à cause de ces alignements qui apparaissent quand on affiche des points au hasard dans une grille 2D. En outre comme tu affiches 32000 points, et que le génétateur a probablement un état interne sur 16bits, il est quasi certain que tu retombes plusieurs fois sur la même séquence de points empêchant de remplir totalement le disk pour avoir la meilleur approximation (problème de cycle court des générateurs modulo-congruentiels).
Une implementation pour forths 16bits ne doit pas être trop compliquée dans la mesure où le générateur ne fait que des décalages à gauche ou a droite et des xor bits à bits facilement adaptable sur des paires adjacentes de cellules 16bits.
