Routines sur les Blocs
Cela rappellera peut-être des souvenirs de l'école primaire à certains, nous allons opter pour une représentation dite notation de position. En base 10, cela signifie que les chiffres sont multipliés par une puissance de 10 en fonction de leur position. Les unités sont multipliée par \(10^0 = 1\), les dizaines par \(10^1 = 10\), les centaines par \(10^2 = 100\), etc...
Pour bien comprendre les routines les plus évoluées, commençons par l'implémentation la plus naïve possible en C/C++ des opérations d'addition et de soustraction. Stockons les chiffres, en base 10, dans un tableau d'entiers et effectuons scrupuleusement les méthodes apprises à l'école:
int Addition(int* somme, int* param1, int* param2, int taille)
{
int retenue = 0;
for (int i = 0; i < taille; i++)
{
int s = param1[i] + param2[i] + retenue;
if (s >= 10)
{
s -= 10;
retenue = 1;
}
else
retenue = 0;
somme[i] = s;
}
return retenue;
};
int Soustraction(int* difference, int* param1, int* param2, int taille)
{
int retenue = 0;
for (int i = 0; i < taille; i++)
{
int d = param1[i] - param2[i] - retenue;
if (d < 0)
{
d += 10;
retenue = 1;
}
else
retenue = 0;
difference[i] = s;
}
return retenue;
};
Sauf qu'ici, nous travaillerons non pas en base 10 mais en base \(2^{64}\) et que, en mémoire, les "unités" précèderont les "dizaines" qui précèderont les "centaines", etc...
Ainsi, nous notons habituellement, en base 10, le nombre 'abcd' comme signifiant: \[ d + c * 10 + b * 100 + a * 1000 \] Nous noterons, en mémoire dans des "cases" numérotées à partir de \(0\), un nombre N: \[ N = Mem[0] + \beta * Mem[1] + \beta^2 * Mem[2] + \beta^3 * Mem[3] + ... \] où \(\beta = 2^{64}\) est la base.
Il est à noter qu'alors les "cases" mémoires, ici notées \(Mem[x]\), ne sont pas des octets mais des mots 64 bits. En effet, en base 10, les chiffres sont représentés par 10 symboles de \(0\) à \(9\) et donc, en base \(2^{64}\), les chiffres, dés lors notés mots, seront \(2^{64}\) symboles de \(0\) à \(2^{64}-1\) stockés comme des entiers non signés en 64 bits, donc sur 8 octets.
Nous voyons donc que l'utilisation mémoire est optimale au sens où la totalité des valeurs d'un mot mémoire (8 octets) sont significatives. De plus, cette notation est unique pourvu que l'on spécifie que le mot de poids fort soit non nul (sauf pour la représentation de 0, nous y reviendrons). Cette représentation, dite little endian car les mots de poids faible précèdent (en mémoire) les mots de poids plus élevé.
Un tel "groupe" de n mots sera, par la suite, appelé un bloc de taille n. Gardez bien en tête cette représentation si vous voulez comprendre tout le code qui suivra car ces blocs seront les briques de nos routines de plus bas niveau implémentant les opérations d'addition, soustraction, multiplication et division.
Concrètement
En C++, un bloc est donc un pointeur sur un entier non signé 64 bits et une taille (aussi un entier non signé). Donc en admettant les déclarations suivantes:
using Word = uint64_t;
struct Bloc
{
Word* data;
uint64_t size;
};
Un bloc représente l'entier naturel:
\[ N = \sum_{k=0}^{size-1} data[k] \; {\beta}^{k} \]
où \(\beta = 2^{64}\).
Cependant, nous n'aurons pas de type Bloc dans le code. Nous aurons à traiter des "morceaux" de Bloc: des parties de nombre entier ou de nombre à virgule flottante. C'est pourquoi nos routines de plus bas niveau prendront les pointeurs et tailles en paramètres, directement. Cette 'struct' précédente étant juste pour donner une représentation mentale... Voyons, enfin, les premières routines réellement utilisées.
Routines de base sur les Blocs
Voici les premières routines de bas niveau en assembleur. Conformément aux conventions en assembleur, j'ai utilisé le mot 'QWORD' (pour Quad Word: 4 mots de 2 octets soit 8 octets) mais ce n'est qu'en commentaire. Les "vraies" déclarations C++ seront un peu plus tard. Notez également qu'aucune vérification n'est faite. Dans un souci de performance, les vérifications sont faites le plus en amont possible. Ces routines traiterons des blocs de tailles identiques et leur préfixe sera donc la lettre "B".
Routines de copie
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- B(n)
;
; void BCopy(QWORD* a, QWORD* b, QWORD n)
align 16
BCopy PROC
lea rcx,[rcx+8*r8] ; rcx <-- adresse de a[n]
lea rdx,[rdx+8*r8] ; rdx <-- adresse de b[n]
neg r8 ; r8 <-- (-n)
align 16
bcopy_loop:
mov rax,[rdx+8*r8] ; rax <-- b[n+r8]
mov [rcx+8*r8],rax ; a[n+r8] <-- rax
inc r8 ; r8++
jnz bcopy_loop ; tant que r8 < 0
ret
BCopy ENDP
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- B(n)
;
; void BCopyR(QWORD* a, QWORD* b, long long n)
align 16
BCopyR PROC
mov rax,[rdx+8*r8-8] ; rax <-- b[r8-1]
mov [rcx+8*r8-8],rax ; a[r8-1] <-- rax
dec r8 ; r8--
jnz BCopyR ; tant que r8 > 0
ret
BCopyR ENDP
La routine BCopy copie un tableau dans l'ordre mémoire croissant. elle pourra donc être utilisée pour effectuer des décalages à droite in place (qui sont des décalages "à gauche" en mémoire) sans risque de corruption. Pour les décalages à gauche (donc "à droite" en mémoire), il faudra utiliser BCopyR qui copie dans l'ordre décroissant et qui ne sera utilisée que pour ces décalages "in place". Car pour copier un tableau dans un autre, la version croissante est plus performante (les processeurs actuels sont optimisés pour les accés incrémentiels à des cases mémoires adjacentes, et pas seulement pour des raisons de cache...).
Les connaisseurs s'étonneront peut-être que je n'utilise pas d'instruction du type "REP MOVSQ". Il s'avère qu'il est plus performant de faire une boucle "conventionnelle" pour des raisons liées à la structure superscalaire des processeurs actuels. D'autre part, la boucle telle qu'elle est faite ici (avec un compteur négatif) permet de conserver le flag CARRY (mais ce n'est pas utile ici) et d'effectuer le test de boucle sur le flag ZERO (plutôt que de faire un CMP suivit d'un jump conditionnel).
Routine de remplissage
Cette routine nous sera utile pour remplir par des zéros, mais aussi par des #FFFFFFFFFFFFFFFF (Base-1). Est-ce nécessaire de la commenter ?
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- x
;
; void BFill(QWORD* a, long long n, QWORD x)
align 16
BFill PROC
lea rcx,[rcx+8*rdx]
neg rdx
align 16
bfill_loop:
mov [rcx+8*rdx],r8
inc rdx
jnz bfill_loop
ret
BFill ENDP
Pour les mêmes raisons que citées plus haut, ce code est plus rapide que REP STOSQ.
Routine de comparaison
Cette routine compare deux entiers de même taille.
;//////////////////////////////////////////////////////////////////////////////
; result <- A(n) comparison B(n)
;
; long BCmp(QWORD* a, QWORD* b, long long n)
align 16
BCmp PROC
sub rcx,8 ; rcx <-- adresse de a[-1]
sub rdx,8 ; rdx <-- adresse de b[-1]
align 16
bcmp_loop:
mov rax,[rcx+8*r8] ; rax <-- a[r8-1]
cmp rax,[rdx+8*r8] ; compare rax et b[r8-1] (carry mis si b[r8-1] > rax)
jnz bcmp_break ; si différents goto break
dec r8 ; r8--
jnz bcmp_loop ; tant que r8 > 0
mov rax,0 ; retourne 0
ret
bcmp_break:
sbb rax,rax ; rax <-- (-1 si carry) (0 sinon) donc (-1) si (b > a)
shl rax,1 ; rax <-- (-2) ou 0
inc rax ; retourne (-1) si (a < b) (+1) si (a > b)
ret
BCmp ENDP
Routines d'addition sur les Blocs
L'algorithme d'addition est directement issu de l'addition que nous avons tous appris à l'école, le passage du décimal à une base quelconque \(B \geq 2\) étant immédiat.
L'opération élémentaire de l'addition en base quelconque \(B \geq 2\) peut être écrite comme suit, pour n allant des mots de poids faible vers les mots de poids fort:
somme[n] = a[n] + b[n] + retenue;
if (somme[n] >= B)
{
somme[n] -= B;
retenue = 1;
}
else
retenue = 0;
Sous l'hypothèses que: \[ 0 \leq a[n] \leq (B-1) \] \[ 0 \leq b[n] \leq (B-1) \] \[ 0 \leq retenue \leq 1 \] on a donc: \[ 0 \leq a+b+retenue \leq (2B-1) \] d'où, si \( B \leq a+b+retenue \leq (2B-1) \) alors nous entrons dans le "si" et, en ôtant \(B\), on trouve: \[ 0 \leq somme[n] \leq (B-1) \] \[ retenue = 1 \] Sinon, c'est que \( 0 \leq a+b+retenue \leq (B-1) \) et donc: \[ 0 \leq somme[n] \leq (B-1) \] \[ retenue = 0 \] Ainsi, en répétant ce processus en allant des mots de poids faibles vers les mots de poids forts et en propageant à chaque fois la retenue précédemment obtenue, nous avons un algorithme valide de calcul de la représentation de la somme de deux représentations de taille arbitraire.
Cependant, nous avons un terme intermédiaire \(0 \leq somme[n] \leq (2B-1)\) qui nécessite un bit de plus que les mots en base \(B\). Si nous choisissons \(2^{64}\) comme base, alors \(somme[n]\) ne pourra pas être représenté dans un registre 64 bits sans dépassement.
Fort heureusement, en assembleur il existe une instruction qui correspond exactement au pseudo-code précédent: l'instruction ADC (pour Addition avec Carry). En C/C++, il faudrait recourir à une base inférieure ou à une implémentation spécifique (et coûteuse).
Additions élémentaires
L'addition se fait du poids faible vers le poids fort avec propagation d'une éventuelle retenue à travers le flag CARRY. On appréciera le fait que le flag CARRY soit conservé avec ce type de boucle.
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- B(n) + C(n)
;
; QWORD BAdd(QWORD* a, QWORD* b, QWORD* c, long long n)
align 16
BAdd PROC
lea rcx,[rcx+8*r9] ; rcx <-- adresse de a[n]
lea rdx,[rdx+8*r9] ; rdx <-- adresse de b[n]
lea r8,[r8+8*r9] ; r8 <-- adresse de c[n]
neg r9 ; r9 <-- (-n)
clc ; carry <-- 0
align 16
badd_loop:
mov rax,[rdx+8*r9] ; rax <-- b[n+r9]
adc rax,[r8+8*r9] ; rax <-- c[n+r9] + b[n+r9] + carry et carry <-- retenue
mov [rcx+8*r9],rax ; a[n+r9] <-- rax
inc r9 ; r9++
jnz badd_loop ; tant que r9 < 0
mov rax,0
rcl rax,1 ; retourne carry
ret
BAdd ENDP
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- A(n) + B(n)
;
; QWORD BAdd_ip(QWORD* a, QWORD* b, long long n)
align 16
BAdd_ip PROC
lea rcx,[rcx+8*r8]
lea rdx,[rdx+8*r8]
neg r8
clc
align 16
baddip_loop:
mov rax,[rcx+8*r8]
adc rax,[rdx+8*r8]
mov [rcx+8*r8],rax
inc r8
jnz baddip_loop
mov rax,0
rcl rax,1
ret
BAdd_ip ENDP
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- B(n) + C(n) + r (0 <= r <= 1)
;
; QWORD BAdc(QWORD* a, QWORD* b, QWORD* c, long long n, QWORD r)
align 16
BAdc PROC
lea rcx,[rcx+8*r9]
lea rdx,[rdx+8*r9]
lea r8,[r8+8*r9]
neg r9
mov rax,[rsp+40] ; rax <-- r
rcr rax,1 ; carry <-- bit 0 de rax
align 16
badc_loop:
mov rax,[rdx+8*r9]
adc rax,[r8+8*r9]
mov [rcx+8*r9],rax
inc r9
jnz badc_loop
mov rax,0
rcl rax,1
ret
BAdc ENDP
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- A(n) + B(n) + r (0 <= r <= 1)
;
; QWORD BAdc_ip(QWORD* a, QWORD* b, long long n, QWORD carry)
align 16
BAdc_ip PROC
lea rcx,[rcx+8*r8]
lea rdx,[rdx+8*r8]
neg r8
rcr r9,1
align 16
badcip_loop:
mov rax,[rcx+8*r8]
adc rax,[rdx+8*r8]
mov [rcx+8*r8],rax
inc r8
jnz badcip_loop
mov rax,0
rcl rax,1
ret
BAdc_ip ENDP
Routines d'incrémentation
Dans le cadre des routines d'addition, je vous propose les fonctions dites "d'incrémentation". Ce sont des fonctions qui ajoutent \(1\) à une représentation. Cela peut paraître trivial mais il n'en est rien. Ces routines peuvent s'obtenir comme des spécialisations de BAdd où le second paramètre ne contient que des mots nuls sauf le premier qui est égal à \(1\), ou bien de BAdc où le second paramètre est nul et la retenue est égale à \(1\).
Le second cas est peut-être le plus facile à comprendre. Dés lors que la retenue est nulle, ça ne devient qu'une simple copie. Et la retenue devient nulle dés lors que nous tombons sur un mot différent de \(B-1\).
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- B(n) + 1
;
; QWORD BInc(QWORD* a, QWORD* b, long long n)
align 16
BInc PROC
lea rcx,[rcx+8*r8] ; rcx <-- adresse de a[n]
lea rdx,[rdx+8*r8] ; rdx <-- adresse de b[n]
neg r8 ; r8 <-- (-n)
mov r9,-1 ; r9 <-- #FFFFFFFFFFFFFFFF
mov r10,0 ; r10 <-- 0
align 16
binc_loop:
mov rax,[rdx+8*r8] ; rax <-- b[n+r8]
cmp rax,r9 ; compare rax avec #FFFFFFFFFFFFFFFF
jnz binc_next ; si différent goto next
mov [rcx+8*r8],r10 ; a[n+r8] <-- 0
inc r8 ; r8++
jnz binc_loop ; tant que r8 < 0
mov rax,1 ; retourne 1
ret
binc_next:
inc rax ; rax <-- b[n+r8] + 1
mov [rcx+8*r8],rax ; a[n+r8] <-- rax
inc r8 ; r8++
jz binc_end ; si r8 == 0 goto end
align 16
binc_loop2: ; boucle de copie
mov rax,[rdx+8*r8]
mov [rcx+8*r8],rax
inc r8
jnz binc_loop2
binc_end:
mov rax,0 ; retourne 0
ret
BInc ENDP
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- A(n) + 1
;
; QWORD BInc_ip(QWORD* a, long long n)
align 16
BInc_ip PROC
lea rcx,[rcx+8*rdx]
neg rdx
mov r8,-1
mov r9,0
align 16
bincip_loop:
mov rax,[rcx+8*rdx]
cmp rax,r8
jnz bincip_next
mov [rcx+8*rdx],r9
inc rdx
jnz bincip_loop
mov rax,1
ret
bincip_next:
inc rax
mov [rcx+8*rdx],rax
mov rax,0
ret
BInc_ip ENDP
Addition de blocs de tailles différentes
Jusqu'ici, les routines présentées n'additionnaient que des blocs de tailles identiques et renvoyaient un bloc de taille égale à la taille des paramètres avec une éventuelle retenue. Ces routines peuvent servir de "briques" à un algorithme de plus haut niveau qui additionnerait deux blocs de tailles différentes et renverrait leur somme ainsi que sa taille.
Voici un exemple de ce que seraient de telles routines en C:
///////////////////////////////////////////////////////////////////////////////
// A <- B(nb) + C(nc)
//
// return the size of A
SIZE_T MAdd(LIMB_T* a, LIMB_T* b, SIZE_T nb, LIMB_T* c, SIZE_T nc)
{
LIMB_T r;
if( nb < nc )
{
r = BAdd(a, b, c, nb);
if( r )
r = BInc(a+nb, c+nb, nc-nb);
else
BCopy(a+nb, c+nb, nc-nb);
if( r )
{
a[nc] = 1;
return nc+1;
}
return nc;
}
/// nb >= nc
r = BAdd(a, b, c, nc);
if( nb > nc )
{
if( r )
r = BInc(a+nc, b+nc, nb-nc);
else
BCopy(a+nc, b+nc, nb-nc);
}
if( r )
{
a[nb] = 1;
return nb+1;
}
return nb;
}
///////////////////////////////////////////////////////////////////////////////
// A(result) <- A(na) + B(nb)
//
// return the new size of A
SIZE_T MAdd_ip(LIMB_T* a, SIZE_T na, LIMB_T* b, SIZE_T nb)
{
LIMB_T r;
if( na < nb )
{
r = BAdd_ip(a, b, na);
if( r )
r = BInc(a+na, b+na, nb-na);
else
BCopy(a+na, b+na, nb-na);
if( r )
{
a[nb] = 1;
return nb+1;
}
return nb;
}
/// na >= nb
r = BAdd_ip(a, b, nb);
if( (r != 0) && (na > nb) )
r = BInc_ip(a+nb, na-nb);
if( r )
{
a[na] = 1;
return na+1;
}
return na;
}
On voit bien que nos "briques" sont vraiment de bas niveau... En réalité, j'ai longtemps utilisé ces routines. Jusqu'à ce que je décide d'écrire directement MAdd et MAdd_ip directement en assembleur, gagnant ainsi beaucoup de temps (lié à l'appel de sous-fonctions et au passage de paramètres).
Ce sont ces fonctions qui seront appelées, la plupart du temps, pour additionner. En voici le code complet:
;//////////////////////////////////////////////////////////////////////////////
; A <- B(nb) + C(nc)
;
; return the size of A
;
; QWORD MAdd(QWORD* a, QWORD* b, QWORD nb, QWORD* c, QWORD nc)
;
align 16
MAdd PROC ; rcx := a rdx := b r8 := nb r9 := c [rsp+40] := nc
mov r10,[rsp+40] ; r10 <- nc
cmp r8,r10
jz madd_equal ; if (nb == nc) goto equal
jc madd_main ; if (nb < nc) goto main
;////////// here: nb > nc
xchg rdx,r9 ; b <-> c
xchg r8,r10 ; nb <-> nc
mov [rsp+40],r10 ; save nb
madd_main: ;////////// here: nb < nc
lea rcx,[rcx+8*r8] ; rcx <- addr a[nb]
lea rdx,[rdx+8*r8] ; rdx <- addr b[nb]
lea r9,[r9+8*r8] ; r9 <- addr c[nb]
sub r10,r8 ; r10 <- nc - nb
neg r8 ; r8 <- (-nb)
clc ; clear CF
align 16
madd_loop1:
mov rax,[rdx+8*r8] ; rax <- b[nb+r8]
adc rax,[r9+8*r8] ; rax <- b[nb+r8] + c[nb+r8] + CF
mov [rcx+8*r8],rax ; a[nb+r8] <- rax
inc r8
jnz madd_loop1 ; while r8 < 0
jnc madd_no_carry
;////////// here CF is set so: increment
lea rcx,[rcx+8*r10] ; rcx <- addr a[nc]
lea r9,[r9+8*r10] ; r9 <- addr c[nc]
neg r10 ; r10 <- -(nc-nb)
mov r11,-1 ; r11 <- (-1)
align 16
madd_loop2:
mov rax,[r9+8*r10] ; rax <- c[nc+r10]
cmp rax,r11
jnz madd_next ; if rax != (-1) goto next
mov [rcx+8*r10],r8 ; a[nc+r10] <- 0
inc r10
jnz madd_loop2 ; while r10 < 0
mov rax,1
mov [rcx],rax ; a[nc] <- 1
mov rax,[rsp+40] ; rax <- nc
inc rax ; return nc+1
ret
madd_next:
inc rax ; rax <- c[nc+r10] + 1
mov [rcx+8*r10],rax ; a[nc+r10] <- rax
inc r10
jz madd_end
align 16
madd_loop3:
mov rax,[r9+8*r10] ; rax <- c[nc+r10]
mov [rcx+8*r10],rax ; a[nc+r10] <- rax
inc r10
jnz madd_loop3 ; while r10 < 0
madd_end:
mov rax,[rsp+40] ; rax <- nc
ret ; return nc
madd_no_carry:
lea rcx,[rcx+8*r10] ; rcx <- addr a[nc]
lea r9,[r9+8*r10] ; r9 <- addr c[nc]
neg r10 ; r10 <- -(nc-nb)
align 16
madd_loop4:
mov rax,[r9+8*r10] ; rax <- c[nc+r10]
mov [rcx+8*r10],rax ; a[nc+r10] <- rax
inc r10
jnz madd_loop4 ; while r10 < 0
mov rax,[rsp+40] ; rax <- nc
ret ; return nc
;////////// nb == nc //////////////////////////////////////////
madd_equal:
lea rcx,[rcx+8*r8] ; rcx <- addr a[nb]
lea rdx,[rdx+8*r8] ; rdx <- addr b[nb]
lea r9,[r9+8*r8] ; r9 <- addr c[nb]
neg r8 ; r8 <- (-nb)
clc
align 16
madd_equ_loop:
mov rax,[rdx+8*r8] ; rax <- b[nb+r8]
adc rax,[r9+8*r8] ; rax <- b[nb+r8] + c[nb+r8] + CF
mov [rcx+8*r8],rax ; a[nb+r8] <- rax
inc r8
jnz madd_equ_loop ; while r8 < 0
jc madd_equ_carry
mov rax,r10 ; rax <- nc
ret ; return nc
madd_equ_carry:
mov rax,1
mov [rcx],rax ; a[nc] = 1
mov rax,r10 ; rax <- nc
inc rax
ret ; return nc+1
MAdd ENDP
;//////////////////////////////////////////////////////////////////////////////
; A <- B(nb) + C(nc)
;
; return the size of A
;
; QWORD MAdd_ip(QWORD* a, QWORD na, QWORD* b, QWORD nb)
;
align 16
MAdd_ip PROC ; rcx := a rdx := na r8 := b r9 := nb
cmp rdx,r9
jz maddip_equal ; if (na == nb) goto equal
jc maddip_m ; if (na < nb) goto m
;////////// here: na > nb
lea rcx,[rcx+8*r9] ; rcx <- addr a[nb]
lea r8,[r8+8*r9] ; r8 <- addr b[nb]
mov r11,rdx ; save na
sub rdx,r9 ; rdx <- na - nb
neg r9 ; r9 <- (-nb)
clc ; clear CF
align 16
maddip_loop11:
mov rax,[rcx+8*r9] ; rax <- a[nb+r9]
adc rax,[r8+8*r9] ; rax <- a[nb+r9] + b[nb+r9] + CF
mov [rcx+8*r9],rax ; a[nb+r9] <- rax
inc r9
jnz maddip_loop11 ; while r9 < 0
jc maddip_carry
mov rax,r11 ; rax <- na
ret ; return na
;////////// here CF is set so: increment
maddip_carry:
lea rcx,[rcx+8*rdx] ; rcx <- addr a[na]
neg rdx ; rdx <- -(na-nb)
mov r10,-1 ; r10 <- (-1)
align 16
maddip_loop22:
mov rax,[rcx+8*rdx] ; rax <- a[na+rdx]
cmp rax,r10
jnz maddip_next1 ; if rax != (-1) goto next1
mov [rcx+8*rdx],r9 ; a[na+rdx] <- 0
inc rdx
jnz maddip_loop22 ; while rdx < 0
mov rax,1
mov [rcx],rax ; a[na] <- 1
mov rax,r11 ; rax <- na
inc rax ; return na+1
ret
maddip_next1:
inc rax ; rax <- a[na+rdx] + 1
mov [rcx+8*rdx],rax ; a[na+rdx] <- rax
mov rax,r11 ; rax <- na
ret ; return na
;////////// here: na < nb
maddip_m:
lea rcx,[rcx+8*rdx] ; rcx <- addr a[na]
lea r8,[r8+8*rdx] ; r8 <- addr b[na]
mov r11,r9 ; save nb
sub r9,rdx ; r9 <- nb - na
neg rdx ; rdx <- (-na)
clc ; clear CF
align 16
maddip_loop1:
mov rax,[rcx+8*rdx] ; rax <- a[na+rdx]
adc rax,[r8+8*rdx] ; rax <- a[na+rdx] + b[na+rdx] + CF
mov [rcx+8*rdx],rax ; a[na+rdx] <- rax
inc rdx
jnz maddip_loop1 ; while rdx < 0
jnc maddip_no_carry
;////////// here CF is set so: increment
lea rcx,[rcx+8*r9] ; rcx <- addr a[nb]
lea r8,[r8+8*r9] ; r8 <- addr b[nb]
neg r9 ; r9 <- -(nb-na)
mov r10,-1 ; r10 <- (-1)
align 16
maddip_loop2:
mov rax,[r8+8*r9] ; rax <- b[nb+r9]
cmp rax,r10
jnz maddip_next ; if rax != (-1) goto next
mov [rcx+8*r9],rdx ; a[nb+r9] <- 0
inc r9
jnz maddip_loop2 ; while r9 < 0
mov rax,1
mov [rcx],rax ; a[nb] <- 1
mov rax,r11 ; rax <- nb
inc rax ; return nb+1
ret
maddip_next:
inc rax ; rax <- b[nb+r9] + 1
mov [rcx+8*r9],rax ; a[nb+r9] <- rax
inc r9
jz maddip_end
align 16
maddip_loop3:
mov rax,[r8+8*r9] ; rax <- b[nb+r9]
mov [rcx+8*r9],rax ; a[nb+r9] <- rax
inc r9
jnz maddip_loop3 ; while r9 < 0
maddip_end:
mov rax,r11 ; rax <- nb
ret ; return nb
maddip_no_carry:
lea rcx,[rcx+8*r9] ; rcx <- addr a[nb]
lea r8,[r8+8*r9] ; r8 <- addr b[nb]
neg r9 ; r10 <- -(nb-na)
align 16
maddip_loop4:
mov rax,[r8+8*r9] ; rax <- b[nb+r9]
mov [rcx+8*r9],rax ; a[nb+r9] <- rax
inc r9
jnz maddip_loop4 ; while r9 < 0
mov rax,r11 ; rax <- nb
ret ; return nb
;////////// na == nb //////////////////////////////////////////
maddip_equal:
lea rcx,[rcx+8*r9] ; rcx <- addr a[nb]
lea r8,[r8+8*r9] ; r8 <- addr b[nb]
neg r9 ; r9 <- (-nb)
clc
align 16
maddip_equ_loop:
mov rax,[rcx+8*r9] ; rax <- a[nb+r9]
adc rax,[r8+8*r9] ; rax <- a[nb+r9] + b[nb+r9] + CF
mov [rcx+8*r9],rax ; a[nb+r9] <- rax
inc r9
jnz maddip_equ_loop ; while r9 < 0
jc maddip_equ_carry
mov rax,rdx ; rax <- na
ret ; return na
maddip_equ_carry:
mov rax,1
mov [rcx],rax ; a[nc] = 1
mov rax,rdx ; rax <- na
inc rax
ret ; return na+1
MAdd_ip ENDP
On voit bien que, dés que la complexité augmente un peu, le code assembleur devient vite assez long... Mais je puis vous assurer que le gain est substantiel et en vaut largement la chandelle. Cela fait partie du prix à payer pour rivaliser avec les meilleures implémentations.
Routines de négation
Ces routines suivent le principe du complément à deux. Le principe est le suivant: en commençant par les mots de poids faible, tant qu'ils sont nuls, on les copie. Et à partir du premier qui n'est pas nul, on fait un 'neg' sur ce mot (générant une retenue) puis on inverse tous les bits des mots suivants jusqu'au bout (correspondant à un 'neg' puis à l'ajout de cette retenue).
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- -B(n)
;
; QWORD BNeg(QWORD* a, QWORD* b, long long n)
align 16
BNeg PROC
lea rcx,[rcx+8*r8] ; rcx <-- adresse de a[n]
lea rdx,[rdx+8*r8] ; rdx <-- adresse de b[n]
neg r8 ; r8 <-- (-n)
align 16
bneg_loop:
mov rax,[rdx+8*r8] ; rax <-- b[n+r8]
or rax,rax
jnz bneg_next ; si rax != 0 goto next
mov [rcx+8*r8],rax ; a[n+r8] <-- rax (qui est nul)
inc r8 ; r8++
jnz bneg_loop ; tant que r8 < 0
mov rax,0 ; retourne 0
ret
bneg_next:
neg rax ; rax <-- -b[n+r8]
mov [rcx+8*r8],rax ; a[n+r8] <-- rax
inc r8 ; r8++
jz bneg_end ; si r8 == 0 goto end
align 16
bneg_loop2: ; boucle de copie en XOR #FFFFFFFFFFFFFFFF
mov rax,[rdx+8*r8]
xor rax,-1
mov [rcx+8*r8],rax
inc r8
jnz bneg_loop2
bneg_end:
mov rax,1 ; retourne 1
ret
BNeg ENDP
;//////////////////////////////////////////////////////////////////////////////
; A(n) <- -B(n)
;
; QWORD BNeg_ip(QWORD* a, long long n)
align 16
BNeg_ip PROC
lea rcx,[rcx+8*rdx]
neg rdx
align 16
bnegip_loop:
mov rax,[rcx+8*rdx]
or rax,rax
jnz bnegip_next
inc rdx
jnz bnegip_loop
mov rax,0
ret
bnegip_next:
neg rax
mov [rcx+8*rdx],rax
inc rdx
jz bnegip_end
align 16
bnegip_loop2:
mov rax,[rcx+8*rdx]
xor rax,-1
mov [rcx+8*rdx],rax
inc rdx
jnz bnegip_loop2
bnegip_end:
mov rax,1
ret
BNeg_ip ENDP
Un entier dit de taille \(n\) sera représenté par un tableau \(T\) de taille \(n\) et d'un bit de signe, noté \(s\): \[ N = {(-1)}^{s} \sum_{k=0}^{n-1} T[k] \; {\beta}^{k} \qquad 0 \leq T[k] < \beta \] Cette représentation est dite little endian car les mots de poids faible précèdent (en mémoire) les mots de poids plus élevé.
Un nombre représenté en virgule flottante (dit plus succintement: un flottant) de taille \(n\) sera également représenté par un tableau de taille \(n\), un exposant (noté \(E\) et représenté par un entier 64 bits signé) et un signe \(s\) : \[ x = {(-1)}^{s} \; {\beta}^{E} \; \sum_{k=0}^{n-1} T[k] \; {\beta}^{k-n} \qquad 0 \leq T[k] < \beta \] Remarquons que ces notations ne sont pas uniques pour un entier ou un flottant donné. Pour avoir l'unicité, il faut ajouter une condition à ces notations:
La représentation d'un entier ou d'un flottant est dite normalisée si le tableau associé, de taille \(n\), vérifie: \( T[n-1] \neq 0 \).
Notons que la représentation de \(0\) pose alors problème, aussi bien en tant qu'entier qu'en tant que flottant. Nous imposerons à la taille d'être strictement supérieure à \(0\) et donc le nombre \(0\) sera le seul à avoir une représentation non normalisé (en entier comme en flottant).
Structures C
struct st_mpInt
{
long long sgn;
long long size;
unsigned long long m[0];
};
typedef struct st_mpInt MpInt;
///////////////////////////////////////
struct st_mpFloat
{
long long sgn;
long long size;
long long exp;
long long prec;
unsigned long long m[0];
};
typedef struct st_mpFloat MpFloat;
Notons tout d'abord l'astuce de mettre un champ m tableau de taille \(0\) afin de le manipuler comme un tableau/pointeur sur la zone mémoire qui suit. Ce n'est pas une pratique recommandée, mais ça évite de faire une allocation supplémentaire et ça assure la contigüité mémoire. Il faudra évidemment bien faire attention à ne pas déborder sur la mémoire allouée.
Le champ sgn servira pour stocker le signe. Un seul bit suffirait pour cette information mais c'est plus simple comme ça. De plus m devra être aligné sur 16 octets.
Le champ size contiendra la taille allouée du tableau m et ne changera jamais tout au long de la "vie" de notre nombre (sauf réallocation).
Le champ exp sera l'exposant du nombre à virgule flottante et prec correspondra à la précision, c'est à dire au nombre de mots effectivement utilisés dans le champ m car dans le code nous aurons besoin de gérer des précisions variables.
Performances
L'objectif du code proposé dans cette section précision arbitraire n'est pas seulement pédagogique. Si vous cherchez des informations sur les méthodes de calcul en précision arbitraire dans un soucis exclusif de compréhension, beaucoup d'autres sites proposent ces explications de manière simple. Ici, une plus grande attention sera portée sur les performances (d'où l'usage d'assembleur) et cela se fera nécessairement au détriment de la clarté et de la simplicité. Mon souhait est surtout de montrer comment il est possible de s'approcher des implémentations les plus rapides du monde.
Dans le chapitre suivant, nous verrons comment implémenter du code assembleur à un projet sous Visual C++ et ensuite, nous attaquerons les routines de plus bas niveau.
Nous n'allons pas coder directement les routines d'addition, de soustraction, etc... sur les différentes représentations de nombres auxquelles nous aurons à faire, c'est à dire, les nombres entiers et les nombres à virgule flottante (que nous appellerons flottants).
Il est bien plus judicieux et performant d'implémenter une couche dite de bas niveau qui va traiter ce que nous appellerons des blocs. Nous allons en effet représenter nos nombres par des tableaux d'entiers non signés 64 bits.