Racine carrée en virgule fixe
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
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.