最終更新: 2021/8/2

特殊関数の精度保証

柏木雅英

1. はじめに

本ライブラリが持つ精度保証付き数値積分solverや精度保証付き常微分方程式solverを用いて、特殊関数の精度保証を行う関数をいくつか試験的に作成している。

内部で数値積分や常微分方程式を解いているため、基本的に遅い。しかし、単純な原理に基づいているため精度保証されていることが明快に理解でき、また計算に使う型を変えることによって高精度計算を行うこともできるので、役に立つ場面もあるのではないかと考えている。現状複素数に対応していない。順次対応させるつもりではあるが、予定は未定である。

2. Gamma関数

gamma.hppの中に、

gamma関数
lgamma関数
digamma関数
trigamma関数

の実装がある。いずれもinterval<T>を入力し、interval<T>を出力する。

詳細は、ガンマ関数の精度保証付き計算メモ(Version: 2021/6/26) に書いた。

なお、gamma関数の極値点を計算するのにdigammaに対するKrawczyk法を行う関係で、呼び出しの初回は非常に時間がかかる。

3. Lobachevsky関数

Λ(θ) = - ∫₀^θ log|2sin(t)|dt

で定義される関数をLobachevsky関数という。 lobachevsky.hppの中に、関数lobachevskyがある。 interval<T>を入力し、interval<T>を出力する。

詳細はまだ書いてない。

4. Airy関数

airy.hppの中に、

Airy_Ai
Airy_Bi
Airy_Ai_d
Airy_Bi_d

の実装がある。それぞれ、Ai関数、Bi関数、Ai関数の微分、Bi関数の微分である。いずれもinterval<T>を入力し、interval<T>を出力する。

アルゴリズムは、単純にAiry関数を定義する初期値問題

f''(t) = t f(t)

を初期値

f(0) = 1/(3^2/3Γ(2/3)), f'(0) = -1/(3^1/3Γ(1/3)) (Aiの場合)
f(0) = 1/(3^1/6Γ(2/3)), f'(0) = 1/(3^1/6Γ(1/3)) (Biの場合)

を使って解くだけ。 2階に直して解くので副産物として微分も得られる。初期値問題を時刻0から単に解いていくだけなので、絶対値の大きい入力に対しては精度も出ないし計算時間もかかる。

5. Bessel関数

bessel.hppの中に、第1種Bessel関数の実装がある。~~整数次のみ~~ 非整数次も試験的に実装。

アルゴリズムは、積分 J_n(x) = 1/π * ∫₀^π cos(x sin t - nt) dt を計算するだけ。
非整数次の場合は、 J_nu(x) = (x/2)^nu / sqrt(π) / gamma(nu + 0.5) * ∫₀^π cos(x * cos(t)) * sin(t)^2nu dt を計算。nu>-0.5が必要。

6. Beta関数

beta.hppにある。引数を計算しやすいところにreduceした後、直接数値積分している。

7. ガウスの超幾何関数

hypergeom.hppにある。ガウスの超幾何関数 ₂F₁(a,b,c;z) を計算する。

~~オイラーの積分表示を直接計算している。そのため、 c > b > 0, |z| < 1 のようなパラメータの制限がある。~~

計算方法は hypergeom.pdf (Version: 2015/11/8) の通り。途中に出てくる区間公比を持つ等比級数の和の公式は、 geoseries.pdf (Version: 2015/11/8) にある。