Forth : La méthode de Monte-Carlo

Dominique · Message par **Dominique** » 01 mai 2016 12:35

@Sam,

Que penses-tu de la solution proposée ici ?

https://groups.google.com/forum/#!topic ... glVTqncYzQ

Elle est facile à mettre en pratique

Code : Tout sélectionner

HEX
65E8 CONSTANT RMULT
: 2VARIABLE <BUILDS SWAP , , DOES> ; 
: 2@ DUP @ SWAP 2+ @ ; 
: 2! DUP >R 2+ ! R> ! ; 3 1 2VARIABLE RLOC
: RNDM RLOC 2@ RMULT U* ROT 0 D+ OVER RLOC 2! ; 
: RANDOM RNDM U* NIP ;

__sam__ · Message par **__sam__** » 01 mai 2016 12:51

Oulà, ca c'est de la vache de code expert avec passage sur la return stack et tout

Bon si je lis à la place le texte de googlegroup, je trouve x[n+1] = x[n] * multiplier + carry mod 2^32. C'est le "+carry" qui change tout et n'en fait pas un générateur modulo-congruentiel. En effet, carry n'est pas une constante, mais la partie haute du cycle précédent.

Je crois que c'est précisément cet algo que j'utilise dans les démos 8bits pour avoir des nombres aléatoires sur 8 bits avec un code riki et très rapide: (extrait de "oh la belle bleue")

Code : Tout sélectionner

***************************************
* Genere un nombre pseudo aleatoire sur
* 8 bit. C'est une version modifiee
* du multiply with carry. La periode
* de la sequence est 31870. Compact
* et rapide, que demander de plus?
***************************************
rnd     ldd     #3*256+249   ! A=X[n] B=mutliplier
        mul                  ! X[n]*multiplier
rnd1    addd    #0           ! + carry
        sta     rnd1+2       ! mise a jour carry
        stb     rnd+1        ! mise a jour X[n+1] pour cycle suivant
        rts                  ! B=resultat

Donc si tu fais la version 16bit elle ne doit pas être mal du tout

Dominique · Message par **Dominique** » 01 mai 2016 13:24

Avec ce nouvel algo je suis allé à la limite : 25748 Rges / 32766 essais -> 3.1432
J'ai fait une copie écran juste avant la fin pour voir la distribution des points.

Je fais tourner maintenant l ancien algo à 32767 essais pour voir à l'image la différence de distribution

Dominique · Message par **Dominique** » 01 mai 2016 14:25

Voici donc comment se présente la distribution des points :

Ancien Algo

Code : Tout sélectionner

0 VARIABLE RND HERE RND ! 
: RANDOM RND @ 31421 * 6927 + DUP RND ! ;

: AncienAlgo32766.jpg (34.63 Kio) Consulté 4836 fois

Nouvel Algo

Code : Tout sélectionner

HEX
65E8 CONSTANT RMULT
: 2VARIABLE <BUILDS SWAP , , DOES> ;
: 2@ DUP @ SWAP 2+ @ ;
: 2! DUP >R 2+ ! R> ! ; 3 1 2VARIABLE RLOC
: RNDM RLOC 2@ RMULT U* ROT 0 D+ OVER RLOC 2! ;

: NouvAlgo32761.jpg (39.18 Kio) Consulté 4836 fois

__sam__ · Message par **__sam__** » 01 mai 2016 16:14

L'un des générateur est visuellement plus aléatoire que l'autre

Pour comparaison, voici une image 200X200 produite par le générateur aléatoire vrai de http://www.random.org (bruit atmosphérque)

(<== l'image change à chaque refresh de page

)

Maintenant l'impact sur la valeur de pi est finalement faible. Ce serait plus prononcé avec un carré avec encore plus de pixels (genre 10000x10000) je pense. En effet, avec un carré de 200x200 et un générateur parfait on couvrirait l'ensemble des pixels assez vite. Le resultat ne sera alors pas trop différent du programme basic qui prend chacun de ces 40000 pixels un à un.

Dominique · Message par **Dominique** » 01 mai 2016 16:38

Mais @Sam je génère 20000*20000 points

Code : Tout sélectionner

