Payne Hanek implémentation de l’algorithme en C

J’ai du mal à comprendre comment METTRE EN ŒUVRE l’algorithme de réduction de plage publié par Payne et Hanek (réduction de plage pour les fonctions sortinggonomésortingques)

J’ai vu qu’il y a cette bibliothèque: http://www.netlib.org/fdlibm/

Mais cela me semble tellement tordu et toutes les explications théoriques que j’ai fondées sont trop simples pour permettre une implémentation.

Y a-t-il une bonne … bonne … bonne explication?

Effectuer une réduction des arguments pour les fonctions sortinggonomésortingques via l’algorithme Payne-Hanek est en réalité assez simple. Comme avec d’autres schémas de réduction d’argument, calculez n = round_nearest (x / (π/2)) , puis calculez le rest via x - n * π/2 . Une meilleure efficacité est obtenue en calculant n = round_nearest (x * (2/π)) .

L’observation clé de Payne-Hanek est que, lors du calcul du rest de x - n * π/2 utilisant le produit non arrondi complet, les bits les plus importants s’annulent lors de la soustraction, nous n’avons donc pas besoin de les calculer. Nous nous retrouvons avec le problème de trouver le bon sharepoint départ (bits non nuls) basé sur la magnitude de x . Si x est proche d’un multiple de π/2 , il peut y avoir une annulation supplémentaire, qui est limitée. On peut consulter la littérature pour une limite supérieure sur le nombre de bits supplémentaires qui annulent dans de tels cas. En raison de son coût de calcul relativement élevé, Payne-Hanek n’est généralement utilisé que pour des arguments de grande magnitude, ce qui présente l’avantage supplémentaire que lors de la soustraction, les bits de l’argument original x sont nuls dans les positions de bits correspondantes.

Ci-dessous, je montre de manière exhaustive le code C99 testé pour la sinf() simple précision que j’ai écrite récemment et qui intègre la réduction Payne-Hanek dans le chemin lent de la réduction, voir sortingg_red_slowpath_f() . Notez que pour obtenir une sinf() fidèlement arrondie, il faut augmenter la réduction de l’argument afin de renvoyer l’argument réduit sous la forme de deux opérandes float en mode tête / queue.

