[インデックス 10473] ファイルの概要

このコミットは、Go言語の math パッケージにおける Cbrt (立方根) 関数のパフォーマンス改善に関するものです。特に、amd64 アーキテクチャでは127 ns/opから105 ns/opへ、386 アーキテクチャでは208 ns/opから169 ns/opへと、処理速度が向上しています。これは、立方根の初期推定方法を最適化することで達成されました。

コミット

commit f3aa54e30de2c2c71a8735c5f61c9c1d93f7cd9f
Author: Charles L. Dorian <cldorian@gmail.com>
Date:   Mon Nov 21 09:56:07 2011 -0500

    math: faster Cbrt
    
    For amd64, from 127 to 105 ns/op; for 386, from 208 to 169 ns/op.
    
    R=rsc, golang-dev
    CC=golang-dev
    https://golang.org/cl/5412056

GitHub上でのコミットページへのリンク

https://github.com/golang/go/commit/f3aa54e30de2c2c71a8735c5f61c9c1d93f7cd9f

元コミット内容

math: faster Cbrt

amd64 では127 ns/opから105 ns/opへ、386 では208 ns/opから169 ns/opへ高速化。

変更の背景

Cbrt 関数は、与えられた数値の立方根を計算する関数です。浮動小数点数の立方根計算は、一般的にニュートン法などの反復計算を用いて近似的に求められます。この反復計算の収束速度は、初期推定値の精度に大きく依存します。初期推定値が真の値に近いほど、必要な反復回数が減り、計算速度が向上します。

このコミットの背景には、Cbrt 関数の既存の実装における初期推定の効率を改善し、より高速な立方根計算を実現するという目的があります。特に、amd64 や 386 といった異なるアーキテクチャでのパフォーマンス向上を目指しています。

前提知識の解説

浮動小数点数表現 (IEEE 754)

Go言語の float64 型は、IEEE 754倍精度浮動小数点数標準に従います。この標準では、浮動小数点数は以下の形式で表現されます。

(-1)^s * m * 2^e

ここで、

s は符号ビット (0または1)
m は仮数 (mantissa または significand)
e は指数 (exponent)

です。

`Frexp` 関数と `Ldexp` 関数

Frexp(x) 関数は、浮動小数点数 x を frac * 2^exp の形式に分解します。ここで、frac は 0.5 <= |frac| < 1.0 の範囲の仮数部、exp は整数指数部です。これにより、浮動小数点数を仮数部と指数部に分離できます。

Ldexp(frac, exp) 関数は、frac * 2^exp を計算し、浮動小数点数を再構築します。これは Frexp の逆操作にあたります。

これらの関数は、浮動小数点数の計算において、仮数部と指数部を個別に操作する際に非常に有用です。特に、数値のスケールを調整したり、特定の範囲に正規化したりする際に利用されます。

立方根の計算と初期推定

立方根 y = x^(1/3) を計算する際、x = f * 2^e と表現すると、 y = (f * 2^e)^(1/3) = f^(1/3) * 2^(e/3) となります。

ここで問題となるのは、e/3 が整数にならない場合があることです。例えば、e が3の倍数でない場合、e/3 は小数になります。これを解決するために、e を3の倍数になるように調整し、その調整を f に反映させることで、f の範囲を適切に保ちつつ、e/3 を整数にすることができます。

初期推定は、反復法（例えばニュートン法）の開始点として使用されます。良い初期推定値は、反復回数を減らし、計算の収束を速めます。一般的に、多項式近似やルックアップテーブルなどが初期推定に用いられます。

技術的詳細

このコミットの技術的詳細は、Cbrt 関数における初期推定のロジック変更にあります。

元のコードでは、Frexp(x) で x を f * 2^e に分解した後、e を3で割った余り m に基づいて f を調整していました。

// 元のコード
// Reduce argument
f, e := Frexp(x)
m := e % 3
if m > 0 {
    m -= 3
    e -= m // e is multiple of 3
}
f = Ldexp(f, m) // 0.125 <= f < 1.0
// Estimate cube root
switch m {
case 0: // 0.5 <= f < 1.0
    f = A1*f + A2 - A3/(A4+f)
case -1: // 0.25 <= f < 0.5
    f = B1*f + B2 - B3/(B4+f)
default: // 0.125 <= f < 0.25
    f = C1*f + C2 - C3/(C4+f)
}

この元のコードでは、f = Ldexp(f, m) の行で f の範囲を 0.125 <= f < 1.0 に正規化していました。そして、その f の範囲に応じて switch m 文で異なる近似多項式（A, B, C の係数を持つ）を使用して初期推定を行っていました。

変更後のコードでは、Frexp(x) の直後に f の範囲を調整するのではなく、switch m 文の各ケース内で f を調整するように変更されています。

// 変更後のコード
// Reduce argument and estimate cube root
f, e := Frexp(x) // 0.5 <= f < 1.0
m := e % 3
if m > 0 {
    m -= 3
    e -= m // e is multiple of 3
}
// Estimate cube root
switch m {
case 0: // 0.5 <= f < 1.0
    f = A1*f + A2 - A3/(A4+f)
case -1:
    f *= 0.5 // 0.25 <= f < 0.5
    f = B1*f + B2 - B3/(B4+f)
default: // m == -2
    f *= 0.25 // 0.125 <= f < 0.25
    f = C1*f + C2 - C3/(C4+f)
}