20000 CONSTANT 20K
: ESSAI 
          20K CHOOSE                     \X= Absisse entre 0 et 20000
          20K CHOOSE                     \ Y=Ordonnée entre 0 et 20000
          2DUP 
          DIAGONALE                      \ X2+Y2 = (20000)2 + (20000)2
          DMINUS                         \ (X2+Y2)*-1
          400M                           \ 400000000 
          D+                             \ 400000000 - (X2+Y2)
          0<                             \  < 0 ?
                                         \ etc
        SWAP DROP 0= IF  1 ROUGE +! 5 PLOT ENDIF ;

Ce qui me manque c'est le nombre d'essais. Je suis sur les non signés pour passer à 65535 essais
Ce qu'on voit sur l'image n'est que la valeur entière de X/100 et n'entre pas dans le calcul. C'est une sorte de mise à l'échelle de l'écran du MO5

Code : Tout sélectionner

: PLOT INKG >R 100 / R> 100 / PSET ;

Daniel · Message par **Daniel** » 01 mai 2016 16:59

Dominique a écrit :Mauvaise manip ! Excuses. Je vais la remettre

Avec la version pi_monte-carloV11.k7, il y a encore un problème (I/O error). C'est peut-être parce que le fichier .k7 a été copié avant d'être mis à jour par dcmoto. Dans l'émulateur la cassette est chargée en mémoire. Le fichier est mis à jour s'il est déchargé, ou rechargé, ou en cas de réinitialisation de l'ordinateur ou de sortie de dcmoto. Sinon il n'est pas modifié.

: ##.png (1.32 Kio) Consulté 4817 fois

__sam__ · Message par **__sam__** » 01 mai 2016 17:11

Dominique a écrit :Mais @Sam je génère 20000*20000 points

Ah oui! Javais pas bien compris le code

Bon il n'y a pas de problèmes alors. 20K x 20K, il y a de quoi avoir une assez bonne approximation.

Ce qui me manque c'est le nombre d'essais. Je suis sur les non signés pour passer à 65535 essais

Passe les compteurs en 32bits directement, ce sera plus simple pour laisser tourner l'algo très très longtemps.

Dominique · Message par **Dominique** » 01 mai 2016 17:18

Hier, j'ai corrigé le fichier de mon premier message et remis en ligne dans ce premier message.
Je viens de télécharger ce zip avec readme et dezippé dans un dossier vide.
Et j'ai pu lancer la K7 normalement à l'instant . Ref de cette K7 : Modifié le 30/04/2016 23H30 . C'est elle ?

edit : Je vais supprimer cette V11 que j'avais faite dans l'urgence et qui n'a plus lieu d'être après la correction de mon 1° message; Merci en tous cas de faire ces vérifications

Dominique · Message par **Dominique** » 01 mai 2016 17:26

__sam__ a écrit : Ah oui! Javais pas bien compris le code

Non c'est moi. Parmi les les défauts du Forth il en est un majeur : la lisibilité.

C'est au programmeur de faire l'effort d'être assez didactique pour expliquer clairement ce qu'il a fait.
Donc je vais tenter de me corriger.

__sam__ a écrit :Passe les compteurs en 32bits directement, ce sera plus simple pour laisser tourner l'algo très très longtemps.

Le problème c est les routines de division et multiplication que j'ai du mal à dominer. ça me prend énormément de temps. Mais je vais y arriver. Tu as raison, il faut que je passe tout en 32 bits. C'est pour bientôt

Daniel · Message par **Daniel** » 01 mai 2016 17:40

Dominique a écrit :Ref de cette K7 : Modifié le 30/04/2016 23H30 . C'est elle ?

J'avais pris la version postée en réponse à mon premier message. Celle-ci est datée du 30/04/2016 à 23h11, elle est incomplète.
Celle du premier message est datée du 30/04/2016 à 23h30 et tout est bon. Merci !

Dominique · Message par **Dominique** » 03 mai 2016 19:42

Depuis samedi je tentais passer aux simples non signés pour obtenir au moins 65536 essais
Mais j'obtenais des résultats farfelus au delà de 33000 qui m'ont amenés à faire des recherches.
Effectivement la procédure Figforth: U/ souffre d'un bug.

http://6502org.wikidot.com/errata-software-figforth