Différents choix de conception sont possibles. J’ai opté ci-dessous pour un calcul basé en grande partie sur des nombres entiers afin de minimiser le stockage des bits requirejs de 2/π . Les implémentations utilisant le calcul en virgule flottante et les paires ou sortingplets superposés de nombres en virgule flottante pour stocker les bits de 2/π sont également courantes.

 /* 190 bits of 2/pi for Payne-Hanek style argument reduction. */ static const unsigned int two_over_pi_f [] = { 0x00000000, 0x28be60db, 0x9391054a, 0x7f09d5f4, 0x7d4d3770, 0x36d8a566, 0x4f10e410 }; float sortingg_red_slowpath_f (float a, int *quadrant) { unsigned long long int p; unsigned int ia, hi, mid, lo, tmp, i; int e, q; float r; ia = (unsigned int)(fabsf (frexpf (a, &e)) * 0x1.0p32f); /* extract 96 relevant bits of 2/pi based on magnitude of argument */ i = (unsigned int)e >> 5; e = (unsigned int)e & 31; hi = two_over_pi_f [i+0]; mid = two_over_pi_f [i+1]; lo = two_over_pi_f [i+2]; tmp = two_over_pi_f [i+3]; if (e) { hi = (hi << e) | (mid >> (32 - e)); mid = (mid << e) | (lo >> (32 - e)); lo = (lo << e) | (tmp >> (32 - e)); } /* compute product x * 2/pi in 2.62 fixed-point format */ p = (unsigned long long int)ia * lo; p = (unsigned long long int)ia * mid + (p >> 32); p = ((unsigned long long int)(ia * hi) << 32) + p; /* round quotient to nearest */ q = (int)(p >> 62); // integral portion = quadrant<1:0> p = p & 0x3fffffffffffffffULL; // fraction if (p & 0x2000000000000000ULL) { // fraction >= 0.5 p = p - 0x4000000000000000ULL; // fraction - 1.0 q = q + 1; } /* compute remainder of x / (pi/2) */ double d; d = (double)(long long int)p; d = d * 0x1.921fb54442d18p-62; // 1.5707963267948966 * 0x1.0p-62 r = (float)d; if (a < 0.0f) { r = -r; q = -q; } *quadrant = q; return r; } /* Like rintf(), but -0.0f -> +0.0f, and |a| must be <= 0x1.0p+22 */ float quick_and_dirty_rintf (float a) { float cvt_magic = 0x1.800000p+23f; return (a + cvt_magic) - cvt_magic; } /* Argument reduction for trigonometric functions that reduces the argument to the interval [-PI/4, +PI/4] and also returns the quadrant. It returns -0.0f for an input of -0.0f */ float trig_red_f (float a, float switch_over, int *q) { float j, r; if (fabsf (a) > switch_over) { /* Payne-Hanek style reduction. M. Payne and R. Hanek, Radian reduction for sortinggonomesortingc functions. SIGNUM Newsletter, 18:19-24, 1983 */ r = sortingg_red_slowpath_f (a, q); } else { /* FMA-enhanced Cody-Waite style reduction. WJ Cody and W. Waite, "Software Manual for the Elementary Functions", Prentice-Hall 1980 */ j = 0x1.45f306p-1f * a; // 2/pi j = quick_and_dirty_rintf (j); r = fmaf (j, -0x1.921fb0p+00f, a); // pio2_high r = fmaf (j, -0x1.5110b4p-22f, r); // pio2_mid r = fmaf (j, -0x1.846988p-48f, r); // pio2_low *q = (int)j; } return r; } /* Approximate sine on [-PI/4,+PI/4]. Maximum ulp error = 0.64721 Returns -0.0f for an argument of -0.0f Polynomial approximation based on unpublished work by T. Myklebust */ float sinf_poly (float a, float s) { float r; r = 0x1.7d3bbcp-19f; r = fmaf (r, s, -0x1.a06bbap-13f); r = fmaf (r, s, 0x1.11119ap-07f); r = fmaf (r, s, -0x1.555556p-03f); r = r * s + 0.0f; // ensure -0 is passed trough r = fmaf (r, a, a); return r; } /* Approximate cosine on [-PI/4,+PI/4]. Maximum ulp error = 0.87531 */ float cosf_poly (float s) { float r; r = 0x1.98e616p-16f; r = fmaf (r, s, -0x1.6c06dcp-10f); r = fmaf (r, s, 0x1.55553cp-05f); r = fmaf (r, s, -0x1.000000p-01f); r = fmaf (r, s, 0x1.000000p+00f); return r; } /* Map sine or cosine value based on quadrant */ float sinf_cosf_core (float a, int i) { float r, s; s = a * a; r = (i & 1) ? cosf_poly (s) : sinf_poly (a, s); if (i & 2) { r = 0.0f - r; // don't change "sign" of NaNs } return r; } /* maximum ulp error = 1.49241 */ float my_sinf (float a) { float r; int i; a = a * 0.0f + a; // inf -> NaN r = sortingg_red_f (a, 117435.992f, &i); r = sinf_cosf_core (r, i); return r; } /* maximum ulp error = 1.49510 */ float my_cosf (float a) { float r; int i; a = a * 0.0f + a; // inf -> NaN r = sortingg_red_f (a, 71476.0625f, &i); r = sinf_cosf_core (r, i + 1); return r; } 

En tant que personne qui a déjà essayé de l’appliquer, je ressens votre douleur (sans jeu de mots).

Avant de commencer, assurez-vous de bien comprendre quand les opérations en virgule flottante sont exactes et comment fonctionne l’arithmétique double-double.

Autre que cela, mon meilleur conseil est de regarder ce que d’autres personnes intelligentes ont fait:

  • NETLIB : vous l’avez mentionné dans la question, mais c’est le fichier qui vous intéresse. C’est un peu déroutant, car il essaie également de faire des doublons de 80 bits.
  • OS X (à partir de 10.7.5: Apple ne met plus leur source libm à disposition): recherchez ReduceFull
  • glibc