Meilleure approximation minimax polynomiale optimisée par la machine pour arctangente sur ?

Pour la mise en œuvre simple et efficace de fonctions mathématiques rapides avec une précision raisonnable, les approximations polynomiales à minimax constituent souvent la méthode de choix. Les approximations minima sont généralement générées avec une variante de l’algorithme de Remez. Divers outils largement disponibles, tels que Maple et Mathematica, disposent d’une fonctionnalité intégrée à cet effet. Les coefficients générés sont généralement calculés en utilisant une arithmétique de haute précision. Il est bien connu que simplement arrondir ces coefficients à la précision de la machine conduit à une précision sous-optimale dans la mise en oeuvre résultante.

Au lieu de cela, on recherche des ensembles de coefficients étroitement liés qui sont exactement représentables sous forme de numéros de machine afin de générer une approximation optimisée pour la machine. Deux documents pertinents sont:

Nicolas Brisebarre, Jean-Michel Muller et Arnaud Tisserand, “Calcul des approximations polynomiales efficaces en machine”, Transactions ACM sur les logiciels mathématiques, vol. 32, n ° 2, juin 2006, p. 236–256.

Nicolas Brisebarre et Sylvain Chevillard, “Approches de L polynomiales efficaces”, 18ème symposium de l’IEEE sur l’arithmétique de l’ordinateur (ARITH-18), Montpellier (France), juin 2007, p. 169-176.

Une implémentation de l’algorithme LLL à partir de ce dernier article est disponible en tant que commande fpminimax() de l’ outil Sollya . Je crois comprendre que tous les algorithmes proposés pour la génération d’approximations optimisées par la machine sont basés sur des méthodes heuristiques et que, par conséquent, on ignore généralement quelle précision peut être obtenue par une approximation optimale. Il n’est pas clair pour moi si la disponibilité de FMA (multiplie-addition fusionné) pour l’évaluation de l’approximation a une influence sur la réponse à cette question. Il me semble naïvement que ce devrait être le cas.

Je suis actuellement en train de regarder une approximation polynomiale simple pour arctangente sur [-1,1] qui est évaluée en arithmétique à simple précision IEEE-754, en utilisant le schéma de Horner et FMA. Voir la fonction atan_poly() dans le code C99 ci-dessous. Pour le moment, faute d’access à une machine Linux, je n’ai pas utilisé Sollya pour générer ces coefficients, mais ma propre heuristique, que l’on pourrait décrire vaguement comme un mélange de recuits les plus raides et simulés (pour éviter de restr bloqué sur des minima locaux). . L’erreur maximale de mon polynôme optimisé pour la machine est très proche de 1 ulp, mais idéalement, j’aimerais que l’erreur ulp maximale soit inférieure à 1 ulp.

Je suis conscient que je pourrais modifier mon calcul pour augmenter la précision, en utilisant par exemple un coefficient d’avance représenté avec une précision supérieure à la précision avec une précision, mais j’aimerais conserver le code tel quel (c’est-à-dire aussi simple que possible). ajuster uniquement les coefficients pour fournir le résultat le plus précis possible.

Un ensemble optimal “prouvé” de coefficients serait idéal, les pointeurs sur la littérature pertinente sont les bienvenus. J’ai fait une recherche documentaire, mais je n’ai trouvé aucun article qui fpminimax() de fpminimax() significative l’état de la technique au-delà du fpminimax() de Sollya, et aucun qui examine le rôle des FMA (le cas échéant) dans ce numéro.

 // max ulp err = 1.03143 float atan_poly (float a) { float r, s; s = a * a; r = 0x1.7ed1ccp-9f; r = fmaf (r, s, -0x1.0c2c08p-6f); r = fmaf (r, s, 0x1.61fdd0p-5f); r = fmaf (r, s, -0x1.3556b2p-4f); r = fmaf (r, s, 0x1.b4e128p-4f); r = fmaf (r, s, -0x1.230ad2p-3f); r = fmaf (r, s, 0x1.9978ecp-3f); r = fmaf (r, s, -0x1.5554dcp-2f); r = r * s; r = fmaf (r, a, a); return r; } // max ulp err = 1.52637 float my_atanf (float a) { float r, t; t = fabsf (a); r = t; if (t > 1.0f) { r = 1.0f / r; } r = atan_poly (r); if (t > 1.0f) { r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, -r); // pi/2 - r } r = copysignf (r, a); return r; } 

La fonction suivante est une implémentation fidèlement arrondie d’ arctan sur [0, 1] :

 float atan_poly (float a) { float s = a * a, u = fmaf(a, -a, 0x1.fde90cp-1f); float r1 = 0x1.74dfb6p-9f; float r2 = fmaf (r1, u, 0x1.3a1c7cp-8f); float r3 = fmaf (r2, s, -0x1.7f24b6p-7f); float r4 = fmaf (r3, u, -0x1.eb3900p-7f); float r5 = fmaf (r4, s, 0x1.1ab95ap-5f); float r6 = fmaf (r5, u, 0x1.80e87cp-5f); float r7 = fmaf (r6, s, -0x1.e71aa4p-4f); float r8 = fmaf (r7, u, -0x1.b81b44p-3f); float r9 = r8 * s; float r10 = fmaf (r9, a, a); return r10; } 

Le harnais de test suivant sera abandonné si la fonction atan_poly ne parvient pas à être arrondie fidèlement sur [1e-16, 1] et affiche “succès” sinon:

 int checkit(float f) { double d = atan(f); float d1 = d, d2 = d; if (d1 < d) d2 = nextafterf(d1, 1.0/0.0); else d1 = nextafterf(d1, -1.0/0.0); float p = atan_poly(f); if (p != d1 && p != d2) return 0; return 1; } int main() { for (float f = 1; f > 1e-16; f = nextafterf(f, -1.0/0.0)) { if (!checkit(f)) abort(); } printf("success\n"); exit(0); } 

