Misez p'tit Optimisez n°94 : les nombres de Motzkin

[Danny]

Sur Casio fx-3800p, 7 secondes pour calculer M(100) (avec l'algo de Gege).

Casio fx-3900p Motzkin.jpg (90.17 Kio) Vu 5934 fois

[C.Ret]

[...]La HP15C Limited Edition donne M(100) en 1.5 seconde.

Avec le code utilisant les deux Cx,y mon HP-15C met un peu plus de 10min pour calculer M(100), c'est long.

Code : Tout sélectionner

LBL A  STO 0   0  STO 1  GTO 1
LBL 0  Cy,x  LSTx  RCL 1  Cy,x  ×  ISG 1  PSE  RCL÷1  +
LBL 1  RCL 0  RCL 1  RCL+1  x≤y?  GTO 0
       R↓  x↔y  RTN

J'ai donc adopté le code proposé par gégé :

Avec le code de gégé qui utilise les registres R1 à R6 mon HP-15C met un peu moins de 2 min pour calculer M(100), c'est bien !

Code : Tout sélectionner

LBL A  STO 1  1  STO 2  STO 5  STO 6  +  1  +  STO 3  1  +  STO 4
LBL 1  RCL 1  RCL 5  RCL+5  x>y?  GTO 2
       RCL 3  RCL÷5  2  -  STO×2  1  STO+5  RCL 4  RCL÷5  2  - STO×2
       RCL 2  STO+6  GTO 1
LBL 2  RCL 6  RTN

Mais cela m'a permis à moi aussi de décomposer les calculs des coefficient binomiaux. A chaque itération, on refait les mêmes multiplications et divisions. L'idée est donc d'utiliser un registre d'accumulation a (R2) dans lequel on mémorise le produit des facteurs provenant de la relation entre n et k. Comme la méthode utilise 2k, on ajoute une fois sur deux la valeur de l'accumulateur a au résultat Mn:
Par exemple, pour M(10) l'itération fait varier R0 de n=10 à 1 et R1 de k=1 à k=5.
Comme n et k interviennent respectivement au numérateur et dénominateur de chaque coefficient, on a les accroissements multiplicatifs suivants dans l'accumulateur a :
a={1 ; 1*10/1=10 ; 10*9/2=45 ; 45*8/2=180 ; 180*7/3=420 ; 420*6/3=840 ; 840*5/4=1050 ; 1050*4/4=1050 ; 1050*3/5=630 ; 630*2/5=252
; 252*1/6=42 }
Et le résultat Mn est obtenu en effectuant la somme M(10) = 1 + 45 + 420 + 1050 + 630 + 42 = 2188.
Ce qui correspond bien à .

Avec le code ci-dessous basé sur l'idée de gégé et utilisant une implémentation n'utilisant que la pile (pour la somme) et les registres R0 R1 et R2 (respectivement n, k et a), la même calculette HP-15C ne met que 1min26" pour afficher le même M(100), c'est un peu mieux !!

Code : Tout sélectionner

LBL A  STO 0  1 x≥y?  RTN  STO 1  STO 2
LBL 1  RCL 2  RCL×0  RCL÷1  ISG 1
LBL 2  DSE 0  GTO 3  R↓  RTN
LBL 3  RCL×0  RCL÷1  STO 2  +  DSE 0  GTO 1  RTN

(ne cherchez pas de GTO 2, l'étiquette LBL 2 peut être remplacer par n'importe quelle instruction, elle est skipée systématiquement par le ISG 1 précèdent)

[Danny]

Avec le code ci-dessous basé sur l'idée de gégé et utilisant une implémentation n'utilisant que la pile (pour la somme) et les registres R0 R1 et R2 (respectivement n, k et a), la même calculette HP-15C ne met que 1min26" pour afficher le même M(100), c'est un peu mieux !!

Code : Tout sélectionner

LBL A  STO 0  1 x≥y?  RTN  STO 1  STO 2
LBL 1  RCL 2  RCL×0  RCL÷1  ISG 1
LBL 2  DSE 0  GTO 3  R↓  RTN
LBL 3  RCL×0  RCL÷1  STO 2  +  DSE 0  GTO 1  RTN

(ne cherchez pas de GTO 2, l'étiquette LBL 2 peut être remplacer par n'importe quelle instruction, elle est skipée systématiquement par le ISG 1 précèdent)

'gnifique !

Et j'avoue que le coup du ISG qui sert juste à incrémenter un registre, je n'y pense jamais... faut que je m'en rappelle !

[C.Ret]

[...]Et j'avoue que le coup du ISG qui sert juste à incrémenter un registre, je n'y pense jamais... faut que je m'en rappelle !

Moi c'est l'inverse, je fais parfois que utiliser les DSE et ISG (respectivement DSZ et ISZ sur les Classic)

Par exemple, cette version qui en abuse avec le même principe de fonctionnement que précédemment, mais plus court:

Code : Tout sélectionner

001- ClrREG  STO 0  1  STO 1
005- LBL 1  STO+2  RCL×0  RCL÷1  DSE 0  ISG 1  CLx  RCL×0  RCL÷1  DSE 0  GTO 1
016- RCL+2

Cette version se lance en tapant n gRTN R/S et affiche le résultat Mn en X: (le registre Y: mémorise la valeur n initiale), les déterminations suivantes peuvent être lancée directement en saisissant n et pressant sur la touche R/S.
Pour M(100) mon HP-15C affiche 7,374155 E 44 au bout de 1min14".

MPO94 - HP15C Motzkin Numbers (v3).gif (16.58 Kio) Vu 5733 fois

[C.Ret]

[...]Je cherche à éviter les conversions. Une bonne idée aussi est de pré-dimensionner les deux vecteurs A et B pour éviter les réarrangement mémoire dans la boucle de calcul.
Malheureusement, il faut à chaque fois inverser l'ordre dans l'un des deux vecteurs. Je bloque !

Débloqué !

Le point bloquant était de vouloir pré-dimensionner les vecteurs. La bonne idée issue des listes était d'ajouter les éléments en début ou en fin de liste.
Il est possible de faire de même avec les vecteurs. Plus besoin de boucle d'inversion, le redimensionnement des deux vecteurs se fait naturellement dans la pile (infinie) du système RPL:

Code : Tout sélectionner

« → n
  « [ 1 ] DUP 1                                           /* Initialise les deux vecteur A:[m0] et B:[1] et le premierterme m1
    IF n 2 ≥ THEN
      2 n FOR i
        ROT ARRY→ DROP i ROLL OVER i 1 + ROLLD i →ARRY    /* Ajoute le nouveau terme à la fin du vecteur A 
        SWAP ROT ARRY→ DROP i →ARRY DUP2 DOT              /* Ajoute l'avant dernier élément de A au début du vecteur B
        DUP2 DOT                                          /* calcule le terme suivant
      NEXT
    END SWAP DROP » »                                     /* Efface vecteur B, donne 1:M(n) et Laisse la vecteur A dans la pile

La pile est initialisée avec [ 1 ] [ 1 ] 1 c'est à dire les vecteurs A:[m0] et B:[1] et le premier terme m1.
Puis le nombre de Motzkin du niveau 1: est ajouté à la fin du vecteur A et l'avant dernier terme du nouveau vecteur A (en réalité l'ancien dernier élément de A) est ajouté au début du vecteur B. Les ajouts aux deux vecteur se faisant en les "étalant" dans la pile par l'instruction ARRY→ et en les regroupant par l'instruction réciproque i →ARRY où i est l'indice de boucle donnant également la taille des vecteurs A et B.

Pour donner une idée des mouvements et jeux de la pile opérationnelle, je donne un aperçu signalétique de deux premières itérations. Pour faciliter la lecture, j'ai repéré les termes m0, m1, m2, etc mais en réalité seules les valeurs numériques sont manipulées par l'HP-28S

(Première Itération i=2)

Code : Tout sélectionner

5:                                                                                     [1]   
4:                              [1]                  [1]                    [1]         m0        [1]       
3:       [m0]        [1]         m1        [1]        m1         [1]         m0         m1         m0        [1]
2:        [1]         m1         m0         m1        m0          m0         m1         m0         m0         m0
1:         m1       [m0]        {1}         m0         2          m1         m0          3         m1    [m0 m1]
             ROT        ARRY→      DROP       i          ROLL       OVER       i 1 +      ROLLD      i →ARRY 
5:
4:                                     [m0 m1]               [m0 m1]               [m0 m1]
3:        [1]        [1]    [m0 m1]         m0    [m0 m1]         m0                [m0 1]    [m0 m1]
2:         m0    [m0 m1]         m0          1         m0          1    [m0 m1]    [m0 m1]     [m0 1]
1:    [m0 m1]         m0        [1]        {1}          1          2     [m0 1]     [m0 1]         m2
             SWAP       ROT        ARRY→      DROP       i          →ARRY      DUP2       DOT        NEXT

(Seconde itération i=3)

Code : Tout sélectionner

6:                                                                                  [m0 1] 
5:                           [m0 1]                [m0 1]                [m0 1]         m0     [m0 1]  
4:                               m2     [m0 1]         m2     [m0 1]         m0         m1         m1  
3:    [m0 m1]     [m0 1]         m0         m2         m0         m0         m1         m2         m0     [m0 1]
2:     [m0 1]         m2         m1         m0         m1         m1         m2         m1         m1         m1
1:         m2    [m0 m1]      { 2 }         m1          3         m2         m1          4         m2 [m0 m1 m2]
             ROT        ARRY→      DROP       i          ROLL       OVER       i 1 +      ROLLD      i →ARRY 
6:
5:                                  [m0 m1 m2]            [m0 m1 m2]
4:                                          m1 [m0 m1 m2]         m1            [m0 m1 m2] 
3:     [m0 1]     [m0 1] [m0 m1 m2]         m0         m1         m0             [m1 m0 1] [m0 m1 m2]
2:         m1 [m0 m1 m2]         m1          1         m0          1 [m0 m1 m2] [m0 m1 m2]  [m1 m0 1]
1: [m0 m1 m2]         m1     [m0 1]        {2}          1          3  [m1 m0 1]  [m1 m0 1]         m3
             SWAP       ROT        ARRY→      DROP       i          →ARRY      DUP2       DOT        NEXT

(et ainsi de suite...)
A la fin des itération, la pile contient les deux vecteurs A et B ainsi que le Nombre de Motzkin souhaité. Les instruction SWAP DROP permettent de ne conserver que le vecteur A contenant les nombres de M(0) à M(n-1) et le résultat souhaité M(n).

Ce qui permet de calculer M(100) en environ 2min24".

C'est loin d'être le plus rapide sur HP-28C/S, le calcul de M(100) prend moins de 3.52sec et M(260) environ 7.96sec avec le code suivant :

Code : Tout sélectionner

« → n « 1 1 1                              // Initialise respectivement Mn j et a 
        n 1 FOR i
            i i 1 - * *                    // Multiplie a par i*(i-1)
            OVER 1 + SWAP ROT / OVER /	   // Divise a par j.(j+1) et laisse (j+1) dans la pile
            ROT OVER + ROT ROT             // Ajoute une copie de a à Mn
        -2 STEP  DROP2 » »                 // Efface j et a pour laisser le résultat Mn

Presque aussi rapide (pour un HP-28S) mais plus court car utilisant la matrice et les les fonctions statistiques:

Code : Tout sélectionner

« CL∑ 1 DUP ∑+                    // Initialise matrice statistique et accumulateur a
  SWAP 1 FOR n
     n DUP 1 - * *                // Multiplie a par n et (n-1)
     N∑ DUP 1 + * /                 // Divise a par k et (k+1)
     DUP ∑+                       // Ajoute copie de a dans matrice statistique
  -2 STEP 
  DROP TOT »                      // Efface dernier a et donne résultat

[C.Ret]

Comme personne ne semble plus utiliser de calculatrice Texas Instruments, je suis obligé :
M(100) approché en 1min38" et 25 pas de programme sur TI-57 LCD

Partition n°3 (c'est à dire 25/32 pas de programme et 3/3 registres R0 R1 et R2 utilisés)

MPO94 Motzkin Numbers TI 57 LCD.gif (35.75 Kio) Vu 5847 fois

On notera l'utilisation de inv DSZ et DSZ. Sur TI-57 LCD, il n'y a pas de fonction ISG, les instructions DSZ et inv DSZ décrémentent toutes les deux le registre R0 jusqu'à atteindre zéro (en réalité elles augmentent ou diminuent la partie entière de R0 d'une unité selon son signe (ou pas si R0 est déjà nul)). La différence entre DSZ et inv DSZ provient de la façon dont se fait le branchement conditionnel à l'issue de la décrémentation; DSZ saute (skipe) le pas suivant si R0 est nul après décrémentation alors que inv DSZ saute le pas suivant si R0 est non-nul après décrémentation.

La version que j'ai posté peut être encore optimisée et il faut aussi corriger un 'bug' que tout le monde, et j'espère tous les utilisateurs de TI, a remarqué ici, lors de l'enchainement de plusieurs déterminations !

[C.Ret]

Sur Casio fx-3800p, 7 secondes pour calculer M(100) (avec l'algo de Gege).

A mince, c'est presque deux fois plus long sur une CASIO fx-602p

MPO94 Motzkin Numbers CASIO fx-602p.gif (17.47 Kio) Vu 5823 fois

[Gilles59]

Bonjour à tous.

C.RET, on peut encore gagner 2 pas dans le programme 602P (ça fait quand même 6,7% lol) :

. Le AC est inutile.
. La séquence ... MR01 x 1 M+01 / MR01 peut s'écrire ... MR01 / 1 M+01 MR01

On peut aussi gagner une petite seconde en plus des 2 pas pour 100 P0 comme suit:

Code : Tout sélectionner

[P0]
 x=0 1 Min00 1 Min01 Min02
LBL0
 x MR00 DSZ GOTO1 GOTO9
LBL1
 x MR00 / MR01 / 1 M+01 MR01 = M+02 DSZ GOTO0
LBL9
 MR02

EDIT : j'avais mal recopié le pgm ;D

On peut aussi gagner 3 pas (soit 10%) comme suit :

Code : Tout sélectionner

[P0]
 x=0 1 Min00 1 Min01 Min02
LBL0
 x MR00 DSZ x=0 GOTO9
 x MR00 / MR01 / 1 M+01 MR01 = M+02 DSZ GOTO0
LBL9
 MR02

Ce qui fait 26 pas (ou plutôt 27 car chaque zone programme utilisée compte pour 1 pas)
Je voulais tester la vitesse sur ma 603P mais "Lowbat!"

On remarque au passage que les 502/602/603P sont très proches de la TI57 niveau programmation, bien plus que les TI58/59
Si on ne gère pas le cas "0", le programme passe à 25 pas (en comptant P0) comme la TI57 :

Code : Tout sélectionner

Min00 1 Min01 Min02
LBL0
 x MR00 DSZ x=0 GOTO9
 / MR01 x MR00 / 1 M+01 MR01 = M+02 DSZ GOTO0
LBL9
 MR02

Soit quasiment le programme de C.RET (qui ne gére pas le cas 0 il me semble?) a 2 détails près :
-> inv DSZ n'existe pas sur 502P&co mais ici DSZ x=0 fait l'affaire
-> La 502P&co est plus souple coté syntaxe. Une séquence genre 2 * 1 M+01 MR01 = Min02 ne la géne pas. ca équivaut en pseudo langage à :
MR01 := MR01 + 1
MR02 := 2 * MR01

[Gilles59]

Avec l'algorithme de Gégé ma HP50g en NewRPL trouve Motzkin(260) :

- en 0,058 sec pour 32 chiffres significatifs. Le même programme en RPL tourne en 6,99s soit 120 fois plus lentement avec 12 chiffres significatifs! (0 sec. sur l'émulateur, qui est plutôt un simulateur en fait ... ou une version du NewRPL compilé avec une cible Windows. AU passage le NewRPL tourne sous Android aussi)
- en 0,092 sec pour 121 chiffres significatifs (0,001 sec sur l'émulateur). Dans de cas les 119 premiers chiffres du résultat sont exacts ;D (en comparant avec C.ret page 2). Ca finit par ....140 au lieu de ...139.

Le NewRpl sur la 50g est plus rapide (*) qu'une Prime et bien plus précise en calcul numérique ;D Par contre la Prime (et la HP50g standard) possèdent un mode exact que n'a pas vraiment le NewRPL.

Je me demande si une autre calculatrice gére aussi facilement une précision paramétrable que la 50g/NewRPL?
(*) en fait le NewRpl sur 50g est _parfois_ plus rapide que la Prime. Sur des calculs « pur » avec de la trigo ou autre elle est environ 2 fois plus lente.

Code : Tout sélectionner

«
  -> n
  «
    1 'a' LSTO 1 'k' LSTO 1 'm' LSTO
    n 2 + 'p' LSTO
    n 3 + 'q' LSTO
    WHILE k 2 * n <= REPEAT
     a p k / 2 * 'a' LSTO
     1 'k' STO+
     a q k / 2 - * DUP ' a' STO
     'm' STO+
    END
    m
  »

260 octets

J'ai signalé à l'auteur du NewRPL que

Code : Tout sélectionner

1 {'a' 'k' 'm'} LSTO 

devrait fonctionner plutôt que le moche

Code : Tout sélectionner

1 'a' LSTO 1 'k' LSTO 1 'm' LSTO

LSTO est un nouvelle commande à comprendre comme Local STO. Cela créé une variable locale visible uniquement dans le programme... ca manquait beaucoup en RPL. Une fois la variable locale créé LSTO et STO sont identiques.

Hors les graphiques et sons, le NewRPL est désormais opérationnel et très stable. La gestion de la carte SD comme un disque est parfaite. la connexion USB avec le PC épatante! La gestion des unités (et conversion d'unités en tout genre) est la meilleure que j'ai jamais vue... Reste quelques petites anomalies ici ou là. La seule chose que j'ose critiquer est la nouvelle interface, moi qui avait tellement l'habitude de l'ancienne, je peste encore contre le nouveau fonctionnement de la touche ALPHA par exemple mais vais finir par m'y faire... je suis totalement bluffé de ce qu'a fait Claudio Lapilli. C'est juste extraordinaire pour qui n'est pas allergique au RPL. J'en ai rêvé, il l'a fait ;D

Ah j'oubiais ! En plus le NewRpl laisse 415Ko de RAM dispo pour les programmes contre 240 environ pour le RPL ;D. Bon presque 20 éditions de ce message j'arrête là !

Pour en savoir plus c'est là https://newrpl.wiki.hpgcc3.org/doku.php?id=start

[Gilles59]

* comment on supprime un message? *