Multiplier deux nombres à la vitesse de l’éclair sur Z80

Même s’ils sont moins suivis que les compos, la Revision propose également des seminars. Ce sont des cours ou des exposés sur un thème précis. C’est souvent technique, on y apprend toujours quelque chose. Ce week-end, un séminaire a attiré mon attention : « Writing/building a 3D geometry engine from scratch », c’est à dire « Ecrire/construire un moteur 3D à partir de zéro ». À partir de zéro ? Il y a sûrement quelques tips à retirer de là. Surtout que le moteur en question tourne sur Vic 20 ! Si le cours en-lui même est trop léger à mon goût (l’auteur survole les formules mathématiques comme le code en lui -même. C’est compréhensible étant donné le sujet abordé et le peu de temps disponible), il reste quelques perles à grapiller. Notamment une phrase concernant une optimisation de multiplication. Il parle d’une méthode basée sur la différence de deux carrés.

Wikipedia est prolixe sur le sujet : Article en français sur la différence de deux carrés, mais seule la version anglaise de l’article contient une formule qui correspond à nos besoins :
xy=((x+y)/2)²-((x-y)/2)².
On va la simplifier pour notre CPC de la manière suivante :
xy=((x+y)²/2²)-((x-y)²/2²)
xy=((x+y)²/4)-((x-y)²/4)

En pratique, la seule difficulté est l’élévation au carré, mais ça peut se faire facilement avec une table en mémoire. Le résultat de la multiplication pourra être un nombre 16 bits : Il suffit que x et y soient supérieurs à 16. 16×16=256… Donc si on veut une multiplication utilisable il est impératif d’avoir un résultat sur deux octets. Mais un problème se pose : La division par 4 prend beaucoup trop de temps sur Z80. Huit nops, c’est énorme. Il faudrait vraiment pouvoir s’en passer.

// Division de HL par quatre
    srl h ; 2 NOPs 0 -> [H] -> CF
    rr l  ; 2 NOPs CF -> [L] -> CF
    srl h ; 2 NOPs 0 -> [H] -> CF
    rr l  ; 2 NOPs CF -> [L] -> CF

Et c’est là qu’intervient Overflow qui lisait par-dessus mon épaule : « J’ai appliqué une formule similaire dans Logon’s Run. Jette un oeil à la routine mul.8sx8s dans le fichier sys.Z80 » dans les sources de la démo. Il n’est pas utile de diviser par 4 à la fin : Il suffit de stocker directement les carrés divisés dans ta table. L’arrondi ne pose pas de problème. De mémoire il n’y a qu’un cas problématique et il se compense de lui-même ».

En conclusion, pour notre Z80 il suffira de constituer une table à 256 entrées de valeurs issues de la formule f(a)=a²/4 et le CPU n’aura qu’à effectuer le processus suivant :

• Je calcule x+y, je regarde la valeur correspondante dans la table et je la stocke quelque part
• Je calcule x-y, je regarde la valeur correspondante dans la table et je la stocke quelque part
• Je soustrais les deux valeurs

Simple. Mais bien entendu le diable se cache dans les détails. Depuis le début de l’article je parle de multiplication de nombre positifs, pas de nombres signés. La routine d’Overflow gère tous les cas et dure plus longtemps que les 23 NOPs qu’on ambitionne d’atteindre. Je propose de rester uniquement sur des nombres positifs, ça nous évitera des complications. Cependant, que se passe t’il si y>x ? dans ma fomule j’ai un « x-y », et le résultat de l’opération est un nombre négatif ! Alors qu’on pensait se débarrasser du problème avec une petite remarque au détour de l’article, le voilà qui réapparait !

On va devoir ajouter des restrictions complémentaires (Je vous rassure, il est tout à fait possible de s’en affranchir au prix d’une pirouette supplémentaire, mais ça coûte en NOPs comme en mémoire). Vous savez probablement comment sont stockés les nombres négatifs signés en Z80 : -1=256-1=255, -2=256-2=254, etc. Généralement on considère que si le bit 7 d’un nombre signé est égal à 1 c’est qu’il s’agit d’un nombre négatif. D’où la plage [-128;127] qu’on utilise habituellement. Ce qui veut dire qu’il faut trouver les valeurs maximales de x et de y pour qu’à la fois (x+y) et (x-y) se retrouvent à coup sûr dans la plage dont on dispose. On arrive à x et y inférieur ou égal à 63. Ainsi x+y=126 si x=y=63, ce qui est inférieur à 127, et x-y=-63 si x=0 et y=63, on est dans les clous. On n’utilise même pas 256 valeurs dans notre table, car on va juste avoir besoin de calculer des valeurs pour -63 à 126. Donc, dans notre table, 0 à 126 et 193 à 255. Avec un gros trou entre 127 et 192. C’est moche. On doit pouvoir faire mieux.