Le problème de l’utilisation de s dans chaque multiplication est que les coefficients du polynôme ne se dégradent pas rapidement. Des entrées proches de 1 entraînent l’annulation de nombres presque égaux, ce qui signifie que vous essayez de trouver un ensemble de coefficients de sorte que l’arrondi accumulé à la fin du calcul se rapproche beaucoup du résidu d’ arctan .

La constante 0x1.fde90cp-1f est un nombre proche de 1 pour lequel (arctan(sqrt(x)) - x) / x^3 est très proche du flottant le plus proche. En d’autres termes, c’est une constante qui entre dans le calcul de u afin que le coefficient cubique soit presque complètement déterminé. (Pour ce programme, le coefficient cubique doit être soit -0x1.b81b44p-3f ou -0x1.b81b42p-3f .)

Alterner les multiplications par s et u a pour effet de réduire l’effet de l’erreur d’arrondi de ri sur r{i+2} d’un facteur d’au plus 1/4, car s*u < 1/4 quel a soit a. Cela donne une marge de manœuvre considérable dans le choix des coefficients du cinquième ordre et au-delà.


J'ai trouvé les coefficients à l'aide de deux programmes:

  • Un programme connecte plusieurs points de test, écrit un système d’inégalités linéaires et calcule les bornes des coefficients à partir de ce système d’inégalités. Notez que, étant donné a , on peut calculer la plage de r8 qui conduit à un résultat fidèlement arrondi. Pour obtenir des inégalités linéaires , j'ai prétendu que r8 serait calculé comme un polynôme dans les float s s et u dans l' arithmétique en nombres réels ; les inégalités linéaires ont contraint ce nombre réel r8 à se situer dans un intervalle. J'ai utilisé la bibliothèque de polyèdres de Parme pour gérer ces systèmes de contraintes.
  • Un autre programme a testé de manière aléatoire des ensembles de coefficients dans certaines gammes, en insérant d’abord un ensemble de points de test, puis tous les float de 1 à 1e-8 en ordre décroissant, et en vérifiant atan_poly produit un arrondi fidèle de atan((double)x) . Si certains x échoué, ils l'ont imprimé et pourquoi.

Pour obtenir des coefficients, j'ai piraté ce premier programme afin de corriger c3 , de calculer les limites de r7 pour chaque sharepoint test, puis d'obtenir les limites de coefficients d'ordre supérieur. Ensuite, je l'ai piraté pour corriger c3 et c5 et obtenir des limites sur les coefficients d'ordre supérieur. Je l'ai fait jusqu'à ce que tous les coefficients de l'ordre le plus élevé, c13 , c15 et c17 , c13 c15 c17 .

J'ai développé l'ensemble des points de test dans le deuxième programme jusqu'à ce qu'il cesse d'imprimer quoi que ce soit ou qu'il imprime "succès". Il me fallait étonnamment peu de points de test pour rejeter presque tous les polynômes erronés - je compte 85 points de test dans le programme.


Ici, je montre une partie de mon travail en sélectionnant les coefficients. Afin d’obtenir un arctan fidèlement arrondi pour mon ensemble initial de points de test en supposant que r1 à r8 soient évalués en arithmétique réelle (et arrondis de manière désagréable mais d’une manière dont je ne me souviens plus), mais r9 et r10 sont évalués en arithmétique float , J'ai besoin:

 -0x1.b81b456625f15p-3 <= c3 <= -0x1.b81b416e22329p-3 -0x1.e71d48d9c2ca4p-4 <= c5 <= -0x1.e71783472f5d1p-4 0x1.80e063cb210f9p-5 <= c7 <= 0x1.80ed6efa0a369p-5 0x1.1a3925ea0c5a9p-5 <= c9 <= 0x1.1b3783f148ed8p-5 -0x1.ec6032f293143p-7 <= c11 <= -0x1.e928025d508p-7 -0x1.8c06e851e2255p-7 <= c13 <= -0x1.732b2d4677028p-7 0x1.2aff33d629371p-8 <= c15 <= 0x1.41e9bc01ae472p-8 0x1.1e22f3192fd1dp-9 <= c17 <= 0x1.d851520a087c2p-9 

En prenant c3 = -0x1.b81b44p-3, en supposant que r8 est également évalué en arithmétique float :

 -0x1.e71df05b5ad56p-4 <= c5 <= -0x1.e7175823ce2a4p-4 0x1.80df529dd8b18p-5 <= c7 <= 0x1.80f00e8da7f58p-5 0x1.1a283503e1a97p-5 <= c9 <= 0x1.1b5ca5beeeefep-5 -0x1.ed2c7cd87f889p-7 <= c11 <= -0x1.e8c17789776cdp-7 -0x1.90759e6defc62p-7 <= c13 <= -0x1.7045e66924732p-7 0x1.27eb51edf324p-8 <= c15 <= 0x1.47cda0bb1f365p-8 0x1.f6c6b51c50b54p-10 <= c17 <= 0x1.003a00ace9a79p-8 

En prenant c5 = -0x1.e71aa4p-4, en supposant que r7 est fait en arithmétique float :

 0x1.80e3dcc972cb3p-5 <= c7 <= 0x1.80ed1cf56977fp-5 0x1.1aa005ff6a6f4p-5 <= c9 <= 0x1.1afce9904742p-5 -0x1.ec7cf2464a893p-7 <= c11 <= -0x1.e9d6f7039db61p-7 -0x1.8a2304daefa26p-7 <= c13 <= -0x1.7a2456ddec8b2p-7 0x1.2e7b48f595544p-8 <= c15 <= 0x1.44437896b7049p-8 0x1.396f76c06de2ep-9 <= c17 <= 0x1.e3bedf4ed606dp-9 

