Float 指数・対数アルゴリズム

概要

指数関数 $\exp(x)$ と対数関数 $\log(x)$ は互いの逆関数であり、 calx Float の数学関数群の基盤を成す。 三角関数・双曲線関数・べき乗などの多くの超越関数がこの 2 つに帰着されるため、 その効率が全体の性能を左右する。

両関数とも、入力精度に応じた多段ディスパッチを行い、 低精度では固定小数点演算 (Q128/Q256) による高速パスを、 高精度では漸近的に最適なアルゴリズムを選択する。

内部計算には RawFloat を使用する。 RawFloat は mpn 配列 (d)・ワード数 (nw)・指数 (exp) のみを持つ軽量構造体であり、 Float オブジェクトの構築・破棄に伴うメモリアロケーションを排除する。 Taylor 級数のループ内で多数の中間値が生成される exp/log にとって、 この設計は性能上不可欠である。

なお、数学定数 $\log 2$, $\log 10$ は log 関数を用いて計算されるが、 log の低精度パス自体は定数に依存しない atanh 級数で実装されており、 循環依存を回避している。

exp のアルゴリズム

$\exp(x)$ は精度に応じて 4 段階にディスパッチされる:

実用範囲ではほぼ常に expDoubling が使用される。 exp_newton は log の計算に exp を必要とするため、 log 側の Halley 反復から間接的に呼び出される場合に限られる。

exp_small / exp_medium

引数削減により $x$ を $[0, \ln 2)$ の範囲に縮小する:

$$x = k \ln 2 + z, \quad k = \lfloor x / \ln 2 \rfloor, \quad z \in [0, \ln 2)$$

$$\exp(x) = 2^k \cdot \exp(z)$$

$\exp(z)$ を Q128 (64 bit 精度) または Q256 (128 bit 精度) の 固定小数点 Taylor 級数で計算する。 $|z| < \ln 2 \approx 0.693$ と小さいため、少ない項数で収束する。 $2^k$ の乗算は指数の加算で実現されるため、コストはゼロである。

引数半減と二乗復元

expDoubling は任意精度の exp を計算する主力アルゴリズムであり、 引数半減 (argument halving) と 二乗復元 (squaring reconstruction) を組み合わせる。

引数半減

$x$ を $K$ 回半減して $z = x / 2^K$ を得る。 $|z|$ は十分小さくなり、Taylor 級数が高速に収束する:

$$\exp(z) = \sum_{k=0}^{N} \frac{z^k}{k!} + O(z^{N+1})$$

半減回数 $K$ は精度 $p$ に応じて選択される。 $K$ を増やすと Taylor の項数 $N$ は減るが、復元の二乗回数が増える。 最適な $K$ は $K \approx \sqrt{p / \log_2 p}$ 程度である。

二乗復元

$\exp(z)$ から $\exp(x)$ を復元する:

$$\exp(x) = \exp(z \cdot 2^K) = \exp(z)^{2^K}$$

$K$ 回の自乗で復元する。RawFloat での二乗復元では mpn::square を使用し、対称性を活用して一般乗算より高速に計算する。

ガードビットの必要性

二乗復元では各ステップで丸め誤差が最大 2 倍に増幅される。 $K$ 回の二乗で誤差は最大 $2^K$ 倍になるため、 $K + \alpha$ ビットのガードビットを追加して計算する。 ここで $\alpha$ は Taylor 級数の打ち切り誤差分のマージンである。

Paterson-Stockmeyer 法

Taylor 級数を Horner 法で評価すると $N$ 回の多倍長乗算を要する。 Paterson-Stockmeyer 法は多項式を $\sqrt{N}$ 次のブロックに分割し、 $O(\sqrt{N})$ 回の多倍長乗算で評価する。

$s = \lceil \sqrt{N} \rceil$ とし、 $p(z) = \sum_{k=0}^{N} a_k z^k$ を以下のように書き直す:

