retour index
Multiplication
Il existe un algorithme très efficace pour la multiplication sur des nombres à
rallonge. C'est l'algorithme de Karatsuba-Ofman [Bibl. 11 entre autres]
Tout nombre binaire (base 2) de m digits peut être représenté sous la
forme
a.2(m-1) + b.2(m-2) + .. + y .21 + z.20 .
Maintenant, augmentons la valeur de la base, jusqu'à 2(m/2) .
Pour simplifier, considérons que m est pair (le raisonnement est identique pour m impair, mais les deux parties ne sont pas de longueur égale).
Notre nombre peut alors se diviser en 2 parties :
abcd..2(m/2) + ..wxyz.20
On se retrouve avec deux parties, dont le degré de la plus forte a diminuer de moitié.
En généralisant, le nombre A de m digits devient Ah * X + Al, avec Ah la partie haute composée de la moitié haute des digits du nombre, et Al la partie basse composée du reste, et X la base équivalente à la moitié du degré du nombre.
Prenons un exemple en décimal: soit le nombre 100200.
Il est équivalent à 1.105 + 0.104 + 0.103 + 2.102 + 0.101 + 0.100
Ou encore : 100.103 + 200.100, == (Ah * X) + Al
100.103 correspond à Ah * X, avec X = 103 et Al représente 200.100 (100 = 1).
103 est la nouvelle base X égale à la moitié du degré du nombre A.
Prenons un deuxième nombre B de même longueur.
On lui applique le même traitement pour obtenir Bh * X et Bl.
Faisons la multiplication A * B, en remplaçant chaque nombre par son équivalent:
R = ((Ah * X) + Al) * (Bh * X) + Bl)
Le développement classique d'une telle relation est:
R = X2 (Ah.Bh) + X(Al.Bh + Ah.Bl) + Al.Bl
On constate que toutes les multiplications à effectuer sont de degré moitié des nombres d'entrée.
Les multiplications par X et X carré s'obtiennent tout simplement en décalant le coefficient concerné du nombre correct de digits vers la gauche, ce qui est très facile avec des registres de processeur.
En itérant le processus , on arrive à des multiplication de degré 1.
Allons plus loin. Reprenons le développement de notre multiplication. Le terme central, en X est composé de 2 multiplications.
Nous avons calculé le 2 termes extrêmes Ah.Bh, et Al.Bl.
Si on ajoute ces 4 termes, on obtient Ah.Bh + Ah.Bl + Al.Bh + Al.Bl.
Une mise en facteur nous donne immédiatement (Ah + Al) * (Bh + Bl).
On peut donc dire que le terme central en X est égal a (Ah + Al) * (Bh + Bl) moins les deux autres termes déjà calculés.
On a gagné une multiplication, remplacée par 2 additions et 2 soustractions.
Il est bien connu que les additions et soustractions sont beaucoup moins coûteuses en temps de calcul sur un processeur, qu'une multiplication.
En résumé, pour Karatsuba, il faut calculer 3 valeurs :
Q0 = Al*Bl ; Q2 = Ah*Bh et Q1 = (Al + Ah)(Bl + Bh)- Q0 - Q2; ;
Le résultat sera égal à la somme de 3 termes :
Q2 décalé de 2 fois la base
Q1 décalé d'une fois la base
Q0
Cette méthode permet de réduire le degré de chaque multiplication.
Gestion des dépassements :
Si on additionne 2 nombres d'une base donnée, on cours le risque d'avoir un résultat supérieur à cette
base. C'est un dépassement (une retenue). Dans un système qui utilise des données de longueur fixes,
le dépassement disparaît, puisqu'il dépasse la longueur permise. Il faut donc effectuer des corrections.
1er cas : Ah + Al est supérieur à la base :
Décomposons le résultat de cette addition en ses 2 parties :
la partie haute donnant le nb de fois la base (h), et la partie basse (l), le dépassement de la base :
0xA0 + 0x82 = 0x122 c'est-à-dire 1 fois la base + 22 :
xh + l, dans lequel h ne peut valoir que 0 ou 1.
Reprenons la multiplication, avec (Bh + Bl) :
(x + l) * (Bh + Bl) == x(Bh + Bl) + l(Bh + Bl)
Le terme l(Bh + Bl) est celui que nous avons obtenu au cours de la multiplication (Ah+Al)(Bh+Bl).
Il suffit donc d'additionner au résultat obtenu, (Bh+ Bl) décalé de la base.
2ème cas : Bh + Bl est supérieur à la base :
Ce cas est identique au précédent, mais cette fois, on additionne (Ah + Al).
3ème cas : Il se peut qu'au cours des additions correctives précédentes avec Q1 on obtienne aussi un dépassement.
Le résultat s'exprime comme suit :
Ce schéma montre clairement ce qu'il convient de faire en cas de dépassement sur Q1 :
La retenue doit être additionnée à Q2l.
Dans le programme qui va suivre, on imbrique additions et soustractions. Une variable supplémentaire sera utilisée pour incrémenter la retenue sur dépassement d'addition, et la décrémenter sur dépassement de soustraction.
Dans le cas d'opérations sur des nombres très longs, l'algorithme est utilisé d'une façon récursive, c'est-à-dire que l'on coupe en deux le mot d'entrée, puis si celui-ci est encore trop long pour le calculateur, on continu la découpe, jusqu'à ce que l'on atteigne la taille de mot manipulée par le calculateur. Plus cette taille sera grande, moins d'itération il y aura.
Le programme présenté ci-dessous fait appel au fait que dans un programme en C, une fonction peut s'appeler elle-même.
La structure d'une fonction pour multiplier 2 octets par 2 octets est la même que celle pour multiplier 4 octets par 4 octets, et ainsi de suite. Seule change la dimension et la positions des blocs de mémoire utilisés.
La fonction ci-dessous est appelée avec 3 paramètres :
byte *p pointe sur le buffer contenant premier opérande (le multiplicande).
byte *t pointe sur le buffer contenant le second opérande (le multiplicateur)
int size est la dimension de chaque buffer à l'entrée de la fonction.
Les deux buffers d'entrée auront la dimension donnée par size.
Ces deux buffers devront se suivre en mémoire, car le résultat est rendu dans le buffer commençant
à p, et de longueur 2 fois size.
Deux fonctions annexes sont appelées par mul(). Il s'agit de adc() et sbc(), qui sont
respectivement l'addition et la soustraction, qui prennent comme paramètres:
2 pointeurs sur les buffers d'entrées, 1 pointeur sur le buffer de résultat, un int donnant le
nombre d'octets à travailler, et un byte pour un éventuel report d'entrée.
Chaque fonction retourne 1 ou 0, selon qu'il y a eu dépassement ou pas.
/************************ mul **************************/
void mul(byte *p, byte *t, int size)
{
int n, res;
byte carry1;
byte carry2;
byte carry, x = 0;
byte temp, nb;
byte *s1 = NULL, *s2, *pr; /* pointeur sur bloc-notes des 2 sommes, et du produit de Q0*/
nb = size / 2;
switch(size)
{
case 1: /* mul8 */
res = *p * *t;
*p = res >> 8;
*(p+1) = (byte)res;
return;
default:
s1 = (byte *)malloc(size);
pr = (byte *)malloc(size);
break;
}
s2 = s1 + nb;
carry2 = adc(t, t+nb, s2, nb, 0);
carry1 = adc(p, p+nb, s1, nb, 0);
memcpy(pr, s1, size); /* copie les 2 sommes pour multiplication */
mul(pr, pr+nb, nb);
/* pr = (Ah+Al)(Bh+Bl) */
if(carry2)
{
x = adc(pr, s1, pr, nb, 0);
if(carry1)
{
x++;
carry = adc(pr, s2, pr, nb, 0);
if(carry)
x++;
}
}
else if(carry1)
{
carry = adc(pr, s2, pr, nb, 0);
if(carry)
x++;
}
for(n = 0; n < nb; n++)
/* on aurait pu mettre cela en fonction externe */
{
temp = *(p+nb + n);
*(p+nb+n) = *(t+n);
*(t+n) = temp; /* swap Al et Bh */
}
mul(p, p+nb, nb); /* p = Ah.Bh */
carry = sbc(pr, p, pr, size, 0);
if(carry)
x--;
mul(t, t+nb, nb); /* t = Al.Bl */
carry = sbc(pr, t, pr, size, 0);
if(carry) /* ici, on a Q1 */
x--;
carry = adc(p + nb, pr, p + nb, size, 0);
res = *(p + (nb - 1)) + x; /* corrige les dépassements antérieurs */
if(carry) /* si dépassement */
res++;
carry = res / 256; /* si dépassement carry = 1, sinon = 0 */
*(p + (nb-1)) = (byte)res; /* range la partie basse */
if(carry)
{
p += nb-2;
do {
(*p)++;
} while (!(*(p--)));
}
if(s1)
{
free(s1);
free(pr);
}
}
/**************************** adc *******************************/
byte adc(byte *p, byte *t, byte *d, int nb, byte carin)
{
byte carry;
int res;
carry = carin;
nb--;
for(; nb >= 0; nb--)
{
res = *(p+nb) + *(t+nb) + carry;
carry = (res > 255) ? 1 : 0;
*(d+nb) = (BYTE)res;
}
return(carry);
}
/***************************** sbc ****************************/
byte sbc(byte *p, byte *t, byte *dest, int nb, byte incar)
{
int res, ressub, n;
byte carry;
carry = incar;
nb--;
for(n = nb; n >= 0; n--)
{
ressub = *(t + n) + carry;
res = *(p + n) - ressub;
carry = (ressub > *(p + n)) ? 1 : 0;
*(dest + n) = (BYTE)res;
}
return(carry);
}
Les bloc notes sont créés dynamiquement, à chaque appel de la fonction, de façon que ces derniers ne s'écrasent pas. La dimension des bloc notes dépend de la dimension des buffers d'entrée.
Le principe de l'algorithme de Karastuba-Ofman est d'utiliser une multiplication de degré moitié des datas d'entrée. Par exemple, pour multiplier 4 octets par 4 octets, il faudra 2 buffers de 2 octets pour les sommes partielles et un buffer de 4 octets pour les multiplications partielles. Cette notion de dimension moitié est utilisée tout le long de la fonction. On défini immédiatement une variable nb, égale à la moitié de la dimension d'entrée.
Les buffers temporaires sont alloués dynamiquement, et supprimés à la sortie.
Deux buffers pour les sommes, s1 et s2, et un pour les produits partiels pr.
Le switch(size) du début permet d'utiliser la même fonction pour la multiplication élémentaire 1 octet par 1 octet.
Ensuite, on effectue l'opération (Ah+Al)(Bh+Bl), en appliquant les corrections d'éventuels report, comme il est expliqué plus haut. La multiplication nécessaire est faite en appelant de nouveau cette même fonction dans laquelle on est actuellement. Ce nouvel appel crée de nouveaux buffers différents.
Puis, on calcul les produits Ah.Bh, et Al.Bl, après un swap de Al et Bh, disposant les bonnes datas aux bons endroits, pour les multiplications.
Après les deux soustractions, on dispose de Q1 dans le buffer pr.
Dans les buffers d'entrée, sont disposés à la suite, AhBh et AlBl. Il suffit d'additionner Q0 à partir de Bh, pour obtenir le résultat souhaité, conformément au schéma ci-dessous
Après les corrections entraînées par les différents reports, on libère la mémoire des deux buffers temporaires.
Application à un processeur 8 bits
Malheureusement, dans un processeur industriel, simple 8 bits, on ne peut pas faire simplement d'allocation dynamique de mémoire. Il va donc falloir revoir le problème avec des buffers fixes.
Ceci nous oblige à abandonner la récursivité sur la même fonction, et écrire autant de fonctions qu'il y aura de multiplications élémentaires
J'ai choisi pour l'exemple, deux processeurs RISC de la firme Microchip. Le premier est un 17C44 qui
possède une instruction de multiplication sur l'octet, donnant un résultat sur 2 octets. C'est donc cette valeur qui nous servira de base. Pour la multiplication 4 octets * 4 octets, il faudra donc développer la multiplication 2 octets * 2 octets.
Pour plus de détails, lire les bibliographies indiquées. D'autre part, la longueur des mots d'adresse
est de 16 bits, permettant une programmation linéaire, sans changement de page. En plus, il possède deux pointeurs
facilitant les manipulations de données.
Le second choisi, 16F876, complique le problème par l'absence d'instruction de multiplication 8*8, de l'obligation de pagination
de la mémoire programme, et n'a qu'un seul pointeur.
Pour ce dernier, le code est un peu plus long, et un peu plus lent, mais fonctionne également correctement, et la
disposition en RAM est différente.
Pour des multiplications de 32 * 32 octets, comme on vient de le dire, le problème va être scindé en morceaux :
multiplication de 16 * 16, puis de 8 * 8, et enfin, de 4 * 4. C'est par cette dernière que nous allons commencer.
Souvenons-nous que l'on doit non seulement réduire le temps d'exécution, mais en plus, utiliser un minimum de mémoire RAM.
Le résultat de la multiplication de 32 octets par 32 octets, est une sortie sur 64 octets.
Plusieurs fonctions doivent être établies:
La multiplication M256 de 32 * 32 appelle la multiplication M128 de 16 * 16 qui appelle la multiplication M64 de 8 * 8, qui appelle la multiplication M32 de 4 * 4, qui appelle la multiplication M16 de 2 * 2, avant de trouver la multiplication de 1 octet * 1 octet. Cette cascade d'appels nous oblige à maintenir des buffers indépendants pour chaque opérations.
Il faut :
2 buffers de 32 octets pour M256
2 buffers de 16 octets pour M128
2 buffers de 8 octets pour M64
2 buffers de 4 octets pour M32
2 buffers de 2 octets pour M16
soit 124 octets. Il reste 4 octets pour un bloc-notes de la multiplication M16.
On a aussi besion de quelques octets pour la gestion des dépassements.
Le 17C44 possède une RAM en 2 pages. La première de 0 à 0xFF, et la seconde de 0x120 à 0x1FF.
On va utiliser la première page, afin d'éviter des incessants changements de page, entre 0x80 et 0xFF.
On se réserve la partie basse pour d'autres opérations.
Pour le 16F876, la répartition des buffers est un peu plus complexe. Chaque page de RAM comporte 80 octets
utilisable facilement, plus une zone d'octets non paginée.
La répartition des buffers se passe de la façon suivante:
De 0x80 à 0xBF, les 2 buffers de M256
De 0xC0 à 0xDF, les 2 buffers de M128
De 0xE0 à 0xEF, les 2 buffers de M64
De 0xF0 à 0xF7, les 2 buffers de M32
De 0xF8 à 0xFB, les 2 buffers de M16
Et enfin, de 0xFC à 0xFF, le bloc-notes des opérations.
On va voir maintenant comment ordonner les opérations afin d'utiliser ces buffers le plus astucieusement possible.
Et comme un dessin vaut mieux qu'un long discours, les schémas joints donnent les plans de
fonctionnement pour la multiplication de 256 bits x 256 bits, résultat sur 512 bits.
Les buffers sont alignés de sorte que le bloc-notes de chaque multiplication correspond aux buffers d'entrée
de l'opération de degré inférieur.
Ainsi, pour la multiplication de 256 * 256, la fonction de départ de la chaîne, les datas d'entrée sont dans
les 2 buffers situés de 0x80 à 0x9F, et de 0xA0 à 0xBF.
Les calculs intermédiaires
arrivent dans les buffers 0xC0 à 0xCF et 0xD0 à 0xDF qui seront les buffers d'entrée de la multiplication
de 128 * 128, et ainsi de suite.
Dans les dessins, les flèches indiquent les mouvements de datas. Il y a des échanges de place (swap), permettant
de positionner les valeurs aux bons endroits pour les opérations suivantes, telles que les décalages, par exemple.
Les buffers sont symbolisés par des rectangles, dans lesquels les chiffres indiqués à droite et
à gauche sont les index des registres (par exemple, 80 indique qu'il s'agit de l'octet d'offset 0x80).
Le signe '+' dans un cercle indique la fonction d'addition ADC, addition avec carry.
De même le signe '-' dans un cercle indique la soustraction SBC, soustraction avec report.
Les listings de programme dépendent du processeur choisi.
Nota
Je ne donne en exemple, que le listing relatif au 16F876.
Mais, pour ceux qui ont quelque chose entre les oreilles, cela ne devrait pas poser de problèmes, d'autant plus qu'il
vaut bien mieux travailler sur son propre code, que d'utiliser un code écrit par quelqu'un d'autre, où très
souvent, les différentes astuces ne sont pas forcément bien comprises.
Bibliographie:
[10] Microchip : le site officiel de Microchip oú l'on peut trouver
toutes les datas sheets de processeurs.
[11] Algorithme de KARATSUBA pour la multiplication rapide
Il suffit de demander à Google "Algorithme de KARATSUBA" pour avoir quantité de site très bien fournis
[12] D.E. Knuth, The art of computer programming 3rd edition vol. 1 pages 13 et 14
retour index