sine cosinus arithmétique de précision étendue modulaire

Question

J'ai vu dans de nombreuses interprétations de sinus/cosinus ce que l'on appelle une arithmétique de précision modulaire étendue. Mais à quoi cela sert-il? Par exemple dans le cephes implemetation, après réduction à la plage [0, pi/4], ils font cette arithmétique de précision modulaire pour améliorer la précision.

le code Ci-dessous:

z = ((x - y * DP1) - y * DP2) - y * DP3; 

où DP1, DP2 et DP3 sont un certain coefficient hardcoded. Comment trouver ces coefficients mathématiquement? J'ai compris le but de "l'arithmétique d'extension modulaire" pour grand num, mais ici quel est son but exact?

Answer 1

Dans le contexte de la réduction des arguments pour les fonctions trigonométriques, ce que vous cherchez à est la réduction de l'argument Cody-Waite, une technique introduite dans le livre: William J. Cody et William Waite, Manuel du logiciel pour les fonctions élémentaires, Prentice-Hall, 1980. L'objectif est d'obtenir, pour les arguments jusqu'à une certaine ampleur, un argument réduit précis, malgré subtractive cancellation dans le calcul intermédiaire. A cette fin, la constante pertinente est représentée avec plus que la précision native, en utilisant une somme de plusieurs nombres de grandeur décroissante (ici: DP1, DP2, DP3), de sorte que tous les produits intermédiaires sauf le moins significatif peuvent être calculé sans erreur d'arrondi. Considérons comme exemple le calcul de sin (113) dans IEEE-754 binary32 (simple précision). La réduction d'argument typique calculerait conceptuellement i=rintf(x/(π/2)); reduced_x = x-i*(π/2). Le numéro binary32 le plus proche de π/2 est 0x1.921fb6p+0. Nous calculons i=72, le produit arrondit à 0x1.c463acp+6, ce qui est proche de l'argument x=0x1.c40000p+6. Pendant la soustraction, certains bits principaux s'annulent, et nous finissons par reduced_x = -0x1.8eb000p-4. Notez les zéros de fin introduits par la renormalisation. Ces bits zéro ne portent aucune information utile. En appliquant une approximation précise à l'argument réduit, sin(x) = -0x1.8e0eeap-4, alors que le vrai résultat est -0x1.8e0e9d39...p-4. Nous nous retrouvons avec une erreur relative importante et une grande erreur ulp.

Nous pouvons y remédier en utilisant une réduction d'argument Cody-Waite en deux étapes. Par exemple, nous pourrions utiliser pio2_hi = 0x1.921f00p+0 et pio2_lo = 0x1.6a8886p-17. Notez les huit bits de fin arrière dans la représentation à simple précision de pio2_hi, ce qui nous permet de multiplier avec tout nombre entier de 8 bits i et toujours le produit i * pio2_hi représentable exactement comme un nombre à simple précision. Lorsque nous calculons ((x - i * pio2_hi) - i * pio2_lo), nous obtenons reduced_x = -0x1.8eafb4p-4, et donc sin(x) = -0x1.8e0e9ep-4, un résultat assez précis. La meilleure façon de diviser la constante en une somme dépendra de l'amplitude de i que nous devons gérer, sur le nombre maximum de bits soumis à une annulation soustractive pour une plage d'arguments donnée (sur la base des multiples entiers proches de π/2 peut obtenir des entiers), et des considérations de performance. Les cas typiques d'utilisation dans la vie réelle impliquent des schémas de réduction Cody-Waite de deux à quatre étapes. La disponibilité de FMA (fusion multiple-add) permet d'utiliser des constantes constitutives avec moins de bits de fin de queue. Voir ce papier: Sylvie Boldo, Marc Daumas, et Ren-Cang Li, "Réduction d'argument formellement vérifiée avec un multiplier fusionné." Transactions IEEE sur les ordinateurs, 58: 1139-1145, 2009. Pour un exemple travaillé utilisant fmaf(), vous souhaiterez peut-être regarder le code dans one of my previous answers.