$$p(z) = \sum_{i=0}^{s-1} \left( \sum_{j=0}^{s-1} a_{is+j} \, z^j \right) (z^s)^i$$

1. $z, z^2, \ldots, z^s$ を事前計算する ($s-1$ 回の乗算)
2. 各ブロックの内側多項式を Horner 法で評価する (係数とのスカラー乗算のみ)
3. 外側多項式を Horner 法で $z^s$ について評価する ($s-1$ 回の乗算)

合計 $O(\sqrt{N})$ 回の多倍長乗算で済むため、 大精度計算で Horner 法に対して大幅な性能改善をもたらす。

log のアルゴリズム

$\log(x)$ は精度に応じて 5 段階にディスパッチされる:

引数の分解

入力 $x > 0$ を仮数と指数に分解する:

$$x = f \cdot 2^k, \quad f \in [1, 2)$$

$$\log(x) = \log(f) + k \cdot \log 2$$

$f \in [1, 2)$ なので $\log(f) \in [0, \ln 2)$ と小さく、 級数の収束が保証される。 $k \cdot \log 2$ は定数キャッシュを用いて計算する。

atanh 級数による低精度パス

$\log(f)$ を atanh 級数で計算する:

$$\log(f) = 2 \operatorname{atanh}\!\left(\frac{f-1}{f+1}\right) = 2 \sum_{k=0}^{\infty} \frac{1}{2k+1} \left(\frac{f-1}{f+1}\right)^{2k+1}$$

$t = (f-1)/(f+1) \in [0, 1/3)$ であるため、$t^2 < 1/9$ と小さく、 級数は急速に収束する。 Q128 (64 bit) / Q256 (128 bit) 固定小数点で直接計算する。

このパスは $\log 2$ などの定数を使用しないため、 定数キャッシュの初期化時にも安全に呼び出せる。

Halley 反復による log

231 -- 3000 ビットの精度範囲では、Halley 反復を用いて $\log(f)$ を計算する。 関数名は log_newton であるが、実際の反復公式は Halley 法である。

反復公式

$$y_{n+1} = y_n + 2 \cdot \frac{f - \exp(y_n)}{f + \exp(y_n)}$$

これは $g(y) = \exp(y) - f = 0$ に対する Halley 反復であり、 3 次収束を持つ。 すなわち、正確なビット数が各ステップで約 3 倍に増加する:

$$53 \to 150 \to 450 \to 1350 \to \cdots$$

初期近似

初期近似は精度範囲に応じて選択される:

• $\leq 450$ bit: std::log(double) で $\sim$53 bit の初期近似を得て、1 -- 2 回の Halley 反復で目標精度に到達する
• $> 450$ bit: Q256 atanh (log_medium) で $\sim$240 bit の初期近似を得て、1 回の Halley 反復で目標精度に到達する

exp の呼び出し

各 Halley ステップで $\exp(y_n)$ を計算する必要がある。 $y_n$ の精度は反復が進むにつれて増加するため、 対応する exp ディスパッチも段階的に変わる:

• 第 1 ステップ ($\sim$53 bit): exp_small (Q128 Taylor)
• 第 2 ステップ ($\sim$150 bit): exp_medium (Q256 Taylor) または expDoubling
• 第 3 ステップ ($\sim$450 bit): expDoubling

各ステップの exp 計算は目標精度ではなく現在のステップの精度で行われるため、 初期のステップは非常に軽い。

AGM 法による log

6000 ビットを超える高精度では、算術幾何平均 (AGM) 法を用いる。 3001 -- 6000 ビットでは素数冪分解法 (log_multiprime) が中間段として使用される。

AGM の定義

$a_0, b_0 > 0$ に対し、以下の反復を定義する:

$$a_{n+1} = \frac{a_n + b_n}{2}, \quad b_{n+1} = \sqrt{a_n \cdot b_n}$$

$a_n$ と $b_n$ は共通の極限 $M(a_0, b_0)$ に 2 次収束する。 各ステップは乗算 1 回 + 平方根 1 回 = $O(M(n))$ であり、 全体で $O(\log p)$ ステップ、合計 $O(M(p) \log p)$ の計算量である。