En prenant c7 = 0x1.80e87cp-5, en supposant que r6 est fait en arithmétique float :

 0x1.1aa86d25bb64fp-5 <= c9 <= 0x1.1aca48cd5caabp-5 -0x1.eb6311f6c29dcp-7 <= c11 <= -0x1.eaedb032dfc0cp-7 -0x1.81438f115cbbp-7 <= c13 <= -0x1.7c9a106629f06p-7 0x1.36d433f81a012p-8 <= c15 <= 0x1.3babb57bb55bap-8 0x1.5cb14e1d4247dp-9 <= c17 <= 0x1.84f1151303aedp-9 

En prenant c9 = 0x1.1ab95ap-5, en supposant que r5 est fait en arithmétique float :

 -0x1.eb51a3b03781dp-7 <= c11 <= -0x1.eb21431536e0dp-7 -0x1.7fcd84700f7cfp-7 <= c13 <= -0x1.7ee38ee4beb65p-7 0x1.390fa00abaaabp-8 <= c15 <= 0x1.3b100a7f5d3cep-8 0x1.6ff147e1fdeb4p-9 <= c17 <= 0x1.7ebfed3ab5f9bp-9 

J'ai choisi un point proche du milieu de la plage pour c11 et choisi au hasard c13 , c15 et c17 .


EDIT: J'ai maintenant automatisé cette procédure. La fonction suivante est également une implémentation fidèlement arrondie d’ arctan sur [0, 1] :

 float c5 = 0x1.997a72p-3; float c7 = -0x1.23176cp-3; float c9 = 0x1.b523c8p-4; float c11 = -0x1.358ff8p-4; float c13 = 0x1.61c5c2p-5; float c15 = -0x1.0b16e2p-6; float c17 = 0x1.7b422p-9; float juffa_poly (float a) { float s = a * a; float r1 = c17; float r2 = fmaf (r1, s, c15); float r3 = fmaf (r2, s, c13); float r4 = fmaf (r3, s, c11); float r5 = fmaf (r4, s, c9); float r6 = fmaf (r5, s, c7); float r7 = fmaf (r6, s, c5); float r8 = fmaf (r7, s, -0x1.5554dap-2f); float r9 = r8 * s; float r10 = fmaf (r9, a, a); return r10; } 

Je trouve surprenant que ce code existe même. Pour les coefficients proches de ceux-ci, vous pouvez obtenir une borne sur la distance entre r10 et la valeur du polynôme évaluée en arithmétique réelle de l'ordre de quelques ulps grâce à la convergence lente de ce polynôme lorsque s est proche de 1 . Je m'attendais à ce que l’erreur d’arrondi se comporte de manière fondamentalement "indomptable" simplement en ajustant les coefficients.

Ce n’est pas une réponse à la question, mais c’est trop long pour tenir dans un commentaire:

Votre question concerne le choix optimal des coefficients C3, C5,…, C17 dans une approximation polynomiale à arctangente où vous avez épinglé C1 à 1 et C2, C4,…, C16 à 0.

Le titre de votre question indique que vous recherchez des approximations sur [-1, 1], et une bonne raison d’épingler les coefficients pairs à 0 est qu’il est suffisant et nécessaire que l’approximation soit exactement une fonction impaire. Le code de votre question “contredit” le titre en appliquant l’approximation polynomiale uniquement sur [0, 1].

Si vous utilisez l’algorithme de Remez pour rechercher les coefficients C2, C3,…, C8 à une approximation polynomiale d’arctangente sur [0, 1] à la place, vous pouvez obtenir quelque chose comme les valeurs ci-dessous:

 #include  #include  float atan_poly (float a) { float r, s; s = a; // s = a * a; r = -3.3507930064626076153585890630056286726807491543578e-2; r = fmaf (r, s, 1.3859776280052980081098065189344699108643282883702e-1); r = fmaf (r, s, -1.8186361916440430105127602496688553126414578766147e-1); r = fmaf (r, s, -1.4583047494913656326643327729704639191810926020847e-2); r = fmaf (r, s, 2.1335202878219865228365738728594741358740655881373e-1); r = fmaf (r, s, -3.6801711826027841250774413728610805847547996647342e-3); r = fmaf (r, s, -3.3289852243978319173749528028057608377028846413080e-1); r = fmaf (r, s, -1.8631479933914856903459844359251948006605218562283e-5); r = fmaf (r, s, 1.2917291732886065585264586294461539492689296797761e-7); r = fmaf (r, a, a); return r; } int main() { for (float x = 0.0f; x < 1.0f; x+=0.1f) printf("x: %f\n%a\n%a\n\n", x, atan_poly(x), atan(x)); } 

Cela a à peu près la même complexité que le code de votre question: le nombre de multiplications est similaire. En regardant ce polynôme, il n'y a aucune raison particulière de vouloir épingler un coefficient à 0. Si on voulait approximer une fonction impaire sur [-1, 1] sans épingler les coefficients pairs, ils seraient automatiquement très petits et sujets à l'absorption, et alors nous voudrions les épingler à 0, mais pour cette approximation sur [0, 1], ils ne le font pas, donc nous n'avons pas à les épingler.

