Implementasi ganda presisi yang cepat dan akurat untuk fungsi gamma tidak lengkap

Apa cara canggih menerapkan fungsi khusus presisi ganda? Saya membutuhkan integral berikut: untuk dan , yang dapat ditulis dalam bentuk fungsi gamma bawah yang tidak lengkap. Inilah implementasi Fortran dan C saya: m=0,1,2,. . . t>

$$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, ...$$t>0$

https://gist.github.com/3764427

yang menggunakan ekspansi seri, meringkas istilah sampai akurasi yang diberikan, dan kemudian menggunakan hubungan rekursi untuk secara efisien mendapatkan nilai untuk lebih rendah . Saya mengujinya dengan baik dan saya mendapatkan akurasi 1e-15 untuk semua nilai parameter yang saya butuhkan, lihat komentar versi Fortran untuk detailnya.$m$

Apakah ada cara yang lebih baik untuk mengimplementasikannya? Berikut ini adalah implementasi fungsi gamma di gfortran:

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

itu menggunakan perkiraan fungsi rasional alih-alih menyimpulkan beberapa seri tak terbatas yang saya lakukan. Saya pikir itu pendekatan yang lebih baik, karena seseorang harus mendapatkan akurasi yang seragam. Apakah ada beberapa cara kanonik untuk mendekati hal-hal ini, atau apakah kita harus mencari algoritma khusus untuk setiap fungsi khusus?

Pembaruan 1 :

Berdasarkan komentar, berikut adalah implementasinya menggunakan SLATEC:

https://gist.github.com/3767621

itu mereproduksi nilai dari fungsi saya sendiri, kira-kira pada tingkat akurasi 1e-15. Namun, saya melihat masalah bahwa untuk t = 1e-6 dan m = 50, istilah sama dengan 1e-303 dan untuk "m" yang lebih tinggi, ia mulai memberikan jawaban yang salah. Fungsi saya tidak memiliki masalah ini, karena saya menggunakan hubungan ekspansi / perulangan rangkaian langsung untuk . Berikut ini contoh nilai yang benar: F$t^{m+{1\over2}}$$F_m$

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

tapi saya tidak bisa mendapatkan ini menggunakan SLATEC karena penyebutnya meledak. Seperti yang Anda lihat, nilai sebenarnya dari bagus dan kecil.$F_m$

Pembaruan 2 :

Untuk menghindari masalah di atas, seseorang dapat menggunakan fungsi dgamit( fungsi Gamma Tricomi yang tidak lengkap), maka F(m, t) = dgamit(m+0.5_dp, t) * gamma(m+0.5_dp) / 2, jadi tidak ada masalah dengan lagi, tapi sayangnya pukulan untuk . Namun ini mungkin cukup tinggi untuk tujuan saya.m ≈ 172 $t$gamma(m+0.5_dp)$m\approx 172$$m$

Integral yang dimaksud juga dikenal sebagai fungsi Boys, setelah ahli kimia Inggris Samuel Francis Boys yang memperkenalkan penggunaannya pada awal 1950-an. Beberapa tahun yang lalu, saya perlu menghitung fungsi ini dalam presisi ganda, secepat mungkin tetapi akurat. Saya berhasil mencapai kesalahan relatif pada urutan di seluruh domain input.$10^{-15}$

Secara umum menguntungkan untuk menggunakan pendekatan yang berbeda untuk argumen kecil dan besar, di mana peralihan optimal antara "besar" dan "kecil" paling baik ditentukan secara eksperimental, dan secara umum merupakan fungsi dari . Untuk kode saya, saya mendefinisikan argumen "kecil" sebagai kondisi yang memuaskan a ≤ m + 1 $m$ .$a \le m + 1{1\over 2}$

Untuk argumen besar, saya menghitung

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

Urutan operasi ini menghindari underflow dini. Karena kita hanya memerlukan fungsi gamma rendah tidak lengkap dari pesanan setengah bilangan bulat di sini daripada fungsi gamma lengkap tidak lengkap yang lebih rendah, menguntungkan dari perspektif kinerja untuk menghitung

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

menggunakan nilai tabulasi dan komputasiΓ(m+$\Gamma \left(m + {1\over 2}\right)$sesuai dengan jawaban ini, sambil dengan hati-hati menghindari masalah pembatalan kurang melalui penggunaan operasi penambahan-penggandaan yang digabungkan. Sebuah optimasi lebih lanjut potensial adalah untuk mengamati bahwa untuk cukup besaryang,γ(m+$\Gamma\left(m + {1\over 2}, a\right)$$a$ke dalam presisi floating-point yang diberikan.$\gamma \left(m + {1\over 2}, a\right) = \Gamma \left(m + {1\over 2}\right)$

Untuk argumen kecil, saya mulai dengan ekspansi seri untuk fungsi gamma rendah lengkap dari

A. Erdelyi, W. Magnus, F. Oberhettinger, dan FG Tricomi, "Fungsi Transendental Tinggi, Vol. 2". New York, NY: McGraw-Hill 1953

dan dimodifikasi untuk menghitung Boys fungsi sebagai berikut (truncating seri ketika istilah ini cukup kecil untuk presisi yang diberikan):$\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)$

Anda dapat melihat Metode Numerik untuk Fungsi Khusus oleh Amparo Gil, Javier Segura, dan Nico M. Temme.

Saya akan melihat buku Abramowicz & Stegun, atau revisi yang lebih baru yang diterbitkan NIST beberapa tahun yang lalu dan itu tersedia secara online, saya percaya. Mereka juga membahas cara-cara untuk mengimplementasikan berbagai hal dengan cara yang stabil.

Tampaknya tidak canggih, tetapi SLATEC di Netlib menawarkan "1.400 tujuan umum matematika dan statistik." Gamma tidak lengkap tersedia di bawah fungsi khusus di sini .

Menerapkan fungsi-fungsi seperti itu memakan waktu dan rawan kesalahan jadi saya tidak akan melakukannya sendiri kecuali benar-benar diperlukan. SLATEC telah ada selama beberapa waktu sekarang dan banyak digunakan, setidaknya berdasarkan jumlah unduhan , jadi saya berharap implementasinya akan matang.