対数への応用

以下の公式で $\log(x)$ を AGM に帰着させる:

$$\log(x) = \frac{\pi}{2 \cdot \operatorname{AGM}(1,\, 4/s)} - m \cdot \log 2$$

ここで $s = x \cdot 2^m$ であり、$m$ は AGM の収束を加速するために $s$ を十分大きくするパラメータである。 $m$ が大きいほど $4/s$ が小さくなり、AGM の初期値 $(1, 4/s)$ の比が大きくなるため、 少ないステップで収束する。

この公式は $\pi$ と $\log 2$ を必要とするが、 これらは thread_local キャッシュに保持されており、 同一精度での再計算は発生しない。

閾値の設定

Halley 反復は 3 次収束だが、各ステップで exp を呼ぶため、 ステップあたりのコストが高い。 AGM 法は各ステップが乗算 + 平方根のみで構成され、 $O(M(p) \log p)$ の漸近計算量を持つ。

実測ベンチマークにより、以下の閾値が設定されている:

• 3000 ビット ($\sim$900 桁): Halley 反復と素数冪分解法 (log_multiprime) の交差点。 250〜900 桁では Halley が素数冪分解法より 30〜40% 高速
• 6000 ビット ($\sim$1800 桁): 素数冪分解法と AGM の交差点

log2, log10 の計算

$\log_2(x)$ と $\log_{10}(x)$ は底変換公式により $\log(x)$ に帰着されるが、 低精度では直接計算の方が高速である。

低精度パスでは $\log_2 e$ や $\log_{10} e$ をコンパイル時定数として Q128/Q256 固定小数点で保持し、atanh 級数の結果と乗算する。 これにより定数キャッシュへの依存を排除し、 キャッシュ初期化中の循環呼び出しを回避する。

高精度パスでは log(x) を計算した後、 thread_local キャッシュに保持された $\log 2$ または $\log 10$ で除算する。 定数は一度計算されると同一精度以下では再計算されない。

設計判断

RawFloat による mpn 直接計算

Taylor 級数のループ内では乗算・除算・加算が毎イテレーション発生する。 Float オブジェクトを使うと、各演算で仮数の Int を構築・破棄するため、 ヒープアロケーションが支配的になる。 RawFloat は d (mpn ポインタ)・nw (ワード数)・exp (指数) のみの構造体であり、mpn 関数を直接呼ぶことでアロケーションゼロの Taylor ループを実現する。

NTT キャッシュ (rf_mul_cached)

Halley 反復では $\exp(y_n)$ の値が分子・分母の両方に現れる。 rf_mul_cached は NTT の forward 変換結果をキャッシュし、 同じオペランドとの乗算を 2 回目以降は逆変換のみで実行する。 これにより Halley 反復 1 ステップあたりの乗算コストを削減する。

effective_bits_ による精度追跡

Float は内部に effective_bits_ を持ち、 値の有効精度を追跡する。 例えば 64 bit 精度の入力に対して 1000 bit 精度の exp を計算しても 結果は 64 bit 精度にしかならない。 effective_bits_ により、低精度入力に対しては 低精度のディスパッチパスが選択され、不要な高精度計算を回避する。

参考文献

• Brent, R. P. and Zimmermann, P. (2010). Modern Computer Arithmetic. Cambridge University Press. Chapter 4: Elementary and Special Function Evaluation.
• Paterson, M. S. and Stockmeyer, L. J. (1973). "On the Number of Nonscalar Multiplications Necessary to Evaluate Polynomials". SIAM Journal on Computing, 2(1), 60--66.
• Borwein, J. M. and Borwein, P. B. (1987). Pi and the AGM. Wiley.
• Muller, J.-M. et al. (2018). Handbook of Floating-Point Arithmetic. 2nd ed. Birkhäuser. Chapter 11: Evaluating Floating-Point Elementary Functions.