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

コミット

このコミット f379ea0b07a28aad1f95abcc5ec26254978c0745 は、Go言語の標準ライブラリ src/lib/math における Log (自然対数)、Exp (指数関数)、Pow (べき乗) 関数の精度向上を目的としたものです。また、テストファイルの整理として test.go が all_test.go にリネームされています。

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

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

元コミット内容

commit f379ea0b07a28aad1f95abcc5ec26254978c0745
Author: Russ Cox <rsc@golang.org>
Date:   Thu Nov 20 10:54:02 2008 -0800

    more accurate Log, Exp, Pow.
    move test.go to alll_test.go.
    
    R=r
    DELTA=1024  (521 added, 425 deleted, 78 changed)
    OCL=19687
    CL=19695
--
 src/lib/math/Makefile                 |  19 ++---
 src/lib/math/{test.go => all_test.go} |  36 ++++----
 src/lib/math/exp.go                   | 154 +++++++++++++++++++++++++++-------
 src/lib/math/log.go                   | 146 ++++++++++++++++++++++++--------
 src/lib/math/pow.go                   | 106 ++++++++++++++---------
 src/lib/math/sin.go                   |   4 +
 6 files changed, 329 insertions(+), 136 deletions(-)

変更の背景

このコミットの主な背景は、Go言語の math パッケージにおける Log、Exp、Pow といった基本的な数学関数の計算精度を向上させることにあります。浮動小数点演算の精度は、科学技術計算、金融、グラフィックスなど、多くの分野で極めて重要です。初期のGo言語の math パッケージは、必ずしも最高の精度を提供しているわけではありませんでした。

コミットメッセージにある "more accurate Log, Exp, Pow." は、これらの関数が以前よりも正確な結果を返すように改善されたことを示しています。特に、exp.go と log.go の変更には、FreeBSDの msun ライブラリ（これはIEEE 754浮動小数点標準に準拠した高品質な数学関数実装で知られています）からのコードとコメントが引用されており、これはより堅牢で正確なアルゴリズムへの移行を示唆しています。

また、test.go から all_test.go へのリネームは、テストコードの整理と、将来的にすべてのテストをこのファイルに集約する意図があったと考えられます。

前提知識の解説

浮動小数点数とIEEE 754標準

コンピュータにおける数値表現の一つで、非常に広い範囲の数値を表現できます。しかし、その性質上、常に厳密な精度を保証するわけではありません。IEEE 754は、浮動小数点数の表現形式、演算方法、丸め規則などを定めた国際標準であり、多くのプログラミング言語やハードウェアで採用されています。この標準に準拠することで、異なるシステム間での浮動小数点演算の互換性と予測可能性が向上します。

ULP (Unit in the Last Place)

浮動小数点数の精度を測る指標の一つです。ある浮動小数点数のULPは、その数と、その数に最も近い表現可能な次の浮動小数点数との間の距離を指します。数学関数の精度が「1 ULP以内」であるとは、計算結果が真の値から1 ULP以上離れていないことを意味し、これは非常に高い精度を示します。

Remesアルゴリズム

関数近似のための数値アルゴリズムの一つで、与えられた区間内で多項式や有理関数を用いて、元の関数を最も均一に近似する（最大誤差を最小化する）方法です。数学ライブラリで高精度な関数を実装する際によく用いられます。

引数削減 (Argument Reduction)

数学関数、特に周期性を持つ関数（例: 三角関数）や、特定の範囲で近似が容易な関数（例: 指数関数、対数関数）において、入力引数をより小さな、扱いやすい範囲に変換する技術です。これにより、近似多項式の適用範囲を限定し、精度を維持しつつ計算効率を高めることができます。例えば、exp(x) の計算では x = k*ln2 + r の形に変換し、exp(x) = 2^k * exp(r) として計算することで、exp(r) の計算をより小さな r の範囲で行うことができます。

sys.ldexp と sys.frexp

Go言語の math パッケージ（または内部の sys パッケージ）で提供される関数で、浮動小数点数を mantissa * 2^exponent の形式に分解・再構築するために使用されます。

• frexp(x): x を mantissa と exponent に分解します。x = mantissa * 2^exponent となり、0.5 <= |mantissa| < 1.0 です。
• ldexp(mantissa, exponent): mantissa * 2^exponent を計算します。 これらの関数は、浮動小数点数の内部表現を操作し、特定の計算を効率的かつ正確に行うために利用されます。

