Discussion memo of https://github.com/ruby/bigdecimal/pull/524 and https://github.com/tompng/bigdecimal/tree/gamma_multipoint_evaluation by Fable 5

gamma_algorithm.md

BigMath.gamma の計算法: 旧実装・新実装・既存手法・その先

任意精度のガンマ関数 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)(小さい方のサイズでブロック化)と数える。

1. 問題の構造

ガンマ関数の任意精度計算の難しさは、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 を考える。

2. 旧実装 (〜v4.1.x): Spouge の近似式

$$\Gamma(z+1) = (z+a)^{z+1/2} e^{-(z+a)} \left[\sqrt{2\pi} + \sum_{k=1}^{a-1} \frac{c_k}{z+k} + \varepsilon\right],\quad c_k = \frac{(-1)^{k-1}}{(k-1)!}(a-k)^{k-1/2}e^{a-k}$$

精度 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 が再現される)。

3. 標準的な既存手法(他ライブラリ)

3.1 Stirling 漸近展開 + 引数シフト + Bernoulli 数 (mpmath / MPFR / Arb)

$$\log\Gamma(x) \sim \left(x-\tfrac12\right)\log x - x + \tfrac12\log 2\pi + \sum_{k\ge1} \frac{B_{2k}}{2k(2k-1)x^{2k-1}}$$

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。

3.2 不完全ガンマ級数 (Brent 1978 / D.M. Smith, TOMS Algorithm 814)

$$\Gamma(a) = \gamma(a,r) + \Gamma(a,r),\qquad \gamma(a,r) = r^a e^{-r}\sum_{n\ge0}\frac{r^n}{a(a+1)\cdots(a+n)}$$

r ≈ p·ln10 に取ると裾 Γ(a,r) ≈ r^{a−1}e^{−r} < 10^{−p}。級数は全項が正(桁落ちが存在しない)で、 必要項数は ~2.72·p·ln10 ≈ 6.3p 項(裾指数は Poisson 型 γ−(1+γ)ln(1+γ) で決まる。 ガウス近似で数えると (1+√2)·p·ln10 ≈ 5.6p と ~13% 過少になるので注意 — 実測で発見)。 項比 r/(a+n) は有理数なので binary splitting / バッチ化が効く。

• 利点: 誤差解析が厳密かつ自明(裾の明示上界 + 正項級数の打ち切り)。Bernoulli 不要
• 弱点: x ≫ p では項数が √(x·p) に発散するため、巨大 x には別の手法が結局必要
• 補足: 多点評価機構(§6.2)の存在下では評価が逆転する — 項比の分子が定数なので √バッチ化との相性が良く、第2クライアントとして実装した結果、full-digit x ∈ [0.5, 3] では Lagrange 多点評価版より 1.4〜2.3× 速い(gamma_multipoint.md §9.3)

3.3 整数 x の厳密階乗 (Arb など)

整数 x はコスト ∝ x·polylog で n! を厳密計算して丸めるのが、x が小さいうちは最速 (GMP の積木は極めて速い)。x が大きくなると Stirling に切り替える。 Arb は dps=10⁵ でこの切り替え境界が x ≈ 2×10⁷ にあり、境界直上 + Bernoulli キャッシュが 冷えている初回はどちらの戦略も最悪になる「谷」が実測で観測できる(§7)。

3.4 有理数の特殊値

Γ(1/3) のような固定有理数は hypergeometric 級数の binary splitting で準線形時間で計算できる (CLN の手法)。一般の x には適用できない。

4. 新実装: b^x/x! のラグランジュ補間

4.1 核心のアイデア

Γ(x) を直接等間隔ノードで多項式補間することは、極と Runge 現象のため破綻する。代わりに

$$f(x) = \frac{b^x}{x!}$$

を補間する。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) が 「小さい有理数」になるため:

1. barycentric 形式の係数列 w_i·f(x_i) が有理数の漸化式で生成できる (比は −b(2l−j+1)/(j·(b−l+j)))→ binary splitting が効く(Spouge には √ があって効かない)
2. 同じ理由で、full 桁 x にはバッチ多項式評価(BSGS)が効く
3. ノード値が階乗そのものなので、階乗の高速計算と相互再帰できる(→ 倍数公式の実用化)

candidates としては x!/b^x, x!/x^x, x!/x^b なども考えられるが、 比が有理数になり、かつ補間点数が最少になるのは b^x/x! である。

4.2 補間の形と誤差

barycentric(第1形式):

$$f(x) \approx \omega(x)\sum_{i}\frac{w_i f(x_i)}{x-x_i},\qquad \omega(x)=\prod_i (x-x_i)$$

誤差はガウス曲線モデルから 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次因子の処理」に帰着する。

4.3 評価戦略(x の桁数で分岐)

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。

4.4 Factorial doubling(倍数公式の実用化)

補間の再構成には (b−l)! が要る。x が巨大だとこれが巨大整数の階乗になる。そこで Legendre 倍数公式

$$n! = \left(\frac{n}{2}\right)!\left(\frac{n-1}{2}\right)!\cdot\frac{2^n}{\sqrt\pi}$$

で半分に落とす。片方は整数の階乗(再帰)、もう片方は半整数の階乗で、 これがまさに上の補間の 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。

4.5 lgamma の遠方だけ Stirling

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つになる。

4.6 実装上の要点(保守のための invariant)

• 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 の自己整合、で行う

5. この方法のどこが嬉しいのか

1. Bernoulli 数が本質パスから消える。ステートレスでキャッシュ不要。 ワンショットの呼び出しで、キャッシュ済み mpmath より短桁 x で速い
2. 近似誤差モデルが gamma 全域で1つ。精度が x に対して一様で、 手法の継ぎ目(Arb の x≈2×10⁷ の谷のような)が存在しない。テスト・保守面積が小さい
3. 項比が有理数という一点から、短桁の準線形(BSM)・フル桁の p²(BSGS)・ 巨大 x の倍数公式が全部導かれる。1つの式 + 評価戦略の分岐、という構造
4. 指数(オーダー)は標準手法と互角以上: Spouge の p²·log p に対し BSM は p·log³p、 BSGS は p²·loglog p(モデル値)

