任意精度のガンマ関数 BigMath.gamma(x, prec) / BigMath.lgamma(x, prec) の計算法についてのまとめ。
このドキュメントは単体で読めるように書いてある(コードや diff を参照しなくても完結する)。
記法: p = 要求精度(10進桁数)。乗算コストは断りのない限り準線形 M(n) = O(n log n) モデルで数え、 大きい数×小さい数の積は M(m, n) = (n/m)·M(m) = O(n·log m)(小さい方のサイズでブロック化)と数える。
ガンマ関数の任意精度計算の難しさは、x の種類によってコスト構造が全く違うこと:
| x の種類 | 例 | 効く手法 |
|---|---|---|
| 整数 | gamma(100000) | 階乗そのもの(積) |
| 桁数の少ない実数 | gamma(1.25) | 級数/補間 + binary splitting |
| フル桁の実数 | gamma(√2) を10万桁 | 級数/補間の全精度演算 ~p 回が不可避 |
| 巨大な引数 | gamma(10¹⁷), lgamma(1e400) | 漸近展開 or 何らかの引数削減 |
x < 0.5 は反射公式 Γ(z)Γ(1−z) = π/sin(πz) で正側に移せるので、以下 x ≥ 0.5 を考える。
精度 p 桁には a ≈ p/log₁₀(2π) ≈ 1.25p 項が必要。
- 全精度の除算(と項ごとの sqrt)が ~1.25p 回 → O(p²·log p)
- 係数 c_k に平方根 (a−k)^(k−1/2) が入るため、項比が有理数にならず binary splitting が効かない。 桁数の少ない x でも p² から逃げられない
- コストが x に対して平坦: gamma(10¹⁷) でも安くならない(旧実装で 71s @1万桁)
- 利点: x > 0 全域で収束する単一の式、実装が短い
実測(旧実装): gamma(1.25, 10⁴) = 49s、gamma(1/3 フル桁, 10⁴) = 56s、gamma(10¹⁷, 10⁴) = 71s。
なお Spouge の式は「Γ をスターリング風にスケールした有理型関数の、留数による部分分数展開の打ち切り」 であり、後述の新実装と同じ「スケールして補間する」哲学の一種と解釈できる (スケーリング関数を (x+a−1)!·e^{−x−a}/(x+a)^{x+1/2} に取り、無限遠での収束条件を課すと Spouge が再現される)。
x が小さければ rising factorial x(x+1)…(x+n−1) で x ≈ 0.4p までシフトしてから適用。 全精度演算は級数 ~0.4p 項 + シフト ~0.4p 回で、演算回数は最少級。業界標準。
弱点は Bernoulli 数の生成コスト:
- 素朴な漸化式は N 個で O(N²) 回の p 桁加算。full 桁の x では N ≈ 0.2p 必要なので O(p³) 級になり全く使えない
- 高速生成(ζ 経由、冪級数反転 + Kronecker 代入)はそれ自体がライブラリ級の仕事
- 実用上はキャッシュ前提: mpmath は gamma(1.25) 2万桁が初回 226s、キャッシュ後 3.7s
さらに原理的な限界がある(§6.1「Bernoulli データ床」): full 桁の x では、 必要な B₂ₖ 列の総データ量が Θ(p²/log p) ビットあり、係数列を実体化する限り読むだけで p²/log p。
r ≈ p·ln10 に取ると裾 Γ(a,r) ≈ r^{a−1}e^{−r} < 10^{−p}。級数は全項が正(桁落ちが存在しない)で、 必要項数は (1+√2)·p·ln10 ≈ 5.6p 項。項比 r/(a+n) は有理数なので binary splitting / バッチ化が効く。
- 利点: 誤差解析が厳密かつ自明(裾の明示上界 + 正項級数の打ち切り)。Bernoulli 不要
- 弱点: x ≫ p では項数が √(x·p) に発散するため、巨大 x には別の手法が結局必要。 また通常帯域でも新実装より全精度演算が 1.5〜2 倍多い
整数 x はコスト ∝ x·polylog で n! を厳密計算して丸めるのが、x が小さいうちは最速 (GMP の積木は極めて速い)。x が大きくなると Stirling に切り替える。 Arb は dps=10⁵ でこの切り替え境界が x ≈ 2×10⁷ にあり、境界直上 + Bernoulli キャッシュが 冷えている初回はどちらの戦略も最悪になる「谷」が実測で観測できる(§7)。
Γ(1/3) のような固定有理数は hypergeometric 級数の binary splitting で準線形時間で計算できる (CLN の手法)。一般の x には適用できない。
Γ(x) を直接等間隔ノードで多項式補間することは、極と Runge 現象のため破綻する。代わりに
を補間する。f は整関数(極なし)で、x = b を中心とするほぼ対称なベル曲線 (スターリング近似より f(b+y) ≈ f(b)·e^{−y²/2b})になり、等間隔整数ノード x_i = b−l, …, b+l での局所補間が高精度に効く。x! = b^x / f(x) で復元する。
スケーリングに b^x を選んだことが全てを決める。f(n+1)/f(n) = b/(n+1) が 「小さい有理数」になるため:
- barycentric 形式の係数列 w_i·f(x_i) が有理数の漸化式で生成できる (比は −b(2l−j+1)/(j·(b−l+j)))→ binary splitting が効く(Spouge には √ があって効かない)
- 同じ理由で、full 桁 x にはバッチ多項式評価(BSGS)が効く
- ノード値が階乗そのものなので、階乗の高速計算と相互再帰できる(→ 倍数公式の実用化)
candidates としては x!/b^x, x!/x^x, x!/x^b なども考えられるが、 比が有理数になり、かつ補間点数が最少になるのは b^x/x! である。
barycentric(第1形式):
誤差はガウス曲線モデルから E ~ (l/2eb)^l と見積もられ、E ≤ 10^{−p} より l ≈ p / log₁₀(2eb/l)(不動点反復2回で解く。初期値 l = p は根の上側から単調収束するので安全側)。
x < 2p のときは x を 2p 以上へシフトしてから補間する(x! = x(x−1)…·(x−shift)! で戻す)。 ノードの正値性の理論的下限は b ≈ 1.36p だが、b = 2p が経験的に総コスト最適。 b = 2p のとき l ≈ p、シフト積 ≈ 2p 項で、全体は「~4p 個の1次因子の処理」に帰着する。
BSM(短桁 x、準線形): x を 10^(小数部桁数) 倍して全体を整数化し、 [部分和, 乗率, 分母] の3つ組を整数の binary splitting で畳む。 中間整数は prec·log₂10 + 64 bit で下位ビットを捨てる(bit-drop)。 捨てた 2^s は後述の exp2 チャネルで無料回収。O(p·log³p)。
BSGS(フル桁 x): log₂(p) 個ずつのバッチで Π(x−k) を整数係数多項式として展開し、 事前計算した x の冪で評価する。合成除算(Ruffini)で Σc_k/(x−k) のバッチも同じ多項式評価に落とす。 バッチを √l でなく log p に取るのは係数サイズ (b+l)^batch の爆発を防ぐため。 O(p²·log log p)(unbalanced 積をブロック化して数えるモデル。実装は小さい係数を schoolbook で掛けるので文字通りには log 一つ悪いが、word パッキングの定数 ~1/word² が 効いて実測は p²。クロスオーバーは 10 の数十〜百乗桁の彼方)。
分岐条件は「BSM の全精度演算換算 l·n_sig/p 個」対「BSGS の l·p ビット」の比較で n_sig·log₂(p) > p。
補間の再構成には (b−l)! が要る。x が巨大だとこれが巨大整数の階乗になる。そこで Legendre 倍数公式
で半分に落とす。片方は整数の階乗(再帰)、もう片方は半整数の階乗で、 これがまさに上の補間の BSM(x = k+0.5、シフト不要、指数が半整数なので b 冪は整数冪×√b)で 準線形時間で求まる。
通常この公式は無意味である — (n/2±0.5)! を求めるにはガンマの近似計算が要り、 それなら n! を直接近似した方が速い(Stirling は x が大きいほど収束が良い)からだ。 半整数階乗が準線形で出るこの方法で初めて実用になる。私の知る限り、 倍数公式を階乗計算の主経路に使っている既存ライブラリはない。
再帰の帳簿は n! = base·small!·(large!)²·2^n/√π の形(large! が再帰で2乗される)を [base_power_part, factorial_power_part, exp2, exp_sqrtpi] として持ち、 2^i 乗は最後にまとめて計算する。gamma(10¹⁷, 10⁴) は log₂(10¹⁷/4p) ≈ 41 段で 1.0s。
lgamma の定義域は指数が無制限(x = 1e400 など)で、doubling のコストは ∝ log x で増える一方、 Stirling は必要項数 N ≈ p·ln10/(2·ln 2πx) が x が大きいほど減る。必ず交差するので、 遠方の Stirling は不可避。切り替え条件を log₁₀x > √p/6 に置くと N < 3√p が保証され、 素朴な O(N²) の Bernoulli 漸化式が漸近的に無害になる(= 高速 Bernoulli 生成は永遠に不要)。
gamma 側は結果が overflow する x ≲ 4×10¹⁷ で定義域が切れるため doubling だけで足り、 gamma は全定義域で近似誤差モデルが1つになる。
- bit-drop の予算: 保持ビット = (10·prec + 192)/3 = prec·log₂10 強 + 64bit。 マージ1段 ~1bit の損失 × 木の深さ log₂(2l) を 64bit が吸収
- exp2 チャネル: bit-drop で捨てた 2^s は、補間の返り値 [base, large!, small!, exp2] の exp2 に負号で乗せ、doubling の既存の 2^exp2 会計に合流させる(冪計算を払わない)。 階乗パラメータ側の fact range 積だけは厳密に計算する(捨てた 2^s が 2^index 乗されて 指数レンジを壊すため)
- BSGS の batch_prod 相殺: x がノード直近のとき batch_prod は激しく桁落ちするが、 同じ計算値が sum の分母と prod の因子の両方に使われるため最終結果では厳密に相殺する。 batch_prod を別々に計算するリファクタはノード直近入力だけを壊す
- factorial_power_part の単調性: 深い再帰ほど b が小さく l が大きいので非減少。 増分更新はこれに依存
- 検証は (i) 恒等式(漸化式・反射・倍数公式・lgamma=log∘gamma)のクロスチェック、 (ii) mpmath との高精度突き合わせ、(iii) 各経路の prec vs prec+300 の自己整合、で行う
- Bernoulli 数が本質パスから消える。ステートレスでキャッシュ不要。 ワンショットの呼び出しで、キャッシュ済み mpmath より短桁 x で速い
- 近似誤差モデルが gamma 全域で1つ。精度が x に対して一様で、 手法の継ぎ目(Arb の x≈2×10⁷ の谷のような)が存在しない。テスト・保守面積が小さい
- 項比が有理数という一点から、短桁の準線形(BSM)・フル桁の p²(BSGS)・ 巨大 x の倍数公式が全部導かれる。1つの式 + 評価戦略の分岐、という構造
- 指数(オーダー)は標準手法と互角以上: Spouge の p²·log p に対し BSM は p·log³p、 BSGS は p²·loglog p(モデル値)
正直な弱点: C + GMP の Arb には定数倍で負ける(乗算プリミティブ 4〜10× + Ruby の 呼び出し/GC 税 3〜5×)。ただし指数は同じで、ワンショット比較では差が縮む。
full 桁の x を Stirling で p 桁計算するとき、B₂ₖ は項の大きさに応じて p_k ≈ p(1−k/N) 桁必要で、
つまり Bernoulli 係数列を float として実体化する手法は、係数を読むだけで Θ(p²/log p)。 生成をどれだけ速くしても超えられない。sub-quadratic への道は Bernoulli を捨てた側(補間・不完全ガンマ系 = 係数が有理漸化式で生成される側)にしかない。
BSGS の p² の由来は「バッチ係数を厳密整数で持つ」制約で batch = log p 止まりなこと。 係数を p 桁 float に落とし、batch = m = √n (n ≈ 4p) に取ると:
- バッチ多項式 D(z) = Π_{j<m}((x−j)−z), N(z) を係数 p 桁 float の積木で構築 — M(m·p)·log m
- z = 0, m, 2m, … の等差数列での多点評価。一般点の remainder tree は float では不安定だが、 等差数列なら Newton(下降階乗)基底変換 + 畳み込みで安定に評価できる。 誤差増幅は 2^O(m) × 値域比 → ガード桁 O(√p·log p) = o(p)
- 階乗型 prefactor は m 回の逐次乗法更新
合計 O(p^1.5·polylog p)。整数の世界での完全な先例が Pollard–Strassen (1976) の n! mod N アルゴリズム(p 桁 float ≈ mod 10^p の類推)。近代的整理は Chudnovsky–Chudnovsky / Bostan–Gaudry–Schost(いずれも exact/modular)。 float 版の安定性解析つき実装は私の知る限り世界に存在しない(新規性がある部分)。
よくある誤解: 「係数に x^m が現れて p^1.5 桁に爆発する」— これは厳密表現と float 表現の混同。 x の magnitude は ~4p なので x^m は指数 ~m·log₁₀(4p) の普通の p 桁 float であり、 問題は表現サイズではなく桁落ち量(上記の通り o(p) で抑えられる見込み)。
留意点:
- 中間形がない: 係数を float 化した瞬間、1点評価では全精度乗算に戻る。 多点評価まで一気に作って初めて利得が出る(オール・オア・ナッシング)
- 前提条件: 準線形の整数乗算(素の CRuby Integer は GMP なしだと Toom3 で本末転倒。 BigDecimal 自前 NTT に Kronecker で載せるのが正道)。メモリは係数列 p^1.5 桁
- 効く領域は「フル桁 x × 数万桁以上」のみ。クロスオーバー予想 10⁴〜10⁵ 桁、10⁵ 桁で ~10倍
- 成立条件: 誤差台帳の文書化・遅い経路を oracle にした敵対的テスト・kill switch の3点セット
なお Γ(x) は x について holonomic でない(線形微分方程式を満たさない)ため、 bit-burst 系の準線形評価は適用できない。p^1.5 が現実的なフロンティアと思われる。
同一マシン (Apple Silicon)。「旧」= Spouge 実装、「新」= 本方式(整数化 BSM 統合後, 2026-07)。
| 計算 | 旧 (Spouge) | 新 | 経路 |
|---|---|---|---|
| gamma(1.25, 10⁴) | 49s | 0.24s | BSM |
| gamma(1.125, 10⁵) | ~6000s (推定) | 2.1s | BSM |
| gamma(1/3 フル桁, 10⁴) | 56s | 8.2s | BSGS |
| gamma(1/3 フル桁, 10⁵) | ~7000s (推定) | 800s | BSGS |
| gamma(10¹⁷, 10⁴) | 71s | 1.0s | doubling |
| gamma(10¹⁷, 10⁵) | ~9000s (推定) | 13s | doubling |
| 計算 | 桁 | mpmath 初回 | mpmath キャッシュ後 | 新 |
|---|---|---|---|---|
| gamma(1.25) | 10⁴ | 26.9s | 0.69s | 0.24s |
| gamma(1.25) | 2×10⁴ | 226s | 3.7s | 0.57s |
| gamma(1/3) | 2×10⁴ | 226s | 3.8s | 33s |
短桁 x ではキャッシュ済み mpmath より速い(ステートレスのまま)。フル桁 x では キャッシュ済み Stirling に負ける(演算数の差 + 実装言語)。
| 計算 | Arb | 新 | 備考 |
|---|---|---|---|
| gamma(1.125) | 0.27s | 2.1s | 乗算プリミティブ差 (GMP は BigDecimal NTT の 4〜10×) + Ruby 税 |
| gamma(√2) 10⁴桁 | 0.1s | 7.8s | フル桁。Arb は rectangular splitting |
| gamma(2×10⁶) | 0.19s | 2.3s | Arb は厳密階乗 (コスト ∝ x) |
| gamma(2×10⁷) | 3.7s (初回) / 1.6s | 3.4s | Arb の戦略境界の谷 + cold Bernoulli では互角 |
| gamma(2×10⁸) | 1.3s | 3.9s | Arb は Stirling (x↑ で軽くなる)、新は doubling (∝ log x) |
指数は互角、定数は「基礎演算 4〜10×」×「Ruby の呼び出し固定費 + GC 3〜5×」の積で説明できる (プロファイル実測: doubling 経路は Integer#* (GMP) 48%、GC は整数化後 ~8%)。
- J.L. Spouge, "Computation of the gamma, digamma, and trigamma functions" (1994)
- R.P. Brent, "Unrestricted algorithms for elementary and special functions" (IFIP 1980) — 不完全ガンマ級数による Γ
- D.M. Smith, "Algorithm 814: Fortran 90 software for floating-point multiple precision arithmetic, gamma and related functions" (ACM TOMS, 2001)
- F. Johansson, "Arbitrary-precision computation of the gamma function" (2021) — 既存手法の包括的サーベイ
- F. Johansson, "Evaluating parametric holonomic sequences using rectangular splitting" (ISSAC 2014)
- J.M. Pollard (1974) / V. Strassen (1976) — n! mod N の √n·polylog 算法(§6.2 の原型)
- A. Bostan, P. Gaudry, É. Schost, "Linear recurrences with polynomial coefficients and application to integer factorization and Cartier–Manin operator" (2007)
- R.P. Brent, P. Zimmermann, "Modern Computer Arithmetic" — M(n) 表記・基数変換・特殊関数の教科書