特殊関数: ガンマ・誤差関数の評価

概要

ガンマ関数族 ($\Gamma$, $\ln\Gamma$, $B$, $\psi$, 不完全ガンマ・ベータ) と誤差関数族 ($\operatorname{erf}$, $\operatorname{erfc}$, $\operatorname{erfcx}$, $\operatorname{erfi}$, Fresnel, Dawson) は、 確率分布・統計・物理・信号処理のいたるところに現れる。 いずれも初等関数では閉じた式に書けず、引数の領域ごとに収束の速い手法へ切り替えるのが評価の要となる。

sangi の特殊関数群は、この「領域分割 + 手法切替」を一貫した設計で実装する:

• 反射公式で左半面を右半面の評価に帰着する ($\Gamma$, $\ln\Gamma$, $\psi$)。
• 引数シフトで小さい引数を漸近展開が効く大きい引数へ持ち上げる ($\psi$, $\psi^{(n)}$, $\ln\Gamma$ の任意精度版)。
• 級数と連分数の切替で、収束の速い側だけを使う (不完全ガンマ・ベータ、$\operatorname{erfcx}$)。
• スケーリングでオーバーフロー/アンダーフローを避ける ($\ln\Gamma$, $\operatorname{erfcx}$)。

native 浮動小数点型 (float / double / long double) では、 標準ライブラリにある関数は std::tgamma, std::lgamma, std::beta, std::erf, std::erfc へ委譲し、 標準にない関数 (ディガンマ・ポリガンマ・不完全ガンマ/ベータ・$\operatorname{erfcx}$・$\operatorname{erfi}$・$\operatorname{erfInv}$・Fresnel・Dawson) は本ライブラリ独自に実装する。 Complex 型および任意精度 Float 型には専用のオーバーロードを用意し、 Lanczos 近似 (複素・倍精度) と Stirling 漸近展開 (任意精度) を使い分ける。

関連 API: SpecialFunctions。

ガンマ関数 (Lanczos と反射公式)

ガンマ関数は階乗の連続化であり、$\Gamma(n+1) = n!$、関数等式 $\Gamma(z+1) = z\,\Gamma(z)$ を満たす:

$$\Gamma(z) = \int_0^\infty t^{z-1} e^{-t}\,dt \qquad (\operatorname{Re}(z) > 0)$$

Lanczos 近似

複素・倍精度版の gamma(Complex<R>) は Lanczos 近似 ($g=7$, 係数 9 個) を用いる。 $\operatorname{Re}(z) \geq 1/2$ で次の形に評価する:

$$\Gamma(z) = \sqrt{2\pi}\;t^{\,z-\frac12}\,e^{-t}\,A(z),\qquad t = z - 1 + g + \tfrac12$$

ここで $A(z)$ は有理関数で、定数係数 $c_0,\dots,c_8$ を用いて

$$A(z) = c_0 + \sum_{k=1}^{8} \frac{c_k}{(z-1)+k}$$

と書ける。実装の係数は $g = 7$、$c_0 = 0.99999999999980993$、$c_1 = 676.5203681218851$、… という標準的な Lanczos の組で、倍精度で約 15 桁の精度を与える。

反射公式 (左半面)

$\operatorname{Re}(z) < 1/2$ では Lanczos が効かないので、反射公式で右半面に帰着する:

$$\Gamma(z)\,\Gamma(1-z) = \frac{\pi}{\sin(\pi z)} \quad\Longrightarrow\quad \Gamma(z) = \frac{\pi}{\sin(\pi z)\,\Gamma(1-z)}$$

$1-z$ は右半面に入るので、そこに Lanczos を適用する。$\sin(\pi z)$ の零点 (= $z$ が整数) が $\Gamma$ の極に対応する。

対数ガンマ $\ln\Gamma$ (オーバーフロー回避)

$\Gamma(z)$ は急速に増大し、$z$ がやや大きいだけで倍精度の範囲を超える。 ベータ関数 $B(a,b)$ や不完全ガンマの前因子のように $\Gamma$ の比・積を取る計算では、 まず 対数ガンマ $\ln\Gamma$ を評価し、最後に指数を取ることで桁あふれを避ける:

$$\ln\Gamma(z) = \left(z-\tfrac12\right)\ln t - t + \tfrac12\ln(2\pi) + \ln A(z)$$

実装では $B(a,b) = \exp\!\big(\ln\Gamma(a) + \ln\Gamma(b) - \ln\Gamma(a+b)\big)$ のように、 積・商を加減算に置き換えて評価する。 左半面では $\ln\Gamma(z) = \ln\pi - \ln\sin(\pi z) - \ln\Gamma(1-z)$ という対数形の反射公式を使う。

Stirling 漸近展開 (任意精度版)

任意精度 Float 型の lnGamma(Complex<Float>, precision) は、固定係数の Lanczos ではなく Stirling 漸近展開を使う:

$$\ln\Gamma(z) = \left(z-\tfrac12\right)\ln z - z + \tfrac12\ln(2\pi) + \sum_{k=1}^{K} \frac{B_{2k}}{2k\,(2k-1)\,z^{2k-1}}$$

$B_{2k}$ はベルヌーイ数で、要求精度に応じて Akiyama-Tanigawa 法で生成する。 漸近展開は $|z|$ が大きいほど速く収束するので、$|z|$ が小さいときは引数シフト $\Gamma(z+1) = z\,\Gamma(z)$ を繰り返して $z$ を大きくしてから展開する。 漸近級数は最良項を過ぎると発散に転じるため、項が増大し始めたら打ち切る (最適打切り) 必要がある。 要求桁に応じてシフト量とベルヌーイ項数を調整する。 $\Gamma$ 自体は $\exp(\ln\Gamma)$ で得る。

ディガンマ・ポリガンマ

ディガンマ関数は対数ガンマの導関数、ポリガンマはその高階導関数:

$$\psi(z) = \frac{d}{dz}\ln\Gamma(z) = \frac{\Gamma'(z)}{\Gamma(z)},\qquad \psi^{(n)}(z) = \frac{d^{\,n+1}}{dz^{\,n+1}}\ln\Gamma(z)$$

$\psi_1 = \psi^{(1)}$ をトリガンマと呼ぶ。これらは標準ライブラリにないので独自に実装する。

引数シフト + 漸近展開

$\psi$ の漸近展開は $z$ が大きいほど精度がよい:

$$\psi(z) \approx \ln z - \frac{1}{2z} - \sum_{k=1}^{K} \frac{B_{2k}}{2k\,z^{2k}}$$

そこで、再帰式 (関数等式の対数微分)

$$\psi(z+1) = \psi(z) + \frac{1}{z}$$

を使って引数を漸近展開が効く大きさまで持ち上げる。 native 版は $z \geq 8$ まで上げてから 7 項のベルヌーイ係数で評価する。 トリガンマも同様に $\psi_1(z) = \psi_1(z+1) + 1/z^2$ でシフトしてから

$$\psi_1(z) \approx \frac{1}{z} + \frac{1}{2z^2} + \sum_{k=1}^{K} \frac{B_{2k}}{z^{2k+1}}$$

を評価する。一般の $\psi^{(n)}$ ($n \geq 2$) は $z \geq 10$ までシフトし、係数に $(-1)^{n+1}$ の符号と階乗因子 $\prod_{j=1}^{n-1}(2k+j)$ を含む一般の漸近展開を用いる (項が相対精度を切るか、項が増大に転じたら打ち切る)。

反射公式 (負の引数)

$\psi$ の左半面は反射公式で右半面に帰着する:

$$\psi(1-z) - \psi(z) = \pi\cot(\pi z),\qquad \psi_1(1-z) + \psi_1(z) = \frac{\pi^2}{\sin^2(\pi z)}$$

$z$ が非正整数のときは極なので NaN を返す。任意精度 Float 版は、ベルヌーイ数を要求桁で生成し、 $\ln\Gamma$ と同じく最適打切りで漸近級数を制御する。

不完全ガンマ・不完全ベータ

正則化不完全ガンマ関数は、ガンマ積分を途中で切ったものを $\Gamma(a)$ で割って $[0,1]$ に正規化したもの:

$$P(a,x) = \frac{\gamma(a,x)}{\Gamma(a)} = \frac{1}{\Gamma(a)}\int_0^x t^{a-1}e^{-t}\,dt,\qquad Q(a,x) = \frac{\Gamma(a,x)}{\Gamma(a)} = 1 - P(a,x)$$

これらはカイ二乗分布・ガンマ分布・ポアソン分布の累積分布などに直結する。

級数と連分数の切替

下側 $P(a,x)$ にはテイラー級数が、上側 $Q(a,x)$ には連分数が向く。 sangi は境界 $x = a+1$ で両者を切り替える:

• $x < a+1$: $P$ をテイラー級数 $$P(a,x) = \frac{x^a e^{-x}}{\Gamma(a)}\sum_{n=0}^{\infty}\frac{x^n}{a(a+1)\cdots(a+n)}$$ で評価し、$Q = 1 - P$ とする。
• $x \geq a+1$: $Q$ を Legendre の連分数 $$Q(a,x) = \frac{x^a e^{-x}}{\Gamma(a)}\cdot \cfrac{1}{x+1-a-\cfrac{1\cdot(1-a)}{x+3-a-\cfrac{2\cdot(2-a)}{x+5-a-\cdots}}}$$ で評価し、$P = 1 - Q$ とする。

前因子 $x^a e^{-x}/\Gamma(a)$ は $\exp(-x + a\ln x - \ln\Gamma(a))$ として対数経由で計算し、桁あふれを避ける。 連分数は modified Lentz 法で評価する: 分母が 0 に近づくと破綻するので、微小量 $\texttt{tiny}$ でクランプしながら $C_n = b_n + a_n/C_{n-1}$, $D_n = 1/(b_n + a_n D_{n-1})$, $f_n = f_{n-1}\,C_n D_n$ を進め、$|C_n D_n - 1|$ が機械精度を切ったら止める。 非正則化版 $\gamma(a,x) = P(a,x)\,\Gamma(a)$, $\Gamma(a,x) = Q(a,x)\,\Gamma(a)$ も提供する。

正則化不完全ベータ $I_x(a,b)$

不完全ベータはベータ分布・スチューデントの $t$ 分布・$F$ 分布の累積に現れる:

$$I_x(a,b) = \frac{1}{B(a,b)}\int_0^x t^{a-1}(1-t)^{b-1}\,dt$$

これも連分数で評価する。連分数は $x$ が小さい側で速く収束するので、 対称関係 $I_x(a,b) = 1 - I_{1-x}(b,a)$ を使い、$x \geq (a+1)/(a+b+2)$ なら $x \to 1-x$, $a \leftrightarrow b$ と入れ替えて収束の速い側に持ち込む。 前因子は $\exp(a\ln x + b\ln(1-x) + \ln\Gamma(a+b) - \ln\Gamma(a) - \ln\Gamma(b))/a$ と対数経由で計算し、 連分数本体は不完全ガンマと同じ modified Lentz 法で評価する (奇数項・偶数項で係数 $d_n$ の形が変わる標準の展開)。

誤差関数 (erf / erfc / erfcx / erfi / erfInv)

誤差関数とその仲間は正規分布・拡散・Voigt 関数などに現れる:

$$\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\,dt,\qquad \operatorname{erfc}(x) = 1 - \operatorname{erf}(x)$$

native 版の $\operatorname{erf}$, $\operatorname{erfc}$ は標準ライブラリ (std::erf, std::erfc) へ委譲する。 標準にない $\operatorname{erfcx}$, $\operatorname{erfi}$, $\operatorname{erfInv}$ と、 複素・任意精度の各オーバーロードを独自に実装する。

スケール化補誤差関数 $\operatorname{erfcx}$

大きな $x$ では $\operatorname{erfc}(x)$ が指数的にアンダーフローし、$e^{x^2}$ がオーバーフローするので、 その積を 1 つの関数として安定に評価する:

$$\operatorname{erfcx}(x) = e^{x^2}\operatorname{erfc}(x) \sim \frac{1}{x\sqrt{\pi}} \sum_{n=0}^{\infty}\frac{(-1)^n (2n-1)!!}{(2x^2)^n} \qquad (x \to +\infty)$$

native 版は領域で手法を切り替える:

• $x > 4$: 上の漸近展開 (項が増大に転じたら打ち切る最適打切り)。
• $0 < x \leq 4$: $e^{x^2}$ がまだ溢れないので直接 $e^{x^2}\operatorname{erfc}(x)$。
• $x < 0$: 反射 $\operatorname{erfcx}(x) = 2e^{x^2} - \operatorname{erfcx}(-x)$。

複素版は $|z| \le$ しきい値で直接計算、$\operatorname{Re}(z) \ge 0$ かつ $|z|$ 大では Laplace 連分数 ($a_n = n/2$) を modified Lentz 法で評価し、 $\operatorname{erfcx}(z) = 1/(\sqrt{\pi}\,g)$ とする。$\operatorname{Re}(z) < 0$ では同じ反射公式を使う。

虚数誤差関数 $\operatorname{erfi}$

虚数誤差関数は

$$\operatorname{erfi}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{t^2}\,dt = -i\,\operatorname{erf}(ix)$$

で定義される奇関数。$\operatorname{erf}(ix)$ は純虚数なので、複素版 $\operatorname{erf}$ に $z = ix$ を渡してその虚部を返す。複素 $\operatorname{erf}$ は、桁落ちを避けるため正項のみを足し上げる 級数 ($_1F_1$ 由来) で全平面に対して収束する。$e^{x^2}$ のオーダーで発散するので、 倍精度では $|x| \gtrsim 27$ で $\pm\infty$ にクランプする。

逆誤差関数 $\operatorname{erfInv}$

$\operatorname{erf}(y) = p$ を満たす $y$ を Newton 反復で求める。導関数は $\dfrac{d}{dy}\operatorname{erf}(y) = \dfrac{2}{\sqrt{\pi}}e^{-y^2}$ なので:

$$y_{n+1} = y_n - \frac{\operatorname{erf}(y_n) - p}{\frac{2}{\sqrt{\pi}}e^{-y_n^2}}$$

初期値は領域で分ける:

• $|p| \leq 0.7$: 一次近似 $y_0 = \frac{\sqrt{\pi}}{2}\,p$。
• $|p| > 0.7$: 裾の近似 $y_0 = \operatorname{sign}(p)\sqrt{-\ln\big((1-|p|)(1+|p|)\big)}$。

数回の反復で収束する。$p = \pm 1$ で $\pm\infty$、$|p| > 1$ で NaN を返す。

Fresnel 積分・Dawson 関数

Fresnel 積分 $C(x), S(x)$

Fresnel 積分は回折・光学・クロソイド曲線に現れる (DLMF と同じ正規化):

$$C(x) = \int_0^x \cos\!\Big(\frac{\pi t^2}{2}\Big)\,dt,\qquad S(x) = \int_0^x \sin\!\Big(\frac{\pi t^2}{2}\Big)\,dt$$

sangi は $|x| \leq 4$ でテイラー級数、$|x| > 4$ で漸近展開に切り替える。 級数は

$$C(x) = \sum_{k=0}^{\infty}\frac{(-1)^k (\pi/2)^{2k}\,x^{4k+1}}{(2k)!\,(4k+1)},\qquad S(x) = \sum_{k=0}^{\infty}\frac{(-1)^k (\pi/2)^{2k+1}\,x^{4k+3}}{(2k+1)!\,(4k+3)}$$

で、$C$ と $S$ を漸化的に同時に積み上げる。$|x| > 4$ では補助関数 $f(x), g(x)$ を用いた漸近形

$$C(x) = \tfrac12 + f(x)\sin\!\Big(\frac{\pi x^2}{2}\Big) - g(x)\cos\!\Big(\frac{\pi x^2}{2}\Big),\qquad S(x) = \tfrac12 - f(x)\cos\!\Big(\frac{\pi x^2}{2}\Big) - g(x)\sin\!\Big(\frac{\pi x^2}{2}\Big)$$

を、最小項での最適打切りで評価する ($f, g$ の項比はそれぞれ $-(4k+1)(4k+3)/(\pi x^2)^2$, $-(4k+3)(4k+5)/(\pi x^2)^2$)。$C, S$ は奇関数なので $x<0$ では符号反転で対応する。 $x \to \pm\infty$ では $C, S \to \pm 1/2$。

Dawson 関数 $D(x)$

Dawson 関数は

$$D(x) = e^{-x^2}\int_0^x e^{t^2}\,dt$$

で定義される奇関数で、$\operatorname{erfi}$ と $D(x) = \frac{\sqrt{\pi}}{2}e^{-x^2}\operatorname{erfi}(x)$ で結ばれる。 $\operatorname{erfi}$ や $e^{x^2}$ と違って有界 (最大 $\approx 0.541$) なので、裾の評価に向く。 sangi は Rybicki のアルゴリズムに沿って領域で切り替える:

• $|x| < 0.2$: テイラー級数 $D(x) = \sum_{n=0}^{\infty}\dfrac{(-1)^n 2^n x^{2n+1}}{(2n+1)!!}$。
• $|x| > 5$: 漸近展開 $D(x) = \dfrac{1}{2x}\sum_{n=0}^{\infty}\dfrac{(2n-1)!!}{(2x^2)^n}$。
• 中間: Simpson 則で $\displaystyle\int_0^x e^{t^2-x^2}\,dt$ を直接数値積分。

奇関数性 $D(-x) = -D(x)$ で負の引数に対応する。

比較表

共通する設計指針は、(1) 反射公式で左半面を右半面に、(2) 引数シフトで小引数を大引数に帰着し、 (3) 級数 (小引数) と連分数・漸近展開 (大引数) を境界で切り替え、 (4) 対数経由・スケーリングで桁あふれを避ける、という 4 点に集約される。 漸近級数は発散級数なので、必ず最良項での最適打切りを伴う。

参考文献

• Abramowitz, M. & Stegun, I. A. (eds.) (1972). Handbook of Mathematical Functions. Dover. (§6 ガンマ関数, §7 誤差関数, 6.3.18 ディガンマ漸近展開)
• NIST Digital Library of Mathematical Functions (DLMF). https://dlmf.nist.gov/. (§5 ガンマ関数, §7 誤差・Fresnel・Dawson, §8 不完全ガンマ・ベータ)
• Lanczos, C. (1964). "A precision approximation of the gamma function". Journal of the SIAM, Series B: Numerical Analysis, 1(1), 86–96.
• Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007). Numerical Recipes. 3rd ed. Cambridge University Press. (不完全ガンマ・ベータの級数/連分数, modified Lentz 法, Dawson の Rybicki 法)
• Cody, W. J. (1969). "Rational Chebyshev approximation for the error function". Mathematics of Computation, 23(107), 631–637.
• Rybicki, G. B. (1989). "Dawson's integral and the sampling theorem". Computers in Physics, 3(2), 85–87.