calx の Float はどのような内部表現を使っていますか？

calx の Float は mantissa × 2^exponent の形式で表現されます。mantissa は多倍長整数 (Int)、exponent は int64_t、符号は別途保持されます。精度は10進桁数で指定し、内部で ceil(precision × 3.32192809488736) によりビット数に変換されます。デフォルト精度は 53 ビット (約16桁) で、thread_local 変数で管理されます。

Float の加減算で桁落ちはどう扱われますか？

ほぼ同じ大きさの数の減算では上位ビットが相殺されて有効桁数が減少します (桁落ち)。calx では effective_bits_ フィールドで各値の有効精度を追跡しており、桁落ちが発生すると effective_bits_ が自動的に減少します。後続の演算ではこの情報を参照して不必要な高精度計算を回避します。

FloatOps 3引数 API はどのような場面で有効ですか？

反復計算やループ内で同じ変数に結果を繰り返し格納する場合に有効です。FloatOps::mul(lhs, rhs, result) のように結果バッファを再利用することで、オブジェクト構築やヒープ確保のオーバーヘッドを排除できます。19桁精度の乗算で約1.8倍、38桁精度で約2.0倍の高速化が得られます。

Float の除算はどのアルゴリズムを使っていますか？

被除数を目標精度分だけ左シフトしてから Int の除算を実行します。Int の除算は 64ワード未満では Schoolbook、64ワード以上では Burnikel-Ziegler アルゴリズムが自動選択されます。除数が1リム (64ビット) の場合は divmod_1 ファストパスにより汎用除算をバイパスします。

Float 四則演算アルゴリズム

概要

Float は多倍長浮動小数点数を表現するクラスであり、内部的には以下の3要素で構成される:

$$\text{value} = \text{mantissa} \times 2^{\text{exponent}}$$

mantissa (Int): 仮数部。多倍長整数として保持される。
exponent (int64_t): 指数部。2の冪乗のスケールを表す。
sign: 符号ビット。仮数部とは独立に管理される。

精度は10進桁数で指定する。デフォルト精度は default_precision_ = 53 (10進約16桁相当) であり、 thread_local 変数として管理される。内部では precisionToBits 関数により10進桁数からビット数への変換が行われる:

$$\text{bits} = \lceil \text{precision} \times 3.32192809488736 \rceil = \lceil \text{precision} \times \log_2 10 \rceil$$

各演算の結果は setResultPrecision() により目標精度に丸められる。丸めモードは RoundingMode 列挙型で指定し、以下の4種類がある:

ToNearest (デフォルト): 最近接偶数丸め (Round-To-Nearest-Even)
Truncation: 0方向への切り捨て
Up: 正の無限大方向への切り上げ
Down: 負の無限大方向への切り下げ

effective_bits_ フィールドは各値の実際の有効精度を追跡する。入力が低精度 (たとえば double から変換された53ビット値) であれば、 effectiveComputePrecision により計算精度が自動的に削減され、不必要な高精度計算を回避する。

API の使い方については Float API リファレンスを参照。

加算・減算

浮動小数点の加減算は、指数の異なる2つのオペランドの仮数部を位置合わせし、整数レベルの加減算を行う操作である。同符号の加算は addUnsigned、異符号の加算 (実質的な減算) は subtractUnsigned に委譲される。

指数の位置合わせ

2つのオペランド $a = m_a \times 2^{e_a}$ と $b = m_b \times 2^{e_b}$ ($e_a \geq e_b$) の加算を考える。指数の差 $\Delta e = e_a - e_b$ を計算し、$m_b$ を $\Delta e$ ビットだけ右シフトして位置合わせする:

$$a + b = (m_a \cdot 2^{\Delta e} + m_b) \times 2^{e_b}$$

実装では、小さいほうの仮数部をシフトするのではなく、大きいほうの仮数部を左シフトして精度を保つ方式を採用している。計算精度は $\max(\text{lhs\_bits}, \text{rhs\_bits}) + \text{guard bits}$ で決定される。

mpn レベルの加減算

位置合わせ後の仮数部は mpn レベルの add / sub で処理される。ワード単位のキャリー/ボロー伝播により、任意長の加減算を効率的に実行する。

rvalue reference による最適化

operator+ と operator- は rvalue reference 版のオーバーロードを提供する。一時オブジェクトが渡された場合、その内部バッファを再利用することでヒープ確保のコストを削減する。

桁落ち (catastrophic cancellation)

ほぼ同じ大きさの数の減算では上位ビットが相殺され、有効桁数が大幅に減少する。 subtractUnsigned では結果の先頭ゼロ数を検出し、 effective_bits_ を対応する分だけ減少させる。この情報は後続の演算で参照され、失われた精度を超える計算を回避する。

乗算

Float の乗算は、仮数部の整数乗算・指数の加算・符号の決定の3ステップで構成される。

仮数部の乗算

仮数部の乗算は Int の乗算にそのまま帰着する。 Int の乗算はオペランドのワード数に応じて Basecase / Karatsuba / Toom-Cook-3 / Toom-Cook-4 / NTT が自動選択される (Int 乗算アルゴリズムを参照)。

指数と符号

$$\text{result.exponent} = \text{lhs.exponent} + \text{rhs.exponent}$$

符号は XOR で決定される (同符号なら正、異符号なら負)。

精度制限

$n$ ビットの仮数部同士の乗算結果は最大 $2n$ ビットになるが、目標精度を超える下位ビットは不要である。 setResultPrecision により目標精度に丸め、余分なビットを切り捨てる。