この変更のポイントは以下の通りです。

f = Ldexp(f, m) の削除: 元のコードでは switch 文に入る前に f を Ldexp(f, m) で調整していましたが、これが削除されました。
switch 文内での f の調整:
- m == -1 の場合 (0.25 <= f < 0.5 の範囲に対応): f *= 0.5 が追加されました。これにより、f の値が適切な範囲に調整されます。
- m == -2 の場合 (0.125 <= f < 0.25 の範囲に対応): f *= 0.25 が追加されました。これにより、f の値が適切な範囲に調整されます。

この変更により、Ldexp 関数を呼び出す回数が減り、また f の調整がより直接的に行われることで、初期推定の計算パスが最適化され、全体的なパフォーマンスが向上したと考えられます。Ldexp は浮動小数点演算を伴うため、その呼び出し回数を減らすことはパフォーマンス改善に直結します。

具体的には、Frexp によって得られる f は常に 0.5 <= f < 1.0 の範囲にあります。e を3の倍数にするために m を調整する際、f の値も 2^m 倍される必要があります。元のコードでは Ldexp(f, m) でこれを行っていましたが、新しいコードでは m の値に応じて f に 0.5 または 0.25 を直接乗算することで、同じ効果をより効率的に実現しています。

コアとなるコードの変更箇所

src/pkg/math/cbrt.go ファイルの Cbrt 関数内で、浮動小数点数の仮数部 f の調整と初期推定のロジックが変更されています。

--- a/src/pkg/math/cbrt.go
+++ b/src/pkg/math/cbrt.go
@@ -45,22 +45,21 @@ func Cbrt(x float64) float64 {
 		x = -x
 		sign = true
 	}
-	// Reduce argument
-	f, e := Frexp(x)
+	// Reduce argument and estimate cube root
+	f, e := Frexp(x) // 0.5 <= f < 1.0
 	m := e % 3
 	if m > 0 {
 		m -= 3
 		e -= m // e is multiple of 3
 	}
-	f = Ldexp(f, m) // 0.125 <= f < 1.0
-
-	// Estimate cube root
 	switch m {
 	case 0: // 0.5 <= f < 1.0
 		f = A1*f + A2 - A3/(A4+f)
-	case -1: // 0.25 <= f < 0.5
+	case -1:
+		f *= 0.5 // 0.25 <= f < 0.5
 		f = B1*f + B2 - B3/(B4+f)
-	default: // 0.125 <= f < 0.25
+	default: // m == -2
+		f *= 0.25 // 0.125 <= f < 0.25
 		f = C1*f + C2 - C3/(C4+f)
 	}
 	y := Ldexp(f, e/3) // e/3 = exponent of cube root

コアとなるコードの解説

変更の核心は、f の初期推定範囲を調整する方法の変更です。

f = Ldexp(f, m) の削除: 元のコードでは、Frexp で得られた f (範囲 0.5 <= f < 1.0) を、m の値に応じて Ldexp(f, m) を使って 0.125 <= f < 1.0 の範囲に調整していました。この Ldexp の呼び出しは、浮動小数点演算を伴うため、コストがかかる可能性があります。
switch m 文内での f の直接調整: 変更後のコードでは、Ldexp(f, m) を削除し、switch m 文の各ケース内で f を直接乗算することで、同じ範囲調整を実現しています。
- m == 0 の場合: f はすでに 0.5 <= f < 1.0 の範囲にあるため、変更は不要です。
- m == -1 の場合: e を3の倍数にするために e から m を引いた結果、f は 2^m 倍される必要があります。m = -1 なので、f は 2^-1 = 0.5 倍される必要があります。f *= 0.5 を行うことで、f の範囲が 0.25 <= f < 0.5 になります。
- m == -2 の場合: 同様に、m = -2 なので、f は 2^-2 = 0.25 倍される必要があります。f *= 0.25 を行うことで、f の範囲が 0.125 <= f < 0.25 になります。

この変更により、Ldexp の呼び出しが不要になり、よりシンプルな乗算で f の範囲調整が行われるため、計算効率が向上し、結果として Cbrt 関数の高速化が実現されました。これは、浮動小数点演算の特性を理解し、より低コストな操作に置き換えることでパフォーマンスを改善する典型的な例と言えます。

参考にした情報源リンク

IEEE 754 浮動小数点数標準に関する情報 (Wikipediaなど): https://ja.wikipedia.org/wiki/IEEE_754
ニュートン法に関する情報 (Wikipediaなど): https://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95
Go言語のコードレビューシステム (Gerrit): https://go-review.googlesource.com/ (コミットメッセージに記載されている https://golang.org/cl/5412056 は、このGerritシステムへのリンクです。)
浮動小数点数の立方根計算に関するアルゴリズムの論文や記事 (一般的な情報源):
- "Fast Inverse Square Root" (Quake III Arenaのコードで有名ですが、浮動小数点数の高速計算の概念が参考になります)
- 数値計算の教科書やオンラインリソース

comemo