Calcul précis du format PDF de la dissortingbution normale standard à l’aide de la bibliothèque mathématique standard C

La fonction de densité de probabilité de la dissortingbution normale standard est définie par e- x 2/2 / √ (2π). Cela peut être rendu directement dans le code C. Un exemple d’implémentation simple précision pourrait être:

float my_normpdff (float a) { return 0x1.988454p-2f * my_expf (-0.5f * a * a); /* 1/sqrt(2*pi) */ } 

Bien que ce code soit exempt de débordement prématuré, il existe un problème de précision car l’erreur commise dans le calcul du 2/2 est amplifiée par l’exponentiation ultérieure. On peut facilement le démontrer avec des tests sur des références plus précises. L’erreur exacte différera en fonction de la précision des implémentations exp() ou expf() utilisées; Pour les fonctions d’exponentiation fidèlement arrondies, on observera généralement une erreur maximale d’environ 2 6 ulps pour IEEE-754 binary32 simple précision, d’environ 2 9 ulps pour IEEE-754 binary64 double précision.

Comment le problème de précision peut-il être traité de manière raisonnablement efficace? Une approche sortingviale consisterait à utiliser un calcul intermédiaire de haute précision, par exemple double calcul double pour la mise en œuvre de float . Mais cette approche ne fonctionne pas pour une double implémentation si l’arithmétique à virgule flottante de précision plus élevée n’est pas facilement disponible et peut être inefficace pour une implémentation float si le double arithmétique est nettement plus coûteux que le calcul float , par exemple sur de nombreux GPU.

La question de la précision soulevée dans la question peut être résolue de manière efficace en utilisant des quantités limitées de calcul double- float ou double- double , facilitées par l’utilisation de l’opération de multiplication-addition fusionnée (FMA).

Cette opération est disponible depuis C99 via les fonctions mathématiques standard fmaf(a,b,c) et fma(a,b,c) qui calculent a * b + c, sans arrondir le produit intermédiaire . Bien que les fonctions correspondent directement aux opérations matérielles rapides sur presque tous les processeurs modernes, elles peuvent utiliser le code d’émulation sur les anciennes plates-formes, auquel cas elles peuvent être très lentes.

Cela permet de calculer le produit avec deux fois la précision normale en utilisant seulement deux opérations, ce qui donne une paire tête: queue de nombres de précision native:

 prod_hi = a * b // head prod_lo = FMA (a, b, -hi) // tail 

Les bits de poids fort du résultat peuvent être transmis à l’exponentiation, tandis que les bits de poids faible sont utilisés pour améliorer la précision du résultat par interpolation linéaire, en tirant parti du fait que e x est sa propre dérivée:

 e = exp (prod_hi) + exp (prod_hi) * prod_lo // exp (a*b) 

Cela nous permet d’éliminer la plupart des erreurs de la mise en œuvre naïve. L’autre source d’erreur de calcul, mineure, est la précision limitée avec laquelle la constante 1 / √ (2π) est représentée. Cela peut être amélioré en utilisant une représentation tête: queue pour la constante qui fournit deux fois la précision native, et en calculant:

 r = FMA (const_hi, x, const_lo * x) // const * x 

Le document suivant souligne que cette technique peut même conduire à une multiplication correctement arrondie pour certaines constantes de précision arbitraire:

Nicolas Brisebarre et Jean-Michel Muller, “Multiplication correctement arrondie par des constantes de précision arbitraires”, IEEE Transactions on Computers, vol. 57, n ° 2, février 2008, p. 165-174

En combinant les deux techniques et en prenant en charge quelques cas complexes impliquant des NaN, nous arrivons à la mise en œuvre de float suivante basée sur la norme binary32 IEEE-754:

 float my_normpdff (float a) { const float RCP_SQRT_2PI_HI = 0x1.988454p-02f; /* 1/sqrt(2*pi), msbs */ const float RCP_SQRT_2PI_LO = -0x1.857936p-27f; /* 1/sqrt(2*pi), lsbs */ float ah, sh, sl, ea; ah = -0.5f * a; sh = a * ah; sl = fmaf (a, ah, 0.5f * a * a); /* don't flip "sign bit" of NaN argument */ ea = expf (sh); if (ea != 0.0f) ea = fmaf (sl, ea, ea); /* avoid creation of NaN */ return fmaf (RCP_SQRT_2PI_HI, ea, RCP_SQRT_2PI_LO * ea); } 

La double implémentation correspondante, basée sur IEEE-754 binary64 , semble presque identique, à l’exception des différentes valeurs constantes utilisées:

 double my_normpdf (double a) { const double RCP_SQRT_2PI_HI = 0x1.9884533d436510p-02; /* 1/sqrt(2*pi), msbs */ const double RCP_SQRT_2PI_LO = -0x1.cbc0d30ebfd150p-56; /* 1/sqrt(2*pi), lsbs */ double ah, sh, sl, ea; ah = -0.5 * a; sh = a * ah; sl = fma (a, ah, 0.5 * a * a); /* don't flip "sign bit" of NaN argument */ ea = exp (sh); if (ea != 0.0) ea = fma (sl, ea, ea); /* avoid creation of NaN */ return fma (RCP_SQRT_2PI_HI, ea, RCP_SQRT_2PI_LO * ea); } 

La précision de ces implémentations dépend de la précision des fonctions mathématiques standard expf() et exp() , respectivement. Lorsque la bibliothèque de mathématiques C fournit des versions fidèlement arrondies de celles-ci, l’erreur maximale de l’une ou l’autre des deux implémentations ci-dessus est généralement inférieure à 2,5 ulps.