Racine carrée en virgule fixe

De Wiki LOGre
Aller à la navigation Aller à la recherche


Projet réalisé par : Edgar Bonet.

Ceci est une tentative d'implémenter la racine carrée en virgule fixe.

Méthodes itératives

On peut calculer la racine carrée par des méthodes d'approximations successives, en partant d'une approximation grossière. Les polynômes d'approximation de faible degré sont des bons candidats pour fournir l'approximation initiale.

Méthode de Héron

Cette méthode a été décrite par Héron d'Alexandrie au Ier siècle ap. J.-C. Elle se base sur le fait que, si un est une approximation de [math]\displaystyle{ \sqrt{x} }[/math], alors x/un en est une autre, et la valeur exacte est la moyenne géométrique de ces deux approximations :

[math]\displaystyle{ \sqrt{x} = \sqrt{u_n \times \frac{x}{u_n}} }[/math]

Aussi, quand deux nombres sont proches, leur moyenne géométrique est très proche de leur moyenne arithmétique. On obtient donc une approximation améliorée en prenant la moyenne arithmétique des deux approximations précédentes :

[math]\displaystyle{ u_{n+1} = \frac{u_n + x/u_n}{2} }[/math]

Cette méthode a une convergence quadratique : le nombre de chiffres corrects est environ doublé à chaque itération. Pour plus de détails, voir Méthode de Héron dans Wikipedia.

Méthode de Newton

On peut aussi utiliser la méthode de Newton pour trouver le zéro (en y) de la fonction

y2 − x = 0

En développant le calcul, il se trouve que cette méthode est identique à la méthode de Héron ci-dessus.

Polynôme d'approximation

Réduction de l'intervalle de travail

En utilisant la relation [math]\displaystyle{ \sqrt{2^{2n}x} = 2^n\sqrt{x} }[/math], on peut se ramener à l'intervalle [0.25, 1] avec des décalages de bits.

Contraintes

On cherche un polynôme qui soit exact aux deux extrémités de l'intervalle :

  • P(1/4) = 1/2
  • P(1) = 1

Le polynôme le plus simple qui satisfait ces contraintes est

  • P1(x) = 1/3 + 2/3 x

À celui-ci on peut ajouter tout polynôme vérifiant P(1/4) = P(1) = 0, c'est à dire

  • P2(x) = (x−1/4)(1−x)Q(x)

où Q(x) est un polynôme quelconque. On cherche donc une solution de la forme

  • P(x) = 1/3 + 2/3 x + (x−1/4)(1−x)Q(x)

Il s'ensuit que Q(x) approxime [math]\displaystyle{ (\sqrt{x} - 1/3 - 2/3\;x)/\big((x-1/4)(1-x)\big) }[/math].

Polynômes optimaux

Pour les petits degrés, il est possible de trouver l'optimum par tâtonnements. Par exemple, pour trouver P(x) de degré 3, on cherche Q(x) de degré 1, et il n'y a donc que deux paramètres à ajuster. À partir du degré 4 en P(x), cette méthode devient pénible. J'ai alors utilisé une variante de l'algorithme de Remez.

Voici les polynômes optimaux pour les degrés 1 à 5 :

degré polynôme erreur max.
1 0.3333333 + 0.6666667 x 4.16670 × 10−2
2 0.2578607 + 1.0440300 x − 0.3018907 x2 5.45518 × 10−3
3 0.2170856 + 1.3158433 x − 0.8046814 x2 + 0.2717525 x3 1.02489 × 10−3
4 0.1908901 + 1.5392302 x − 1.4475870 x2 + 1.0217778 x3 − 0.3043110 x4 2.25604 × 10−4
5 0.1722984 + 1.7336691 x − 2.2038443 x2 + 2.3952900 x3 − 1.4780286 x4 + 0.3806153 x5 5.43654 × 10−5

À titre de comparaison, pour un format 16 bits 0.16, la résolution numérique est de 2−16 ≈ 1,53 × 10−5.

Erreur

Erreur d’approximation : [math]\displaystyle{ P_4(x) - \sqrt{x} }[/math]

Ces approximations ont la particularité d'avoir une erreur qui oscille, avec tous les extrema qui ont la même valeur absolue. Ceci est similaire au comportement des polynômes de meilleure approximation, avec deux différences qui sont dues aux contraintes imposées :

  • l'erreur du polynôme de degré n présente n extrema, contre n+2 pour le polynôme de meilleure approximation ;
  • l'erreur s'annule par construction aux extrémités de l'intervalle.

La figure ci-contre montre l'erreur du polynôme de degré 4.

Implémentations

Polynôme de degré 4

