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

10

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^{2 m} e^{- t u^{2}} d 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}}}$ dla

m = 0, 1, 2, . . .

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

t > 0

$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 dla. 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 . Może to jednak być wystarczająco wysoka dla moich celów. $t$ gamma(m+0.5_dp) $m\approx 172$ $m$

efficiency accuracy special-functions

— Ondřej Čertík
źródło

2

Po co kodować własną funkcję? GSL, cephes i SLATEC to wszystko implementują.

— Geoff Oxberry

Zaktualizowałem pytanie, dlaczego nie używam SLATEC.

— Ondřej Čertík

@ OndřejČertík Chyba odkryłeś błąd! Głosowałem na twoje pytanie!

— Ali

Ali --- to nie jest błąd w SLATEC, ale faktem, że faktycznie muszę podzielić

przez

γ (z, x)

$\gamma(z, x)$

w celu uzyskania wartości

. Zatem metoda numeryczna, która działa dla

może nie działać tak dobrze dla

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, przepraszam, mój błąd, nie sprawdziłem kodu przed opublikowaniem komentarza.

— Ali

9

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 $m$ . $a \le m + 1{1\over 2}$

W przypadku dużych argumentów obliczam

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

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

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

$\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 $\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żego, $\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)$

F_{m} (a) = \frac{1}{2} \frac{1}{m + \frac{1}{2}} \exp (- a) (1 + \sum_{n = 1}^{\infty} \frac{a^{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)$ aby obliczyć wszystkie pozostałe wartości funkcji.

— njuffa
źródło

Dzięki @njuffa za świetną odpowiedź. Jeśli stworzysz kod dla tego open source, myślę, że byłby bardzo użyteczny dla wielu osób.

— Ondřej Čertík

1

Obecnie implementacja CUDA opisanego algorytmu jest dostępna do bezpłatnego pobrania ze strony internetowej dewelopera NVIDIA (wymaga bezpłatnej rejestracji jako programista CUDA, zatwierdzenie zwykle w ciągu jednego dnia roboczego). Kod jest objęty licencją BSD, która powinna być kompatybilna z każdym projektem.

— njuffa

5

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

— John D. Cook
źródło

To świetna książka, dziękuję za podpowiedź!

— Ondřej Čertík

4

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.

— Wolfgang Bangerth
źródło

Korzystałem z tego: dlmf.nist.gov/8 , podczas jego wdrażania, ale to prawdopodobnie inny zasób. Rozdział 5 w Recepturach numerycznych również zawiera interesujące informacje, ale dotyczy tylko funkcji jednej zmiennej.

— Ondřej Čertík

Nie sądzę, że znajdziesz coś znacznie nowszego niż ich wzmianka z 2001 roku; SLATEC będzie starszy.

— Geoff Oxberry

1

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.

— Ali
źródło