Small Buffer Optimization (SBO)

結果の仮数部が2ワード (128ビット) 以下の場合、 Int の SBO により仮数部がスタック上に保持され、ヒープ確保が発生しない。 double 精度 (53ビット) 同士の乗算では、結果が106ビット $\leq$ 128ビットとなるため SBO の恩恵を受ける。

除算

Float の除算は、被除数の仮数部を十分に左シフトしてから Int の除算を実行することで、目標精度の商を得る。

基本手順

$a \div b$ を計算する場合、被除数 $a$ の仮数部 $m_a$ を目標精度 + ガードビット分だけ左シフトする:

$$q = \lfloor m_a \cdot 2^{\text{shift}} \mathbin{/} m_b \rfloor$$

指数は次のように調整される:

$$\text{result.exponent} = e_a - e_b - \text{shift}$$

この手法により、浮動小数点除算を整数除算に帰着できる。 Int の除算は 64 ワード未満では Schoolbook (div_basecase)、 64 ワード以上では Burnikel-Ziegler が自動選択される (Int 除算アルゴリズムを参照)。

1リム除数のファストパス

除数の仮数部が1リム (64ビット) に収まる場合は、 divmod_1 により汎用除算をバイパスする。 divmod_1 は単一ワードの除算を反復するだけであり、 Burnikel-Ziegler の再帰オーバーヘッドや一時バッファの確保が不要になる。整数リテラルや double 由来の小さな除数ではこのパスが頻繁に使われる。

Newton 逆数反復との関係

calx の Float 除算は直接的に Newton 逆数反復を使うわけではない。 Newton 逆数反復 (mu_div_qr) は Int 側の大精度除算で Burnikel-Ziegler の内部で間接的に利用される。具体的には、被除数と除数のワード数の比が大きい不均衡除算 (unbalanced division) において、Newton 法で除数の逆数近似を求め、乗算に帰着する動的閾値方式が適用される。

3引数 FloatOps API

FloatOps は Float の四則演算を3引数形式で提供するクラスである。通常の operator+ 等は新しいオブジェクトを構築して返すが、 3引数版は既存のバッファに結果を直接書き込むことで、オブジェクト構築とヒープ確保のオーバーヘッドを排除する。

// 通常の演算子版
Float c = a * b;        // 一時オブジェクト構築 + コピー/ムーブ

// FloatOps 3引数版
FloatOps::mul(a, b, c); // c のバッファに直接書き込み

FloatOps::mul

乗算の3引数版は、result.mantissa_ のワード配列に mpn 乗算の結果を直接書き込む。通常の operator* ではスタック上の一時バッファに乗算結果を格納してから結果オブジェクトにコピーする必要があるが、 3引数版ではこのコピーが完全に排除される。

FloatOps::add

加算の3引数版は、指数差 ($\text{exp\_diff} \neq 0$) を含む一般的なケースも直接実装している。 BUF_LIMIT = 128 のスタックバッファを活用し、 128ワード (約 $2{,}400$ 桁) 以下の仮数部ではヒープ確保を完全に回避する。

FloatOps::div

除算の3引数版は、結果バッファの再利用に加え、除数が1リムの場合の divmod_1 ファストパスも備える。

3引数 API による改善は、小精度の演算ほど顕著である。オブジェクト構築やヒープ確保のコストが演算本体に対して相対的に大きいためである。

丸めと精度

setResultPrecision

すべての四則演算の最後に setResultPrecision() が呼ばれ、仮数部を目標ビット数に丸める。具体的には:

現在の仮数部のビット数が目標を超えていれば、超過分を右シフトする
切り捨てられるビット列を検査し、丸めモードに従って丸め方向を決定する
必要に応じて仮数部に1を加算 (round-up) する
指数をシフト量だけ加算して補正する

Round-To-Nearest-Even

デフォルトの丸めモードは最近接偶数丸めである。切り捨てビットがちょうど中間値 (midpoint) の場合、結果の最下位ビットが偶数になる方向に丸める。これは IEEE 754 の既定丸めモードと同一であり、統計的な偏り (rounding bias) を最小化する。

round-up の最適化

丸め上げが必要な場合、mantissa_ + Int(1) と記述すると Int(1) の一時オブジェクト構築と多倍長加算が発生する。 calx では IntOps::addDelta(mantissa_, 1) を用いて最下位ワードへの即値加算とキャリー伝播のみで処理する。これにより一時オブジェクトの構築コストが完全に回避される。

normalize

normalize() は仮数部の先頭ゼロを除去し、指数を対応する分だけ調整する。加減算で上位ビットが相殺された場合や、丸め後に先頭ゼロが生じた場合に呼ばれる。 clz (count leading zeros) 命令を用いて先頭ゼロのワード数とビット数を高速に検出する。

effective_bits_ による精度追跡

effective_bits_ は各 Float 値が実際に保持する有効精度を表す。たとえば double から変換された値は effective_bits_ = 53 となる。 effectiveComputePrecision は2つのオペランドの effective_bits_ から実際に必要な計算精度を算出する:

$$\text{compute\_precision} = \min(\text{target\_precision},\; \max(\text{lhs.effective\_bits\_},\; \text{rhs.effective\_bits\_}) + \text{guard})$$

これにより、たとえば10万桁精度が設定されていても、入力が double 由来であれば53ビット精度で計算され、無駄な高精度計算を回避する。結果の effective_bits_ は演算の種類に応じて適切に伝播される。