たった160行で書く誤差関数の実数多倍長ルーチン

はじめに

　数学ライブラリでは初等関数の実装は標準的で、アルゴリズムも簡単なものばかりである。しかし特殊関数となると、実装は処理系依存であることが多い。
　誤差関数やガンマ関数はそんな話題にありながら数学ライブラリとしてもよく使われる関数である。
　技術ノートを共有するコミュニティでは、とても一般的ながら、同時にあまりにも特殊であるため提供者のほうが稀であると言える。
　この記事では、すでにアンダーフローに達しているなど特殊な収束条件を見極めるだけで、参考文献[1]を頼りに実装すると、たった160行程度で書けてしまった誤差関数の多倍長ルーチンを提供する。

目次

• 定義
  • 誤差関数
  • 相補誤差関数
• 実装
  • 誤差関数
  • 相補誤差関数
• まとめ
• 参考文献

定義

誤差関数

　誤差関数$\textrm{erf}(x)$は、次式で定義される。

$\begin{align}
\textrm{erf}(x) & = & \displaystyle\frac{1}{\sqrt{\pi}}\int_{-x}^{x}e^{-t^2}dt, \
& = & \displaystyle\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^2}dt.
\end{align}$

　ただし、$\textrm{erf}(0)=0$、$\textrm{erf}(\infty)=1$、 $\textrm{erf}(-\infty)=-1$ である。

相補誤差関数

　また相補誤差関数$\textrm{erfc}(x)$は、次式で定義される。

$$\begin{align}
\textrm{erfc}(x) & = 1 - \textrm{erf} (x) & = & \displaystyle{1-}\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^2}dt, \
& & = & \displaystyle\frac{2}{\sqrt{\pi}}\int_{x}^{\infty}e^{-t^2}dt.
\end{align}$$

　ただし、$\textrm{erfc}(0)=1$、$\textrm{erfc}(\infty)=0$、$\textrm{erfc}(-\infty)=2$である。

実装

誤差関数

　参考文献[1]では、引数$| x |$が小さいときに$x=0$の級数展開を、やや大きいときに$x=\infty$の級数展開を、その範囲を超えるならば連分数展開を使うことを示している。

$\begin{align}
\textrm{erf}(x) & = & \displaystyle\frac{2}{\sqrt{\pi}}\sum_{n=0}^\infty\frac{(-1)^n x^{2n+1}}{n!(2n+1)} \
& = & \displaystyle\frac{e^{-x^2}}{\sqrt{\pi}}\sum_{n=0}^\infty\frac{n!\ (2x)^{2n+1}}{(2n+1)!} \
& = & \displaystyle{1+}\frac{-e^{-x^2}/\sqrt{\pi}}{x+\displaystyle\frac{1}{2x+\displaystyle\frac{2}{x+\displaystyle\frac{3}{2x+\displaystyle\frac{4}{x+\cdots}}}}}
\end{align}$

　ここで$n!$は階乗関数で、2つ目の級数展開では係数として2つあることに注意する。

　さて、実際の実装だが、収束条件は公式のBigDecimalにあるBigMathのモジュール関数と同様の定義で持っていける。$|x|=0$の級数展開では、以下のように定義できた。