技術的詳細

このコミットでは、exp.go と log.go の実装が大幅に変更されています。

exp.go の変更点

• アルゴリズムの変更: 以前のシンプルな多項式近似から、FreeBSDの msun ライブラリの e_exp.c に基づく、より洗練されたアルゴリズムに移行しています。この新しいアルゴリズムは、以下のステップで構成されます。
  1. 引数削減: x = k*ln2 + r の形式に x を分解します。ここで |r| <= 0.5*ln2 となります。r は hi - lo の形式で表現され、精度が向上します。
  2. exp(r) の近似: exp(r) を区間 [0, 0.34658] で特殊な有理関数 R(r**2) を用いて近似します。R(z) ~ 2.0 + P1*z + P2*z**2 + ... + P5*z**5 (ここで z=r*r) の形式で近似し、exp(r) = 1 + 2*r / (R - r) または exp(r) = 1 + r + r*R1(r) / (2 - R1(r)) (より高精度な形式) を用いて計算します。
  3. スケールバック: exp(x) = 2^k * exp(r) を用いて最終結果を導出します。
• 定数の追加: Ln2, HalfLn2, Ln2Hi, Ln2Lo, Log2e, P1 から P5 までの多項式係数、Overflow, Underflow, NearZero といった定数が追加され、より正確な計算と特殊ケースのハンドリングを可能にしています。
• 特殊ケースのハンドリング: NaN, ±Inf, 0 の引数に対する挙動が明示的に定義され、IEEE 754標準に準拠した振る舞いを保証しています。特に、x が Overflow や Underflow の閾値を超えた場合の ±Inf や 0 の返却が追加されています。
• 精度目標: コメントには「エラーは常に1 ULP未満」と記載されており、非常に高い精度を目指していることがわかります。

log.go の変更点

• アルゴリズムの変更: こちらもFreeBSDの msun ライブラリの e_log.c に基づく、より高精度なアルゴリズムに移行しています。
  1. 引数削減: x = 2^k * (1+f) の形式に x を分解します。ここで sqrt(2)/2 < 1+f < sqrt(2) となります。
  2. log(1+f) の近似: s = f/(2+f) を定義し、log(1+f) = 2s + 2/3 s**3 + ... の展開を利用します。R(z) を多項式で近似し、log(1+f) = f - (hfsq - s*(hfsq+R)) (ここで hfsq = f*f/2) の形式で計算することで、高精度を実現しています。
  3. 最終結果の導出: log(x) = k*ln2 + log(1+f) を計算します。ln2 は ln2_hi + ln2_lo に分割され、精度を維持しています。
• 定数の追加: Ln2Hi, Ln2Lo, Lg1 から Lg7 までの多項式係数、Two54, TwoM20, TwoM1022, Sqrt2 といった定数が追加されています。
• 特殊ケースのハンドリング: NaN, ±Inf, 0 の引数に対する挙動が明示的に定義され、IEEE 754標準に準拠した振る舞いを保証しています。

pow.go の変更点

• 特殊ケースのハンドリングの強化: y=0, y=1, x=0 (y>0, y<0), y=0.5, y=-0.5 といった特殊なべき乗のケースが明示的に処理されるようになりました。これにより、これらの一般的なケースでの計算が最適化され、精度と効率が向上します。
• 負の底の処理: x < 0 かつ y が整数でない場合の NaN の返却が追加され、数学的に未定義なケースが適切に処理されるようになりました。
• Exp(y * Log(x)) へのフォールバック: y が非常に大きい場合 (yi >= 1<<63) には、Exp(y * Log(x)) を用いた計算にフォールバックするロジックが追加されています。これは、直接的なべき乗計算がオーバーフローや精度問題を引き起こす可能性があるため、対数と指数関数に変換して計算する一般的な手法です。
• バイナリ分割による計算: ans *= x^yi の部分で、sys.frexp を用いて x を x1 * 2^xe に分解し、yi のビットをチェックしながら x1 を繰り返し二乗していくことで、効率的に整数べき乗を計算しています。これにより、大きな整数べき乗も正確に計算できるようになります。

all_test.go (旧 test.go) の変更点

