Szybka i dokładna implementacja niepełnej funkcji gamma o podwójnej precyzji

Jaki jest najnowocześniejszy sposób wdrażania funkcji specjalnych podwójnej precyzji? Potrzebuję następującej całki:

$$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}}}$$ dla$m=0, 1, 2, ...$oraz$t>0$, które można zapisać w kategoriach niższej niepełnej funkcji gamma. Oto moja implementacja Fortran i C.

https://gist.github.com/3764427

który korzysta z interpretacji szeregowej, sumuje warunki do danej dokładności, a następnie używa relacji rekurencyjnych w celu skutecznego uzyskania wartości dla niższego $m$ . Testowałem to dobrze i uzyskuję dokładność 1e-15 dla wszystkich wartości parametrów, których potrzebuję, zobacz komentarze do wersji Fortran, aby uzyskać szczegółowe informacje.

Czy istnieje lepszy sposób na jego wdrożenie? Oto implementacja funkcji gamma w gfortran:

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

wykorzystuje przybliżenie funkcji racjonalnej zamiast sumowania niektórych nieskończonych szeregów, które robię. Myślę, że to lepsze podejście, ponieważ należy uzyskać jednolitą dokładność. Czy istnieje jakiś kanoniczny sposób podejścia do tych rzeczy, czy też trzeba opracować specjalny algorytm dla każdej funkcji specjalnej?

Aktualizacja 1 :

W oparciu o komentarze, oto implementacja wykorzystująca SLATEC:

https://gist.github.com/3767621

odtwarza wartości z mojej własnej funkcji, mniej więcej na poziomie dokładności 1e-15. Zauważyłem jednak problem, że dla t = 1e-6 im = 50, wyrażenie równa się 1e-303, a dla większego „m” zaczyna po prostu dawać nieprawidłowe odpowiedzi. Moja funkcja nie ma tego problemu, ponieważ używam relacji rozszerzania / powtarzania serii bezpośrednio dlaFm. Oto przykład poprawnej wartości:$t^{m+{1\over2}}$$F_m$

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

ale nie mogę tego uzyskać za pomocą SLATEC, ponieważ mianownik się wysadza. Jak widać, rzeczywista wartość jest ładna i mała.$F_m$

Aktualizacja 2 :

Aby uniknąć powyższego problemu, można użyć funkcji dgamit(niekompletna funkcja gamma Tricomi) F(m, t) = dgamit(m+0.5_dp, t) * gamma(m+0.5_dp) / 2, więc nie ma już problemu z , ale niestety wysadza się dla m ≈ 172 . Może to jednak być wystarczająco wysoka m dla moich celów.$t$gamma(m+0.5_dp)$m\approx 172$$m$

Całka, o której mowa, znana jest również jako funkcja Boys, od brytyjskiego chemika Samuela Francisa Boysa, który wprowadził jej stosowanie na początku lat 50. XX wieku. Kilka lat temu musiałem obliczyć tę funkcję z podwójną precyzją, tak szybko, jak to możliwe, ale dokładnie. Udało mi się osiągnąć względny błąd rzędu w całej domenie wejściowej.$10^{-15}$

Zasadniczo korzystne jest stosowanie różnych przybliżeń dla małych i dużych argumentów, przy czym optymalne przełączanie między „dużym” i „małym” najlepiej określa się eksperymentalnie i ogólnie jest funkcją . W moim kodzie zdefiniowałem „małe” argumenty jako spełniające warunek a ≤ m + 1 $m$ .$a \le m + 1{1\over 2}$

W przypadku dużych argumentów obliczam

$$\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)}$$

Ta kolejność operacji pozwala uniknąć przedwczesnego niedomiaru. Ponieważ potrzebujemy tutaj tylko niższej niepełnej funkcji gamma rzędów o wartościach całkowitych, a nie w pełni ogólnej niższej niepełnej funkcji gamma, z punktu widzenia wydajności korzystne jest obliczenie

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

przy użyciu tabelarycznych wartości i obliczeniaΓ(m+$\Gamma \left(m + {1\over 2}\right)$zgodnie z tą odpowiedzią, przy jednoczesnym ostrożnym unikaniu problemu subtraktywnego anulowania poprzez zastosowanie połączonej operacji wielokrotnego dodawania. Potencjalną dalszą optymalizacją jest zaobserwowanie, że dla wystarczająco dużegoa,γ(m+$\Gamma\left(m + {1\over 2}, a\right)$$a$w ramach określonej precyzji zmiennoprzecinkowej.$\gamma \left(m + {1\over 2}, a\right) = \Gamma \left(m + {1\over 2}\right)$

W przypadku małych argumentów zacząłem od rozszerzenia serii dla niższej niepełnej funkcji gamma od

A. Erdelyi, W. Magnus, F. Oberhettinger i FG Tricomi, „Higher Transcendental Functions, Vol. 2”. Nowy Jork, NY: McGraw-Hill 1953

i zmodyfikowałem go, aby obliczyć funkcję Boys w następujący sposób (obcięcie serię gdy termin jest dostatecznie mała dla danej dokładności)$\mathrm{F}_{m}(a)$

$$\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)$ aby obliczyć wszystkie pozostałe wartości funkcji.

Możesz rzucić okiem na metody numeryczne dla funkcji specjalnych autorstwa Amparo Gil, Javiera Segury i Nico M. Temme.

Rzuciłbym okiem na książkę Abramowicza i Steguna lub na nowszą wersję, którą NIST opublikował kilka lat temu i wydaje mi się, że jest dostępna online. Dyskutują także o sposobach wdrażania rzeczy w stabilny sposób.

Nie wydaje się to być najnowocześniejszym, ale SLATEC w Netlib oferuje „1400 ogólnych procedur matematycznych i statystycznych”. Niekompletna wartość gamma jest dostępna tutaj w ramach funkcji specjalnych .

Wdrożenie takich funkcji jest czasochłonne i podatne na błędy, więc nie zrobiłbym tego sam, chyba że absolutnie konieczne. SLATEC istnieje już od dłuższego czasu i jest szeroko stosowany, przynajmniej w oparciu o liczbę pobrań , więc oczekiwałbym, że implementacja będzie dojrzała.