Misez Rapide–Accélérez n°2 : La suite des nombres de Hamming

[C.Ret]

Tout cela est très interessant.

Comme je ne comprenais pas bien le principe de la méthode géométrique, je me suis mis à explorer cette voie à ma manière.

J'obtiens un algorithme fort proche de celui proposé par gégé, sans toute fois passer par les coordonnées logarithmes.
Certes les logarithmes interviennent, mais pas dans le calcul des nombres, uniquement dans la détermination des limites de chaque "tétraèdres".

Ce qui fait, que je peux déterminer les nombres de Hamming au-delà de la précision de ma calculette, car le code n'utilise que les coordonnées géométriques (x,y,z) dans le repère (Ox>,(Oy>, (Oz> correspondant respectivement aux puissances de 2,3 et 5.

Comme pour la méthode volumétrique, je cherche à dénombrer les nombres de Hamming par intervalles successifs ] (2^k)/2 , 2^k ]

Pour k=0 le seul point possible est (0,0,0) c'est à dire H(1)=1
Pour k=1 le seul point possible est (1,0,0) c'est à dire H(2)=2
Pour k=2 deux points sont dénombrés (2,0,0) et (0,1,0) c'est à dire 4 et 3.

En fait, les points possibles pour le tétraèdre de niveau général k sont les points de coordonnées positives (ou nulles) vérifiant :
y*log(3) + z*log(5) <= x*log(2)

Ici, la fonction log peut être naturelle, décimale ou d'une base quelconque. Cela revient à 'quantifier' la vitesse à laquelle les sous-branches des n-multiples se développent.

Ce qui est intéressant, c'est que si un point (x,y,z) appartient au tétraèdre de niveau (k-1) alors le point (x+1,y,z) appartient au tétraèdre de niveau k, car les points se projettent orthogonalement sur le plan (Ox>

Par exemple, 3 et 4 appartiennent au tétraèdre k=2. On a donc les points 6 et 8 qui appartiennent au tétraèdre k=3
En effet, les nombres de Hamming de l'intervalle ]4,8] sont dans l'ordre 5 6 8

Dans la copie d'écran suivante, les nombres de Hamming déduits du tétraèdre précédant par multiplication sont indiqués entre parenthèses:

Code : Tout sélectionner

run
? 20
TETRA:                      tot MEM:
2^ 0 :  1                   t 1   1:
2^ 1 : (2)                  t 2
2^ 2 : (4)                         
        3                   t 4   2:
2^ 3 : (8)   5                    3:
       (6)                  t 7
2^ 4 :(16) (10)                   3:
      (12)  15                    4:
        9                   t 12  5:
2^ 5 :(32) (20)  25               6:
      (24) (30)                        
      (18)                             
       27                   t 19  7:
2^ 6 :(64) (40) (50)
      (48) (60)                   
      (36)  45                    8:
      (54)                  t 27

LIST MEM 8:
 64  60  54  50  48  45  40  36
                             ^
t27 t26 t25 t24 t23 t22 t21 t20

RESULTAT :
       2  2  0
H   = 2 .3 .5  = 36                 0'  1.4s
 20

ready.

Dans la première partie de l'algorithme, comme dans celui présenté par gégé, on dénombre les Nombres de Hamming présents dans chaque intervalle ](2^k)/2,2^k]. On en retient la position géographique (x,y,z) en mémorisant leurs coordonnées.

Pour cela je parcours les y et z possibles qui vérifient y*log(3)+z*log(5)<=x*log(2) - J'ai mis des log dans chaque terme, comme cela on utilise indifféremment log népérien, décimaux ou particuliers.

Le compteur n% est incrémenté pour chaque 'nouveau' nombre, c'est à dire une combinaison (y,z) qui n'existait pas dans le tétraèdre précédant; n% est donc utilisé comme indice mémoire étant donné qu'une fois qu'un couple de coordonnées (x,y,z) entre dans le tétraèdre 2^k, les points (1+x,y,z), (2+x,y,z), etc... seront également dans les tétraèdres suivants k+1, k+2, etc...

D'ailleurs, le code n'a pas de variable k%, c'est directement x% qui est utilisée.

Le total t% donne le nombre de Nombres de Hamming depuis le premier tétraèdre. On parcourt donc séquentiellement les tétraères jusqu'à obtenir un total t% au moins égal à l'indice recherché.

A l'issue de cette première partie, on obtient donc la liste des n% nombres de Hamming de l'intervalle ] (2^k)/2,2^k ]
Pour en déduire le r-ième, il est nécessaire de les trier.
C'est à dire de chainer la liste dans l'ordre décroissant. On sait que le premier nombre est le plus grand car il correspond à 2^k. Donc , il suffit d'insérer dans la liste les suivants afin qu'il forment une liste ordonnée.
C'est à cela que sert la valeur du pointeur a%() qui indique l'indice mémoire du Nombre de Hamming immédiatement inférieur (ou contient zéro pour marquer la fin de la liste).
Pour gagner du temps, la valeur (même approchée) du nombre est mémorisé dans le tableau N() ce qui facilite les nombreux tests pour le tri.

La dernière partie se contente de dérouler cette liste chainée ordonnée afin d'atteindre le nombre recherché d'indice r.

Code : Tout sélectionner

run
? 10
       2  1  0
h   = 2 .3 .5  = 12                    0'  .7"
 10

       2  2  0
h   = 2 .3 .5  = 36                    0' 1.4"
 20

        9  1  0
h    = 2 .3 .5  = 1536                 0' 5.1"
 100

         14  0  5
h     = 2  .3 .5  = 51200000           0'50.9"
 1000

          5  10  16
h      = 2 .3  .5   = 2.88325196e+17  12'59.5"
 10000

ready.

Et le listing :

Code : Tout sélectionner

10 rem ---------------------------- initiate
20 l2=log(2):l3=log(3):l5=log(5)
30 m%=52*36/2:dim n(m%),n%(52,36),x%(m%),y%(m%),z%(m%),a%(m%)
40 x%=-1:n%=0:t%=0:input r% : t0=timer
100 rem --------------------------- scan tetra and build list
110 do
120 :  x%=x%+1
130 :  for y=0 to x%*l2/l3
140 :  :  for z=0 to (x%*l2-y*l3)/l5
150 :  :  :  if n%(y,z)=0 then n%=n%+1:n%(y,z)=n%:x%(n%)=x%:y%(n%)=y:z%(n%)=z
160 :  :  next z
170 :  next y
180 :  t%=t%+n%
190 loop while t%<r%
200 rem --------------------------- sort list of h() in ] 2^x-1 , 2^x ]
210 n(1)=2^x%
220 for i=2 to n%:n(i)=int(2^(x%-x%(i)))*int(3^y%(i))*int(5^z%(i)):a%=1
230 :  do while n(i)<n(a%(a%)):a%=a%(a%): loop :a%(i)=a%(a%):a%(a%)=i
240 next i
300 rem --------------------------- display result
310 a%=1:do while r%<t%:t%=t%-1:a%=a%(a%): loop
320 i=a%:print "[dwn]h[dwn][lft]";t%;"[up]= 2[up][lft]";x%-x%(i)"[dwn][lft]
    .3[up][lft]"y%(i)"[dwn][lft].5[up][lft]"z%(i)"[dwn]="n(i);
400 t1=int((timer-t0)/6):t1%=t1/600:t1=t1-600*t1%:print tab(35)t1%".'"t1/10".s[dwn]"

ready.
N()    Nombre de hamming pour le tri (ligne 200-240) puis affichage (ligne 320)
a%() Pointeur de liste chainée désignant le nombre immédiatement inférieur (ou zéro si dernier de la liste)

x%() y%() z%()  coorodnnées du nombre de hamming dans le repère Oxyz correspondant aux expansion x2 x3 et x5

n%( y , z ) marqueur pour découverte des nouveaux points - ce tableau contient l'indice du nombre. Mais celui-ci n'est pas utilisé (sauf que les points non encore découverts ont un indice nul). On pourrait remplacer ce tableau par les fonctions graphiques PSET RPOINT d'un Pocket en allument le pixel (y,z) et en en lisant la valeur pour le test de la ligne 150.
  La dimenson du tableau permet en théorie d'atteindre 3^52 ou 5^36, c'est à dire des nombres de Hamming au-delà de 1.E+24

[C.Ret]

La nuit porte conseil, et de toute façon j'ai passé l'âge de faire nuit blanche, et j'ai repris le code comme suit en modifiant le calcul du modulo pour le calcul de A mais pas pour la position du pointeur.
Du coup en affichage, je peux utiliser la valeur du pointeur la plus faible (toujours C1) et donc afficher la taille du buffer utile pour atteindre le calcul en cours et en plus je gagne quelques dixièmes de secondes (ce qui m'étonne un peu d'ailleurs) :

Code : Tout sélectionner

10 CLEAR:DIM H(500):N=1:A=2:B=3:C=5:INPUT R:TIMER=0
20 H(O MOD 500)=N:O=O+1:IF A=N THEN A1=A1+1:A=2*H(A1 MOD500)
30 IF B=N THEN B1=B1+1:B=3*H(B1 MOD500)
40 IF C=N THEN C1=C1+1:C=5*H(C1 MOD500)
50 IF A<B AND A<C THEN N=A ELSE IF B<C THEN N=B ELSE N=C
60 IF O<R THEN 20
70 PRINT H((O-1)MOD500);TIMER;O-C1+1

Par contre, avec les limites imposées par la taille du buffer, il peut être intéressant de voir jusqu'où on peut aller avec les différents ordipoches en fonction de la taille mémoire disponible et de la précision interne de calcul en entiers.

Cette histoire du pointeur ^c5 produisant les multiples de 5 qui reste toujours à la traine permet dans le cas du tableau circulaire (les indices sont modulo pour tourner en rond) d'arrêter les calculs lors qu'il est rattrapé par le pointeur ^a2 produisant les multiples de 2 plus nombreux.

Le pointeur ^a2 fini toujours par rattraper ^c5; si l'on veut aller vers de plus grands Nombres de Hamming, il faut pouvoir augmenter la taille du tableau.
Une autre solution est de ne pas tourner en rond, mais de récupérer la mémoire en ne conservant pas les valeurs qui sont derrière le pointeur ^c5.


J'ai donc eut l'idée de modifier mon algorithme et d'utiliser l'instruction SUB pour décaler la liste des valeurs de Hamming à chaque fois que le pointeur ^c1 qui n'existe plus, la première valeur de la liste est systématiquement celle qui sert à calculer le prochain multiple de 5. Il faut juste penser à faire reculer ^a2 et ^b3. On économise ainsi pas mal de mémoire et on peut aller un peu plus loin avant d'arriver à la limite mémoire.

Code : Tout sélectionner

EXPORT Ham1(r)
BEGIN
//    a b et c  pointeurs pour les multiples de 2 3 et 5 (mémorisés respectivement dans ha hb et hc 
//       m        minimum		n         compteur de boucle	L1        liste globale  'BUFFER'
LOCAL a:=1,ha:=2,b:=1,hb:=3,hc:=5,m,n;
{1}▶L1;
FOR n FROM 2 TO t DO
   MIN(ha,hb,hc)▶m; m▶L1(0);
   IF m==ha THEN                                   1+a▶a;       2*L1(a)▶ha END;
   IF m==hb THEN                                   1+b▶b;       3*L1(b)▶hb END;
   IF m==hc THEN SUB(L1,2,n+1)▶L1; a-1▶a; b-1▶b; 5*L1(1)▶hc END;
END;
RETURN L1(0);
END;

J'ajoute mes résultats à ceux de caloubugs :

Code : Tout sélectionner

CASIO Z1-GR                          Hewlett Packard HP Prime
Rang        Hamming  Buffer   Temps  Rang           Hamming  Buffer   Temps
  50            243      28     2"7    50               243      27    34ms
 100           1536      46     5"9   100              1536      45    60ms
 200          16200      74    12"0   200             16200      73   122ms
 300          82944     100    18"4   300             82994      99   179ms
 400         311040     121    24"9   400            311040     120   242ms
 500         937500     141    31"5   500            937500     140   309ms
 600        2460375     161    38"1   600           2460375     160   376ms
 700        5989240     179    44"7   700           5989240     178   440ms
 800       12754584     197    51"4   800          12754584     196   515ms
 900       26244000     213    58"1   900          26244000     212   591ms
