特殊関数: ゼータ・超幾何・その他

概要

特殊関数の評価で繰り返し現れるのは、級数和・漸近展開・解析接続の三つの道具立てである。 sangi の special 名前空間は、これらを関数ごとに使い分けて評価する。

その中心にあるのがゼータ関数族と超幾何関数族で、多くの特殊関数を統一的に表す枠組みになっている。 Lerch 超越 $\Phi(z,s,a) = \sum_{k=0}^{\infty} z^k/(k+a)^s$ は、$z=1$ で Hurwitz ゼータ $\zeta(s,a)$、$z=-1$ で交代和 (Dirichlet $\beta$)、$z$ を残せば多重対数 $\mathrm{Li}_s(z)$ をすべて含む。 一方、一般化超幾何 $_pF_q$ は、初等関数・Bessel 関数・合流型関数を冪級数の特別な場合として束ねる。

このページでは、ゼータ族・超幾何族を軸に、そこから派生する指数積分・対数積分・三角積分、さらに Lambert W・Owen's T・Coulomb 波動関数・Fermi-Dirac/Debye/輸送積分・角運動量結合係数 (Wigner 記号) までを、実装した評価アルゴリズムに即して順に解説する。

関連 API: SpecialFunctions。 対応型は概ね float / double / long double (ネイティブ浮動小数点)、一部は任意精度 Float および複素 Complex<double> / Complex<Float> に対応する。

ゼータ関数族

提供する関数: Riemann $\zeta(s)$、Hurwitz $\zeta(s,a)$、Dirichlet $\eta(s)$、Clausen 関数 $\mathrm{Cl}_n(\theta)$、Lerch 超越 $\Phi(z,s,a)$。

Dirichlet η を経由した Riemann ζ

$\zeta(s)$ の素朴な定義級数 $\sum_{n\ge 1} n^{-s}$ は $\mathrm{Re}(s) \le 1$ で発散し、$\mathrm{Re}(s) > 1$ でも収束が遅い。 そこで交代ゼータ (Dirichlet $\eta$) を経由する:

$$\eta(s) = \sum_{k=1}^{\infty} \frac{(-1)^{k-1}}{k^s} = \bigl(1 - 2^{1-s}\bigr)\,\zeta(s), \qquad \zeta(s) = \frac{\eta(s)}{1 - 2^{1-s}}.$$

$\eta(s)$ は交代級数なので加速が効く。複素ネイティブ型では Hasse 級数 (Euler 変換) を使う:

$$\eta(s) = \sum_{n=0}^{N} \frac{1}{2^{n+1}} \sum_{k=0}^{n} (-1)^k \binom{n}{k}\,(k+1)^{-s}.$$

誤差は項あたり約 1 ビット ($2^{-N}$) で減るので、double 相当なら $N=50$ で約 15 桁が出る。 任意精度 Complex<Float> 版では Euler 変換の項数を作業精度のビット数に比例させて取り、相対 $\varepsilon$ で打ち切る。 実数 double 版は std::riemann_zeta / std::pow に委譲しつつ、$\eta(1)=\ln 2$ など極を避けた特別値を直接返す。

関数等式による解析接続

$\mathrm{Re}(s) < 0$ では、関数等式で収束域 $\mathrm{Re}(1-s) > 1$ に折り返す:

$$\zeta(s) = 2^{s}\,\pi^{s-1}\,\sin\!\left(\frac{\pi s}{2}\right) \Gamma(1-s)\,\zeta(1-s).$$

右辺の $\zeta(1-s)$ は収束域にあるので $\eta$ 経由で評価でき、$\Gamma(1-s)$ はガンマ関数の実装を呼ぶ。 負の偶数 $s=-2n$ では $\sin(\pi s/2)=0$ から自明零点 $\zeta(-2n)=0$ を直接返す。

Hurwitz ζ と Euler-Maclaurin 総和

Hurwitz ゼータ $\zeta(s,a) = \sum_{k=0}^{\infty}(a+k)^{-s}$ は、最初の $N$ 項を直接和し、残りをEuler-Maclaurin 総和で近似する:

$$\zeta(s,a) \approx \sum_{k=0}^{N-1}(a+k)^{-s} + \frac{(a+N)^{1-s}}{s-1} + \frac{(a+N)^{-s}}{2} + \sum_{j=1}^{M} \frac{B_{2j}}{(2j)!}\,(s)_{2j-1}\,(a+N)^{-(s+2j-1)},$$