• ファイル名変更: test.go から all_test.go へのリネーム。
• Tolerance, Close, VeryClose 関数の導入:
  • Close(a,b float64) bool 関数が Tolerance(a,b,e float64) bool を呼び出すように変更され、許容誤差 e を明示的に指定できるようになりました。
  • 新たに VeryClose(a,b float64) bool 関数が追加され、より厳しい許容誤差 4e-16 (これは倍精度浮動小数点数の精度限界に近い値) で比較を行うようになりました。
• テストの精度向上: TestAsin, TestAtan, TestExp, TestSinh, TestSqrt, TestTanh, TestHypot などのテストにおいて、Close の代わりに VeryClose を使用するように変更されています。これにより、これらの関数のテストがより厳密になり、精度向上が検証されるようになりました。
• TestLog のテストケースに math.Log(10) のテストが追加され、Ln10 との比較が行われるようになりました。

Makefile の変更点

• src/lib/math/Makefile では、exp.o のビルド順序が O2 から O1 に移動し、pow.o が O3 から O2 に移動しています。これは、依存関係の整理と、コンパイル時の最適化順序の調整を示唆しています。
• O4 というカテゴリが削除され、ビルドプロセスが簡素化されています。

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

src/lib/math/exp.go

// The original C code, the long comment, and the constants
// below are from FreeBSD's /usr/src/lib/msun/src/e_exp.c
// and came with this notice.  The go code is a simplified
// version of the original C.
// ... (中略：FreeBSD msun e_exp.cからの詳細なアルゴリズム説明) ...

export const (
	Ln2				= 0.693147180559945309417232121458176568;
	HalfLn2			= 0.346573590279972654708616060729088284;

	Ln2Hi	= 6.93147180369123816490e-01;
	Ln2Lo	= 1.90821492927058770002e-10;
	Log2e	= 1.44269504088896338700e+00;

	P1   =  1.66666666666666019037e-01; /* 0x3FC55555; 0x5555553E */
	P2   = -2.77777777770155933842e-03; /* 0xBF66C16C; 0x16BEBD93 */
	P3   =  6.61375632143793436117e-05; /* 0x3F11566A; 0xAF25DE2C */
	P4   = -1.65339022054652515390e-06; /* 0xBEBBBD41; 0xC5D26BF1 */
	P5   =  4.13813679705723846039e-08; /* 0x3E663769; 0x72BEA4D0 */

	Overflow	= 7.09782712893383973096e+02;
	Underflow	= -7.45133219101941108420e+02;
	NearZero	= 1.0/(1<<28);		// 2^-28
)

export func Exp(x float64) float64 {
	// special cases
	switch {
	case sys.isNaN(x) || sys.isInf(x, 1):
		return x;
	case sys.isInf(x, -1):
		return 0;
	case x > Overflow:
		return sys.Inf(1);
	case x < Underflow:
		return 0;
	case -NearZero < x && x < NearZero:
		return 1;
	}

	// reduce; computed as r = hi - lo for extra precision.
	var k int;
	switch {
	case x < 0:
		k = int(Log2e*x - 0.5);
	case x > 0:
		k = int(Log2e*x + 0.5);
	}
	hi := x - float64(k)*Ln2Hi;
	lo := float64(k)*Ln2Lo;
	r := hi - lo;

	// compute
	t := r * r;
	c := r - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
	y := 1 - ((lo - (r*c)/(2-c)) - hi);
	// TODO(rsc): make sure sys.ldexp can handle boundary k
	return sys.ldexp(y, k);
}

src/lib/math/log.go

// The original C code, the long comment, and the constants
// below are from FreeBSD's /usr/src/lib/msun/src/e_log.c
// and came with this notice.  The go code is a simpler
// version of the original C.
// ... (中略：FreeBSD msun e_log.cからの詳細なアルゴリズム説明) ...