def erf_ser_zero(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  n = BigDecimal.double_fig + prec

  zero = BigDecimal(0)
  one = BigDecimal(1)
  two = BigDecimal(2)

  x2 = x * x
  f = one
  r = one
  i = one
  d = x
  y = zero

  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    d = x.div(f * r, m)
    y = y + d
    x = -x.mult(x2, n)
    f *= i
    r += two
    i += one
  end
  (2 / BigMath.PI(prec).sqrt(prec)).mult(y.round(prec), prec)
end

　変数nとして任意精度にBigDecimal#double_fig(16くらい)を足し上げた精度数に注目しよう。実際に例えば16くらいの任意精度として計算しても、ちょうど足りなかったり少し切り上げるのが速すぎたりするため、このように余分な精度の数を使っている。また、収束を早めるべく、イテレータ変数としてmを別途用意し、切り上げる。

　$|x|=\infty$の級数展開では、以下のように定義できた。

def erf_ser_inf(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  n = BigDecimal.double_fig + prec

  s = BigMath.exp(-x * x, prec)
  return x.negative? ? BigDecimal(-1) : BigDecimal(1) if s.zero?
  s = s.div(BigMath.PI(prec).sqrt(prec), prec)

  zero = BigDecimal(0)
  one = BigDecimal(1)
  two = BigDecimal(2)

  x = two * x
  x2 = x * x
  f = one
  r = one
  i = zero
  i2 = one
  d = x
  y = zero

  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    d = (f * x).div(r, m)
    y = y + d
    x = x.mult(x2, n)
    f *= (i += one)
    r *= (i2 += one) * (i2 += one)
  end
  s.mult(y, prec)
end

　定義上、$|x|$は無限から数えるため、計算式で出てくる指数関数$e^{-x^2}$の解はどうやっても$0$に近くなる。この性質を利用し、係数$-e^{-x^2}/\sqrt{\pi}$を先に計算させ、この時点で$0$ならば微分せず解を送出する。

　次に連分数展開だが、値を求めるのに

$b_0+\displaystyle\frac{c_1}{b_1+}\frac{c_2}{b_2+}\frac{c_3}{b_3+}\frac{c_4}{b_4+}\cdots\frac{c_n}{b_n}$

　の$n$を減らす(ascend)書き方は、収束条件を必要とする多倍長演算では$n$を加算させた(decend)ほうが適している。
　Rubyで$n$を減らす書き方なら、以下のように書く。

f = 0
while n > 0
  f = c[n] / (b[n] + f); n = n.pred
end
return f + b[0]

　ここでは奥村先生のアルゴリズム事典[2]に従いつつ、前回の値と今回の値が同値のときにイテレーションを抜けるというルーチンを組んでみた。

c = (分数の分子の初期値)
b = (分数の係数の初期値の次)
p1 = one;  q1 = zero;  p2 = (分数の係数の初期値);  q2 = one;
i = one
loop do
  prev = p2
  t = p1 * c + p2 * b;  p1 = p2;  p2 = t;
  t = q1 * c + q2 * b;  q1 = q2;  q2 = t;
  if q2.nonzero?
    p1 = p1.div(q2, n);  q1 = q1.div(q2, n);  p2 = p2.div(q2, n);  q2 = one;
  end
  i += one
  c = (分数の分子の次数)
  b = (分数の係数の次数) 
  break unless q2.nonzero? && prev != p2
end

　ここではp2が返却値である。
　なお、加算する場合はこのようなアルゴリズムなのでテンプレート化も便利そうではあるが、収束条件が異なるためケースバイケースであると言えよう。

　実際に定義すると以下のようになった。

def erf_cf(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  n = BigDecimal.double_fig + prec

  zero = BigDecimal(0)
  one = BigDecimal(1)

  c = -BigMath.exp(-x * x, prec)
  return x < 0 ? BigDecimal(-1) : BigDecimal(1) if c.round(prec).zero?
  c /= BigMath.PI(n).sqrt(n); b = x
  p1 = one;  q1 = zero;  p2 = one;  q2 = one;
  i = one
  loop do
    prev = p2
    t = p1 * c + p2 * b;  p1 = p2;  p2 = t;
    t = q1 * c + q2 * b;  q1 = q2;  q2 = t;
    if q2.nonzero?
      p1 = p1.div(q2, n);  q1 = q1.div(q2, n);  p2 = p2.div(q2, n);  q2 = one;
    end
    i += one
    c = i - one
    b = i % 2 == 1 ? x : 2 * x;
    break unless q2.nonzero? && prev != p2
  end
  (x < 0 ? -p2 : p2).round(prec)
end

　これらを参考文献に従って条件分岐する。すると以下のようなプログラムになる。

def erf_branch(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  return BigDecimal::NAN if x != x
  return BigDecimal(x < 0 ? -1 : 1) if x.infinite?
  absx = x.abs
  branch_condition_a = (BigDecimal(prec + 21) * BigMath.log(10, prec)).sqrt(prec).div(5, prec)
  if absx <= branch_condition_a
    y = erf_ser_zero(absx, prec)
    y = -y if x.negative?
    return y
  end
  branch_condition_b = (((prec + 6) * BigMath.log(10, prec)) ** 0.6r).div(4.1r, prec)
  if absx <= branch_condition_b
    y = erf_ser_inf(absx, prec)
    y = -y if x.negative?
    return y
  else
    y = erf_cf(absx, prec)
    y = -y if x.negative?
    return y
  end
end

　branch_condition_*が分岐する条件となっている。参考文献では$\log(b)$として$b$は処理系の基数であると記されているため、値は10としている。

相補誤差関数

　相補誤差関数では$x$が正として大きい場合に連分数展開し、そうでない場合は1-erf(x)を計算する。

def erfc_branch(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  return BigDecimal::NAN if x != x
  return BigDecimal(x < 0 ? 2 : 0) if x.infinite?
  branch_condition_a = (((prec + 6) * BigMath.log(10, prec)) ** 0.7).div(9, prec)
  if x <= branch_condition_a
    y = 1 - erf_branch(x, prec)
    return y;
  else
    y = erfc_cf(x, prec)
    return y;
  end
end

　相補誤差関数の連分数展開は初期値が異なるだけで、誤差関数のアルゴリズムと同じである。

$$\textrm{erfc}(x) = \displaystyle\frac{e^{-x^2}/\sqrt{\pi}}{x+\displaystyle\frac{1}{2x+\displaystyle\frac{2}{x+\displaystyle\frac{3}{2x+\displaystyle\frac{4}{x+\cdots}}}}}
$$

　以下のように定義できた。

def erfc_cf(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  n = BigDecimal.double_fig + prec

  zero = BigDecimal(0)
  one = BigDecimal(1)

  c = BigMath.exp(-x * x, prec) / BigMath.PI(prec).sqrt(prec)
  return BigDecimal(0) if c.round(prec).zero?
  b = x
  p1 = one;  q1 = zero;  p2 = zero;  q2 = one;
  i = one
  loop do
    prev = p2
    t = p1 * c + p2 * b;  p1 = p2;  p2 = t;
    t = q1 * c + q2 * b;  q1 = q2;  q2 = t;
    if q2.nonzero?
      p1 = p1.div(q2, n);  q1 = q1.div(q2, n);  p2 = p2.div(q2, n);  q2 = one;
    end
    i += one
    c = i - one
    b = i % 2 == 1 ? x : 2 * x;
    break unless q2.nonzero? && prev != p2
  end
  p2.round(prec)
end

まとめ

　面白いほど実装がシンプルにまとまったのがお分かりいただけただろうか。恐らくガンマ関数の実装もそこまで難しくはないだろう。
　テストして使ってみると、十進ならではの振る舞いを確認できる。例えば$x$が$0$より離れるに付き、指定した精度分は位数を引き合わせてしまう。

erfc_branch(BigDecimal('5.0'), 20)
#=> 0.153745979e-11 # 足りない
erfc_branch(BigDecimal('5.0'), 40)
#=> 0.15374597944280348501883434854e-11 # 精度を上げると増える
Math.erfc(5)
#=> 1.5374597944280351e-12

　さて、これまでのルーチンをまとめたのが以下である。

# frozen_string_literal: true

require 'bigdecimal/math'

def erf_branch(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  return BigDecimal::NAN if x != x
  return BigDecimal(x < 0 ? -1 : 1) if x.infinite?
  absx = x.abs
  branch_condition_a = (BigDecimal(prec + 21) * BigMath.log(10, prec)).sqrt(prec).div(5, prec)
  if absx <= branch_condition_a
    y = erf_ser_zero(absx, prec)
    y = -y if x.negative?
    return y
  end
  branch_condition_b = (((prec + 6) * BigMath.log(10, prec)) ** 0.6r).div(4.1r, prec)
  if absx <= branch_condition_b
    y = erf_ser_inf(absx, prec)
    y = -y if x.negative?
    return y
  else
    y = erf_cf(absx, prec)
    y = -y if x.negative?
    return y
  end
end

def erf_ser_zero(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  n = BigDecimal.double_fig + prec

  zero = BigDecimal(0)
  one = BigDecimal(1)
  two = BigDecimal(2)

  x2 = x * x
  f = one
  r = one
  i = one
  d = x
  y = zero

  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    d = x.div(f * r, m)
    y = y + d
    x = -x.mult(x2, n)
    f *= i
    r += two
    i += one
  end
  (2 / BigMath.PI(prec).sqrt(prec)).mult(y.round(prec), prec)
end

def erf_ser_inf(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  n = BigDecimal.double_fig + prec

  s = BigMath.exp(-x * x, prec)
  return x.negative? ? BigDecimal(-1) : BigDecimal(1) if s.zero?
  s = s.div(BigMath.PI(prec).sqrt(prec), prec)

  zero = BigDecimal(0)
  one = BigDecimal(1)
  two = BigDecimal(2)

  x = two * x
  x2 = x * x
  f = one
  r = one
  i = zero
  i2 = one
  d = x
  y = zero

  while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    d = (f * x).div(r, m)
    y = y + d
    x = x.mult(x2, n)
    f *= (i += one)
    r *= (i2 += one) * (i2 += one)
  end
  s.mult(y, prec)
end

def erf_cf(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  n = BigDecimal.double_fig + prec

  zero = BigDecimal(0)
  one = BigDecimal(1)

  c = -BigMath.exp(-x * x, prec)
  return x < 0 ? BigDecimal(-1) : BigDecimal(1) if c.round(prec).zero?
  c /= BigMath.PI(n).sqrt(n); b = x
  p1 = one;  q1 = zero;  p2 = one;  q2 = one;
  i = one
  loop do
    prev = p2
    t = p1 * c + p2 * b;  p1 = p2;  p2 = t;
    t = q1 * c + q2 * b;  q1 = q2;  q2 = t;
    if q2.nonzero?
      p1 = p1.div(q2, n);  q1 = q1.div(q2, n);  p2 = p2.div(q2, n);  q2 = one;
    end
    i += one
    c = i - one
    b = i % 2 == 1 ? x : 2 * x;
    break unless q2.nonzero? && prev != p2
  end
  (x < 0 ? -p2 : p2).round(prec)
end

def erfc_branch(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  return BigDecimal::NAN if x != x
  return BigDecimal(x < 0 ? 2 : 0) if x.infinite?
  branch_condition_a = (((prec + 6) * BigMath.log(10, prec)) ** 0.7).div(9, prec)
  if x <= branch_condition_a
    y = 1 - erf_branch(x, prec)
    return y;
  else
    y = erfc_cf(x, prec)
    return y;
  end
end

def erfc_cf(x, prec)
  raise TypeError, "precision must be in Integer" unless prec.class == Integer
  raise RangeError, "Zero or negative precision" if prec <= 0
  n = BigDecimal.double_fig + prec

  zero = BigDecimal(0)
  one = BigDecimal(1)

  c = BigMath.exp(-x * x, prec) / BigMath.PI(prec).sqrt(prec)
  return BigDecimal(0) if c.round(prec).zero?
  b = x
  p1 = one;  q1 = zero;  p2 = zero;  q2 = one;
  i = one
  loop do
    prev = p2
    t = p1 * c + p2 * b;  p1 = p2;  p2 = t;
    t = q1 * c + q2 * b;  q1 = q2;  q2 = t;
    if q2.nonzero?
      p1 = p1.div(q2, n);  q1 = q1.div(q2, n);  p2 = p2.div(q2, n);  q2 = one;
    end
    i += one
    c = i - one
    b = i % 2 == 1 ? x : 2 * x;
    break unless q2.nonzero? && prev != p2
  end
  p2.round(prec)
end

参考文献

Multiple-Precision Exponential Integral and Related Functions - David M. Smith (Loyola Marymount University)
C言語による標準アルゴリズム事典 - 奥村晴彦