Calculs en virgule fixe

De Wiki LOGre
Aller à : navigation, rechercher

Projet réalisé par : Edgar Bonet.

Le but de ce projet est de réaliser des calculs sur des petits processeurs, tels les AVR qu'on trouve dans les Arduinos, sans être pénalisés par la lenteur et la lourdeur du code en virgule flottante. Cette page comporte une introduction et une routine de multiplication. La suite, qui est dans une autre page, est le calcul de fonctions trigonométriques.


Introduction à la virgule fixe

D'abord une introduction pour ceux qui ne sont pas familiers de la virgule fixe. Les calculs en virgule flottante sont très lents sur un processeur dépourvu de FPU, qui plus est travaillant en 8 bits. Un moyen souvent efficace pour accélérer les choses consiste à travailler en virgule fixe. Le principe est de multiplier la grandeur qu'on veut manipuler par un facteur constant (une sorte de changement d'unités), de sorte que celle-ci puisse être stockée avec une résolution satisfaisante dans un entier. Les nombres à virgule fixe ne sont alors, somme toute, que des entiers.

Exemples en décimal

L'exemple classique est le calcul en unités monétaires. Pour faire les comptes de ce qui rentre et sort de son compte en banque, on peut représenter les montants en euros avec des nombres ayant jusqu'à deux chiffres après la virgule. Mais si on représente les montants en centimes, alors on n'a plus que des entiers à manipuler. Il est toujours possible d'ajouter ou enlever les virgules dans les entrées/sorties, et donner l'illusion que les calculs se font avec des nombres à virgule.

Sinfix-odometer.jpg

Un autre exemple, peut-être plus instructif, est l'odomètre de voiture de la photo ci-contre. Cet odomètre compte une distance en kilomètres avec un chiffre après la virgule. Pourtant, il est architecturalement identique à un compteur d'entiers. Il n'est un compteur « à virgule » que par la vertu d'une convention qui veut qu'en le lisant, on imagine une virgule implicite avant le chiffre qui est sur fond blanc. Selon cette convention, cet odomètre affiche 653,6 km. Mais on peut tout à fait ignorer cette convention, considérer qu'il s'agit d'un compteur d'hectomètres, et lire un entier : 6536 hm. Les deux lectures sont aussi légitimes l'une que l'autre. Les leçons à retenir de cet odomètre, et qui se généralisent aux calculs en binaire, sont :

  • architecturalement, c'est un compteur d'entiers ;
  • la virgule ne fait pas partie du compteur lui-même, elle est uniquement dans l'interprétation qu'on fait de la valeur affichée ;
  • la valeur affichée a plusieurs interprétations légitimes, selon l'unité qu'on lui attribue et la position (ou absence) de la virgule implicite.

En binaire

Pour stocker un nombre x en virgule fixe, on le multiplie par une puissance de deux (2n) et on stocke le résultat, arrondi, dans un entier. On peut alors considérer que les n bits de poids faible de l'entier stocké représentent la partie fractionnaire de x. Le format de virgule fixe utilisé est alors défini par le choix de n, ainsi que type d'entier utilisé pour stocker le produit 2nx. Et comme il semble que dans ce domaine chacun a sa notation pour désigner ces formats, voici donc la mienne : je note les formats à virgule fixe sw,n et uw,m où :

  • s désigne un format signé et u (comme unsigned) un non-signé ;
  • w est la largeur du nombre en bits ;
  • n est le nombre de bits après la virgule.

À noter qu'on rencontre souvent aussi les notations m.n et simplement .n, où m est le nombre de bits avant la virgule, et n le nombre de bits après. Mais ce notations sont ambigües : elles ne précisent pas si le format est signé ou non et, quand il est signé, le bit de poids fort (bit de signe) est parfois compté dans m, mais pas toujours.

Un exemple : pour représenter 5.375 au format s8,4, on calcule 5.375 × 2⁴ = 86, et on le range ce résultat dans un int8_t. La séquence de bits (01010110) sera alors lue « 0101.0110 », c'est à dire :

position 3 2 1 0 −1 −2 −3 −4
bit 0 1 0 1 0 1 1 0
poids −8 4 2 1 1/2 1/4 1/8 1/16

Ici le bit de poids fort a un poids négatif car c'est aussi le bit de signe (représentation en complément à deux). Si on utilisait le format le format u8,4, le bit de poids fort aurait un poids +8.

Par la suite je m'intéresse principalement aux formats 16 bits.


Multiplication en virgule fixe

Je vais approximer la fonction sinus par un polynôme, c'est à dire une séquence d'additions et de multiplications. L'addition ne pose pas de problème en virgule fixe : on fait si nécessaire des décalages de bits pour mettre les termes au même format, et une fois au même format on additionne comme des entiers.

La multiplication, par contre, demande un traitement spécifique. En principe, quand on multiplie deux nombres, le produit a autant de bits que les deux facteurs réunis. Et quand on multiplie un nombre au format uw,n par un autre au format uw′,n′, on obtient un résultat au format u(w+w′),(n+n′). Mais il n'est pas tenable d'avoir une taille de données qui grossit à chaque multiplication. Aussi, on voudra généralement réduire la taille du résultat en abandonnant des bits de poids faible. La difficulté est qu'en langage C, ce sont les bits de poids fort qui sont automatiquement perdus !

Multiplication en C

Lorsqu'on multiplie, en C, deux uint16_t, on ne récupère que les 16 bits de poids faible du résultat, c'est à dire le résultat modulo 216. Pour avoir le résultat complet en 32 bits, il est tentant d'affecter le produit à un uint32_t :

uint16_t x, y;
uint32_t z = x * y;

mais cela ne produit pas le résultat escompté. Le compilateur interprète cette multiplication comme une multiplication 16 bits. Il produit donc un résultat tronqué à 16 bits, qu'il étend ensuite à 32 bits en ajoutant des zéros à gauche. Si on veut le résultat complet, il faut convertir les deux opérandes en 32 bits avant d'effectuer la multiplication. Seul l'un des opérandes a besoin d'être converti explicitement, la conversion étant alors implicite pour le deuxième :

uint16_t x, y;
uint32_t z = (uint32_t) x * y;

Et si on peut se contenter d'une précision de 16 bits pour le résultat, alors on voudra probablement les bits de poids fort, qui sont les plus significatifs. On tronque alors le résultat par un décalage à droite. Ou mieux, on l'arrondit au plus proche en ajoutant, avant le décalage, la moitié du poids du bit de poids faible du résultat final. Ce qui nous donne cette fonction pour la multiplication 16 bits :

/*
 * Fixed point multiplication.
 *
 * Multiply two unsigned FP numbers in u16,16 (unsigned 0.16) format.
 * Returns result in the same format.
 * Rounds to nearest, ties rounded up.
 */
static uint16_t mul_fix_u16(uint16_t x, uint16_t y)
{
    return ((uint32_t) x * y + 0x8000) >> 16;
}

À noter que, bien que le commentaire précise que les opérandes doivent être au format u16,16, la fonction convient à tous les formats u16 : si x est au format u16,n et y au format u16,n′, alors le résultat sera au format u16,(n+n′-16).

En assembleur

La multiplication ci-dessus marche parfaitement, mais elle est inefficace sur un processeur 8 bits. Les AVR équipés d'un « enhanced CPU core » disposent tous d'un multiplieur câblé (8 bits × 8 bits → 16 bits), et de l'instruction mul qui va avec. Pour multiplier des nombres plus grands, il faut décomposer chaque facteur en octets et multiplier chaque octet d'un facteur par chaque octet de l'autre, puis additionner les résultats en tenant compte des décalages. Il faut donc en principe 16 instructions mul, et un bon paquet d'additions (add et adc), pour calculer un produit 32 bits × 32 bits → 64 bits. En pratique, comme le compilateur ne garde que les 32 bits de poids faible, il lui suffit de 10 mul.

Mais 10 mul, à deux cycles pièce sans compter les additions, c'est déjà trop. Le problème initial était de multiplier des facteurs de 16 bits, et pour cela 4 mul suffisent. Le code C produit donc, en plus des 4 mul utiles, 6 autres mul qui donnent zéro, car ils concernent l'extension à 32 bits des opérandes. Le compilateur (gcc -Os) ne sait malheureusement pas optimiser ce cas, et il est donc utile de réécrire cette fonction en assembleur. Alors j'ai commencé par chercher sur le Web où j'ai trouvé quelques documents utiles :

  • Arduino – AVR GCC multiplication : un billet du blog de Norbert Požár, qui explique le problème des multiplications en C et donne du code asm inline pour plusieurs formats de multiplication ;
  • Binary hardware multiplication in AVR Assembler : un chapitre du tutoriel d'assembleur AVR de Gerhard Schmidt, qui détaille bien le code assembleur et donne des routines pour des multiplications 16×8→24, 16×16→32 et 16×24→40 bits.
  • 16 X 16 bit multiply to a 32 bit result : une discussion sur le forum d'AVR freaks avec, vers le milieu de la page, un post de bldrcowboy qui donne du code asm inline pour des multiplications 16×16→32 et 8×16→24.

Je me suis servi de ces documents pour comprendre les bases de l'assembleur AVR, et comment on fait les multiplications. Puis j'ai écrit ma propre routine, qui doit ressembler à certaines de celles citées plus haut, mais qui a eu pour moi l'avantage de me faire apprendre des choses.

On commence par « poser » la multiplication, comme on le ferait pour une multiplication en décimal sur papier, mais en base 256 au lieu de base 10. Chaque opérande et le résultat se décomposent donc en deux octets : x = xH:xL, où l'indice H désigne le poids fort et L le poids faible. De même pour y et z (le produit). Il faut tenir compte du fait que chaque produit d'octets occupe deux octets. Cela donne :

                xH   xL
            ×   yH   yL
          -------------
               (xL * yL)
          (xL * yH)
          (xH * yL)
   + (xH * yH)
   --------------------
      (produit 32 bits)
            +  0x80
   --------------------
      zH   zL

La dernière addition (0x80) est pour ajuster l'arrondi, et on ne récupère que les deux octets de poids fort du résultat final.

Ensuite, l'affectation des registres :

  • x est stocké dans la paire de registres r25:r24 (poids fort dans r25), car d'après les conventions d'appel d'avr-gcc, c'est dans cette paire de registres que la fonction va récupérer son premier argument ;
  • y va dans r23:r22 pour la même raison : ces registres contiennent le deuxième argument ;
  • le résultat des additions sera calculé dans r21:r20:r19, qui sont disponibles (caller-saved) ;
  • le résultat final devra être remis dans r25:r24, où la valeur de retour de la fonction est attendue ;
  • le registre r18 sera mis à zéro, car il faudra additionner des zéros pour propager les retenues entre les colonnes.

En tenant compte de ces affectations, en ajoutant explicitement le registre r18 (à zéro) qui devra être additionné, et en réorganisant les lignes, on peut reposer l'opération ainsi :

               r25  r24
           ×   r23  r22
     ------------------
     (r23*r25) 0x80
     r18  r18  (r22*r24)
     r18  (r22*r25)
   + r18  (r23*r24)
   --------------------
     r21  r20  r19

À noter qu'on est obligés de calculer le troisième octet de la somme (r19), même si celui-ci ne contribue au résultat final que des retenues propagées. Le quatrième octet, en revanche, ne peut pas propager de retenue. L'opération s'écrit alors en pseudo-code comme ceci :

   r21:r20      = r23 * r25
   r19          = 0x80
   r18          = 0
   r21:r20:r19 += r18:r18:high(r22 * r24)
   r21:r20:r19 += r18:(r22 * r25)
   r21:r20:r19 += r18:(r23 * r24)
   r25:r24      = r21:r20

L'implémentation utilise les instructions suivantes du jeu d'instructions AVR :

  • mul ra, rb (2 cycles): r1:r0 ← ra * rb, le résultat va toujours dans r1:r0, on n'a pas le choix ;
  • add ra, rb (1 cycle): ra ← ra + rb, additionne deux registres sans retenue entrante, la retenue sortante est sauvegardée dans le registre d'état ;
  • adc ra, rb (1 cycle): ra ← ra + rb + carry, additionne deux registres plus la retenue de la dernière opération ;
  • clr ra (1 cycle): ra ← 0, c'est un alias de eor ra, ra (ou exclusif) qui met à zéro un registre ;
  • ldi ra, const (1 cycle): ra ← const, chargement d'un registre par une constante (valeur « immédiate ») ;
  • movw ra, rb (1 cycle): r(a+1):ra ← r(b+1):rb, copie d'une paire de registres ;
  • ret (4 cycles): retour de fonction.

Finalement, le code assembleur :

mul_fix_u16:
    mul  r25, r23  ; r1:r0 = r25 * r23
    movw r20, r0   ; r21:r20 = r1:r0
    ldi  r19, 0x80 ; r19 = 0x80
    clr  r18       ; r18 = 0
    mul  r24, r22  ; r1:r0 = r24 * r22
    add  r19, r1   ; r19 += r1
    adc  r20, r18  ; r20 += r18 + carry
    adc  r21, r18  ; r21 += r18 + carry
    mul  r25, r22  ; r1:r0 = r25 * r22
    add  r19, r0   ; r19 += r0
    adc  r20, r1   ; r20 += r1 + carry
    adc  r21, r18  ; r21 += r18 + carry
    mul  r24, r23  ; r1:r0 = r24 * r23
    add  r19, r0   ; r19 += r0
    adc  r20, r1   ; r20 += r1 + carry
    adc  r21, r18  ; r21 += r18 + carry
    movw r24, r20  ; r25:r24 = r21:20
    clr  r1        ; r1 = 0
    ret            ; return r25:r24

L'instruction clr r1 vers la fin est là en vertu d'une convention qui veut que ce registre soit à zéro, et soit remis à zéro par tout code qui le modifie. Cette fonction prend 38 octets de flash et s'exécute en 26 cycles. En pratique elle sera réécrite en assembleur inline, ce qui donne à gcc plus de latitude d'optimisation pour le code qui va autour.

Avec cette routine de multiplication, on peut évaluer de façon efficace des polynômes qui approximent les fonctions usuelles :