const (
	Ln2Hi = 6.93147180369123816490e-01;	/* 3fe62e42 fee00000 */
	Ln2Lo = 1.90821492927058770002e-10;	/* 3dea39ef 35793c76 */
	Lg1 = 6.666666666666735130e-01;  /* 3FE55555 55555593 */
	Lg2 = 3.999999999940941908e-01;  /* 3FD99999 9997FA04 */
	Lg3 = 2.857142874366239149e-01;  /* 3FD24924 94229359 */
	Lg4 = 2.222219843214978396e-01;  /* 3FCC71C5 1D8E78AF */
	Lg5 = 1.818357216161805012e-01;  /* 3FC74664 96CB03DE */
	Lg6 = 1.531383769920937332e-01;  /* 3FC39A09 D078C69F */
	Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */

	Two54 = 1<<54;				// 2^54
	TwoM20 = 1.0/(1<<20);		// 2^-20
	TwoM1022 = 2.2250738585072014e-308;	// 2^-1022
	Sqrt2 = 1.41421356237309504880168872420969808;
)

export func Log(x float64) float64 {
	// special cases
	switch {
	case sys.isNaN(x) || sys.isInf(x, 1):
		return x;
	case x < 0:
		return sys.NaN();
	case x == 0:
		return sys.Inf(-1);
	}

	// reduce
	f1, ki := sys.frexp(x);
	if f1 < Sqrt2/2 {
		f1 *= 2;
		ki--;
	}
	f := f1 - 1;
	k := float64(ki);

	// compute
	s := f/(2+f);
	s2 := s*s;
	s4 := s2*s2;
	t1 := s2*(Lg1 + s4*(Lg3 + s4*(Lg5 + s4*Lg7)));
	t2 := s4*(Lg2 + s4*(Lg4 + s4*Lg6));
	R :=  t1 + t2;
	hfsq := 0.5*f*f;
	return k*Ln2Hi - ((hfsq-(s*(hfsq+R)+k*Ln2Lo)) - f);
}

src/lib/math/pow.go

export func Pow(x, y float64) float64 {
	// TODO: x or y NaN, ±Inf, maybe ±0.
	switch {
	case y == 0:
		return 1;
	case y == 1:
		return x;
	case x == 0 && y > 0:
		return 0;
	case x == 0 && y < 0:
		return sys.Inf(1);
	case y == 0.5:
		return Sqrt(x);
	case y == -0.5:
		return 1 / Sqrt(x);
	}

	absy := y;
	flip := false;
	if absy < 0 {
		absy = -absy;
		flip = true;
	}
	yi, yf := sys.modf(absy);
	if yf != 0 && x < 0 {
		return sys.NaN();
	}
	if yi >= 1<<63 {
		return Exp(y * Log(x));
	}

	ans := float64(1);

	// ans *= x^yf
	if yf != 0 {
		if yf > 0.5 {
			yf--;
			yi++;
		}
		ans = Exp(yf * Log(x));
	}

	// ans *= x^yi
	// by multiplying in successive squarings
	// of x according to bits of yi.
	// accumulate powers of two into exp.
	// will still have to do ans *= 2^exp later.
	x1, xe := sys.frexp(x);
	exp := 0;
	if i := int64(yi); i != 0 {
		for {
			if i&1 == 1 {
				ans *= x1;
				exp += xe;
			}
			i >>= 1;
			if i == 0 {
				break;
			}
			x1 *= x1;
			xe <<= 1;
			if x1 < .5 {
				x1 += x1;
				xe--;
			}
		}
	}

	// ans *= 2^exp
	// if flip { ans = 1 / ans }
	// but in the opposite order
	if flip {
		ans = 1 / ans;
		exp = -exp;
	}
	return sys.ldexp(ans, exp);
}

src/lib/math/all_test.go (旧 test.go)

func Tolerance(a,b,e float64) bool {
	d := a-b;
	if d < 0 {
		d = -d;
	}

	if a != 0 {
		e = e*a;
		if e < 0 {
			e = -e;
		}
	}
	return d < e;
}
func Close(a,b float64) bool {
	return Tolerance(a, b, 1e-14);
}
func VeryClose(a,b float64) bool {
	return Tolerance(a, b, 4e-16);
}

// ... (各テスト関数でCloseからVeryCloseへの変更) ...

export func TestLog(t *testing.T) {
	for i := 0; i < len(vf); i++ {
		a := math.Fabs(vf[i]);
		// 変更前: if f := math.Log(a); !Close(log[i], f) {
		// 変更後: if f := math.Log(a); log[i] != f {
		if f := math.Log(a); log[i] != f {
			t.Errorf("math.Log(%g) = %g, want %g\n", a, f, log[i]);
		}
	}
	const Ln10 = 2.30258509299404568401799145468436421;
	if f := math.Log(10); f != Ln10 {
		t.Errorf("math.Log(%g) = %g, want %g\n", 10, f, Ln10);
	}
}

