Implémentation double précision rapide et précise de la fonction gamma incomplète

10

Quelle est la manière la plus moderne de mettre en œuvre des fonctions spéciales double précision? J'ai besoin de l'intégrale suivante: pour et , qui peut être écrit en termes de fonction gamma incomplète inférieure. Voici mon implémentation Fortran et C:

F_{m} (t) = \int_{0}^{1} u^{2 m} e^{- t u^{2}} ré u = \frac{γ (m + \frac{1}{2}, t)}{2 t^{m + \frac{1}{2}}}

$F_m(t) = \int_0^1 u^{2m} e^{-tu^2} d u = {\gamma(m+{1\over 2}, t)\over 2 t^{m+{1\over 2}}}$

m = 0, 1, 2, . . .

$m=0, 1, 2, ...$

t > 0

$t>0$

https://gist.github.com/3764427

qui utilise l'expansion en série, résume les termes jusqu'à la précision donnée, puis utilise les relations de récursivité pour obtenir efficacement des valeurs pour inférieur . Je l'ai bien testé et j'obtiens une précision de 1e-15 pour toutes les valeurs de paramètres dont j'ai besoin, voir les commentaires de la version Fortran pour plus de détails. $m$

Existe-t-il un meilleur moyen de le mettre en œuvre? Voici une implémentation de la fonction gamma dans gfortran:

https://github.com/mirrors/gcc/blob/master/libgfortran/intrinsics/c99_functions.c#L1781

il utilise l'approximation de fonction rationnelle au lieu de résumer une série infinie que je fais. Je pense que c'est une meilleure approche, car il faut obtenir une précision uniforme. Existe-t-il une manière canonique d'aborder ces choses, ou faut-il trouver un algorithme spécial pour chaque fonction spéciale?

Mise à jour 1 :

Sur la base des commentaires, voici l'implémentation utilisant SLATEC:

https://gist.github.com/3767621

il reproduit les valeurs de ma propre fonction, à peu près au niveau de précision 1e-15. Cependant, j'ai remarqué un problème: pour t = 1e-6 et m = 50, le terme est égal à 1e-303 et pour un "m" supérieur, il commence simplement à donner des réponses incorrectes. Ma fonction n'a pas ce problème, car j'utilise une série de relations expansion / récurrence directement pour . Voici un exemple de valeur correcte: $t^{m+{1\over2}}$ $F_m$

$F_{100}$ (1e-6)=4.97511945200351715E-003 ,

mais je ne peux pas obtenir cela en utilisant SLATEC parce que le dénominateur explose. Comme vous pouvez le voir, la valeur réelle de est belle et petite. $F_m$

Mise à jour 2 :

Pour éviter le problème ci-dessus, on peut utiliser la fonction dgamit(la fonction Gamma incomplète de Tricomi), donc F(m, t) = dgamit(m+0.5_dp, t) * gamma(m+0.5_dp) / 2, il n'y a plus de problème avec , mais malheureusement le souffle pour . Cela pourrait cependant être assez élevé pour mes besoins. $t$ gamma(m+0.5_dp) $m\approx 172$ $m$

efficiency accuracy special-functions Ondřej Čertík
la source

2

Pourquoi coder votre propre fonction? GSL, cephes et SLATEC l'implémentent tous.

Geoff Oxberry

J'ai mis à jour la question de savoir pourquoi je n'utilise pas SLATEC.

Ondřej Čertík

@ OndřejČertík Vous avez apparemment découvert un bug! A voté pour votre question!

Ali

Ali --- ce n'est pas un bug dans SLATEC, mais dans le fait, j'ai en fait besoin de diviser le

par

γ (z, x)

$\gamma(z, x)$

afin d'obtenir une valeur pour

. Ainsi, la méthode numérique qui fonctionne pour

pourrait ne pas fonctionner aussi bien pour

.

t^{m + \frac{1}{2}}

$t^{m+{1\over2}}$

F_{m} (t)

$F_m(t)$

γ (z, x)

$\gamma(z, x)$

F_{m} (t)

$F_m(t)$

Ondřej Čertík

@ OndřejČertík OK, désolé, mon erreur, je n'ai pas vérifié votre code avant de faire mon commentaire.

Ali

9

L'intégrale en question est également connue sous le nom de fonction Boys, d'après le chimiste britannique Samuel Francis Boys qui a introduit son utilisation au début des années 1950. Il y a quelques années, j'avais besoin de calculer cette fonction en double précision, le plus rapidement possible mais avec précision. J'ai réussi à obtenir une erreur relative de l'ordre de dans le domaine d'entrée entier. $10^{-15}$

Il est généralement avantageux d'utiliser différentes approximations pour les petits et les grands arguments, où la commutation optimale entre "grand" et "petit" est mieux déterminée expérimentalement, et est en général fonction de . Pour mon code, j'ai défini les "petits" arguments comme ceux satisfaisant la condition $m$ . $a \le m + 1{1\over 2}$

Pour les gros arguments, je calcule

F_{m} (a) = \frac{1}{2} γ (m + \frac{1}{2}, a) \times p \times p, p = a^{- \frac{1}{2} (m + \frac{1}{2})}

$\mathrm{F}_m(a) = {1\over 2}\gamma\left(m + {1\over 2}, a\right) \times p \times p, \space \space p = a^{-{1\over 2}\left(m+ {1\over 2}\right)}$