ここで $B_{2j}$ は Bernoulli 数、$(s)_{2j-1}$ は上昇階乗である。 積分項 $(a+N)^{1-s}/(s-1)$ と中点補正 $(a+N)^{-s}/2$ に Bernoulli 補正級数が続く。 引数シフト $N$ で $a+N$ を十分大きくすると補正級数が速く収束する。

double 版は経験的に最適な $N=15$・$M=8$ を固定する ($M>8$ は Bernoulli 数の発散で精度が劣化するため)。 任意精度 Complex<Float> 版は、作業精度に対し $a+N \gtrsim 0.4\,\mathrm{wp}$ を満たすよう $N$ を取り、補正項数の上限を $M \approx \pi(a+N)$ として、項の大きさが相対 $\varepsilon$ を下回った時点で早期に打ち切る。 Bernoulli 数 $B_{2k}$ はガンマ実装と共有する。 $a=1$ なら $\zeta(s,1)=\zeta(s)$ に帰着し、$s=1$ は極なので NaN を返す。

Clausen 関数

Clausen 関数 $\mathrm{Cl}_2(\theta) = \sum_{k=1}^{\infty}\sin(k\theta)/k^2 = -\int_0^\theta \ln\!\bigl|2\sin(t/2)\bigr|\,dt$ は、周期性と反対称性で引数を $[0,\pi]$ に畳んでから評価する。

• $\theta$ が小さいとき: 対数を含む小角展開 $\mathrm{Cl}_2(\theta) = \theta(1-\ln\theta) + \theta^3/36 - \theta^5/3600 + \cdots$
• 一般の $\theta$: Fourier 級数 $\sum \sin(k\theta)/k^2$ を Kahan 補正和で直接加算

一般化版 $\mathrm{Cl}_n(\theta)$ は、$n$ が偶数なら $\sum\cos(k\theta)/k^n$、奇数なら $\sum\sin(k\theta)/k^n$ を相対 $\varepsilon$ まで和する。

Lerch 超越と統一表現

Lerch 超越 $\Phi(z,s,a) = \sum_{k=0}^{\infty} z^k/(k+a)^s$ は、多重対数・Hurwitz ゼータ・Riemann ゼータ・Dirichlet $\eta$ を一つの枠組みに統一する:

$$\mathrm{Li}_s(z) = z\,\Phi(z,s,1),\quad \zeta(s,a) = \Phi(1,s,a),\quad \zeta(s) = \Phi(1,s,1),\quad \eta(s) = -\Phi(-1,s,1).$$

実装は引数の値で分岐する: $z=0$ なら $a^{-s}$、$s=0$ なら $1/(1-z)$、$z=1$ なら Hurwitz $\zeta(s,a)$ ($s>1$ のみ有限)、$z=-1$ なら Hasse 加速した交代和、$|z|<1$ なら直接 Taylor 級数。 $|z|>1$ ($z \neq \pm 1$) の解析接続は将来拡張で、現状は NaN を返す。

超幾何関数族

提供する関数: Gauss $_2F_1(a,b;c;z)$、合流型 $_1F_1(a;b;z)$ (Kummer M)、合流極限 $_0F_1(;b;z)$、一般化 $_pF_q$、Meijer G、Mittag-Leffler $E_{\alpha,\beta}(z)$。

冪級数とその収束域

一般化超幾何関数は冪級数で定義される:

$$_pF_q(a_1,\ldots,a_p;\,b_1,\ldots,b_q;\,z) = \sum_{k=0}^{\infty} \frac{(a_1)_k \cdots (a_p)_k}{(b_1)_k \cdots (b_q)_k}\,\frac{z^k}{k!}.$$

隣接項の比 $t_{k+1}/t_k = z\,\prod_i(a_i+k)\big/\bigl[(k+1)\prod_j(b_j+k)\bigr]$ を漸化的に掛けて積み上げるので、Pochhammer 記号を陽に計算する必要はない。収束域は次のように決まる:

• $p \le q$: 任意の $z$ で収束 (整関数)。$_0F_1$, $_1F_1$ がこれにあたる
• $p = q+1$: $|z| < 1$ で収束、$z=1$ ではパラメータ依存。$_2F_1$ がこれにあたる
• $p > q+1$: 形式的 (漸近) 級数 — 本実装は扱わず NaN を返す