Fig a donné un correctif et j'ai fini par trouver un U/ (que je teste sous le nom DU/ ici)
Procédure actuelle de notre Forth MO5

U/ : effectue la division non signée avec reste
d'un double non signé par un simple précision
(ud1 u1 . . . u2 reste , u3 quotient)

Code : Tout sélectionner

                     cf_uslash
 4A5C 4A5E                    FDB c_uslash
                      c_uslash
 4A5E 8DAA                    BSR    mod5              
 4A60 0EAC                    JMP    PULL  

...  ...  ...  ... 

                     mod5              
 4A0A DD4F                    STD    l454F               
 4A0C 8E0020                  LDX    #$0020             
 4A0F CC0000                  LDD    #$0000             
 4A12 DD51                    STD    l4551
                      mod3               
 4A14 0850                    ASL    l4550               
 4A16 094F                    ROL    l454F               
 4A18 094E                    ROL    l454E               
 4A1A 094D                    ROL    l454D               
 4A1C 59                      ROLB                      
 4A1D 49                      ROLA                      
 4A1E 0952                    ROL    l4552               
 4A20 0951                    ROL    l4551               
 4A22 10A3C4                  CMPD   ,U                 
 4A25 2504                    BCS    mod2              
 4A27 A3C4                    SUBD   ,U                 
 4A29 0C50                    INC    l4550
                      mod2               
 4A2B 301F                    LEAX   -$01,X             
 4A2D 26E5                    BNE    mod3              
 4A2F ED44                    STD    $04,U              
 4A31 DC4D                    LDD    l454D               
 4A33 EDC4                    STD    ,U                 
 4A35 DC4F                    LDD    l454F               
 4A37 ED42                    STD    $02,U              
 4A39 39                      RTS   
                      mod6                    
 4A3A CC0000                  LDD    #$0000             
 4A3D A302                    SUBD   $02,X              
 4A3F ED02                    STD    $02,X              
 4A41 CC0000                  LDD    #$0000             
 4A44 E201                    SBCB   $01,X              
 4A46 A284                    SBCA   ,X                 
 4A48 ED84                    STD    ,X                 
 4A4A 39                      RTS

Procédure trouvée sur le Net

Code : Tout sélectionner

HEX CREATE DU/ EC42 , AE44 , AF42 , ED44 , 6843 , 6942 , 8E C, 0010 , 6945 , 
           6944 , EC44 , A3C4 , 1CFE , 2B04 , ED44 , 1A01 , 6943 , 6942 , 
           301F , 26E8 , 3342 , 0EB6 , SMUDGE

ce qui donne

Code : Tout sélectionner

69C7 EC42       LDD    $02,U               6
69C9 AE44       LDX    $04,U               6
69CB AF42       STX    $02,U               6
69CD ED44       STD    $04,U               6
69CF 6843       ASL    $03,U               7
69D1 6942       ROL    $02,U               7
69D3 8E0010     LDX    #$0010              3
69D6 6945       ROL    $05,U               7
69D8 6944       ROL    $04,U               7
69DA EC44       LDD    $04,U               6
69DC A3C4       SUBD   ,U                  6
69DE 1CFE       ANDCC  #$FE                3
69E0 2B04       BMI    $69E6               3
69E2 ED44       STD    $04,U               6
69E4 1A01       ORCC   #$01                3
69E6 6943       ROL    $03,U               7
69E8 6942       ROL    $02,U               7
69EA 301F       LEAX   -$01,X              5
69EC 26E8       BNE    $69D6               3
69EE 3342       LEAU   $02,U               5
69F0 0EB6       JMP    /$B6                3

: DU.jpeg (3.57 Kio) Consulté 4784 fois

Si le Forth du MO5 avait été un grand succès commercial
peut être les auteurs auraient-ils fait une procédure corrective comme les Américains.
Y a-t-il lieu de le faire dans notre version en recompilant, peut être en
mettant une procédure patch sous forme de note dans un TXT qui accompagnerait la K7 ?
Est il possible de corriger le bug sans déplacer les adresses ?

Je laisse aux experts parce que je n'ai pas les compétences en LM 6809

__sam__ · Message par **__sam__** » 04 mai 2016 11:20