Cela aurait pu être meilleur ou pire que le polynôme étrange de votre question. Il s'avère que c'est pire (voir ci-dessous). Cette application rapide de LolRemez 0.2 (code inclus au bas de la question) semble toutefois être suffisante pour poser la question du choix des coefficients. Je serais en particulier curieux de savoir ce qui se passera si vous soumettez les coefficients de cette réponse à la même étape d'optimisation du «mélange de recuits les plus prononcés et simulés» que vous avez appliquée pour obtenir les coefficients de votre question.

Donc, pour résumer cette remarque postée en tant que réponse, êtes-vous sûr de rechercher les coefficients optimaux C3, C5,…, C17? Il me semble que vous recherchez la meilleure séquence d'opérations à virgule flottante simple précision qui produisent une approximation fidèle de l'arctangente et que cette approximation ne doit pas nécessairement être la forme de Horner d'un polynôme de degré 17 impair.

 x: 0,000000
 0x0p + 0
 0x0p + 0

 x: 0,100000
 0x1.983e2cp-4
 0x1.983e28938f9ecp-4

 x: 0,200000
 0x1.94442p-3
 0x1.94441ff1e8882p-3

 x: 0,300000
 0x1.2a73a6p-2
 0x1.2a73a71dcec16p-2

 x: 0,400000
 0x1.85a37ap-2
 0x1.85a3770ebe7aep-2

 x: 0,500000
 0x1.dac67p-2
 0x1.dac670561bb5p-2

 x: 0,600000
 0x1.14b1dcp-1
 0x1.14b1ddf627649p-1

 x: 0,700000
 0x1.38b116p-1
 0x1.38b113eaa384ep-1

 x: 0,800000
 0x1.5977a8p-1
 0x1.5977a686e0ffbp-1

 x: 0,900000
 0x1.773388p-1
 0x1.77338c44f8faep-1

C'est le code que j'ai lié à LolRemez 0.2 afin d'optimiser la précision relative d'une approximation polynomiale de degré 9 d'arctangente sur [0, 1]:

 #include "lol/math/real.h" #include "lol/math/remez.h" using lol::real; using lol::RemezSolver; real f(real const &y) { return (atan(y) - y) / y; } real g(real const &y) { return re (atan(y) / y); } int main(int argc, char **argv) { RemezSolver<8, real> solver; solver.Run("1e-1000", 1.0, f, g, 50); return 0; } 

J’ai réfléchi aux différentes idées que j’ai reçues dans les commentaires et ai également mené quelques expériences basées sur ces commentaires. En fin de compte, j’ai décidé qu’une recherche heuristique raffinée constituait la meilleure voie à suivre. J’ai maintenant réussi à réduire l’erreur maximale pour atanf_poly() à 1.01036 ulps, avec seulement trois arguments dépassant l’objective déclaré d’une erreur de 1 ulp:

 ulp = -1.00829 @ |a| = 9.80738342e-001 0x1.f62356p-1 (3f7b11ab) ulp = -1.01036 @ |a| = 9.87551928e-001 0x1.f9a068p-1 (3f7cd034) ulp = 1.00050 @ |a| = 9.99375939e-001 0x1.ffae34p-1 (3f7fd71a) 

En se basant sur la manière de générer l’approximation améliorée, rien ne garantit qu’il s’agisse d’une approximation optimale; pas de percée scientifique ici. Comme l’erreur ulp de la solution actuelle n’est pas encore parfaitement équilibrée et que la recherche continue de fournir de meilleures approximations (même si les intervalles de temps augmentent de façon exponentielle), j’imagine qu’une borne d’erreur de 1 ulp est réalisable, mais en même temps, semble être très proche de la meilleure approximation optimisée par la machine déjà.

La meilleure qualité de la nouvelle approximation est le résultat d’un processus de recherche raffiné. J’ai observé que toutes les plus grandes erreurs ulp dans le polynôme se produisent près de l’unité, disons dans [0.75,1.0] pour être conservateur. Cela permet de rechercher rapidement des ensembles de coefficients intéressants dont l’erreur maximale est inférieure à une valeur liée, disons 1,08 ulps. Je peux ensuite tester en détail et de manière exhaustive tous les ensembles de coefficients dans un hyper-cône choisi de manière heuristique et ancré à ce point. Cette deuxième étape recherche l’erreur primaire ulp comme objective principal et le pourcentage maximal de résultats correctement arrondis comme objective secondaire. En utilisant ce processus en deux étapes sur les quatre cœurs de mon processeur, j’ai pu accélérer considérablement le processus de recherche: j’ai pu vérifier environ 2 21 ensembles de coefficients.

Sur la base de la plage de chaque coefficient parmi toutes les solutions “proches”, j’estime maintenant que l’espace de recherche utile total pour ce problème d’approximation est> = 2 24 ensembles de coefficients plutôt que le nombre plus optimiste de 2 20 j jeté auparavant. Cela semble être un problème possible à résoudre pour quelqu’un qui est très patient ou qui dispose de beaucoup de puissance de calcul.

Mon code mis à jour est le suivant:

 // max ulp err = 1.01036 float atanf_poly (float a) { float r, s; s = a * a; r = 0x1.7ed22cp-9f; r = fmaf (r, s, -0x1.0c2c2ep-6f); r = fmaf (r, s, 0x1.61fdf6p-5f); r = fmaf (r, s, -0x1.3556b4p-4f); r = fmaf (r, s, 0x1.b4e12ep-4f); r = fmaf (r, s, -0x1.230ae0p-3f); r = fmaf (r, s, 0x1.9978eep-3f); r = fmaf (r, s, -0x1.5554dap-2f); r = r * s; r = fmaf (r, a, a); return r; } // max ulp err = 1.51871 float my_atanf (float a) { float r, t; t = fabsf (a); r = t; if (t > 1.0f) { r = 1.0f / r; } r = atanf_poly (r); if (t > 1.0f) { r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, -r); // pi/2 - r } r = copysignf (r, a); return r; } 