Considérons que a est la valeur maximale pour x comme pour y. Pour utiliser tous les octets de la table, il faudrait que (x+y) soit juste inférieur à la valeur de (0-y) en notation signée Z80. Cela revient à dire que 2a doit être inférieure à 256-a, donc que 3*a<256
256/3=85.33
85*2=170
-85 en notation signée c’est 256-85=171
Conclusion : En utilisant l’intégralité des 256 valeurs possibles pour la table, on pourra multiplier deux nombres positifs de 0 à 85 sans problème.

Sans plus attendre, voici la macro en LUA pour SJASM+ qui génèrera la table pour nous :

    macro MATH_PUT_MULTAB2
        LUA ALLPASS
            _pl("    align 256")
            _pl("MULTAB2:")
            print(string.format("Début MULTAB2 : %04x",_c("$")))

            for x = 0,170,1
            do
                sj.add_byte(math.floor((x*x/4)) % 256)
            end
            for x = -85,-1,1
            do
                sj.add_byte(math.floor((x*x/4)) % 256)
            end
            for x = 0,170,1
            do
                sj.add_byte(math.floor((x*x/4)) >> 8)
            end
            for x = -85,-1,1
            do
                sj.add_byte(math.floor((x*x/4)) >> 8)
            end

            print(string.format("Fin MULTAB2 : %04x",_c("$-1")))
        ENDLUA
    endm

Une petite astuce que j’ai oubliée de vous expliquer : Quand vous avez une table de valeurs 16 bits, il est préférable de stocker d’abord tous les poids faibles dans les 256 premiers octets, puis tous les poids forts dans les 256 suivants. C’est contre intuitif au premier abord : On pourrait stocker les valeurs 16 bits à la suite, c’est a priori moins compliqué à faire. Mais il faut toujours penser à notre Z80, bien embêté qu’il est lorsqu’il s’agit de manipuler des valeurs 16 bits !

Voici ce que pourrait donner une routine de lecture sur une table de données 16 bits stockée linéairement :

    ld hl,matable
    ld a,mon_offset_dans_ma_table
// On va devoir faire hl=hl+a*2
    add l    ; 1 NOP 
    ld l,a   ; 1 NOP
    ld a,0   ; 2 NOPs
    add h    ; 1 NOP
    ld h,a   ; 1 NOP
// Là on a hl=hl+a
    add l    ; 1 NOP 
    ld l,a   ; 1 NOP
    ld a,0   ; 2 NOPs
    add h    ; 1 NOP
    ld h,a   ; 1 NOP
// On a enfin hl=hl+a*2. 12 NOPs pour le calcul d'adresse !
    ld e,(hl); 2 NOPs
    inc hl   ; 2 NOPs
    ld d,(hl); 2 NOPs
// 18 NOPS pour lire la valeur dans la table

Avec le stockage des 256 poids faibles suivis des 256 poids forts :

    ld h,matable>>8 ; 2 NOPs - Il faut que la table soit alignée sur un multiple de 256
    ld l,a          ; 1 NOP - HL contient désormais l'adresse du poids faible
    ld e,(hl)       ; 2 NOPs
    inc h           ; 1 NOP - HL contient l'adresse du poids fort
    ld d,(hl)       ; 2 NOPs
// 8 NOPs pour lire la valeur dans la table !
    

Je pense qu’on a fait le plus compliqué. Le reste c’est un peu d’assembleur en utilisant toutes les astuces qu’on a pu voir au cours de cet article. C’est parti !

;=========================
; HL=DxE si 0<=D<=85 et 0<=E<=85
;-------------------------
; Durée : 23 NOPs
;-------------------------
; Sortie :
; HL = DxE
; A,DE,HL modifiés
;=========================
    macro MUL_D_E
      ld h,MULTAB2>>8   ; 2
      ld a,d            ; 1
      sub e             ; 1 => A = x-y
      ld l,a            ; 1 => HL pointe vers le poids faible de (x-y)²/4
      add e             ; 1
      add e             ; 1 => A = x-y +y+y = x+y
      ld e,(hl)         ; 2
      inc h             ; 1
      ld d,(hl)         ; 2 => DE = (x+y)²/4
      ld l,a            ; 1 => HL pointe vers le poids FORT de (x+y)²/4
      ld a,(hl)         ; 2
      dec h             ; 1
      ld l,(hl)         ; 2
      ld h,a            ; 1 => HL = (x+y)²/4
      sbc hl,de         ; 4 Attention carry => HL=x*y
    endm

On pourrait gagner un NOP supplémentaire si on changeait l’entête de la fonction en « HL=AxE », ce qui nous permettrait de se débarrasser de l’opcode « ld a,d ». C’est du détail. Dans tous les cas cette routine sera à adapter à vos besoins spécifiques. À noter, elle préserve également BC. En remplaçant les registres utilisés pour pourriez en faire une version qui utilise BC mais préserve DE.

Pour finir, voici le source complet de la multiplication Z80 en 23 NOPs, ainsi qu’un DSK contenant une version exécutable ainsi qu’un petit programme en Basic qui teste la multiplication avec toutes les valeurs possibles de x et de y afin de vérifier que tout fonctionne comme il faut.

Merci à Overflow et Longshot pour les tips au cours de la rédaction de cet article !