正直な弱点: C + GMP の Arb には定数倍で負ける(乗算プリミティブ 4〜10× + Ruby の 呼び出し/GC 税 3〜5×)。ただし指数は同じで、ワンショット比較では差が縮む。

6. より速いオーダーの可能性 (future work)

6.1 Bernoulli データ床

full 桁の x を Stirling で p 桁計算するとき、B₂ₖ は項の大きさに応じて p_k ≈ p(1−k/N) 桁必要で、

$$\sum_k p_k \approx \frac{p\cdot N}{2} \approx \frac{p^2}{4\log_{10} p}$$

つまり Bernoulli 係数列を float として実体化する手法は、係数を読むだけで Θ(p²/log p)。 生成をどれだけ速くしても超えられない。sub-quadratic への道は Bernoulli を捨てた側(補間・不完全ガンマ系 = 係数が有理漸化式で生成される側)にしかない。

6.2 p^1.5·polylog への道(多点評価)

→ 実験実装済み (branch gamma_multipoint_evaluation)。実測は当初見積りより良く、 crossover ≈ 2000 桁、5万桁で既存 BSGS の 4.2×(BGS 値領域エンジン、O(p^1.5·log p))、 全テストで既存実装と厳密一致。BigMath.gamma へ配線済み。 設計と実測の詳細は gamma_multipoint.md。以下は元の理論見積り (安定性の見立ては実測で裏付けられた。予想と違ったのは、crossover が予想 10⁴〜10⁵ より 大幅に低かったことと、高速多点評価 (remainder tree) 自体は Horner の machine-word 定数に 25万桁付近まで勝てないこと)。

BSGS の p² の由来は「バッチ係数を厳密整数で持つ」制約で batch = log p 止まりなこと。 係数を p 桁 float に落とし、batch = m = √n (n ≈ 4p) に取ると:

1. バッチ多項式 D(z) = Π_{j<m}((x−j)−z), N(z) を係数 p 桁 float の積木で構築 — M(m·p)·log m
2. z = 0, m, 2m, … の等差数列での多点評価。一般点の remainder tree は float では不安定だが、 等差数列なら Newton(下降階乗)基底変換 + 畳み込みで安定に評価できる。 誤差増幅は 2^O(m) × 値域比 → ガード桁 O(√p·log p) = o(p)
3. 階乗型 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 で本末転倒。 実装は Ruby Integer/GMP に Kronecker で載せた — BigDecimal 自前 NTT は最大桁数上限があり この用途には不適)。メモリは係数列 p^1.5 桁
• 効く領域は「フル桁 x」のみ。クロスオーバーは実測 ≈ 2×10³ 桁(当初予想 10⁴〜10⁵ より良い)
• 成立条件: 誤差台帳の文書化・遅い経路を oracle にした敵対的テスト・kill switch の3点セット (guard 式は実測で確定済み、既存 BSGS が oracle、kill switch は Multipoint.enabled。 BigMath.gamma へ配線済み: GMP 検出 + prec ≥ 3000 の full-digit x で自動使用)

なお Γ(x) は x について holonomic でない(線形微分方程式を満たさない)ため、 bit-burst 系の準線形評価は適用できない。p^1.5 が現実的なフロンティアと思われる。

7. 実測まとめ

同一マシン (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 (gmpy backend) との比較

計算	桁	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 (FLINT, python-flint) との比較 (dps=10⁵)

計算	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%)。

多点評価版(§6.2 の実験実装)vs 既存 BSGS

フル桁 x = √2、結果は全ケース厳密一致。詳細は gamma_multipoint.md。

桁数	多点評価版	BSGS	比
2×10³	0.28s	0.28s	1.0× (crossover)
10⁴	4.1s	7.9s	1.9×
5×10⁴	49s	207s	4.2×

(5×10⁴ 以遠の数値は BGS 値領域エンジン。) フル桁 x の p² は実用域でも実際に破れる、が現時点の結論 (gamma(1/3) 10万桁の 800s も同率なら ~180s 見込み)。

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) 表記・基数変換・特殊関数の教科書

gamma_calculation.md

gamma 計算の議論メモ (2026-07, Claude とのセッション)

このブランチ (gamma_lagrange) の実装検証と、その先の話 (n^1.5 化など) を再開するためのメモ。 自己完結版(発表・merge 後の参照用)は gamma_algorithm.md — 旧実装/新実装/既存手法/p^1.5 の可能性/実測まとめを diff なしで読める形にしてある。こちらはセッション時系列の生メモ。

2026-07-06 更新: bsm_prod の Integer 化 + BSM の統合 (gamma_lagrange_n_plus_half 削除、10^f スケール一般化、exp2 チャネル) を実装済み。 gamma(1.125, 1e5): 4.1→2.1s / doubling 各点 5〜15% 改善 / BSGS 不変 / 精度は全経路で満額を再検証済み。 注意: rake compile が lib を tmp/*/stage/lib にコピーするため、lib 編集後の測定は -Ilib を先に置くか再 compile しないと旧コードを測る(実際に一度踏んだ)。

1. 実装の検証状況

検証済み (すべてパス):

• 数式の再導出: barycentric 係数漸化式 / 再構成式 / 倍数公式の帳簿 (large! が再帰で2乗になり 2^i 指数で積み上がる構造) / l 見積りの不動点反復 (初期値 l=prec から根の上側に単調収束するので過大評価側で安全)
• 恒等式クロスチェック: 漸化式 gamma(x+1)=x·gamma(x)、反射 gamma(1/4)gamma(3/4)=π√2、倍数公式、lgamma=log(gamma)、符号
• mpmath との 1000 桁突き合わせ: 1/3, 0.3, 1e8+1/3, 5+1e-40, -987.654321, 0.5001, lgamma(1e18), lgamma(1e400)
• パス別高精度セルフチェック: BSM gamma(1.25, 20000) 誤差指数 -20000 (満額)、BSGS gamma(1/3, 5000) -4999 (ulp 級、prec 990 でも -989 でドリフトなし)、doubling gamma(1e15+0.5, 5000) -5000

再現方法: mpmath (venv に pip install) で高精度参照値を出して突き合わせ。恒等式チェックは 「gamma(x+1, p) vs x·gamma(x, p+10)」「gamma(x, p) vs gamma(x, p+300) の truncate 比較」の形。 比較時の除算は / でなく .div(x, prec) を使うこと (/ は巨大 exponent で NoMemoryError — 有効桁から結果精度が決まる仕様由来、既知)。

2. 明文化した invariant (コメント済み、壊すと静かに死ぬ)

• batch_prod 同一値相殺 (BSGS): x がノード直近のとき batch_prod は激しく桁落ちするが、同じ計算値を sum の除算と prod の乗算の両方に使うので最終 prod·sum では厳密相殺。batch_prod を再計算・並べ替えするリファクタはノード直近入力だけを壊す。回帰テスト: gamma(5 + 1e-40) (test_bigmath.rb)
• factorial_power_part の単調性 (integer_factorial_parameter): 深い再帰レベルほど b が小さく l が大きいので非減少。fact_y の増分更新はこれに依存
• bit 落としマージン (n_plus_half): 保持ビット (prec·10+192)/3 = prec·log2(10) 超 + 64 bits。マージ1段 ~1 bit 損失 × 木の深さ log2(2l) を吸収
• 精度予算の分散ヒューリスティック一覧: EXTRA_PREC=16 / l の +10 / b = 2·prec (経験的最適、正の下限は ~1.36·prec) / internal_xn_prec = prec + log10(b+l)·batch / 上記 +64 bits。どれかを「改善」するとテスト帯域 (200–1200桁) では通るが大 prec で末尾が欠ける、が典型的な壊れ方

3. 計算量の整理

• 表記の流儀: bit complexity を M(n) でパラメータ化し、unbalanced 積は M(m,n) = (n/m)·M(m) = n·log(m) と数える (Brent–Zimmermann 流)。この前提で BSM O(p·log³p) / BSGS O(p²·loglog p)
• 実装の文字通りの漸近: 小オペランド (係数 ~log²p 桁) は BigDecimal の schoolbook パス (NTT 閾値 450 words 未満) なので O(p²·log²p)、ただし係数 1/(word桁数)² で実測は p²。schoolbook の方が実サイズでは速いための選択であり、主張は理論値でよい (全ライブラリ共通の慣習)
• Spouge との名目比較: 実装同士だと 1 log 悪い (p²log²p vs p²logp) が、クロスオーバーは p ~ 10^150 桁の彼方
• ベルヌーイ実体化ルートの下限: Stirling で full-digit x を p 桁計算するとき B_2k は p_k ≈ p(1−k/N) 桁必要で Σp_k ≈ p²/(4·log10 p)。係数列を「読むだけ」で Θ(p²/log p) — 生成をいくら高速化しても超えられない。full-digit x の sub-quadratic はベルヌーイを捨てた側 (補完系/不完全ガンマ系) にしか存在しない。この観察が n^1.5 化の動機の核

4. n^1.5·polylog への拡張 — 実験実装済み (2026-07-27, branch: gamma_multipoint_evaluation)

実装: lib/bigdecimal/math/gamma_multipoint.rb + 検証 gamma_mp_check.rb。 BSM の [sum_num, mult_num, den] triple を「バッチオフセット z の多項式」に持ち上げ、 固定小数点係数の Kronecker 詰め込み (Ruby Integer/GMP) で triple 木を構築、 z = 0, m, 2m, ... で評価 (現状は per-point Horner; 高速基底変換は未実装)。 E(z) は F2(z) から導出 (E = F2·(x−A−z)/(B·I)、B·I は厳密整数) して near-node 相殺 invariant を維持。 ノード数 n1 = m² (m 奇数で符号維持) に範囲を非対称拡張、係数比の一般形 ρ(t) = −b(n1−t)/(t(A+t))。

結果: prec 100〜50000 で既存実装と一致 (50k で agree=exact、他は ulp 級)。 crossover ≈ 1e4 (prec 2000: 0.6×, 1e4: 1.1×, 2e4: 1.5×, 5e4: 1.9×)。 プロファイル@1e4: Integer#* 43% / Kronecker pack・unpack ~25% / BigDecimal 14%。 guard = 4·m·(log2(n1)+4)+256 bits — 半減させると prec=2000 で 33 桁欠ける → ほぼ適正、過剰包装ではない。

タスク① (2026-07-28, commit 437d0445): remainder-tree 多点評価を実装。 falling-factorial モジュライ (係数 ~m·log2(m) bit の厳密整数) の subproduct 木 + 逆級数 (Newton 反復) を木にメモ化して f01/f2 で共有 + ワイド被除数は R = t^count mod M_root (小係数) でブロック結合してから1回だけ除算。 Kronecker スロット幅を maxA+maxB に修正 (2·max だった) — 小×大係数積が全体で高速化。 結果: Horner と全テストで厳密一致。ただし 実測 crossover は m ≈ 700 (≈ 25万桁) — それ未満では Horner の machine-word 定数が勝つため eval_mode デフォルトは :auto (m > 700 で fast)。5e4桁時点の内訳: 全体 124s / triple木 68s / eval(fast) 32s vs Horner 13s。 教訓: 唯一の p² 項である Horner は定数が小さく (~13%@5e4)、remainder tree の 定数 (ノード毎 2 conv + pack/unpack) は Ruby では重い。p² 項の除去は「≥ 25万桁で効く 漸近的保険」であり、実用域の本丸は triple 木 (4 mult/merge) と pack/unpack 定数。

タスク②③ (2026-07-28, commit 2a9ae04a): ② pack を borrow 折り込みの単一 hex-join に (正負分割 pack + 巨大減算を排除)、 unpack をスロット毎 borrow 伝播に (巨大バイアス定数の生成+加算を排除)。fast 1.18× / horner 1.11×。 ③ guard=0 での実損失を実測: loss = 2.9〜3.3 · m · bl(n1) bits (prec 300〜10000 で安定、 fast/horner 完全同一、near-node でも同一 = dyn range 支配・remainder tree の丸めは測定不能・相殺 invariant 有効)。 式を 4·m·bl(n1)+256 (マージン ~20%) に確定、根拠をコメント化。 到達点: gamma(√2) 1e4桁 6.1s (BSGS 8.5s の 1.4×) / 5e4桁 91.9s (BSGS 207s の 2.25×)、全て BSGS と厳密一致。

triple 木の削減 (2026-07-28, commit 2ea0de59): barycentric pair tree に置換。 鍵は分解 den_j = L_j·B_j·I_j, num_j = −b·L_{j−1}·G_j (L のみ W-bit、B/I/G は小整数係数) と、 num の L が den の L の添字ずらしであること。級数分子は F01 = Σ_j (Ω/L_j)·w_j (Ω = ΠL_i、 w_j は小係数のみ) の barycentric 形になり、ノード [Ω, Φ, BI, GX] (BI/GX は厳密整数多項式) の マージでワイド×ワイド積は Φ_A·Ω_C と Ω_A·Φ_C の2本だけ (旧: deg-3d triple の4本)。 j=0 項は根で付加し、その Ω·BI 部分がちょうど F2 なので再利用。 効果: 1e4 6.1→4.8s / 5e4 91.9→69.9s (BSGS 207s の 2.96×)、crossover ≈ 2000桁 まで低下。 全テストで BSGS と厳密一致を維持。

到達点まとめ (vs BSGS): 2000桁 1.0× / 5000桁 1.5× / 1e4 1.8× / 5e4 3.0×。 2026-07-29: BGS 値領域エンジン追加 (76c68b4c、詳細 gamma_multipoint.md §6)。 shift of evaluation values で木と評価が消え O(p^1.5·log p)。損失実測 1.4·S·bl (外挿増幅なし)、 engine=:auto (≥8000桁で :values)。5e4: 49.1s = BSGS 比 4.2×。 2026-07-31: 統合フェーズ完了 (d89d1fd3: shift 積も値領域、478774c7: 係数エンジン引退)。 最終形 = 値領域単一エンジン 337 行、全体一様に O(p^1.5·log p)、guard 則1つ。 5e4 = 45.3s (BSGS 比 4.6×)、実測局所指数 1.53。3000〜8000 桁帯は引退エンジン比 ~1.3× 譲歩。 同日: 汎用層切り出し (batch_value_tables を葉 block 化) + 第2クライアント incgamma (incgamma_mp_check.rb, 9d796205)。Lagrange 多点評価版より 1.4〜2.3× 速い (5e4: 32.5s vs 44.7s)。 損失則 0.26〜0.50·S·bl。副産物: 項数見積り 5.6p はガウス裾の誤用で、正しくは ~6.3p (loss 測定が検出)。full-digit の [0.5, O(p)] 帯は incgamma 系が優位という含意 — lane 再編は保留。 ④ dispatch 配線済み (2026-07-28, commit 5ec0788c): gamma_lagrange が full-digit x かつ prec ≥ Multipoint.min_prec (3000) で多点評価へ。GMP 判定は defined?(Integer::GMP_VERSION)、 kill switch は Multipoint.enabled。反射・lgamma・doubling は gamma_lagrange 経由で自動対応 (テスト用ラッパーの委譲が反射後の引数を BSGS に落としていた罠はこれで構造的に消滅)。

以下は実装前の設計分析 (経緯用にそのまま残す):

骨子: BSGS の batch 係数を厳密整数 (batch=log p 止まりの原因) から p桁 float に落とし、 batch = m = √n として batch 多項式 D(z)=Π((x−j)−z), N(z) を積木で構築、 z = 0, m, 2m, ... の等差数列多点評価で全 batch を一括評価する。Pollard–Strassen (n! mod N) の float 版。

• 「x^m が p^1.5 桁に爆発」(gemini の反論) は誤り: 桁数と指数の混同。x の magnitude は ~4p なので x^m は指数 ~m·log10(4p) の普通の p桁 float。切り詰めで正しい
• 真の論点は2つ:
  • (a) 展開評価の桁落ち: ~m·log10(4p) + q 桁 (q = −log10|x−最近接ノード|)。ガード O(√p·log p)、q は同一値相殺トリックで処理
  • (b) fast multipoint evaluation の float 安定性: 一般点集合の remainder tree は不安定で使えないが、評価点が等差数列なので Newton (下降階乗) 基底変換 + 畳み込み (Aho–Steiglitz–Ullman 系) が使える。増幅 2^O(m) × 値域比 → ガード O(√p·log p)。ここの誤差解析が本体で、文献に存在しない (新規性)
• コスト: 積木・多点評価 M(m·p)·log m ≈ p^1.5·log²p、prefactor 逐次更新 p^1.5·log p、結合 m·M(p)。計 O(p^1.5·polylog)
• 中間形なし: float 係数にした瞬間 1 点評価では full-prec 乗算に戻るので、多点評価まで一気に作って初めて利得。オール・オア・ナッシング
• 実装の前提条件: 準線形整数乗算。素の CRuby Integer は Toom3 で M(p^1.5)=p^2.2 になり本末転倒 → BigDecimal 自前 NTT に Kronecker で載せる (仮数配列を固定小数点ベクタに流用、C 側に多項式乗算入口)。メモリは係数列 p^1.5 桁 (1e5桁で数百MB)
• 効く領域: full-digit x × p ≳ 数万桁のみ。クロスオーバー予想 1e4–1e5、1e5 で ~10倍、以後 √p/polylog で開く
• メンテ可能性の条件 (これがないなら作らない方がよい): ①誤差台帳文書 (段階ごとのガード桁と根拠) ②遅いパスを oracle にした敵対的テスト (ノード直近・整数±1e-k・full-digit) ③fast path の kill switch
• 最初の一歩: GMP 環境の素 Integer で m=数十, p=数千の玩具プロトタイプを書き、ガード桁の実測カーブを取る (黒なら早期撤退)
• 文献アンカー: Pollard–Strassen 1976 / Chudnovsky–Chudnovsky BSGS / Bostan–Gaudry–Schost (いずれも exact/modular)。Johansson rectangular splitting (ISSAC 2014) は非スカラー乗算回数を減らすだけで bit cost は落とさない。float 版安定性解析つきの実装・論文は見当たらない

5. Arb (python-flint) 比較の実測 (2026-07, このマシン)

• 乗算単体 BigDecimal NTT vs Ruby Integer (GMP): 10.8× / 5.5× / 3.6× @ 1e4/1e5/1e6 桁 (32bit DECDIG の packing 密度 + asm 差込み。Toom 追加で埋まるのは一部)
• プロファイル: BSM (1.125, 1e5) は BigDecimal#mult 55% + add 10% + GC 25% (小 mult ~1e6 回の呼び出し固定費と GC が主犯)。doubling (2e6, 1e5) は Integer#* 48% (GMP) + BigDecimal 25% + GC 13%
• Arb の整数 gamma の正体 (x スイープ @ dps 1e5): x ≤ 1e7 でコスト ∝ x (厳密階乗パス: 0.08→0.9s for 1e6→1e7)、x ≥ 4e7 で ~1.3–1.5s (Stirling、x↑で微減)。x ≈ 2e7 が crossover で、初回はベルヌーイ生成 ~2s を追加で払う (3.67s→2回目 1.62s)。「なぜか arb が遅いエッジケース」= 両戦略の谷間 × cold cache
• BigMath doubling は log x スケール: 2.64 / 3.7 / 4.61s @ 2e6/2e7/2e8。x=2e8 では warm Arb 1.27s に 3.6× 負ける (Stirling は x↑ で安くなるため)。doubling が構造的に勝つのは「ベルヌーイを払えない一発勝負」の状況
• 総合: 差 = 基礎演算 4–10× × Ruby層 (呼び出し固定費+GC) 3–5× × アルゴリズム定数 1.5–3×。指数は互角
• 改善メニュー (効果測定済み):
  • 整数 x の厳密階乗パス → 実施済み相当: int_bsm_prod 化で直接パスが Integer になり吸収 (2026-07-06)
  • integer_factorial 直接パスの Integer 化 → 実施済み (同上)
  • 残ノブ: doubling 閾値 4 * prec → 正しい形は n·log(n) < K·prec(実測で K を決める。32@1e5 実験は K′≈210 相当)
  • 残ノブ: BSM/BSGS crossover 定数(BSM が整数化で速くなったので BSGS 側に動くはず)
  • BSM/BSGS の merge を C カーネル化 or destructive 演算で alloc/GC 削減
  • 64bit DECDIG 化・mid-size Toom (full-digit の 60× のうち 2–4× ぶん、n^1.5 化とは独立に効く)
  • BSM 残余の主役は power(e, prec) の exp/log 系(統合後プロファイルより)

gamma_multipoint.md

多点評価版 gamma(実験実装)— 設計と実測

full-digit x 向けの Lagrange 補間和を、√サイズのバッチ多項式 + 等差数列多点評価で O(PREC^1.5 · polylog) にする実験実装のまとめ。理論的背景(なぜ sub-quadratic は Bernoulli を捨てた側にしかないか、Pollard–Strassen との関係)は gamma_algorithm.md §6 参照。 本文書は実装の設計判断・実測・教訓に絞り、単体で読めるように書く。

• branch: gamma_multipoint_evaluation
• files: lib/bigdecimal/math/gamma_multipoint.rb(本体)、gamma_mp_check.rb(検証: acc / bench / debug、MP_EVAL=fast|horner / MP_ENGINE=coeff|values で強制)
• commits: 6067f05d(初版)→ 437d0445(remainder-tree 評価)→ 2a9ae04a(pack/unpack・guard)→ 2ea0de59(barycentric pair tree)→ 5ec0788c(dispatch 配線)→ 76c68b4c(値領域エンジン)→ d89d1fd3(shift 積も値領域)→ 478774c7(係数エンジン引退・統合)
• 最終形は値領域エンジン単一(337行)。係数エンジン(§2 pair tree + §5 評価2モード)は BSGS↔値領域の crossover(~2500桁)が min_prec = 3000 の内側に入ったため役割を失い引退。 §2・§5 は経緯の記録として残す。3000〜8000桁帯で引退エンジン比 ~1.3× の譲歩と引き換えに、 エンジン1つ・誤差則1つ・コードパス1つになった
• BigMath.gamma に配線済み (commit 5ec0788c): Gamma.gamma_lagrange が full-digit x かつ prec ≥ Multipoint.min_prec (既定 3000、実測 crossover ~2000 に余裕) のとき多点評価へ dispatch。 反射 (x < 0.5)・lgamma・factorial doubling はすべて gamma_lagrange 経由なので自動的に恩恵を受ける
• GMP 判定は defined?(Integer::GMP_VERSION) (CRuby が GMP 付きビルドで定義する定数)。 非 GMP (例: ruby.wasm) では乗算が Toom-Cook で本パイプラインは BSGS より悪化するため自動 off。 Multipoint.enabled = false が kill switch
• 配線前の教訓: テスト用ラッパーが x < 0.5 を丸ごと既存実装へ委譲していて、反射後の full-digit 引数が BSGS に落ちて O(p²) に戻る事故があった (gamma(√2/3) ベンチで発覚)。dispatch を gamma_lagrange 一箇所に置くことでこの種の隙間を構造的に排除した

記法: p = 要求10進精度、prec2 = p + 16(作業精度)、l ≈ p(補間片側ノード数)、 m = バッチサイズ = バッチ数 ≈ √(2l)、n1 = m²(総ノード数)、W = keep(固定小数点の保持ビット数)。

1. パイプライン全景

既存 BSGS と同じ shift(x を 2·prec 以上へ)・b = round(x−1)・l の設定から始める。

1. ノード一般化: A = b − l から n1 = m² 個の連続整数ノード(対称範囲 b±l をわずかに上へ拡張。 m は奇数 — barycentric 再構成の符号 (−1)^(n1−1) を正に保つため)。 係数比は一般形 ρ(t) = −b·(n1−t)/(t·(A+t))、t はノードの0起点添字
2. バッチ多項式化: バッチ k はノード t = km .. km+m−1。BSM の [sum_num, mult_num, den] の構造をバッチ内オフセット z の多項式に持ち上げ、1 組の多項式族(F01, F2)で全バッチを表す
3. 多項式演算: 固定小数点係数(単一スケール)+ Kronecker 詰め込み(Ruby Integer / GMP)
4. 評価: z = 0, m, 2m, … の等差数列で F01, F2(と shift 積の多項式)を評価
5. 結合: バッチごとに C_k(係数比の累積。厳密な整数比 ×m 個で逐次更新)を掛けて sum に加算、E_k を prod に乗算 — 全精度演算は O(m) 回
6. 復元: base = b^(x−A) / (prod·sum)、gamma = base · A! · (n1−1)!(階乗は既存機構)

2. barycentric pair tree(係数エンジンの最終形, 2ea0de59 — 478774c7 で引退)

葉の因子分解が構造の核:

den_j = L_j · B_j · I_j          num_j = −b · L_{j−1} · G_j
L_j = (x−A−j) − z   ← これだけが full-precision (W-bit) 係数
B_j = z + j,  I_j = A + j + z,  G_j = n1 − j − z   ← 小整数係数

num 側の L が den 側の L の添字ずらしなので、prefix-product 和 Σ_j Π_{i≤j} num_i · Π_{i>j} den_i の L 部分は常に「全積 Ω = Π L_i から 1 本抜き」になり、 級数分子は barycentric 形に落ちる:

$$F_{01} = \sum_{j=0}^{K} \frac{\Omega(z)}{L_j(z)} \cdot w_j(z), \qquad K = m-1,\quad w_j = (-b)^j \prod_{i\le j} G_i \prod_{i>j} B_i I_i$$

w_j は小係数因子のみ。木のノードを [Ω, Φ, BI, GX] とする (BI = Π B_i I_i、GX = (−b)^size · Π G_{i+1}、後2者は厳密整数多項式):

Ω_P = Ω_A · Ω_C
Φ_P = (Φ_A · Ω_C) · BI_C + (Ω_A · Φ_C) · GX_A

ワイド×ワイド積は次数 d の Ω に対する 2 本だけ(旧実装は次数 3d の triple 同士 × 4 本)。 ×小係数の積はスロット幅がほぼ半分で済む。j = 0 項(級数先頭の 1)は根で付加し、 その Ω·BI 部分がちょうど F2 = Π den なので再利用できる(F01 = F2 + L_0·Φ·GX_0)。

これは「ノード値の比が小さい有理数になるスケーリングを選んだ」という元の設計利得が もう一段深く効いた形 — full-precision 部分と有理数部分の分離が多項式レベルでも可能だった。

3. 数値表現と安定性

• 固定小数点・単一スケール [coeffs, exp2]: 係数列全体で 2^exp2 を共有。小さい係数は 相対精度が落ちるが、値への寄与も小さいので無害(絶対誤差モデル)。乗算のたび keep bits へ bit-drop
• keep = prec2·log₂10 + 64 + guard、guard = 4·m·bl(n1) + 256 根拠(実測): guard = 0 に強制したときの損失が 2.9〜3.3 · m · bl(n1) bits で prec 300〜10⁴ にわたり安定。しかも評価モード(Horner / remainder tree)間で完全同一、 near-node 入力でも同一。⇒ ①誤差はバッチ間の値の dynamic range に完全支配 ②remainder tree の追加丸めは測定不能 ③相殺 invariant が有効、の3点が実測で裏付けられた
• near-node の相殺 invariant: prod 側の E_k = Π(x−A−km−j) を独立に計算せず、 E_k = F2_k · (x−A−km) / (B_k·I_k)(B_k, I_k は厳密整数)として sum 側と同一の F2_k 計算値から導出する。x がノードに 10^-q 接近しても、桁落ちした同じ値が分子分母で 厳密相殺する(BSGS の batch_prod 再利用と同じ設計)

4. Kronecker 層(Ruby Integer / GMP)

• pack: 係数を固定幅 hex 文字列にして join → to_i(16)(to_s(16)/to_i(16) は線形)。 負係数は borrow を次スロットへ折り込んで単一 pack(正負分割 pack + 巨大減算を排除)
• unpack: to_s(16) → スロット分割 → スロット毎の borrow 伝播(≥ 2^(W−1) を負と解釈)。 巨大バイアス定数の生成・加算を排除
• スロット幅 = bitsA + bitsB + log₂(出力長) + 2 — 両オペランドの max を別々に取る。 2·max にすると小係数×大係数の積が2倍幅で走る(発見前は実際そうなっていて、修正だけで大きく効いた)
• 次数 16 未満は naive 畳み込みにフォールバック

5. 等差数列多点評価(2モード — 478774c7 で引退)

• Horner(デフォルト): 点ごとの Horner。パイプライン唯一の p² 項だが、乗数が小さい整数 (acc × z は W-bit 整数の線形パス)なので定数が machine-word 級に小さい
• remainder tree(:fast): subproduct モジュライは falling-factorial 型で係数が ~d·log₂(count) bit の厳密整数(除算のスケールが良い)。逆級数(Newton 反復)は木に メモ化して F01 / F2 / shift 積の評価で共有。幅広い被除数は R = t^count mod M_root (小係数)でブロック結合してから 1 回だけ deg < 2·count の除算
• 実測 crossover: m ≈ 700(p ≈ 25万桁) — それ未満では Horner の定数が勝つため eval_mode = :auto(閾値超で fast)。値は両モードで厳密一致

6. 値領域エンジン(BGS shift of evaluation values, engine = :values, 76c68b4c)

係数を一切持たない第2エンジン。バッチ遷移の 2×2 行列 P_s(z) = Π_{t=z+1..z+s} [[den_t, 0], [num_t, num_t]] の3成分 (D, N, M) を z = u·s (u = 0..3s) 上の値テーブルで保持し、P_2s(z) = P_s(z)·P_s(z+s) で倍化する。 倍化に必要な u ≤ 12s+3 への延長が shift of evaluation values:

$$Q(a+k) = \frac{\Delta_k}{d!} \sum_{i=0}^{d} Q(i),(-1)^{d-i}\binom{d}{i}\frac{1}{a+k-i}$$

重み C(d,i) は厳密整数(除算なしで畳み込みへ)、カーネル 1/(a+k−i) は小整数の 固定小数点逆数、Δ_k = Π(a+k−j) は厳密整数の逐次更新 — つまり1回の畳み込み + 線形後処理。木も評価ステップも存在せず、総コストは倍化の幾何級数 O(p^1.5·log p)(係数エンジンより log 1本少ない)。

• n1 = S·G+1、S = 2^κ ≈ √(2l)。S 偶数なので barycentric 符号は常に正 (係数エンジンの「m 奇数」制約が消える)
• C_k 更新は M_k/D_k(テーブル値そのもの)— 厳密整数比の再計算も不要
• 相殺 invariant: prod 側は sum と同じ D_k 計算値から E_k = D_k/BI_k で導出(§4 と同じ契約)
• 安定性(実測): guard=0 の損失 = 1.35〜1.49 · S · bl(n1) bits(prec 300〜10⁴、 near-node 同一)。恐れていた値シフトの外挿増幅 2^O(deg) は現れない — 多項式自身が サンプル窓の外で重みと同率に成長するため、誤差は引き続き dyn range 支配。 guard = 2·S·(bl+4)+256(~1.9× マージン)で、係数エンジン(4·m·bl+…)より小さい
• crossover(対 :coeff): ~7000 桁 → 係数エンジン引退後は対 BSGS で ~2500 桁となり、 min_prec = 3000 の内側に収まったためこれが唯一のエンジン(478774c7)
• shift 積も同じ倍化で値領域化(d89d1fd3、1成分・次数 s のテーブル)。パイプライン全体が 一様に O(p^1.5·log p)(残っていた係数領域由来の log² 成分が消滅)

7. 実測(Apple Silicon / GMP 6.3.0 / x = √2 full-digit)

vs 既存 BSGS(全ケースで結果は厳密一致):

p	multipoint	BSGS	比
2,000	0.28s	0.28s	1.0×(crossover ≈ 2000桁)
5,000	1.37s	2.03s	1.5×
10,000	4.4s	7.9s	1.8×
50,000	69.9s	207s	3.0×

精度: sweep(√2、1/3、near-node 7+10^(−p/2)/3、0.6+ε × prec 100〜2000 × 両評価モード) すべて誤差指数 −p または −(p−1)(既存実装と同じ ulp 級)。

5万桁での改善履歴(Horner 経路): 104〜108s(初版)→ 105s(スロット幅修正)→ 94s(borrow pack/unpack)→ 92s(guard 締め)→ 69.9s(pair tree)

エンジン間比較(engine 強制、agree はすべて exact):

p	:values	:coeff	BSGS	:values vs BSGS
2,000	0.33s	0.26s	0.29s	—
5,000	1.84s	1.36s	2.10s	1.1×
10,000	4.1s	4.5s	7.9s	1.9×
50,000	49.1s	65.2s	207s	4.2×

統合後の最終値(shift 値領域化込み、単一エンジン): 1e4 = 4.0s、5e4 = 45.3s(BSGS 比 4.6×)。 実測局所指数(1e4→5e4): BSGS 2.03 → 係数エンジン 1.66 → 値領域 1.53。

8. 教訓

1. 安定性は理論見積り通りだった。「等差数列点上の float remainder tree は安定」が go/no-go の本丸だったが、実測は追加丸めが検出すらできないレベル。損失は loss ≈ 3·m·bl(n1) bits という綺麗なスケーリングに乗り、guard 式をこの実測から確定できた
2. 定数の攻防が本体。漸近的に優れた remainder tree が単純な Horner に実用域で勝てない (crossover 25万桁)のは、schoolbook vs NTT で見た構図の再現。実際に効いたのは 「アルゴリズムを差し替える」より「ワイド積の本数と幅を減らす」構造変更 (pair tree 1.3×、スロット幅・borrow pack 1.2×)で、合計 1.55×
3. 相殺 invariant は評価アルゴリズム非依存の形で設計に残せる(E を F2 から導出する、 という契約にしておけば、評価を Horner から remainder tree に替えても不変)
4. Ruby 税が軽い計算だった: ノード数 O(m log m) の粗粒度演算なので、Θ(p) 個の葉を Ruby で回す BSM/BSGS と違いプリミティブ速度がほぼそのまま出る

9. アルゴリズム族の中での位置づけと再利用性

9.1 これは BGS の事例である

[Ω, Φ] ペアの merge は、部分和の状態遷移を上三角 2×2 行列 [[num_j, 0], [*, den_j]] の積として畳んでいることと等価で、つまり本実装は 「多項式係数の行列階乗」— Chudnovsky–Chudnovsky / Bostan–Gaudry–Schost (BGS) の 問題クラスそのものである(係数領域 = 積木 + 多点評価、の教科書構成メンバー)。 新規なのは機構ではなく、exact/modular の世界の機構を固定小数点 float に移植して 安定性を実測したことと、適用先(任意精度 Γ)の方。

家族の知見が指した次の一手 BGS の "shift of evaluation values"(値領域)は §6 として 実装済み: 予告通り木と評価二択が消えて log が1本減り、8000 桁以遠で係数エンジンを 逆転(5万桁 1.33×)。値シフトの float 安定性(未踏とされた部分)も実測で決着し、 外挿増幅は dyn range の下に沈むことが分かった。 なお exact 実装(FLINT 等)が使う FFT 変換の再利用は Ruby Integer 経由では不可能で、 構造的に ~2-3× をテーブルに残している(transform-blind な実装)。

9.2 機構が要求する性質(適用条件)

和 Σ_n Π_{k≤n} ρ(k) の項比 ρ(k) = P(k)/Q(k) が添字 k の有理式で、係数が 「小さい整数 + 少数の p-bit パラメータ」であること。一般化は次数 s の多項式係数 漸化式(s×s 行列階乗)。Lagrange 補間であることは要求に含まれない (gamma の barycentric 和がこの形の一例だっただけ)。

適用可否の決定的な境界は p-bit 数が「引数」か「パラメータ」か:

p-bit 数の位置	最良既知	例
holonomic ODE の引数(評価点)	bit-burst(準線形)— 本機構の出番なし	erf(x), exp/log/atan の full-digit 点
漸化式のパラメータ(係数)	本機構(p^1.5·polylog)	Γ(x), γ(a,r) の a, pFq のパラメータ, Σ1/(x+k)(digamma 還元)
どちらでもない(指数に入る等)	対象外	ζ(s) full-digit s(項比 (n/(n+1))^s が有理式でない)

9.3 不完全ガンマ級数 — 第2クライアントとして実装・実測済み (9d796205)

incgamma_mp_check.rb: γ(a,r) = r^a·e^{−r}·(1/a)·(1 + Σ_j Π_{i≤j} r/(a+i))、 full-digit a ∈ [0.5, 3]。層の汎用化: batch_value_tables は葉 [den_t, num_t] を block で受ける形になり(9d796205)、クライアント差分は葉の定義+組立の ~100 行。

分子が定数なので 2×2 行列が退化: M_s = r^s は厳密スカラー、テーブルは D, N の 2本・次数 s(gamma クライアントは 3本・次数 3s)— 倍化が項あたり ~3× 軽い。

実測(x = √2、全て厳密一致): Lagrange 多点評価版より速い — 5000桁 2.27× / 1e4 1.90× / 2e4 1.91× / 5e4 1.37×(0.77s vs 1.74s @5000、32.5s vs 44.7s @5e4)。 損失則 0.26〜0.50·S·bl(正項級数+細いテーブルで gamma クライアントの 1.4 より小さい)。

副産物2つ:

1. 項数見積りの訂正: 従来の (1+√2)·r ≈ 5.6p はガウス裾近似の誤用で ~13% 過少 (c ~ 1.4r は Gaussian 領域外)。真の Poisson 裾 γ−(1+γ)ln(1+γ) から γ* ≈ 1.72、 正しくは ~2.72·r ≈ 6.3p 項。過少分は「精度の一定割合 (~28%) の不足」として loss 測定に現れ、そこから発見された(誤差則の実測が仕様バグを釣り上げた実例)
2. 固定小数点テーブルの罠: 基底テーブルを小さいマンティッサ(exp 0)で作ると、 延長との concat で大きい方の exp に丸められ整数精度まで切り捨てられる。 定数でも必ず fixed-point スケールに正規化してから使う

含意: full-digit x の [0.5, ~O(p)] 帯は「incgamma + (必要なら値領域 rising factorial)」の方が Lagrange 多点評価より速く単純。倍数公式(巨大 x)だけが Lagrange 固有の資産として残る。 本線 gamma の lane 再編成をするかは今後の判断(実験スクリプトのまま留め置き)。

10. 残課題

• 10⁵〜10⁶ 桁での測定とメモリプロファイル(係数列は ~p^1.5 桁; 10⁵ で数百MB、10⁶ で GB 級)
• 定数の余地: pack/unpack は依然全体の ~2割、(Φ·Ω)×小係数積の出力幅、 fast 評価の被除数スライス waste など
• guard 式の S 大域(> 10³)での再確認(現状 300〜10⁴ の実測に基づく)
• BGS 値領域版の float 安定性検証 → 実装・実測済み(§6)
• shift 積の値領域化・係数エンジンとの統合 → 完了(d89d1fd3, 478774c7)
• 機構の gamma 非依存部分の切り出し(§9.2 の適用クラス向け汎用層)と incgamma を2番目のクライアントにした比較実験(§9.3)— 値領域エンジンなら incgamma の「分子が定数」の利点はさらに大きい(M テーブルがスカラー列になる)