Mise à jour (après avoir réexaminé la question deux ans et demi plus tard)

En utilisant le projet de publication de T. Myklebust comme sharepoint départ, j’ai trouvé l’approximation arctangente de [-1,1] qui comporte la plus petite erreur avec une erreur maximale de 0,94528 ulp.

 /* Based on: Tor Myklebust, "Computing accurate Horner form approximations to special functions in finite precision arithmetic", arXiv:1508.03211, August 2015. maximum ulp err = 0.94528 */ float atanf_poly (float a) { float r, s; s = a * a; r = 0x1.6d2086p-9f; // 2.78569828e-3 r = fmaf (r, s, -0x1.03f2ecp-6f); // -1.58660226e-2 r = fmaf (r, s, 0x1.5beebap-5f); // 4.24722321e-2 r = fmaf (r, s, -0x1.33194ep-4f); // -7.49753043e-2 r = fmaf (r, s, 0x1.b403a8p-4f); // 1.06448799e-1 r = fmaf (r, s, -0x1.22f5c2p-3f); // -1.42070308e-1 r = fmaf (r, s, 0x1.997748p-3f); // 1.99934542e-1 r = fmaf (r, s, -0x1.5554d8p-2f); // -3.33331466e-1 r = r * s; r = fmaf (r, a, a); return r; } 

Ce n’est pas une réponse, mais un commentaire étendu aussi.

Les processeurs Intel récents et certains futurs processeurs AMD ont AVX2. Sous Linux, recherchez l’indicateur avx2 dans /proc/cpuinfo pour savoir si votre processeur les prend en charge.

AVX2 est une extension qui nous permet de construire et de calculer à l’aide de vecteurs 256 bits – par exemple, huit nombres à simple précision ou quatre nombres à double précision – au lieu de simples scalaires. Il inclut le support FMA3, ce qui signifie multiplie-addition fusionné pour de tels vecteurs. En résumé, AVX2 nous permet d’évaluer huit polynomes en parallèle, à peu près au même moment que nous en évaluons un seul à l’aide d’opérations scalaires.

La fonction error8() parsing un ensemble de coefficients en utilisant des valeurs prédéfinies de x , en les comparant aux valeurs précalculées de atan(x) , et renvoie l’erreur dans les ULP (en dessous et au-dessus du résultat souhaité séparément), ainsi que le nombre de résultats. qui correspondent exactement à la valeur en virgule flottante souhaitée. Celles-ci ne sont pas nécessaires pour simplement tester si un ensemble de coefficients est supérieur au meilleur ensemble actuellement connu, mais permettent à différentes stratégies de tester les coefficients. (Fondamentalement, l’erreur maximale dans les ULP forme une surface, et nous essayons de trouver le point le plus bas sur cette surface; connaître la “hauteur” de la surface à chaque point nous permet de deviner de façon éclairée la direction à suivre – – comment changer les coefficients.)

Il y a quatre tables précalculées utilisées: known_x pour les arguments, known_f pour les résultats à simple précision correctement arrondis, known_a pour la valeur «précise» en double précision (j’espère juste que la bibliothèque atan() est suffisamment précise pour cela – – mais il ne faut pas s’y fier sans vérifier!), et known_m pour redimensionner la différence de double précision par rapport aux ULP. Étant donné une plage d’arguments souhaitée, la fonction atan() à l’aide de la fonction library atan() . (Il s’appuie également sur les formats à virgule flottante IEEE-754 et l’ordre identique des octets flottants et des nombres entiers, mais c’est le cas pour les CPU sur lesquelles ce code est exécuté.)

Notez que les known_x , known_f et known_a pourraient être stockés dans des fichiers binarys; le contenu de known_m est sortingvialement dérivé de known_a . Utiliser la bibliothèque atan() sans vérifier que ce n’est pas une bonne idée – mais comme les miens correspondent aux résultats de njuffa, je n’ai pas pris la peine de chercher une meilleure référence atan() .

Pour simplifier, voici le code sous la forme d’un exemple de programme:

 #define _POSIX_C_SOURCE 200809L #include  #include  #include  #include  #include  #include  /** poly8() - Compute eight polynomials in parallel. * @x - the arguments * @c - the coefficients. * * The first coefficients are for degree 17, the second * for degree 15, and so on, down to degree 3. * * The comstackr should vectorize the expression using vfmaddXXXps * given an AVX2-capable CPU; for example, Intel Haswell, * Broadwell, Haswell E, Broadwell E, Skylake, or Cannonlake; * or AMD Excavator CPUs. Tested on Intel Core i5-4200U. * * Using GCC-4.8.2 and * gcc -O2 -march=core-avx2 -mtune=generic * this code produces assembly (AT&T syntax) * vmulps %ymm0, %ymm0, %ymm2 * vmovaps (%rdi), %ymm1 * vmovaps %ymm0, %ymm3 * vfmadd213ps 32(%rdi), %ymm2, %ymm1 * vfmadd213ps 64(%rdi), %ymm2, %ymm1 * vfmadd213ps 96(%rdi), %ymm2, %ymm1 * vfmadd213ps 128(%rdi), %ymm2, %ymm1 * vfmadd213ps 160(%rdi), %ymm2, %ymm1 * vfmadd213ps 192(%rdi), %ymm2, %ymm1 * vfmadd213ps 224(%rdi), %ymm2, %ymm1 * vmulps %ymm2, %ymm1, %ymm0 * vfmadd132ps %ymm3, %ymm3, %ymm0 * ret * if you omit the 'static inline'. */ static inline __v8sf poly8(const __v8sf x, const __v8sf *const c) { const __v8sf xx = x * x; return (((((((c[0]*xx + c[1])*xx + c[2])*xx + c[3])*xx + c[4])*xx + c[5])*xx + c[6])*xx + c[7])*xx*x + x; } /** error8() - Calculate maximum error in ULPs * @x - the arguments * @co - { C17, C15, C13, C11, C9, C7, C5, C3 } * @f - the correctly rounded results in single precision * @a - the expected results in double precision * @m - 16777216.0 raised to the same power of two as @a normalized * @n - number of vectors to test * @max_under - pointer to store the maximum underflow (negative, in ULPs) to * @max_over - pointer to store the maximum overflow (positive, in ULPs) to * Returns the number of correctly rounded float results, 0..8*n. */ size_t error8(const __v8sf *const x, const float *const co, const __v8sf *const f, const __v4df *const a, const __v4df *const m, const size_t n, float *const max_under, float *const max_over) { const __v8sf c[8] = { { co[0], co[0], co[0], co[0], co[0], co[0], co[0], co[0] }, { co[1], co[1], co[1], co[1], co[1], co[1], co[1], co[1] }, { co[2], co[2], co[2], co[2], co[2], co[2], co[2], co[2] }, { co[3], co[3], co[3], co[3], co[3], co[3], co[3], co[3] }, { co[4], co[4], co[4], co[4], co[4], co[4], co[4], co[4] }, { co[5], co[5], co[5], co[5], co[5], co[5], co[5], co[5] }, { co[6], co[6], co[6], co[6], co[6], co[6], co[6], co[6] }, { co[7], co[7], co[7], co[7], co[7], co[7], co[7], co[7] } }; __v4df min = { 0.0, 0.0, 0.0, 0.0 }; __v4df max = { 0.0, 0.0, 0.0, 0.0 }; __v8si eqs = { 0, 0, 0, 0, 0, 0, 0, 0 }; size_t i; for (i = 0; i < n; i++) { const __v8sf v = poly8(x[i], c); const __v4df d0 = { v[0], v[1], v[2], v[3] }; const __v4df d1 = { v[4], v[5], v[6], v[7] }; const __v4df err0 = (d0 - a[2*i+0]) * m[2*i+0]; const __v4df err1 = (d1 - a[2*i+1]) * m[2*i+1]; eqs -= (__v8si)_mm256_cmp_ps(v, f[i], _CMP_EQ_OQ); min = _mm256_min_pd(min, err0); max = _mm256_max_pd(max, err1); min = _mm256_min_pd(min, err1); max = _mm256_max_pd(max, err0); } if (max_under) { if (min[0] > min[1]) min[0] = min[1]; if (min[0] > min[2]) min[0] = min[2]; if (min[0] > min[3]) min[0] = min[3]; *max_under = min[0]; } if (max_over) { if (max[0] < max[1]) max[0] = max[1]; if (max[0] < max[2]) max[0] = max[2]; if (max[0] < max[3]) max[0] = max[3]; *max_over = max[0]; } return (size_t)((unsigned int)eqs[0]) + (size_t)((unsigned int)eqs[1]) + (size_t)((unsigned int)eqs[2]) + (size_t)((unsigned int)eqs[3]) + (size_t)((unsigned int)eqs[4]) + (size_t)((unsigned int)eqs[5]) + (size_t)((unsigned int)eqs[6]) + (size_t)((unsigned int)eqs[7]); } /** precalculate() - Allocate and precalculate tables for error8(). * @x0 - First argument to precalculate * @x1 - Last argument to precalculate * @xptr - Pointer to a __v8sf pointer for the arguments * @fptr - Pointer to a __v8sf pointer for the correctly rounded results * @aptr - Pointer to a __v4df pointer for the comparison results * @mptr - Pointer to a __v4df pointer for the difference multipliers * Returns the vector count if successful, * 0 with errno set otherwise. */ size_t precalculate(const float x0, const float x1, __v8sf **const xptr, __v8sf **const fptr, __v4df **const aptr, __v4df **const mptr) { const size_t align = 64; unsigned int i0, i1; size_t n, i, sbytes, dbytes; __v8sf *x = NULL; __v8sf *f = NULL; __v4df *a = NULL; __v4df *m = NULL; if (!xptr || !fptr || !aptr || !mptr) { errno = EINVAL; return (size_t)0; } memcpy(&i0, &x0, sizeof i0); memcpy(&i1, &x1, sizeof i1); i0 ^= (i0 & 0x80000000U) ? 0xFFFFFFFFU : 0x80000000U; i1 ^= (i1 & 0x80000000U) ? 0xFFFFFFFFU : 0x80000000U; if (i1 > i0) n = (((size_t)i1 - (size_t)i0) | (size_t)7) + (size_t)1; else if (i0 > i1) n = (((size_t)i0 - (size_t)i1) | (size_t)7) + (size_t)1; else { errno = EINVAL; return (size_t)0; } sbytes = n * sizeof (float); if (sbytes % align) sbytes += align - (sbytes % align); dbytes = n * sizeof (double); if (dbytes % align) dbytes += align - (dbytes % align); if (posix_memalign((void **)&x, align, sbytes)) { errno = ENOMEM; return (size_t)0; } if (posix_memalign((void **)&f, align, sbytes)) { free(x); errno = ENOMEM; return (size_t)0; } if (posix_memalign((void **)&a, align, dbytes)) { free(f); free(x); errno = ENOMEM; return (size_t)0; } if (posix_memalign((void **)&m, align, dbytes)) { free(a); free(f); free(x); errno = ENOMEM; return (size_t)0; } if (x1 > x0) { float *const xp = (float *)x; float curr = x0; for (i = 0; i < n; i++) { xp[i] = curr; curr = nextafterf(curr, HUGE_VALF); } i = n; while (i-->0 && xp[i] > x1) xp[i] = x1; } else { float *const xp = (float *)x; float curr = x0; for (i = 0; i < n; i++) { xp[i] = curr; curr = nextafterf(curr, -HUGE_VALF); } i = n; while (i-->0 && xp[i] < x1) xp[i] = x1; } { const float *const xp = (const float *)x; float *const fp = (float *)f; double *const ap = (double *)a; double *const mp = (double *)m; for (i = 0; i < n; i++) { const float curr = xp[i]; int temp; fp[i] = atanf(curr); ap[i] = atan((double)curr); (void)frexp(ap[i], &temp); mp[i] = ldexp(16777216.0, temp); } } *xptr = x; *fptr = f; *aptr = a; *mptr = m; errno = 0; return n/8; } static int parse_range(const char *const str, float *const range) { float fmin, fmax; char dummy; if (sscanf(str, " %f %f %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %f:%f %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %f,%f %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %f/%f %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %ff %ff %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %ff:%ff %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %ff,%ff %c", &fmin, &fmax, &dummy) == 2 || sscanf(str, " %ff/%ff %c", &fmin, &fmax, &dummy) == 2) { if (range) { range[0] = fmin; range[1] = fmax; } return 0; } if (sscanf(str, " %f %c", &fmin, &dummy) == 1 || sscanf(str, " %ff %c", &fmin, &dummy) == 1) { if (range) { range[0] = fmin; range[1] = fmin; } return 0; } return errno = ENOENT; } static int fix_range(float *const f) { if (f && f[0] > f[1]) { const float tmp = f[0]; f[0] = f[1]; f[1] = tmp; } return f && isfinite(f[0]) && isfinite(f[1]) && (f[1] >= f[0]); } static const char *f2s(char *const buffer, const size_t size, const float value, const char *const invalid) { char format[32]; float parsed; int decimals, length; for (decimals = 0; decimals <= 16; decimals++) { length = snprintf(format, sizeof format, "%%.%df", decimals); if (length < 1 || length >= (int)sizeof format) break; length = snprintf(buffer, size, format, value); if (length < 1 || length >= (int)size) break; if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value) return buffer; decimals++; } for (decimals = 0; decimals <= 16; decimals++) { length = snprintf(format, sizeof format, "%%.%dg", decimals); if (length < 1 || length >= (int)sizeof format) break; length = snprintf(buffer, size, format, value); if (length < 1 || length >= (int)size) break; if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value) return buffer; decimals++; } length = snprintf(buffer, size, "%a", value); if (length < 1 || length >= (int)size) return invalid; if (sscanf(buffer, "%f", &parsed) == 1 && parsed == value) return buffer; return invalid; } int main(int argc, char *argv[]) { float xrange[2] = { 0.75f, 1.00f }; float c17range[2], c15range[2], c13range[2], c11range[2]; float c9range[2], c7range[2], c5range[2], c3range[2]; float c[8]; __v8sf *known_x; __v8sf *known_f; __v4df *known_a; __v4df *known_m; size_t known_n; if (argc != 10 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) { fprintf(stderr, "\n"); fprintf(stderr, "Usage: %s [ -h | --help ]\n", argv[0]); fprintf(stderr, " %s C17 C15 C13 C11 C9 C7 C5 C3 x\n", argv[0]); fprintf(stderr, "\n"); fprintf(stderr, "Each of the coefficients can be a constant or a range,\n"); fprintf(stderr, "for example 0.25 or 0.75:1. x must be a non-empty range.\n"); fprintf(stderr, "\n"); return EXIT_FAILURE; } if (parse_range(argv[1], c17range) || !fix_range(c17range)) { fprintf(stderr, "%s: Invalid C17 range or constant.\n", argv[1]); return EXIT_FAILURE; } if (parse_range(argv[2], c15range) || !fix_range(c15range)) { fprintf(stderr, "%s: Invalid C15 range or constant.\n", argv[2]); return EXIT_FAILURE; } if (parse_range(argv[3], c13range) || !fix_range(c13range)) { fprintf(stderr, "%s: Invalid C13 range or constant.\n", argv[3]); return EXIT_FAILURE; } if (parse_range(argv[4], c11range) || !fix_range(c11range)) { fprintf(stderr, "%s: Invalid C11 range or constant.\n", argv[4]); return EXIT_FAILURE; } if (parse_range(argv[5], c9range) || !fix_range(c9range)) { fprintf(stderr, "%s: Invalid C9 range or constant.\n", argv[5]); return EXIT_FAILURE; } if (parse_range(argv[6], c7range) || !fix_range(c7range)) { fprintf(stderr, "%s: Invalid C7 range or constant.\n", argv[6]); return EXIT_FAILURE; } if (parse_range(argv[7], c5range) || !fix_range(c5range)) { fprintf(stderr, "%s: Invalid C5 range or constant.\n", argv[7]); return EXIT_FAILURE; } if (parse_range(argv[8], c3range) || !fix_range(c3range)) { fprintf(stderr, "%s: Invalid C3 range or constant.\n", argv[8]); return EXIT_FAILURE; } if (parse_range(argv[9], xrange) || xrange[0] == xrange[1] || !isfinite(xrange[0]) || !isfinite(xrange[1])) { fprintf(stderr, "%s: Invalid x range.\n", argv[9]); return EXIT_FAILURE; } known_n = precalculate(xrange[0], xrange[1], &known_x, &known_f, &known_a, &known_m); if (!known_n) { if (errno == ENOMEM) fprintf(stderr, "Not enough memory for precalculated tables.\n"); else fprintf(stderr, "Invalid (empty) x range.\n"); return EXIT_FAILURE; } fprintf(stderr, "Precalculated %lu arctangents to compare to.\n", 8UL * (unsigned long)known_n); fprintf(stderr, "\nC17 C15 C13 C11 C9 C7 C5 C3 max-ulps-under max-ulps-above correctly-rounded percentage cycles\n"); fflush(stderr); { const double percent = 12.5 / (double)known_n; size_t rounded; char c17buffer[64], c15buffer[64], c13buffer[64], c11buffer[64]; char c9buffer[64], c7buffer[64], c5buffer[64], c3buffer[64]; char minbuffer[64], maxbuffer[64]; float minulps, maxulps; unsigned long tsc_start, tsc_stop; for (c[0] = c17range[0]; c[0] <= c17range[1]; c[0] = nextafterf(c[0], HUGE_VALF)) for (c[1] = c15range[0]; c[1] <= c15range[1]; c[1] = nextafterf(c[1], HUGE_VALF)) for (c[2] = c13range[0]; c[2] <= c13range[1]; c[2] = nextafterf(c[2], HUGE_VALF)) for (c[3] = c11range[0]; c[3] <= c11range[1]; c[3] = nextafterf(c[3], HUGE_VALF)) for (c[4] = c9range[0]; c[4] <= c9range[1]; c[4] = nextafterf(c[4], HUGE_VALF)) for (c[5] = c7range[0]; c[5] <= c7range[1]; c[5] = nextafterf(c[5], HUGE_VALF)) for (c[6] = c5range[0]; c[6] <= c5range[1]; c[6] = nextafterf(c[6], HUGE_VALF)) for (c[7] = c3range[0]; c[7] <= c3range[1]; c[7] = nextafterf(c[7], HUGE_VALF)) { tsc_start = __builtin_ia32_rdtsc(); rounded = error8(known_x, c, known_f, known_a, known_m, known_n, &minulps, &maxulps); tsc_stop = __builtin_ia32_rdtsc(); printf("%-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %-13s %lu %.3f %lu\n", f2s(c17buffer, sizeof c17buffer, c[0], "?"), f2s(c15buffer, sizeof c15buffer, c[1], "?"), f2s(c13buffer, sizeof c13buffer, c[2], "?"), f2s(c11buffer, sizeof c11buffer, c[3], "?"), f2s(c9buffer, sizeof c9buffer, c[4], "?"), f2s(c7buffer, sizeof c7buffer, c[5], "?"), f2s(c5buffer, sizeof c5buffer, c[6], "?"), f2s(c3buffer, sizeof c3buffer, c[7], "?"), f2s(minbuffer, sizeof minbuffer, minulps, "?"), f2s(maxbuffer, sizeof maxbuffer, maxulps, "?"), rounded, (double)rounded * percent, (unsigned long)(tsc_stop - tsc_start)); fflush(stdout); } } return EXIT_SUCCESS; } 

The code does comstack using GCC-4.8.2 on Linux, but might have to be modified for other comstackrs and/or OSes. (I'd be happy to include/accept edits fixing those, though. I just don't have Windows or ICC myself so I could check.)

To comstack this, I recommend

 gcc -Wall -O3 -fomit-frame-pointer -march=native -mtune=native example.c -lm -o example 

Run without arguments to see usage; ou

 ./example 0x1.7ed24ap-9f -0x1.0c2c12p-6f 0x1.61fdd2p-5f -0x1.3556b0p-4f 0x1.b4e138p-4f -0x1.230ae2p-3f 0x1.9978eep-3f -0x1.5554dap-2f 0.75:1 

to check what it reports for njuffa's coefficient set, compared against standard C library atan() function, with all possible x in [0.75, 1] considered.

Instead of a fixed coefficient, you can also use min:max to define a range to scan (scanning all unique single-precision floating-point values). Each possible combination of the coefficients is tested.

Because I prefer decimal notation, but need to keep the values exact, I use the f2s() function to display the floating-point values. It is a simple brute-force helper function, that uses the shortest formatting that yields the same value when parsed back to float.

Par exemple,

 ./example 0x1.7ed248p-9f:0x1.7ed24cp-9f -0x1.0c2c10p-6f:-0x1.0c2c14p-6f 0x1.61fdd0p-5f:0x1.61fdd4p-5f -0x1.3556aep-4f:-0x1.3556b2p-4f 0x1.b4e136p-4f:0x1.b4e13ap-4f -0x1.230ae0p-3f:-0x1.230ae4p-3f 0x1.9978ecp-3f:0x1.9978f0p-3f -0x1.5554d8p-2f:-0x1.5554dcp-2f 0.75:1 

computes all the 6561 (3 8 ) coefficient combinations ±1 ULP around njuffa's set for x in [0.75, 1]. (Indeed, it shows that decreasing C 17 by 1 ULP to 0x1.7ed248p-9f yields the exact same results.)

(That run took 90 seconds on Core i5-4200U at 2.6 GHz -- pretty much in line in my estimate of 30 coefficient sets per second per GHz per core. While this code is not threaded, the key functions are thread-safe, so threading should not be too difficult. This Core i5-4200U is a laptop, and gets pretty hot even when stressing just one core, so I didn't bother.)

(I consider the above code to be in public domain, or CC0-licensed where public domain dedication is not possible. In fact, I'm not sure if it is creative enough to be copyrightable at all. Anyway, feel free to use it anywhere in any way you wish, as long as you don't blame me if it breaks.)

Des questions? Enhancements? Edits to fix Linux/GCC'isms are welcome!