1000       51200000     230  1'04"8  1000          51200000     229   654ms
1500      859963392     304  1'38"8  1500         859963392     303   1"037
2000     8062156800     369  2'13"0  2000        8062156800     368   1"435
2500  5,36870912e10     431  2'47"5  2500       53687091200     430   1"865
                                    10000  2.88325195312e17    3420  19"073
                                    15000  1.23695058125e20    6419  48"853

Comme on peut le voir, la taille du buffer est la même et la stratégie ne fait rien gagner sauf peut-être un peu moins de calculs (il n'y a pas les tests ou les MOD des indices à réaliser. Par contre, je ne connais pas l'efficacité du 'ramasse miettes' de l'instruction SUB et s'il gère convenablement le fractionnement de la mémoire ?!

De toute façon, la quantité de mémoire augmente trop rapidement; entre 2500 à 10000 la taille de l'indice est multiplié par quatre et il faut presque 8x plus de mémoire. On ne pourra pas atteindre 20000 avant de dépasser la limite des 10'000 éléments dans le buffer (les listes sur HP Prime ne peuvent dépasser cette limite).


Je vais essayer l'autre méthode, celle qui découpe l'ensemble des Nombres de Hamming en fines tranches par puissances croissantes de 2^x.

Cette méthode est organisée en trois parties :

1. - énumération des Nombres de Hamming par tranche [ 2^x 2.2^x [ jusqu'à parcourir entièrement la tranche contenant le Nombre Recherché.
   Au minium, à chaque tranche on double la quantité de nombres répertoriés sans accroitre significativement l'espace mémoire nécessaire car les doubles des nombre de la tranche précédente se retrouvent systématiquement dans celle-ci. On a juste à ajouter des multiples de 3 ou de 5 qui apparaissent plus ou moins sporadiquement. par ailleurs, je viens de découvrir que ces nouveaux nombres n'apparaissent qu' à la périphérie de la tranche, jamais dans le tronc central. J'ai donc pu considérablement optimiser l'algorithme qui n'a plus à balayer chaque tranche en y et z, mais uniquement le long de l'axe des y. les z possibles se déduisant maintenant par le calcul.
2. -tri des nombres énumérés: l'avantage est que l'ordre des nombres reste le même d'une tranche à l'autre; On peut donc, à l'aide de structures de données appropriées effectuer cette opération au fur et à mesure de l'énumération. Ou ne le faire que sur une sous-partie ce qui peut faire gagner (ou perdre ) du temps. Dans l'algorithme que j'ai donné, tous les nombres sont triés après l'énumération par un tri peu efficace à bulles. On peut aussi comme le fait gege utiliser une méthode plus efficace basée sur une liste doublement chainée afin de profiter des propriétés des nombres triés.
   Sur la prime et afin d'économiser la mémoire, j'utiliserai la très efficace instruction de tri des listes SORT qui permet maintenant de donner une clef de tri. Cette nouvelle fonctionnalité va me permettre de conserver pour chaque nombre les puissances de 2 3 et 5 qui le définissent avec précision (pour avoir toutes les décimales) et de trier à l'aide d'une clef calculée à partir des logarithmes qui donneront, je l'espère, une valeur approchée suffisante pour un tri rapide, fiable et efficace.
3. -extraction du nombre recherché dans la liste triée: En tenant à jour le nombre d'éléments en mémoire et à partir de la différence entre le total de valeurs énumérées (et mémorisée) et la cible (l'argument de la fonction), il est facile de pointer sur le nombre demandé. Le nombre sera définit par le vecteur [ a b c ] formé des puissances respectives de 2 3 et 5 qui le définisse parfaitement. Ce vecteur, permettra, en mode CAS, d'afficher tous les chiffres en utilisant un entier long.

Code : Tout sélectionner

EXPORT Ham2(r)
BEGIN
  LOCAL h:={},l2:=LN(2),l3:=LN(3),l5:=LN(5),lb:=l2/l5,lc:=l3/l5;
  LOCAL n:=0,k:=0,x:=−1,y,z,e;    e:=MAKELIST(−1,y,1,300);
  REPEAT                                               // --------------- Enumération  tranche  [ 2^x ... 2.2^x [
    1+x▶x;
    FOR y FROM 0 TO x*la DO                  //  ---------------  Scan selon l'axe y (multiples de 3 ) 
      FLOOR(x*lb-y*lc)▶z; 
      IF e(1+y)<z THEN z▶e(1+y); 1+n▶n; {[x,y,z],y*l3+z*l5-x*l2}▶h(n) END;
    END;
    k+n▶k;                                            //  --------------- Compteur totalisateur
  UNTIL k≥r;
  SORT(h,2)▶h;                                   // ------------------ Tri (clef de tri valeur log en position 2 )
  RETURN [x-h(n+r-k,1,1),h(n+r-k,1,2),h(n+r-k,1,3)];
END;

Et le code pour le mode CAS:

Code : Tout sélectionner

#cas
HamCas(k):=
BEGIN 
  LOCAL t,h;
  t:=TICKS;
  h:=Ham2(k);
  { 2^h(1)*3^h(2)*5^h(3), ETIME(t)  }
END;
#end

Bon, autant prévenir tout de suite, ça dépote !

MRA2 - HP Prime Full speed all digits upto 900'000.png (15.66 Kio) Vu 2117 fois

Seul regret, je ne peux pas déterminer H(1000000). Mais j'obtiens toutes les valeurs exactes (entiers long) des Nombres de Hamming jusqu'à H(900000) en moins d'une minute dix secondes.

En tout cas , j'ai bien fait de ne pas implémenter une structure de données complexe permettant de pré-trier les valeurs énumérées, car les 9'975 valeurs énumérées pendant les 1.05min sont triées par SORT en moins de 756ms; Alors qu'avec l'algorithme implémenté sur mon vieux (mais fidèle) Commodore C128 passe 4 fois plus de temps à trier qu'à énumérer. Je suspecte fort l'HP Prime de trier (au moins en partie par 'chucks') les valeurs au fur et à mesure de la construction des listes, ce qui d'ailleurs expliquerait la limite du nombre d'élément qui serait en fait la limite des tables d'adressages rapides des 'chucks'.

Code : Tout sélectionner

Rang                                  Hamming           Buffer   Temps  Buffer   Temps1  NbEnum     Temps2
    50                                 243                  28     2"7      27    34ms       14      25ms
   100                                 1536                 46     5"9      45    60ms       23      39ms
   200                                 16200                74    12"0      73   122ms       34      52ms
   300                                 82944               100    18"4      99   179ms       48      68ms
   400                                311040               121    24"9     120   242ms       60      80ms
   500                                937500               141    31"5     140   309ms       65      84ms
   600                                2460375              161    38"1     160   376ms       78      93ms
   700                                5989240              179    44"7     178   440ms       85      97ms
   800                               12754584              197    51"4     196   515ms       91     102ms
   900                               26244000              213    58"1     212   591ms       98     111ms
  1000                               51200000              230  1'04"8     229   654ms      107     119ms
  1500                               859963392             304  1'38"8     303   1"037      138     156ms
  2000                              8062156800             369  2'13"0     368   1"435      166     179ms
  2500                              53687091200            431  2'47"5     430   1"865      195     200ms
 10000                           288325195312500000                       3420  19"073      504     553ms
 15000                          123695058124800000000                     6419  48"853      646     737ms
 20000                         15424418419015680000000                                      783     959ms
 50000                     2379126835648701693226200858624                                 1440     2"227
100000                 290142196707511001929482240000000000000                             2294     4"626
500000     1962938367679548095642112423564462631020433036610484123229980468750             6754    30"813
750000 9844116074857088479400896102109753641824661421606444941755765750304761446400        8831    50"432
900000 632306706993546185913146473467350006103515625000000000000000000000000000000000000   9975  1'03"598
900001 632373585725680579031685776908592159413202566649380479277264004972544000000000000   9975  1'03"193

Il manque quand un buzzer sur cette HP Prime, s'aurait été bien pratique qu'elle 'bip' quand elle a fit ou rencontre une difficulté !