コアとなるコードの解説

exp.go と log.go の精度向上

これらのファイルでは、FreeBSDの msun ライブラリから導入されたアルゴリズムが核となっています。これは、単なる多項式近似ではなく、引数削減、高精度な定数の使用（Ln2Hi と Ln2Lo のように ln2 を2つの浮動小数点数に分割する手法など）、そしてRemesアルゴリズムによって導出された最適化された多項式係数（P1〜P5、Lg1〜Lg7）を組み合わせることで、IEEE 754倍精度浮動小数点数の精度限界に近い「1 ULP未満」のエラーを実現しています。

特に注目すべきは、exp 関数における y := 1 - ((lo - (r*c)/(2-c)) - hi) のような計算式です。これは、浮動小数点演算の丸め誤差を最小限に抑えるための工夫であり、中間結果の精度を保ちながら最終的な高精度を実現しています。

pow.go の堅牢性向上

Pow 関数は、x^y の計算において、様々な特殊ケース（y=0、y=1、x=0、y=±0.5など）を明示的に処理することで、これらの一般的な入力に対する正確性と効率を向上させています。 負の底 x と非整数の指数 y の組み合わせは数学的に複素数を返すため、実数演算を行う math パッケージでは NaN を返すのが適切です。このコミットではその挙動が追加されています。 また、Exp(y * Log(x)) へのフォールバックは、x^y = exp(y * log(x)) という数学的恒等式を利用したもので、直接計算が困難な大きな指数や、浮動小数点数の特性上、対数と指数に変換した方が精度を保ちやすい場合に用いられます。 整数べき乗の計算に sys.frexp とビット演算を用いた「逐次二乗法」のようなアプローチを採用している点も、効率的かつ正確な計算のための重要な変更です。

all_test.go による厳密な検証

Tolerance、Close、VeryClose 関数の導入は、テストの厳密性を大幅に向上させています。特に VeryClose で使用される 4e-16 という許容誤差は、倍精度浮動小数点数の精度限界（約15〜17桁の10進数精度）に非常に近い値であり、これにより Log、Exp、Pow などの関数の精度向上が実際に達成されていることを、より高い信頼性で検証できるようになりました。TestLog で math.Log(10) が Ln10 と厳密に比較されるようになったのも、特定の重要な定数に対する精度保証を強化する意図が見られます。

関連リンク

• IEEE 754 浮動小数点数標準: https://ja.wikipedia.org/wiki/IEEE_754%E6%B5%AE%E5%8B%95%E5%B0%8F%E6%95%B0%E7%82%B9%E6%A8%99%E6%BA%96
• ULP (Unit in the Last Place): https://en.wikipedia.org/wiki/Unit_in_the_last_place
• Remesアルゴリズム: https://en.wikipedia.org/wiki/Remez_algorithm

参考にした情報源リンク

• FreeBSD msun ライブラリのソースコード (特に e_exp.c と e_log.c):
  • e_exp.c: https://github.com/freebsd/freebsd-src/blob/main/lib/msun/src/e_exp.c
  • e_log.c: https://github.com/freebsd/freebsd-src/blob/main/lib/msun/src/e_log.c
• Hart, Cheney, Lawson の数値計算に関する書籍 (引用されている係数の出典):
  • "Computer Approximations" by John F. Hart, E. W. Cheney, Charles L. Lawson, Hans J. Maehly, Charles K. Mesztenyi, William L. Rice, Henry C. Thacher Jr., Christopher Witzgall. (この書籍は数値近似に関する古典的な文献であり、多くの数学ライブラリの基礎となっています。)
• Go言語の math パッケージのドキュメント: https://pkg.go.dev/math
• Go言語の sys パッケージ (内部パッケージのため直接参照は難しいが、frexp や ldexp の挙動に関する情報): https://pkg.go.dev/runtime/internal/sys (これはGoの内部パッケージであり、直接利用することは推奨されません。しかし、その機能は math パッケージを通じて利用されています。)
• 浮動小数点数の精度に関する一般的な情報源 (例: Wikipediaの「倍精度浮動小数点数」など)