Calculs précis et rapides pour les nombres à virgule flottante en utilisant l'exemple de la fonction sinus. Partie 2: libm

Je continue la série d'articles sur le travail avec la virgule flottante. Dans le premier article, j'ai donné une petite introduction mathématique et montré la manière la plus simple et la plus évidente de calculer le sinus avec des exemples de programmes avec différents «pièges». L'article d'aujourd'hui aura un style légèrement différent. Il n'y aura pas de pratique ici, mais nous approfondirons les mathématiques et entrerons dans le saint des saints - le code de la bibliothèque standard. Je donnerai également une réponse à la question à la fin du premier article. Alors allons-y.

Encore quelques maths

De toute évidence, plus la série de Taylor diminue rapidement en valeur absolue, moins il faut de termes pour atteindre la précision requise. Et donc, il semble, le résultat sera plus précis (ci-dessous, cela sera discuté plus en détail). À titre de comparaison, par exemple, prenons un terme du septième degré de la série de Taylor ($$$\frac{x^7}{7!}$) à $$$x=1.0$ et $$$x=0.1$... Les valeurs de l'expression seront$$$1.198\cdot 10^{-4}$ et $$$1.198\cdot 10^{-11}$respectivement. Une énorme différence, n'est-ce pas? Essayons donc de trouver un moyen de réduire la limite supérieure de l'intervalle de calcul de la fonction sinus.

Expansion de la série autour de valeurs données

Pour comprendre cette méthode, il faut revenir à la première année de l'institut et rappeler les définitions de la série Taylor ( wiki ). En un mot: connaissant la fonction et ses dérivées à un moment donné, vous pouvez trouver les valeurs de la fonction au voisinage de ce point en développant une série de Taylor. Pour la fonction sinus, cela signifie ce qui suit

$$

$$\sin(x_0+\Delta x) \approx \sin(x_0)+\sin'(x_0)\frac{\Delta x}{1!}+\sin''(x_0)\frac{\Delta x^2}{2!}+\sin'''(x_0)\frac{\Delta x^3}{3!}+\cdots$$

Que nous apporte cette approche d'un point de vue pratique? Imaginez que nous ayons un intervalle de$$$0$ avant $$$\pi/2$... Choisissons 10 points linéairement distribués sur cet intervalle (le choix n'est pas optimal):$$$x_0=0$, $$$x_1=\pi/20$, $$$x_2=2\pi/20$, $$$x_i=i\cdot \pi/20$... Pour chaque point, calculez la plaque avec le sinus et ses dérivées à ce point. Vous pouvez maintenant modifier la fonction afin que lors de l'obtention de la valeur$$$x$ la fonction prend la valeur la plus proche $$$x_i$ et dispose dans une rangée autour du point $$$x_i$, pas autour de zéro ($$$\Delta x = x - x_i$).

Utilisation des transformations trigonométriques

Si nous remontons encore plus loin, vers les classes supérieures de l'école, alors nous pouvons rappeler une formule très importante:

$$

$$\sin(x_0+\Delta x)=\sin(x_0)\cos(\Delta x)+\cos(x_0)\sin(\Delta x)$$

Et puis tout est le même que dans le paragraphe précédent. Nous sélectionnons des points dans l'intervalle, calculons le sinus et le cosinus pour eux, et lors de l'appel de la fonction sinus, nous cherchons le plus proche et, en utilisant la formule ci-dessus, calculons le sinus en utilisant une petite valeur$$$\Delta x$...

Pensez à laquelle de ces deux méthodes est préférable de choisir, mais pour l'instant, nous allons passer des mathématiques aux calculs pratiques.

Propriété de distribution de la multiplication dans le monde en virgule flottante

J'ai dû demander conseil à Internet, comment ça s'appelle $$$a\cdot(b+c)=a\cdot b + a\cdot c$... Cela s'avère être une propriété distributive. Revenons à la question que j'ai posée à la fin de la première partie. À savoir, pourquoi des expressions mathématiquement équivalentes$$$a\cdot(b+c)$ et $$$a\cdot b + a\cdot c$peut donner des résultats différents dans les calculs en virgule flottante? La manière la plus simple d'illustrer ceci est avec un exemple. Prenons un système hypothétique qui fonctionne avec des nombres à virgule flottante au format décimal avec 4 chiffres de précision. Faisons semblant que$$$b=1.000\text{E}0$, $$$c=1.234\text{E-}2$, et $$$a=5.678\text{E-}2$... Prenons d'abord une expression entre parenthèses et calculons-la étape par étape, en nous rappelant d'arrondir à chaque étape:

1)$$$b+c=1.000\text{E}0 + 1.234\text{E-}2=1.012\text{E}0$

2) $$$a\cdot(b+c)=5.678\text{E-}2\times 1.012\text{E}0=5.746\text{E-}2$

Réponse reçue $$$y=5.746\text{E-}2$

Calculons maintenant la deuxième expression de la même manière étape par étape:

$$$a\cdot b + a\cdot c=5.678\text{E-}2\times 1.000\text{E}0 + 5.678\text{E-}2\times 1.234\text{E-}2=5.678\text{E-}2 + 7.007\text{E-}4=5.748\text{E-}2$

Réponse reçue $$$y=5.748\text{E-}2$

La vraie réponse est 0,0574806652.