Pour utiliser ce polynôme, il faut d'abord multiplier l'argument par 4 autant de fois que nécessaire pour arriver dans l'intervalle [0.25, 1[. Ceci se fait en décalant de deux bits vers la gauche (x <<= 2). Une fois le polynôme évalué, le résultat doit être divisé par deux (y >>= 1) autant de fois que l'argument a été multiplié par 4. Il faut bien sûr commencer par s'assurer que l'argument est non nul. Voici le code: 

#include "fixed_point.h"

/* Max error = 2.56e-4. */
uint16_t sqrt_fix(uint16_t x)
{
    if (x == 0) return 0;
    uint8_t n;
    for (n = 0; x < 0x4000; n++) x <<= 2;  // scale the operand
    uint16_t y;
    y = FIXED(0.30435, 15);
    y = FIXED(1.02173, 15) - mul_fix_u16(x, y);
    y = FIXED(1.44754, 15) - mul_fix_u16(x, y);
    y = FIXED(1.53925, 15) - mul_fix_u16(x, y);
    y = FIXED(0.19089, 16) + (mul_fix_u16(x, y) << 1);
    while (n--) y >>= 1;                   // scale the result
    return y;
}

À cause des erreurs d'arrondi, le polynôme qui est optimal sur des réels ne l'est pas tout à fait sur les nombres à virgule fixe. J'ai donc ajusté légèrement les coefficients pour qu'ils soient optimaux, d'où les légères différences entres ces coefficients et ceux du tableau précédent. L'erreur maximale est 2.56 × 10−4, soit seulement 13% de plus que le polynôme optimal sur des réels.

Le temps d'exécution dépend de l'argument. Pour x ∈ [0.25, 1[, il est de 179 cycles, soit 11,19 µs environ à 16 MHz. Pour des x plus petits, il augmente de 17 cycles fois le nombre de fois que l'argument a dû être multiplié par 4, jusqu'à un maximum de 298 cycles (≈ 18,63 µs). Pour x = 0 ça prend 27 cycles.

x cycles µs @ 16 MHz
0 27 1.6875
[1/65536, 1/16384[ 298 18.625
[1/16384, 1/4096[ 281 17.5625
[1/4096, 1/1024[ 264 16.5
[1/1024, 1/256[ 247 15.4375
[1/256, 1/64[ 230 14.375
[1/64, 1/16[ 213 13.3125
[1/16, 1/4[ 196 12.25
[1/4, 1[ 179 11.1875
moyenne 184.7 11.54

Polynôme de degré 4 + itération de Héron

Le polynôme de degré 4 est plus précis que nécessaire si on veut appliquer une itération de Héron. Mais comme je l'ai déjà implémenté, c'est plus simple de partir de celui-ci. Cela permet aussi de voir l'effet de l'amélioration itérative à la fois sur la précision et le temps d'exécution. 

#include "fixed_point.h"

/* Max error = 7.67e-6. */
uint16_t sqrt_fix(uint16_t x)
{
    if (x == 0) return 0;
    uint8_t n;
    if (x == 0xffff) return 0xffff;
    uint16_t x0 = x;
    for (n = 0; x < 0x4000; n++) x <<= 2;  // scale the operand
    uint16_t y;
    y = FIXED(0.30435, 15);
    y = FIXED(1.02173, 15) - mul_fix_u16(x, y);
    y = FIXED(1.44754, 15) - mul_fix_u16(x, y);
    y = FIXED(1.53925, 15) - mul_fix_u16(x, y);
    y = FIXED(0.19089, 16) + (mul_fix_u16(x, y) << 1);
    while (n--) y >>= 1;                   // scale the result
    y = (y + ((uint32_t)x0 << 16) / y + 1) / 2;
    return y;
}

L'erreur maximale est 7.67 × 10−6, soit seulement 0.502 fois le pas de discrétisation. On est très près de la limite imposée par le format de virgule fixe choisi.

Le temps d'exécution est beaucoup plus irrégulier que dans le cas précédent. Il est minimal pour x = 0 et maximal pour x = 3.

x cycles µs @ 16 MHz
min (0) 46 2.875
max (3) 937 58.5625
moyenne 828.1 51.76

On peut constater que, si la précision est presque parfaite, cette simple itération coûte beaucoup plus de temps de calcul (environ 3.5 fois plus) que l'évaluation du polynôme. J'attribue ça à l'inefficacité de la division qui, contrairement à la multiplication, ne bénéficie pas de support matériel.

Conclusion provisoire : Il vaut mieux adapter le degré du polynôme à la précision voulue que d'essayer l'amélioration itérative.

Récupérée de « https://wiki.logre.eu/index.php?title=Racine_carrée_en_virgule_fixe&oldid=13460 »