De ce que je comprends du code ASM original:

"mod6" calcule le négatif du nombre 32bits pointé par X. Cette routine n'est pas utilisée.

"'mod5" est supposé calculer la division d'un 32bits par un 16bit. Si je me souviens bien le sommet de la pile et dans D et le reste sont des offsets par rapport au registre U. Je m'attendrais donc à ce qu'on lise au moins ,U (fait dans le CMPD), mais aussi 2,U et 4,U mais ces deux derniers accès mémoire ne sont jamais lus. Ca me perturbe. Ca ressemble juste à une division 16bits par 16bits. Il y a un truc que je pige pas dans "mod5".

La version "DU/" me semble plus logique car on accède bien à 2,U et 4,U. Par contre le ANDCC #$FE et ORCC #1 sont des portages direct des instructions CLC et SEC du 6502 (positionner la carry à 1 ou 0). Ce que ca fait c'est inverser l'état de la carry. Or on peut faire une division avec la carry inversée. Il suffit de tout complémenter à la fin. C'est ce que fait la division 16bits par 16bits de http://www.logicielsmoto.com/phpBB/view ... 9&start=60 il faudrait juste l'adapter pour du 32bits par 16bits.

Par ailleurs qu'entends tu par "Est il possible de corriger le bug sans déplacer les adresses ?" tu veux dire reprogrammer la zone du cf_uslash pour qu'elle fasse la même chose ? tout dépend de l'occupation mémoire. Si la 2eme version occupe moins d'octets ce sera joueable. Cependant je trouve le code de "mod5" trop bizarre car il ne travaille pas sur 2,U ou 4,U. Est tu certain que c'est la bonne adresse et pas celle de "/" ?

Dominique · Message par **Dominique** » 04 mai 2016 12:27

SAM, je crois que tu as raison ! Erreur de ma part : L'adresse du BSR Est BSR 4A04 :
Je redonne le code d ici peu

Dominique · Message par **Dominique** » 04 mai 2016 12:44

Effectivement j'avais inversé A4 et 4A

La définition de U/ est

Code : Tout sélectionner

                     ;****************
                      ;* LE MOT U/
                      ;****************
                      w_uslash
 4A57 82                      FCB $82 
 4A58 55                      FCC "U"
 4A59 AF                      FCB $80+'/
 4A5A 4A4B                    FDB w_ustar
                  cf_uslash
4A5C 4A5E                    FDB c_uslash
                      c_uslash
4A5E 8DA4                    BSR    mod7             
4A60 0EAC                    JMP    PULL

Ce qui donne ... Ou on retrouve la lecture de la pile de données par LDD $02,U et LDD $04,U

Code : Tout sélectionner

                     mod7              
 4A04 EC42                    LDD    $02,U              
 4A06 DD4D                    STD    l454D               
 4A08 EC44                    LDD    $04,U
                      mod5              
 4A0A DD4F                    STD    l454F               
 4A0C 8E0020                  LDX    #$0020             
 4A0F CC0000                  LDD    #$0000             
 4A12 DD51                    STD    l4551
                      mod3               
 4A14 0850                    ASL    l4550               
 4A16 094F                    ROL    l454F               
 4A18 094E                    ROL    l454E               
 4A1A 094D                    ROL    l454D               
 4A1C 59                      ROLB                      
 4A1D 49                      ROLA                      
 4A1E 0952                    ROL    l4552               
 4A20 0951                    ROL    l4551               
 4A22 10A3C4                  CMPD   ,U                 
 4A25 2504                    BCS    mod2              
 4A27 A3C4                    SUBD   ,U                 
 4A29 0C50                    INC    l4550
                      mod2               
 4A2B 301F                    LEAX   -$01,X             
 4A2D 26E5                    BNE    mod3              
 4A2F ED44                    STD    $04,U              
 4A31 DC4D                    LDD    l454D               
 4A33 EDC4                    STD    ,U                 
 4A35 DC4F                    LDD    l454F               
 4A37 ED42                    STD    $02,U              
 4A39 39                      RTS

Ceci dit la version corrigée ne marche pas pour les valeurs non signées élevées; C à desesperer....

forum.system-cfg.com

Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo

Re: Forth : La méthode de Monte-Carlo