Comme vous pouvez le voir, la réponse obtenue dans le second cas est beaucoup plus proche du vrai que dans le premier. Si nous expliquons cela sur les doigts, imaginez que lorsque dans le premier cas, nous ajoutons à 1,0 le nombre$$$c=1.234\text{E-}2$nous jetons simplement les deux derniers chiffres. Ils ne sont plus. Dans le second cas, le rejet intervient à la toute fin, après la multiplication. Ceux. dans le second cas, la ou les opérations de multiplication sont plus précises.



Il semble que vous puissiez terminer là-dessus, mais regardez de plus près la première méthode et dites-moi quel sera le résultat du calcul$$$b+c-b$... Et ... nous avons un moyen d'arrondir les nombres à virgule flottante! Ne manquez pas cet exemple. Donnez-vous le temps de le comprendre. L'arrondi des nombres sera utilisé de manière très intensive par nous plus tard dans cet article et les suivants.



Remarquons une autre caractéristique de cette expression. Imaginez que la précision à 4 chiffres de la variable ne nous soit pas suffisante. Que faire? Et ici, nous avons déjà la réponse - pour représenter le nombre sous la forme$$$b+c$et enregistrez-le en mémoire sous forme de somme de deux chiffres. Et, en conséquence, effectuez les opérations (par exemple la multiplication) séparément pour les deux termes. Cette technique est décrite plus en détail dans l'article Ajout de deux nombres à virgule flottante sans perte de précision .



Dans l'article précédent, j'ai également écrit que la méthode$$$a\cdot(b+c)$il y a une caractéristique désagréable. Et c'est comme suit. Nombre$$$c$ toujours tronqué au dernier chiffre significatif d'un nombre $$$b$... Cela signifie que quel que soit le nombre$$$c$, si un $$$b+c \neq b$, alors une erreur dans le dernier signe est toujours possible même pour les petits $$$c$... Cela n'est pas autorisé dans l'approche du chapitre suivant.

Comment cela fonctionne en utilisant la bibliothèque GNU comme exemple

Comment c'est? Avez-vous choisi laquelle des deux méthodes décrites au début de l'article avez-vous choisie pour le calcul précis du sinus? Quelle que soit la méthode choisie, les deux sont correctes. De plus, ils sont absolument identiques. Croyez-moi, vérifiez-le. Ci-dessous, j'utiliserai des formules scolaires. Ils sont plus faciles à expliquer.



Armé des connaissances acquises dans l'article précédent et dans cet article, vous pouvez facilement comprendre le code de la bibliothèque standard. Ouvrons le fichier s_sin.c et trouvons la fonction __sin là - bas :



son code est assez simple. Il est facile de comprendre qu'il appelle un ensemble de fonctions différent en fonction des limites de la variable d'entrée. Dans cet article, nous aborderons la section de code 218-224 pour les angles 2 ^ -26 <| x | <0.855469. Vous pouvez voir que dans cette section du code la fonction do_sin (x, 0) est appelée. Nous allons nous attarder sur cette fonction plus en détail:

1. , dx=0 .
2. 129-130 , abs(x)<0.126, .. x , . , , , .
3. 136-137. , . x 2 . u x. , 0.345678. u=0.34, 0.005678.
4. 140-142. ( s ) ( c ) x . , cos(x)=1-c, 1.0, (. ), .
5. 143. u. , u=0.34 34. sin(u)=sn+ssn, cos(u)=cs+ccs. sn cs — «» u, ssn ccs — .
6. 144-145. sin(u+x)=(sn+ssn)*(1-c)+(cs+ccs)*s. , , 144-145. — .

En fait, je n'ai décrit que la partie la plus simple du calcul du sinus de cette manière. Il reste encore beaucoup de mathématiques. Par exemple, comment calculez-vous la taille d'un tableau et les éléments qu'il contient? D'où viennent les nombres magiques 0,126 et 0,855469? Quand couper le calcul par le nombre de Taylor? Corrections des coefficients de la série de Taylor pour affiner le résultat.



Tout cela, bien sûr, est intéressant, mais, objectivement, la méthode présentée présente de nombreux inconvénients: il est nécessaire de calculer simultanément le sinus (s) et le cosinus (c), ce qui nécessite deux fois plus de calculs de la série de Taylor ¹ . La multiplication par des valeurs tabulaires, comme nous pouvons le voir, n'est pas non plus gratuite. De plus, stocker une table de 3520 octets en RAM n'est bien sûr pas un problème, mais y accéder (même dans le cache) peut être coûteux.



Par conséquent, dans la partie suivante, nous allons essayer de nous débarrasser de la plaque et de calculer le sinus dans l'intervalle [0,126, 0,855469] directement, mais plus précisément que dans le premier chapitre.



Avant de terminer - une question d'esprit vif. Le grand nombre dans cet exemple est 52776558133248 = 3 * 2 ⁴⁴ . D'où vient un tel nombre, pas, par exemple, 2 ⁴⁵ ? Je formulerai la question plus précisément. Pourquoi le nombre 3 * 2 ^{N est-il} optimal lors de l'arrondissement des nombres , et non, par exemple, 2 ^{N + 1} ? Une autre question, quel N devriez-vous choisir pour arrondir un nombre à un entier?



¹ Il est à noter qu'un avantage significatif de cette approche peut apparaître lorsque le sinus et le cosinus sont calculés simultanément à partir du même angle. La deuxième fonction peut être calculée presque gratuitement.