はじめに

オイラーの定数

　数学上では$e$、$\pi$、$\imath$、$\log 2$などの重要な定数がある。また指数積分などで現れる定数も重要である。それがオイラーの定数$\gamma\approx0.5772156649015328606\ldots$である。
　発見者はその名から分かるとおりオイラーである。当初オイラーは$C$といった記号を用いていたが、ロレンツォ・マスケローニが$\gamma$を記号表記したためこれが知られ、相まってオイラー・マスケローニ定数とも呼ばれるようになった。

　数値計算でも利用頻度が高めだが、では自前で実装せよといえば、専門的見地の範ちゅうである。実際にQiitaでこれを投稿したのは、筆者くらいのようである。(2025年6月時点)
　プログラム例は載っていても、玄人っぽく紹介していて全く真値があっていないというサイトもあるが、筆者にはそんなことは滅多にない。

目次

本稿
- 定義
- 積分での表示
- ブレント・マクミラン型公式と数値計算

本稿

定義

　定義式は

$$\gamma=\displaystyle\lim_{n\to\infty}(\sum_{k=0}^n\frac{1}{k}-\log n)$$

　とする。この$\frac{1}{k}$は調和級数 harmonic series と呼ばれる。収束は非常に遅い発散級数である。調和級数の$\frac{1}{n}$までの部分和と$\log n$の差

$$a_n=1+\frac{1}{2}+\frac{1}{3}+\cdots+\frac{1}{n}-\log n$$

　を見てみると、

$$\begin{array}{rcl}
a_3 & = & 0.7347210\ldots \newline
a_{10} & = & 0.6263831\ldots \newline
a_{100} & = & 0.5822073\ldots \newline
a_{1000} & = & 0.5777155\ldots \newline
a_{1000000} & = & 0.5772161\ldots
\end{array}$$

　$a$が$1000000$の項でも僅か5桁しか正しく表示されていないのが分かる。
　数値計算でさえあまりにも膨大だが、オイラーは15桁まで求めている。どのようにしてオイラーの時代で15桁も求めたのかと、未だに話題になることがある。

　この定義式は調和級数を変数とした

$$\gamma=\displaystyle\lim_{n\to\infty}(H_n-\log n)$$

　とすることがある。また$H_n$を調和数 harmonic number と呼ぶことがある。これは収束の遅い調和級数を代替するためのものとして数学者たちの試行錯誤に趣がある。

積分での表示

　この定数はもともと先述するように数列の極限で定義された。これは積分によっても表される。一例では以下の表示式

$$\gamma=-\int_0^\infty e^{-t} \log t\ dt$$

　と書き表される。また$t=e^{-x}$とおけば、置換積分により

$$\log x=\log\log\frac{1}{t},\ dx=-\frac{dt}{e^{-x}}$$

　であり

$$\gamma=-\int_0^1 \log\log\frac{1}{t}\ dt$$

　を得る。

　この定数は様々な公式が発見されている。例えば以下のような表示である。

$$\gamma=\int_1^\infty(\frac{1}{[t]}-\frac{1}{t})\ dt$$
$$\gamma=1-\int_1^\infty\frac{t-[t]}{t^2}\ dt$$

　ここで$[x]$はガウス記号と呼ばれ、ある整数の最大を意味する。数値計算では床関数がこれに該当する。

ブレント・マクミラン型公式と数値計算

　ブレント・マクミラン型公式による数値計算を示す。ブレント・マクミラン型はベッセル関数を用いて高速に収束することで知られる。
　定義式は

$$\displaystyle\lim_{n\to\infty}[\sum_{k=0}^\infty\frac{(\frac{n^k}{k!})^2 H_k}{\sum_{k=0}^\infty(\frac{n^k}{k!})^2}-\ln n]$$

　で与えられる(1980)。当初、公式は以下のような式

$$\gamma=\frac{S_0(2n)-K_0(2n)}{I_0(2n)}-\ln(n),$$

　で与えられた。
　ここで $n$は整数として$n > 0$であり、自由に決めてよい。多倍長では精度になる。$K_0(x)$と$I_0(x)$はベッセル関数である。そして$S_0(x)$は

$$S_0(x)=\displaystyle\sum_{k=0}^\infty\frac{H_k}{(k!)^2}(\frac{x}{2})^{2k}$$

　とする。この公式は微分方程式・差分方程式・級数展開・漸近展開の4つのアルゴリズムが混在しており、定義式をそのまま使うのは難しい。ここで改良型にボールウェイン・ベイリー公式がある。
　公式では、$B_k=\frac{B_{k-1} n^2}{k^2}$、$A_k=\frac{1}{k}(\frac{A_{k-1}n^2}{k}+B_k)$を各々項としつつ、$U_k=U_{k-1}+A_k$、$V_k=V_{k-1}+B_k$として足し上げていく。開始点は$A_0=-\ln n$、$B_0=1$、$U_0=A_0$、$V_0=1$とする。収束したとするならばイテレーションを抜け$U/V$で解を求める(2003)。

　下記プログラム例は、上手く四則演算のメソッドを駆使し、素早くかつ正確に真値を導出するよう工夫している。ちなみにプログラムを別のものに書き換えると真値にならず、${10}^{-30}$あたりで情報落ちが生じる。

require 'bigdecimal/math'

def EulerGamma_borwein_bailey(prec)
  raise TypeError, "precision must be an Integer" unless prec.class == Integer
  raise ArgumentError, "Zero or negative precision" if prec <= 0
  zero = BigDecimal(0)
  one = BigDecimal(1)
  n = prec < BigDecimal.double_fig ? BigDecimal.double_fig : prec + 2
  nn = BigDecimal(n)
  a = u = -BigMath.log(n, n)
  b = v = one
  k = zero
  while a.nonzero? && ((m = n - (u.exponent - a.exponent).abs) > 0)
    m = BigDecimal.double_fig if m < BigDecimal.double_fig
    k = k + one
    b = (b * nn * nn).div(k * k, m)
    a = one.div(k, m) * ((a * nn * nn).div(k, m).add(b, m))
    u = u + a
    v = v + b
  end
  u.div(v, prec)
end

　実行結果を以下に。

EulerGamma_borwein_bailey(50)
#=> 0.57721566490153286060651209008240243104215933593992e0

　Pythonで真値を得たいなら、参考サイトを辿っていくとそのプログラムがある。

参考文献

美しい数学を描く$\pi$，$e$，とオイラーの定数$\gamma$ | 若原龍彦 | 講談社エディトリアル
A bound for the error term in the Brent-McMillan algorithm | Richard P. Brent and Fredrik Johansson

参考サイト

Euler-Mascheroni Constant - Wolfram Mathworld
What is the fastest/most efficient algorithm for estimating Euler's Constant γ
? - Mathematica Stack Exchange