calx の Float sqrt は correct rounding をどのように実現していますか？

整数平方根 sqrtRem で根と余りを正確に計算し、余りが非ゼロなら根の LSB を 1 にセットする sticky bit 方式を用います。この情報により IEEE 754 準拠の Round-To-Nearest-Even を含む全丸めモードで correct rounding が保証されます。

Float sqrt はどのように整数平方根に帰着させていますか？

仮数 m と指数 e を持つ浮動小数点数 m·2^e に対し、仮数を 2(p+1) ビットにスケーリングし、指数を偶数に調整します。√(m·2^(2k)) = √m · 2^k なので、整数 sqrtRem で √m を計算すれば浮動小数点の平方根が得られます。

Float sqrt の低精度ファストパスとは何ですか？

64 bit 精度以下では sqrt_newton_1 (std::sqrt 初期近似 + Newton 1 回 + _udiv128/_umul128 intrinsic) を、128 bit 精度以下では sqrt_small_2 (dc_sqrtrem n=2 直接呼び出し) を使用します。Float オブジェクトの構築を可能な限りバイパスし、CPU intrinsic のみで計算します。

calx の Float sqrt と MPFR の性能比はどの程度ですか？

MPFR 4.2.1 との比較で、50 桁 7.1x、100 桁 5.1x、500 桁 3.3x (calx/MPFR 比、MPFR が高速) です。低精度で差が大きいのは Float オブジェクト構築のオーバーヘッドが原因であり、大精度では Zimmermann アルゴリズムの漸近性能が効いて比率が改善します。

Float 平方根アルゴリズム

概要

Float::sqrt は、任意精度の浮動小数点平方根を 整数平方根 (sqrtRem) への帰着により計算する。余りの情報から sticky bit を構成し、 IEEE 754 準拠の correct rounding を全丸めモードで実現する。

$$\sqrt{m \cdot 2^e} = \sqrt{m'} \cdot 2^{e'/2}$$

ここで $m'$ は $2(p+1)$ ビットにスケーリングされた仮数、$e'$ は偶数に調整された指数である。整数の sqrtRem が返す根と余りから、浮動小数点の正しく丸められた平方根を構成する。

入力精度に応じて 3 段階のディスパッチを行う:

精度	仮数部	関数	方式
$\leq 64$ bit	1 ワード	`sqrt_newton_1`	std::sqrt + Newton 1 回 + intrinsic
$\leq 128$ bit	$\leq 2$ ワード	`sqrt_small_2`	dc_sqrtrem (n=2) 直接呼び出し
$> 128$ bit	任意	`sqrt_core`	Int::sqrtRem への帰着

API の使い方については Float API リファレンスを参照。整数平方根の詳細 (sqrtrem1, sqrtrem2, Zimmermann の dc_sqrtrem) は Int 平方根アルゴリズムを参照。

整数平方根への帰着

Float の平方根を整数演算に帰着させる手順を示す。すべてのディスパッチパスで共通の原理である。

スケーリング

浮動小数点数 $x = m \cdot 2^e$ (仮数 $m$、指数 $e$) に対し、目標精度 $p$ ビットの平方根を求める。

根のビット数の決定: $\text{target\_root\_bits} = p + 1$ (丸め用に 1 ビット余分)
仮数のスケーリング: $m$ を $2(p+1)$ ビットの整数 $I$ にスケーリングする。 $\text{base\_shift} = 2(p+1) - \text{bitLength}(m)$ ビットだけ $m$ を左シフトする。
指数の偶数化: 調整後の指数 $e' = e - \text{base\_shift}$ が奇数ならば、さらに 1 ビット追加シフトして偶数に揃える: $$e' \bmod 2 = 0$$
整数平方根: $s = \lfloor \sqrt{I} \rfloor$、$r = I - s^2$ を計算
結果の指数: $e_{\text{result}} = e' / 2$

