\( \def\vector#1{\boldsymbol{#1}} \)

ガンマ関数と階乗

Takami Torao
  • このエントリーをはてなブックマークに追加

ガンマ関数 (gamma function) は 0 または負の整数以外の複素数に対して以下の積分で定義される特殊関数。\[ \Gamma(z) = \int_0^\infty t^{z-1} e^{-t} dt \] また漸化式 \(\Gamma(z+1) = z \Gamma(z)\) でも表すことができる。正の整数 \(n\) に対して階乗関数 (factorial function) \(\Gamma(n+1) = n!\) となる。

アルゴリズム

ガンマ関数

Lanczos の近似法より実数 \(z \gt 0\) において以下の近似が成り立つ。\[ \begin{align*} \Gamma(z+1) & = \sqrt{2\pi} \left(z + g + \frac{1}{2}\right)^{z + \frac{1}{2}} e^{-(z+g+\frac{1}{2})} A_g(z) \\ A_g(z) & = \frac{1}{2}p_0(g) + p_1(g)\frac{z}{z+1} + p_2(g)\frac{z(z-1)}{(z+1)(z+2)}+\cdots \end{align*} \] ここで \(g\) は制約 \({\rm Re}(z) \gt \frac{1}{2}\) を満たせば任意に選択できる定数である。\(p_k(g)\) は少々複雑で: \[ p_k(g) = \sum_{i=0}^k C(2k+1,2i+1) \frac{\sqrt{2}}{\pi} \left(i-\frac{1}{2}\right)! \left(i+g+\frac{1}{2}\right)^{-\left(i+\frac{1}{2}\right)} e^{i+g+\frac{1}{2}} \] \(C(i,j)\) は Chebyshev 多項式の係数行列における \((i,j)\) 番目を意味する。

\(C(i,j)\) 自体の計算式は以下の通り難しくはないが、\(C(5,5)\) より上を求めるには再帰が多すぎて StackOverflow が発生する。右の実装は Scala のトランポリンを使用して再帰を書き換えたものだがそれでも OutOfMemoryError が発生する。

def c(i:Int, j:Int):Int = (i, j) match {
  case (1, 1) => 1
  case (2, 2) => 1
  case (i, 1) => - c(i-2, 1)
  case (i, j) if i == j => 2 * c(i-1, j-1)
  case (i, j) => 2 * c(i-1, j-1) - c(i-2, j)
}
import scala.util.control.TailCalls._
def c(i:Int, j:Int):TailRec[Int] = (i, j) match {
  case (1, 1) => done(1)
  case (2, 2) => done(1)
  case (i, 1) => tailcall{ c(i-2, 1) }.map{ x => -x }
  case (i, j) if i == j => tailcall{ c(i-1, j-1) }.map{ x => 2 * x }
  case (i, j) => tailcall{ c(i-1, j-1) }.flatMap{ x => tailcall{ c(i-2, j) }.map(y => 2 * x - y) }
}

このままの数式で係数 \(p_k(g)\) を算出するのは難しいと思われるため先人の計算結果を利用する。級数 \(A_g(z)\) を以下のように置き換えた場合: \[ A_g(z) = c_0 + \sum_{k=1}^N \frac{c_k}{z + k} + \epsilon \] 適切な \(g\) (通常は小さな整数) を選択すれば \(N = 5\sim 10\) 程度で \(|\epsilon|\) が一般的な単精度/倍精度浮動小数点での実数演算に十分な精度まで収束する。例えば \(g=5\) とすると定数 \(c_0\) は 1.0 に近い値をとり \(N=6\) の近傍で \(|\epsilon| < 2\times 10^{-10}\) まで収束する。

対数ガンマ関数の実装

数値演算で扱うことのできるガンマ関数の限界について考える。引数を整数 \(n\) としたとき \(\Gamma(n)\) が倍精度浮動小数点の最大値 \(1.8\times 10^{308}\) を超えるのは \(\Gamma(172)=1.2\times 10^{309}\) である。つまり (小数点以下のどこまでかは省略するが) 実数演算で扱うことの出来るガンマ関数は \(0 \lt z \lt 172\) の範囲である (同様に単精度浮動小数点の場合は \(0 \lt z \lt 36\))。このような小さな値でオーバーフローを起こすことは数値演算を使う上でいささか都合が悪い。

しかし、数式にガンマ関数が出現する時は分数の分子/分母で互いに打ち消し合うケースがしばしば見られる。このためコンピュータ上での数値演算は以下のように対数ガンマを求めた後に減算して指数演算できた方が有用である。\[ \begin{align*} \frac{\Gamma(a)}{\Gamma(b)} & = {\rm exp}\left\{\log \frac{\Gamma(a)}{\Gamma(b)} \right\} \\ & = {\rm exp}\left(\log \Gamma(a) - \log \Gamma(b) \right) \end{align*} \] 通常のガンマ関数は対数ガンマ関数の結果を指数関数に適用するだけなので難しくはないだろう。ここで対数ガンマ関数の実装を行うために式を整理する。\[ \begin{align*} \log \Gamma(z) & \simeq \log \frac{\Gamma(z+1;g=5,N=6)}{z} \\ & = - \left(z+5.5\right) + \left(z+0.5\right) \log \left(z + 5.5 \right) \\ & \,\,\, + \log\left\{\frac{\sqrt{2\pi}}{z}\left(c_0 + \sum_{k=1}^6\frac{c_k}{z+k}\right)\right\} \end{align*} \] ここで係数 \(c_k\) は下記のコード上に記した定数 COF である。その他の係数を試したい場合は Coefficients for the Lanczos Approximation to the Gamma Function などを参照。

val COF = Seq(
  1.000000000190015, 76.18009172947146, -86.50532032941677,
  24.01409824083091, -1.231739572450155, 0.1208650973866179e-2,
  -0.5395239384953e-5
)
val SQRT2PI = 2.5066282746310005

def gammaln(x:Double):Double = {
  val x55 = x + 5.5
  var a = COF(0)
  for(k <- 1 until COF.length){
    a += COF(k) / (x + k)
  }
  - x55 + (x + 0.5) * math.log(x55) + math.log(SQRT2PI * a / x)
}

定数 \(\sqrt{2\pi}\) を数値リテラルで記述しているのは精度のためである。コード上で math.sqrt(2*math.Pi) を計算すると誤差のため最後の桁が 2 となり、結果として gammaln(1) := -4.440892098500626E-16 のような精度低下を引き起こす。数値演算ではこのように可読性を優先してコード上で定数を算出すると精度低下につながることがあるので注意したい。

Scala の REPL での実行結果は以下の通り。

scala> gammaln(0)
res0: Double = Infinity

scala> gammaln(1)
res1: Double = 0.0

scala> gammaln(3)
res2: Double = 0.6931471805599443

scala> gammaln(1.5)
res3: Double = -0.12078223763524987

scala> gammaln(Long.MaxValue)
res4: Double = 3.9354535028702885E20

Long の最大値においても値が発散せず Wolfram|Alpha での計算結果 \(3.93545\times 10^{20}\) とも一致していることが分かる。

ガンマ関数の実装

ガンマ関数の実装は先の対数ガンマ関数の結果を単純に指数演算するだけである。

def gamma(x:Double):Double = math.exp(gammaln(x))

先述の通りこの関数 gamma() は \(171 \sim 172\) でオーバーフローを起こすことに注意したい。

scala> Stream.from(1).find(x => gamma(x).isInfinite).get
res0: Int = 172

scala> gamma(172)
res1: Double = Infinity

scala> gamma(171)
res2: Double = 7.257415616277517E306

計算結果を \(x \in (0,4]\) の範囲でプロットしたものが下記のグラフである。

階乗関数

階乗関数の実装

ガンマ関数は自然数 \(n\) に対して \(\Gamma(n+1)=n!\) が成り立つ。やはり階乗も倍精度浮動小数点の精度限界から \(n=172-1\) でオーバーフローを起こすが、アプリケーションが \(0 \leq n \leq 171-1\) の範囲でしか利用しないことが分かっているのであれば通常の階乗関数として利用するのもよいだろう。

ここで倍精度浮動小数点で結果がオーバーフローしない階乗はたかだか 171 個しかないことを考えると、あらかじめ有効な階乗をすべて計算しメモリ上にキャッシュしておいても問題はないだろう。さらに、その算出に近似のガンマ関数を使用せず、多倍長精度整数 BigInt (Java では BigInteger) で直接的に正確な階乗を計算しても性能にはほぼ影響がない。以下はその実装である。

// BigInt で階乗を計算するストリーム (単純に for 文でも良い)
def factorials(n:Int = 0, pf:BigInt = 1):Stream[BigInt] = {
  def _factl(n:Int, pf:BigInt):Stream[BigInt] = {
    val f:BigInt = if(n <= 1) 1 else pf * n
    f #:: _factl(n+1, f)
  }
  _factl(n, pf)
}

// 配列は不変ではないが添え字参照が高速なため
lazy val PMLPF = factorials().takeWhile(! _.toDouble.isInfinity).toArray
lazy val PF = PMLPF.map(_.toDouble)

// 倍精度浮動小数点版 (事前計算分を超えたら無限大を返す)
def factorial(n:Int):Double = if(n < PF.length) PF(n) else Double.PositiveInfinity

// 多倍長精度整数版 (事前計算分を超えたら続きから逐次計算する)
def mlpfactorial(n:Int):BigInt = if(n < PMLPF.length) PMLPF(n) else {
  factorials(PMLPF.length, PMLPF.last).drop(n - PMLPF.length).head
}

この関数 factorial() は事前に計算済みの値を配列から参照するだけなので \(O(1)\) と高速である。n=172-1 以上で倍精度浮動小数点の限界を超えオーバーフローする。mlpfactorial() は多倍長精度整数版でメモリと CPU の許す限り正確な階乗を計算する。

scala> (factorial(0), factorial(1), factorial(2), factorial(3))
res27: (Double, Double, Double, Double) = (1.0, 1.0, 2.0, 6.0)

scala> (factorial(170), factorial(171))
res28: (Double, Double) = (7.257415615307999E306, Infinity)

scala> (factorial(100), gamma(100+1))
res29: (Double, Double) = (9.332621544394415E157, 9.332621545372994E157)

scala> math.abs(factorial(100) - gamma(100+1))/factorial(100)
res30: Double = 1.0485576290199364E-10

scala> mlpfactorial(171)
res31: BigInt = 12410180702176678234248405241031039926166055775016931...

以下は \(n=0\sim 170\) 範囲での近似ガンマ関数 gamma(n+1) と多倍長精度整数を使用した正確な階乗 factorial(n) との誤差を示すグラフである。

最も大きな部分で概ね \(1.3\times 10^{-10}\) 程度の誤差を持つことが分かる。Double の有効数字が 15 桁であることを考えれば多倍長精度整数を使って算出する方が有意な精度向上をもたらすといえる。

対数階乗関数の実装

階乗もガンマ関数と同様にしばしば分数の形で現れることから対数階乗関数としての実装が必要だろう。

上に実装した factorial() は精度が良く高速であるため、その結果がオーバーフローしない範囲であれば対数を計算する方法が良いだろう。それより大きな値については gammaln() の結果を使用する。また先に示した通り gammaln() は Long の最大値でもオーバーフローしないことから引数は Long 型とする。

ここで Long のまま n+1 すると Long.MaxValue + 1 → Long.MinValue となり不正な結果をもたらすことにも注意。なお Double で正確に表現できる整数が 53bit までであるため 53bit 値より大きな n は桁落ちにより計算的等価 n! == (n±α)! となる可能性がある。

def factorialln(n:Long):Double = if(n < PF.length){
  math.log(factorial(n.toInt))
} else {
  gammaln(n.toDouble+1)
}

factorialln() による \(\log n! = \log \Gamma(n+1)\) のグラフを以下に示す。\(n\) の増加と共にほぼ線形に上昇していることが分かる。