いずれの場合も、$a_i$ のどれかが非正整数なら有限和の多項式になり、項数を打ち切って全 $z$ で評価できる。 $b_j$ が非正整数 (ガンマの極) かつ対応する $a_i$ が先に終端しないなら極として NaN を返す。

変換公式による収束域の拡大

$_2F_1$ は $|z| < 1$ でしか級数収束しないため、収束の悪い領域では引数を写す:

$$\text{Pfaff 変換:}\quad {}_2F_1(a,b;c;z) = (1-z)^{-a}\,{}_2F_1\!\left(a,\,c-b;\,c;\,\frac{z}{z-1}\right).$$

$z < -0.5$ のとき $z/(z-1)$ は原点に近づくので、写してから級数和すると収束が速い。 合流型 $_1F_1$ では Kummer 変換 $M(a,b,z) = e^{z}\,M(b-a,b,-z)$ を $z<0$ で使い、桁落ちを避ける。 $_0F_1$ は全 $z$ で収束する整関数で、$_0F_1(;1;-z^2/4)=J_0(z)$、$_0F_1(;1/2;-z^2/4)=\cos z$ といった関係をもつ。

Meijer G (部分実装)

Meijer G 関数 $G^{m,n}_{p,q}$ は超幾何族をさらに一般化する。本実装は留数表現を使った部分的なもので、$n=0$ かつ $a$ が空、$m \le q$、各 $b_j$ が相異なる単純極という設定に限る。 この設定では G は $b_h$ での単純極における留数和として、係数 (ガンマ関数の積) と $z^{b_h}$ と $_0F_{q-1}$ の積になる。$b_j$ 間の差が整数になる共鳴 (対数項が現れる) や $n>0$ の一般形は将来拡張で、NaN を返す。

Mittag-Leffler 関数

Mittag-Leffler 関数は指数関数を分数階に一般化したもので、分数階微分方程式の解に現れる:

$$E_{\alpha,\beta}(z) = \sum_{k=0}^{\infty} \frac{z^k}{\Gamma(\alpha k + \beta)}, \qquad E_{\alpha}(z) = E_{\alpha,1}(z).$$

$|z|$ が小さいときは Taylor 級数を逆ガンマ $1/\Gamma(\alpha k+\beta)$ ごとに項を積む (ガンマの極は $1/\Gamma=0$ として処理)。 $|z|$ が閾値を超えると Taylor が遅く不安定になるため、$z>0$ では Wiman 漸近、$z<0$ では指数項を除いた減衰漸近に切り替える:

$$E_{\alpha,\beta}(z) \approx \frac{1}{\alpha}\,z^{(1-\beta)/\alpha}\,e^{z^{1/\alpha}} - \sum_{k=1}^{N} \frac{z^{-k}}{\Gamma(\beta-\alpha k)}\quad(z>0).$$

漸近級数は発散するので、項の大きさが増加に転じた時点で打ち切る。 $\alpha=1,\beta=1$ で $E_{1,1}(z)=e^z$、$\alpha=2,\beta=1,z\ge 0$ で $E_{2,1}(z)=\cosh\sqrt{z}$ など、初等関数に帰着する特別値は直接返す。$\alpha \in (0,2]$ に限り、それ以外は例外を投げる。

指数積分・対数積分・三角積分

提供する関数: 指数積分 $\mathrm{Ei}(x)$、対数積分 $\mathrm{li}(x)$、一般化指数積分 $E_n(x)$、正弦積分 $\mathrm{Si}(x)$、余弦積分 $\mathrm{Ci}(x)$、双曲線版 $\mathrm{Shi},\mathrm{Chi}$、dilogarithm $\mathrm{Li}_2(x)$。

級数と連分数の使い分け

指数積分は原点近傍の収束級数を使う:

$$\mathrm{Ei}(x) = \gamma + \ln|x| + \sum_{n=1}^{\infty} \frac{x^n}{n\cdot n!},$$

ここで $\gamma$ は Euler-Mascheroni 定数である。対数積分は $\mathrm{li}(x) = \mathrm{Ei}(\ln x)$ で実装する。

一般化指数積分 $E_n(x) = \int_1^\infty e^{-xt}/t^n\,dt$ は、引数の大きさで手法を切り替える:

• $x>1$ または $n>10$: 修正 Lentz 法による連分数 $E_n(x) = e^{-x}\,\mathrm{CF}$。大きい $x$ で速く安定
• それ以外: まず級数で $E_1(x)$ を求め、漸化式 $n\,E_{n+1}(x) = e^{-x} - x\,E_n(x)$ で上方に押し上げる

三角積分

正弦積分・余弦積分は Taylor 級数で評価する:

$$\mathrm{Si}(x) = \sum_{k=0}^{\infty} \frac{(-1)^k\,x^{2k+1}}{(2k+1)\,(2k+1)!}, \qquad \mathrm{Ci}(x) = \gamma + \ln|x| + \sum_{k=1}^{\infty} \frac{(-1)^k\,x^{2k}}{2k\,(2k)!}.$$

項は漸化的に更新し、相対 $\varepsilon$ で打ち切る。$\mathrm{Si}$ は奇関数で $\mathrm{Si}(\infty)=\pi/2$、$\mathrm{Ci}$ は $x>0$ のみ実数値をとる。 双曲線版 $\mathrm{Shi},\mathrm{Chi}$ は交代符号を除いた級数で評価し、$|x|>20$ では $\mathrm{Ei}$ 経由 ($\mathrm{Shi}=(\mathrm{Ei}(x)-\mathrm{Ei}(-x))/2$ 等) に切り替える。

dilogarithm

dilogarithm $\mathrm{Li}_2(x) = \sum_{n\ge 1} x^n/n^2$ は、$|x|\le 0.5$ で直接 Taylor、$0.5 < |x| \le 1$ で反射公式を使って収束を確保する:

$$\mathrm{Li}_2(x) = -\mathrm{Li}_2(1-x) + \frac{\pi^2}{6} - \ln x\,\ln(1-x).$$

複素版はさらに $|z|>1$ で反転公式 $\mathrm{Li}_2(z) = -\mathrm{Li}_2(1/z) - \pi^2/6 - (\ln(-z))^2/2$ を併用し、反射で無限再帰しないよう $|1-z| < |z|$ を確認してから適用する。

その他の特殊関数

Lambert W (Halley 反復・分枝)

Lambert W は $W(x)\,e^{W(x)} = x$ を満たす関数で、Halley 反復 (三次収束) で解く:

$$w_{n+1} = w_n - \frac{w_n e^{w_n} - x}{e^{w_n}(w_n+1) - \dfrac{(w_n+2)(w_n e^{w_n} - x)}{2w_n+2}}.$$

主枝 $W_0$ ($x \ge -1/e$、$W_0 \ge -1$) と第二枝 $W_{-1}$ ($-1/e \le x < 0$、$W_{-1} \le -1$) を別々に扱う。 収束を速めるため、初期推定を引数の領域ごとに変える:

• $-1/e$ 近傍: $p=\sqrt{2(ex+1)}$ を使った級数 $W \approx -1 + p - p^2/3 + \cdots$
• 小さい $x$: Fritsch 近似
• 大きい $x$: 漸近形 $W \approx \ln x - \ln(\ln x)$

複素引数 $W_0,W_{-1}$ も同じ Halley 反復で評価する。任意精度版は double で初期推定を作り、Halley 反復で目標精度に詰める。

Owen's T

Owen's T 関数 $T(h,a) = \dfrac{1}{2\pi}\displaystyle\int_0^a \dfrac{e^{-h^2(1+t^2)/2}}{1+t^2}\,dt$ は、二変量正規分布の象限確率や非心 $t$・歪正規分布の CDF に使う。

• ネイティブ浮動小数点: 10 点 Gauss-Legendre 求積で内部積分を直接計算
• 任意精度 Float: 被積分関数 $e^{-h^2 t^2/2}/(1+t^2)$ の Cauchy 積級数 $\sum d_n t^{2n}$ ($d_n = a_n - d_{n-1}$) を展開

$a>1$ では Owen の変換 $T(h,a) = \tfrac{1}{2}[\Phi(h)+\Phi(ah)] - \Phi(h)\Phi(ah) - T(ah,1/a)$ で $a \le 1$ の評価に帰着させる。

Coulomb 波動関数

正則 Coulomb 波動関数 $F_L(\eta,\rho)$ は、規格化定数 $C_L(\eta)$ と冪級数で評価する:

$$F_L(\eta,\rho) = C_L(\eta)\,\rho^{L+1}\sum_{k\ge 0} a_k\,\rho^k,\qquad a_{k+1} = \frac{2\eta a_k - a_{k-1}}{(k+1)(k+2L+2)}.$$

規格化定数 $C_L(\eta)$ は Sommerfeld パラメータ $2\pi\eta/(e^{2\pi\eta}-1)$ から漸化的に作る。 非正則関数 $G_L(\eta,\rho)$ は、中程度以上の $\rho$ では漸近位相 (Coulomb 位相シフト) と Wronskian 関係 $F_L G_L' - F_L' G_L = 1$ から構成し、小さい $\rho$ では特異なため NaN を返す。

Fermi-Dirac 積分

完全 Fermi-Dirac 積分 $F_s(x) = \dfrac{1}{\Gamma(s+1)}\displaystyle\int_0^\infty \dfrac{t^s}{e^{t-x}+1}\,dt$ は、半導体物理で電子密度に現れる。

• $x>20$: Sommerfeld 展開 $F_s(x) \approx \dfrac{x^{s+1}}{(s+1)\Gamma(s+1)}\bigl[1 + \tfrac{s(s+1)\pi^2}{6x^2} + \cdots\bigr]$
• それ以外: 16 点 Gauss-Laguerre 求積

$F_{1/2},F_{-1/2},F_{3/2}$ の専用入口を用意している。

Debye 関数・輸送積分

Debye 関数 $D_n(x) = \dfrac{n}{x^n}\displaystyle\int_0^x \dfrac{t^n}{e^t-1}\,dt$ は固体の比熱 ($D_3$) などに、輸送積分 $J_n(x) = \displaystyle\int_0^x \dfrac{t^n e^t}{(e^t-1)^2}\,dt$ は輸送現象に現れる。 どちらも 15 点 Gauss-Legendre 求積で内部積分を計算し、原点近傍の被積分関数の特異 ($t^{n-1}$ や $t^{n-2}$ への漸近) を別扱いして可積分性を保つ。 シンクロトロン関数 $F(x),G(x)$ は変形 Bessel $K_\nu$ の漸近・積分近似を経由する。

角運動量結合係数 (Wigner 記号)

Wigner 3j/6j/9j 記号と Clebsch-Gordan 係数は、角運動量の結合に現れる。 半整数の角運動量は $2j$ を整数で持つ ($\texttt{two\_j} = 2j$) ことで丸め誤差を避ける。

• 3j: 三角条件・選択則を確認したうえで Racah 公式を評価する。階乗は $\ln$ 階乗 (lgamma) で扱い、桁あふれを避けて $\exp(\ln\text{前因子} - \ln\text{項})$ で和をとる
• Clebsch-Gordan: 位相因子と $\sqrt{2J+1}$ を掛けた 3j に帰着
• 6j: 4 つの三角条件と Racah 和公式
• 9j: 中間変数 $t$ について 6j の三重積の和 $\sum_t (-1)^{2t}(2t+1)\,\{\cdots\}\{\cdots\}\{\cdots\}$ に展開

比較表

大づかみには、ゼータ族は級数加速 + 関数等式、超幾何族は冪級数 + 変換公式、積分で定義される関数 (Owen's T・Fermi-Dirac・Debye・輸送) はGauss 求積、方程式で定義される関数 (Lambert W) はHalley 反復、というのが基本方針である。

参考文献

• NIST Digital Library of Mathematical Functions (DLMF). https://dlmf.nist.gov/ — §16 (一般化超幾何・Meijer G), §25 (ゼータ・多重対数・Lerch).
• Abramowitz, M. & Stegun, I. A. (1964). Handbook of Mathematical Functions. National Bureau of Standards.
• Borwein, P. (2000). "An efficient algorithm for the Riemann zeta function". Constructive, Experimental, and Nonlinear Analysis, CMS Conf. Proc. 27, 29–34.
• Corless, R. M., Gonnet, G. H., Hare, D. E. G., Jeffrey, D. J. & Knuth, D. E. (1996). "On the Lambert W function". Advances in Computational Mathematics, 5(1), 329–359.
• Gorenflo, R., Kilbas, A. A., Mainardi, F. & Rogosin, S. V. (2014). Mittag-Leffler Functions, Related Topics and Applications. Springer.
• Patefield, M. & Tandy, D. (2000). "Fast and accurate calculation of Owen's T function". Journal of Statistical Software, 5(5), 1–25.
• Edmonds, A. R. (1957). Angular Momentum in Quantum Mechanics. Princeton University Press (Wigner 記号・Racah 公式).