$$\sqrt{m \cdot 2^{e'}} = \sqrt{I} \cdot 2^{e'/2} \approx s \cdot 2^{e'/2}$$

偶数指数の必要性

$\sqrt{2^{2k}} = 2^k$ は整数だが、$\sqrt{2^{2k+1}} = 2^k \sqrt{2}$ は無理数である。指数を偶数に揃えることで、平方根を整数の世界で完結させる。奇数指数の場合は仮数を 2 倍 ($+1$ ビットシフト) して指数を 1 減らすことで対処する。

Correct Rounding (正確な丸め)

浮動小数点平方根の丸めを正確に行うには、「真の平方根が $s$ と $s+1$ の厳密にどちら側にあるか」を知る必要がある。 calx では sticky bit 方式でこれを実現する。

Sticky bit の仕組み

整数 sqrtRem が返す余り $r = I - s^2$ を利用する:

$r = 0$: $\sqrt{I} = s$ (正確)。根の LSB はそのまま
$r \neq 0$: $s < \sqrt{I} < s + 1$。根の LSB を 1 にセット (sticky bit)

$$s_{\text{sticky}} = \begin{cases} s & \text{if } r = 0 \\ s \mathbin{|} 1 & \text{if } r \neq 0 \end{cases}$$

この情報が後段の setResultPrecision() (丸め処理) に伝搬し、 midpoint (丸め境界上の値) を正しく識別する。丸め処理で $p+1$ ビットの根を $p$ ビットに丸める際、 sticky bit が 1 であれば「真値は midpoint の上」と判定でき、 Round-To-Nearest-Even で正しい方向に丸められる。

丸めモードとの関係

sticky bit により、以下の丸めモードすべてで correct rounding が保証される:

Round-To-Nearest-Even (デフォルト): midpoint で偶数方向に丸める。sticky bit が midpoint の識別に必要
Truncation: 常に切り捨て。sticky bit は影響しない
Round-Up / Round-Down: 余りの有無で丸め方向が決まる

この方式は MPFR が採用している手法と同じであり、 IEEE 754 の要件を満たす。

sqrt_newton_1 ($\leq 64$ bit 精度)

1 ワード仮数の平方根を、Float の抽象化を完全にバイパスして計算する。 Int オブジェクトの構築・破棄を回避し、CPU intrinsic のみで全計算を行う。

128 ビットへのスケーリング: 仮数 $m$ と指数 $e$ から 128 ビット整数 $M$ (M_hi:M_lo) を構成する。 desired_bits = 126 とし、パリティ調整で $+1$ されても 128 ビットに収まる。
初期近似: std::sqrt(double(M)) で $\sim$53 ビット精度の初期近似 $s_0$ を得る。 $M$ は 128 ビットだが double の仮数部は 53 ビットのため、 $s_0$ は上位ビットのみの近似である。
Newton 反復 1 回: $$s_1 = \frac{s_0 + \lfloor M / s_0 \rfloor}{2}$$ _udiv128(M_hi, M_lo, s0, &rem) で 128 ÷ 64 ビット除算を実行する。オーバーフロー安全な平均: (s0 >> 1) + (q >> 1) + ((s0 & q) & 1)。この 1 回の反復で精度が $\sim$53 → $\sim$64 ビットに倍増する。
±1 補正: _umul128(s1, s1, &sq_hi) で $s_1^2$ を計算し、 $s_1^2 \leq M < (s_1 + 1)^2$ を検証する。条件を満たさない場合は $s_1$ を $\pm 1$ 調整する。
Sticky bit: $s_1^2 \neq M$ ならば $s_1$ の LSB を 1 にセットする。

結果を Int(s1) から直接 Float を構成し、 setResultPrecision() で丸める。 SBO (Small Buffer Optimization) により、1 ワードの Int はヒープ確保なしで構成される。

sqrt_small_2 ($\leq 128$ bit 精度)

1–2 ワード仮数の平方根を、mpn レベルの dc_sqrtrem で直接計算する。 Int 型の除算や乗算のオーバーヘッドを回避し、固定サイズのスタックバッファで完結する。

4 ワードバッファへのスケーリング: 仮数を 256 ビット (4 ワード) のスタックバッファ scaled[6] に配置する。 dc_sqrtrem は np[n..n+2l-1] の領域を使うため、 2 ワード余分に確保する。
正規化: dc_sqrtrem は scaled[3] >= 2^{62} を要求する。 std::countl_zero() で先頭ゼロを測定し、偶数ビットのみ左シフトする (奇数シフトは平方根の結果を半端にするため)。
dc_sqrtrem 呼び出し: dc_sqrtrem(sp, scaled, n=2, scratch) を直接呼び出す。結果の根は sp[0..1] (2 ワード)、余りは scaled[0..1] に格納される。キャリーフラグも余りの一部として返される。
逆正規化: 正規化で $t$ ビット左シフトした場合、根を $t/2$ ビット右シフトして復元する ($\sqrt{x \cdot 2^t} = \sqrt{x} \cdot 2^{t/2}$)。
Sticky bit: キャリーまたは余り (scaled[0], scaled[1]) が非ゼロなら sp[0] |= 1 でセットする。

Int::fromRawWordsPreNormalized で 2 ワードから Int を構築し (SBO: ヒープ確保なし)、 Float を構成して setResultPrecision() で丸める。

sqrt_core (汎用)

128 ビットを超える任意精度の一般ルーチン。整数帰着の手順に忠実に従い、 IntSqrt::sqrtRem() で整数平方根を計算する。

precision_bits = precisionToBits(precision), target_root_bits = precision_bits + 1
仮数 $m$ を desired_bits = 2 * target_root_bits ビットに左シフト
指数を偶数に調整
IntSqrt::sqrtRem(scaled, remainder) で $s = \lfloor\sqrt{I}\rfloor$ と余りを計算
余りが非ゼロなら s |= Int(1) (sticky bit)
Float(s, result_exponent, false) を構成し setResultPrecision で丸め

内部では入力サイズに応じて Zimmermann の dc_sqrtrem が再帰的に呼び出される。 $n$ ワードの整数平方根は $O(M(n))$ で計算され ($M(n)$ は $n$ ワード乗算のコスト)、 Newton 反復による除算と同じ漸近計算量である。

関連記事: Int 平方根アルゴリズム — sqrtrem1・sqrtrem2・Zimmermann の dc_sqrtrem・isSquare の詳細

参考文献

Zimmermann, P. (1999). "Karatsuba Square Root". INRIA Research Report RR-3805.
Brent, R. P. and Zimmermann, P. (2010). Modern Computer Arithmetic. Cambridge University Press. Chapter 3.5: Square Root.
Muller, J.-M. et al. (2018). Handbook of Floating-Point Arithmetic. 2nd ed. Birkhäuser. Chapter 8: Square Root.
MPFR Documentation. https://www.mpfr.org/mpfr-current/mpfr.html