Cet ordre d'opérations évite un débordement prématuré. Comme nous n'avons besoin que de la fonction gamma incomplète inférieure des ordres demi-entiers ici plutôt que d'une fonction gamma incomplète inférieure entièrement générale, il est avantageux du point de vue des performances de calculer

γ (m + \frac{1}{2}, a) = Γ (m + \frac{1}{2}) - Γ (m + \frac{1}{2}, a)

$\gamma \left(m + {1\over 2}, a\right) = \Gamma \left(m + {1\over 2}\right) - \Gamma\left(m + {1\over 2}, a\right)$

en utilisant des valeurs tabulées de et calcul $\Gamma \left(m + {1\over 2}\right)$ selon cette réponse, tout en évitant soigneusement le problème de l'annulation soustractive par l'utilisation d'une opération multipliée-ajoutée fusionnée. Une optimisation supplémentaire potentielle consiste à observer que pour,suffisamment grand $\Gamma\left(m + {1\over 2}, a\right)$ $a$ à l'intérieur d'une précision en virgule flottante donnée. $\gamma \left(m + {1\over 2}, a\right) = \Gamma \left(m + {1\over 2}\right)$

Pour les petits arguments, j'ai commencé par une expansion en série pour la fonction gamma incomplète inférieure de

A. Erdelyi, W. Magnus, F. Oberhettinger et FG Tricomi, "Fonctions transcendantales supérieures, Vol. 2". New York, NY: McGraw-Hill 1953

et l'a modifié pour calculer la fonction de garçons comme suit (tronquer la série lorsque le terme est suffisamment petit pour une précision donnée): $\mathrm{F}_{m}(a)$

F_{m} (une) = \frac{1}{2} \frac{1}{m + \frac{1}{2}} \exp (- une) (1 + \sum_{n = 1}^{\infty} \frac{{une}^{n}}{(1 + m + \frac{1}{2}) \times . . . \times (n + m + \frac{1}{2})})

$\mathrm{F}_{m}(a) = {1\over 2}\frac{1}{m + {1\over 2}}\exp(-a)\left(1+\sum_{n=1}^{\infty}\frac{a^{n}}{(1 + m + {1\over 2}) \times\space ...\space \times (n + m + {1\over 2})}\right)$

$m = 0, 1, 2, 3$ $\mathrm{F}_{0}(a) = \sqrt{\frac{\pi}{4a}}\mathrm{erf}\left(\sqrt{a}\right)$ $\mathrm{erf}$ ERFerferff

$m = 1, 2, 3$ $a \lt {2{1\over 2}}$ $\mathrm{F}_{m}(a) = \frac{1}{2a}\left(\left(2m-1\right)\mathrm{F}_{m-1}(a) - \exp(-a)\right)$

$a$ $m$ $m$ $\mathrm{F}_{m-1} = \frac{1}{2m-1} \left(2a \space \mathrm{F}_{m}(a) + \exp\left(-a\right)\right)$

njuffa
la source

Merci @njuffa pour la bonne réponse. Si vous créez votre code pour cette source ouverte, je pense qu'il serait très utile à beaucoup de gens.

Ondřej Čertík

1

Actuellement, une implémentation CUDA de l'algorithme décrit est disponible en téléchargement gratuit sur le site Web des développeurs de NVIDIA (nécessite une inscription gratuite en tant que développeur CUDA, approbation généralement dans un délai d'un jour ouvrable). Le code est sous licence BSD, qui devrait être compatible avec à peu près n'importe quel type de projet.

njuffa

5

Vous pouvez jeter un œil aux méthodes numériques pour les fonctions spéciales d'Amparo Gil, Javier Segura et Nico M. Temme.

John D. Cook
la source

Ceci est un excellent livre, merci pour le conseil!

Ondřej Čertík

4

Je voudrais jeter un œil au livre d'Abramowicz & Stegun, ou à la nouvelle révision que le NIST a publiée il y a quelques années et qui est disponible en ligne, je crois. Ils discutent également des moyens de mettre en œuvre les choses de manière stable.

Wolfgang Bangerth
la source

J'utilisais ceci: dlmf.nist.gov/8 , lors de sa mise en œuvre, mais c'est probablement une autre ressource. Le chapitre 5 de Recettes numériques contient également des informations intéressantes, mais uniquement applicables aux fonctions d'une variable.

Ondřej Čertík

Je ne pense pas que vous trouverez quelque chose de beaucoup plus récent que leur référence de 2001; SLATEC sera plus ancien que cela.

Geoff Oxberry

1

Il ne semble pas être à la pointe de la technologie, mais SLATEC chez Netlib propose «1400 routines mathématiques et statistiques à usage général». Le Gamma incomplet est disponible sous les fonctions spéciales ici .

L'implémentation de telles fonctions prend du temps et est sujette aux erreurs, donc je ne le ferais pas moi-même sauf en cas d'absolue nécessité. SLATEC existe depuis un certain temps maintenant et est largement utilisé, au moins en fonction du nombre de téléchargements , donc je m'attends à ce que l'implémentation soit mature.

Ali
la source

Implémentation double précision rapide et précise de la fonction gamma incomplète

Réponses: