Float 数学定数アルゴリズム

概要

calx の Float クラスは以下の 7 つの数学定数を任意精度で提供する:

共通の設計方針は以下の 3 点である:

• thread_local キャッシュ: 同じ精度の再計算を回避する。mutex 不要
• Binary Splitting (BS): 級数を整数演算で再帰分割し、最終段で 1 回だけ Float 除算を行う。計算量は $O(M(n) \cdot \log N)$
• 循環依存の回避: 定数計算は基本四則演算 + sqrt のみを使用し、exp/log 等の超越関数に依存しない

$\pi$ — Chudnovsky 級数

calx は Chudnovsky brothers (1988) の級数を Binary Splitting で評価する。 AGM (Arithmetic-Geometric Mean) 法と比較して、 Binary Splitting + FFT 乗算との相性が良く実測で高速であるため、 Chudnovsky 級数を採用した。

公式

$$\pi = 426880 \sqrt{10005} \cdot \frac{Q}{T}$$

ここで $T/Q$ は以下の級数で与えられる:

$$\frac{T}{Q} = \sum_{k=0}^{N} (-1)^k \cdot \frac{(6k)! \cdot (A + Bk)}{(3k)! \cdot (k!)^3 \cdot C^{3k}}$$

パラメータは $A = 13591409$, $B = 545140134$, $C = 640320$ である。 1 項あたり約 14.18 桁の収束速度を持ち、既知の $\pi$ 級数の中で最速クラスである。

Binary Splitting の構造

なぜ Binary Splitting が必要か

Chudnovsky 級数は有限項で打ち切る ($N$ 項)。素朴に前から順に足しても答えは得られるが、 毎回の加算が多倍長 Float 演算になり、$N$ 項で $O(N \cdot M(n))$ と重い。 Binary Splitting は中間計算をすべて整数 (Int) で行い、 最後の 1 回だけ Float 除算する手法である。 これにより計算量が $O(M(n) \cdot \log N)$ に改善される。

級数の各項を整数関数で表現する

Chudnovsky 級数の第 $k$ 項を分子・分母に分解すると、以下の 3 つの整数関数になる:

$p(k)$ は「前の項から次の項へ進むときに分子に掛かる因子」であり、 $q(k)$ は「前の項から次の項へ進むときに分母に掛かる因子」である。 つまり、級数の隣接項の比 $a_{k}/a_{k-1}$ は $p(k)/(q(k))$ に符号と線形項を掛けたものになっている。

3 つ組 $(P,\; Q,\; T)$ の意味

ある区間 $[a, b)$ の部分和を整数で保持するために、 以下の 3 つの量を定義する:

• $P_{a,b} = \prod_{k=a}^{b-1} p(k)$ — 区間内の分子因子の積
• $Q_{a,b} = \prod_{k=a}^{b-1} q(k)$ — 区間内の分母因子の積
• $T_{a,b}$ — 区間内の部分和 (分母 $Q_{a,b}$ を通分した状態の分子)

$T_{a,b}$ が表しているのは、$Q_{a,b}$ を共通の分母としたときの級数の部分和の分子である。 つまり:

$$\frac{T_{a,b}}{Q_{a,b}} = \sum_{k=a}^{b-1} \frac{t(k)}{q(a) \cdot q(a+1) \cdots q(k)}$$

この等式が成り立つように $T$ を定義する。

再帰マージの導出

区間 $[a, b)$ を中間点 $m = \lfloor (a+b)/2 \rfloor$ で 2 つに分割する。 左半分 $[a, m)$ と右半分 $[m, b)$ から全体を復元する公式は:

$$P_{a,b} = P_{a,m} \cdot P_{m,b}$$

$$Q_{a,b} = Q_{a,m} \cdot Q_{m,b}$$

$$T_{a,b} = T_{a,m} \cdot Q_{m,b} + P_{a,m} \cdot T_{m,b}$$

$P$ と $Q$ はそれぞれ積なので、分割しても単に掛け合わせるだけである。 $T$ のマージが核心で、以下のように理解できる:

• $T_{a,m} \cdot Q_{m,b}$: 左半分の和を、右半分の分母で通分する。 左半分は分母 $Q_{a,m}$ で計算されているが、全体の分母は $Q_{a,m} \cdot Q_{m,b}$ なので、 $T_{a,m}$ に $Q_{m,b}$ を掛けて分母を揃える
• $P_{a,m} \cdot T_{m,b}$: 右半分の和を、左半分の分子因子で補正する。 右半分の各項は「第 $m$ 項から始まる」ので、 $p(a) \cdots p(m-1)$ つまり $P_{a,m}$ が分子に掛かる

葉ノード (1 項の場合)

$b - a = 1$ のとき (区間に 1 項だけ)、再帰の底に到達する。 第 $k$ 項の値を直接計算する:

• $k = 0$ のとき: $P = 1$, $Q = 1$, $T = A = 13591409$
• $k \geq 1$ のとき: $P = p(k)$, $Q = q(k)$, $T = t(k) = (-1)^k \cdot p(k) \cdot (A + Bk)$

具体例: $N = 4$ 項の場合

4 項の級数を以下のように分割する:

              [0, 4)
             /      \
         [0, 2)    [2, 4)
         /   \      /   \
      [0,1) [1,2) [2,3) [3,4)
        ↓     ↓     ↓     ↓
       k=0   k=1   k=2   k=3

まず葉ノードで各 $k$ の $(P, Q, T)$ を計算する。 次にボトムアップでマージしていく。 最終的に得られる $(P_{0,4},\; Q_{0,4},\; T_{0,4})$ から:

$$\pi = 426880 \sqrt{10005} \cdot \frac{Q_{0,4}}{T_{0,4}}$$

を 1 回の Float 除算で得る。$P_{0,4}$ は最終結果には不要なので、 最上位のマージでは $P$ の計算を省略する。

最適化

• 4-way 再帰アンローリング: 通常の 2-way 分割 (左右) を 2 段展開し、4 つの部分問題を 1 関数内で処理する。 関数呼び出しのオーバーヘッドを削減し、並列度を 2 から 4 に拡大する
• 並列 Binary Splitting: BS_PARALLEL_THRESHOLD = 1000 以上の項数で、 左右の部分問題を別スレッドで計算する。 1M 桁の場合、約 71K 項が必要で上位約 6 レベルが並列化される

$e$ — 階乗級数

Napier の定数 $e$ は以下の級数で計算する:

$$e = \sum_{k=0}^{N} \frac{1}{k!}$$

Binary Splitting の構造

この級数では $p(k) = 1$ (常に 1) なので $P$ の追跡は不要であり、 $(Q,\; T)$ の 2 つ組のみで再帰を行う:

$$q(k) = k + 1$$

$$Q_{a,b} = Q_{a,m} \cdot Q_{m,b}$$

$$T_{a,b} = T_{a,m} \cdot Q_{m,b} + T_{m,b}$$

最終的に $e = 1 + T_{0,N} / Q_{0,N}$ を 1 回の Float 除算で得る。

最適化

factorial_bs は $\pi$ の Chudnovsky 級数と同様に 4-way アンローリングと並列マージを適用している。 $P$ を省略できる分、乗算回数が少なく $\pi$ より軽い計算である。

$\log 2$, $\log 10$ — atanh 級数

対数定数は log 関数を呼び出さずに、 atanh の級数展開を Binary Splitting で評価して計算する。 これにより循環依存を完全に回避する。

公式

逆双曲正接関数の級数展開は以下のとおりである:

$$\operatorname{atanh}\!\left(\frac{1}{n}\right) = \frac{1}{n} \sum_{k=0}^{N} \frac{1}{(2k+1) \cdot n^{2k}}$$

これと $\log$ の関係式 $\log\!\left(\frac{1+x}{1-x}\right) = 2 \operatorname{atanh}(x)$ から、 以下の Machin-like 公式を導く:

$$\log 2 = 2 \operatorname{atanh}\!\left(\frac{1}{3}\right)$$

$$\log 10 = 6 \operatorname{atanh}\!\left(\frac{1}{3}\right) + 2 \operatorname{atanh}\!\left(\frac{1}{9}\right)$$

Binary Splitting の構造

$\operatorname{atanh}(1/n)$ の級数に対して:

$$p(k) = 2k + 1$$

$$q(k) = (2k + 3) \cdot n^2$$

$(P,\; Q,\; T)$ を再帰的にマージし、最終段で 1 回の Float 除算を行う。 $\log 10$ では 2 つの atanh 級数を独立に評価して足し合わせる。

$\gamma$ — Brent-McMillan 法

Euler-Mascheroni 定数 $\gamma$ は収束の遅い極限で定義される:

$$\gamma = \lim_{n \to \infty} \left( \sum_{k=1}^{n} \frac{1}{k} - \ln n \right) \approx 0.5772156649\ldots$$

Brent & McMillan (1980) は修正 Bessel 関数を用いた高速収束公式を与えた。

アルゴリズム

パラメータ $n$ を十分大きく取り ($n \approx p / \log 2$、$p$ は目標ビット数)、 以下の 2 つの級数を Binary Splitting で評価する:

• $I = I_0(2n)$: 第一種修正 Bessel 関数の級数
• $S$: $I$ に調和級数の重みを付けた級数

最終的に

$$\gamma = \frac{S}{I} - \log(n)$$

で $\gamma$ を得る。$\log(n)$ は $n$ が整数であるから Float::log() で計算できる ($\gamma$ への循環依存は生じない)。

カタラン定数 $G$

カタラン定数は以下の級数で定義される:

$$G = \sum_{k=0}^{\infty} \frac{(-1)^k}{(2k+1)^2} = 1 - \frac{1}{9} + \frac{1}{25} - \frac{1}{49} + \cdots \approx 0.9159655941\ldots$$

この級数は収束が遅いため、直接使用しない。 代わりに Euler の高速収束公式を Binary Splitting で評価する。

高速収束公式

Euler の公式は中心二項係数を含む級数であり、 1 項あたりの収束が指数的に速い。 Binary Splitting により $(P,\; Q,\; T)$ を再帰的にマージし、 最終段で 1 回の Float 除算で結果を得る。

他の定数と同様、四則演算のみで計算され、超越関数には依存しない。

$\sqrt{2}$

$\sqrt{2}$ は級数展開を使わず、Float::sqrt(Float(2)) で直接計算する。

内部では精度に応じてディスパッチされる:

• 低精度 ($\leq 64$ bit): sqrt_newton_1 — std::sqrt + Newton 1 回 + intrinsic
• 高精度: Zimmermann の dc_sqrtrem — 整数平方根への帰着

詳細は Float 平方根アルゴリズムおよび Int 平方根アルゴリズムを参照。

thread_local キャッシュ

定数の再計算を避けるため、各定数に thread_local なキャッシュを設ける。

キャッシュの構造

各定数に対して以下の 2 つの変数を thread_local で保持する:

thread_local Float cached_pi;
thread_local int  cached_pi_prec = 0;

キャッシュの動作

1. 要求精度 ≤ キャッシュ精度: キャッシュ済みの値を要求精度に丸めて返す
2. 要求精度 > キャッシュ精度: 級数を再評価し、キャッシュを更新する

thread_local の利点

• mutex 不要: 各スレッドが独立したキャッシュを持つため、スレッド間の競合が発生しない
• ロックフリー: キャッシュの読み書きにアトミック操作や排他制御が不要
• トレードオフ: スレッド間でキャッシュを共有しないため、複数スレッドが同じ定数を要求した場合は各スレッドで独立に計算される。ただし、定数計算は通常 1 回のみ呼び出されるため、このオーバーヘッドは小さい

Binary Splitting の計算量

Binary Splitting は級数 $\sum_{k=0}^{N} a_k / b_k$ の効率的な評価手法であり、 calx のすべての定数計算 ($\sqrt{2}$ を除く) の基盤となっている。

原理

$N$ 項の級数を再帰的に半分ずつ分割し、各レベルで整数乗算を行う:

1. 区間 $[0, N)$ を $[0, m)$ と $[m, N)$ に分割 ($m = \lfloor N/2 \rfloor$)
2. 左右の部分和を整数の組 $(P, Q, T)$ として再帰計算
3. 整数乗算と加算でマージ: $T_{0,N} = T_{0,m} \cdot Q_{m,N} + P_{0,m} \cdot T_{m,N}$
4. 再帰の葉 (1 項) では $P$, $Q$, $T$ を定数項から直接計算
5. 最終段で 1 回だけ Float 除算: $\text{result} = T_{0,N} / Q_{0,N}$

計算量解析

再帰の深さは $O(\log N)$ であり、各レベルで $O(M(n))$ の整数乗算を行う ($M(n)$ は $n$ ビット整数の乗算コスト)。 したがって全体の計算量は:

$$O(M(n) \cdot \log N)$$

浮動小数点除算は最終段の 1 回のみであり、 ナイーブに級数を Float で累積する方式の $O(N \cdot M(n))$ と比較して大幅に効率的である。

整数演算に留まる利点

Binary Splitting の核心は、中間計算をすべて整数で行う点にある。 浮動小数点演算では各ステップで丸め誤差が蓄積するが、 整数演算は正確であるため中間結果に誤差がない。 最終段の 1 回の除算でのみ丸めが発生するため、 精度管理が単純で correct rounding を実現しやすい。

参考文献

• Chudnovsky, D. V. and Chudnovsky, G. V. (1988). "Approximations and Complex Multiplication According to Ramanujan". In: Ramanujan Revisited. Academic Press, pp. 375–472.
• Brent, R. P. and McMillan, E. M. (1980). "Some New Algorithms for High-Precision Computation of Euler's Constant". Mathematics of Computation, 34(149), pp. 305–312.
• Haible, B. and Papanikolaou, T. (1998). "Fast Multiprecision Evaluation of Series of Rational Numbers". Algorithmic Number Theory (ANTS-III). LNCS 1423, pp. 338–350.
• Brent, R. P. and Zimmermann, P. (2010). Modern Computer Arithmetic. Cambridge University Press. Chapter 4.8: Binary Splitting.
• Borwein, J. M. and Borwein, P. B. (1987). Pi and the AGM. Wiley-Interscience.