Keyboard shortcuts

Press ← or → to navigate between chapters

Press S or / to search in the book

Press ? to show this help

Press Esc to hide this help

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

このコミットは、Go言語の初期段階におけるbignumパッケージのRational型（有理数）の算術演算をテストするために追加されたtest/hilbert.goファイルに関するものです。このテストは、ヒルベルト行列とその逆行列の積が単位行列になることを検証することで、有理数演算の正確性を確認する「レクリエーション的なプログラミング演習」として位置づけられています。

コミット

commit ce164403dab6d5f493ce155ad206769a39bc34e6
Author: Robert Griesemer <gri@golang.org>
Date:   Thu Nov 6 14:23:49 2008 -0800

    A recreational programming exercise:
    
    Multiplication of a Hilbert matrix with its inverse using
    Bignum.Rationals as a test case for rational arithmetic.
    
    R=r
    OCL=18706
    CL=18706
---
 test/hilbert.go | 167 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 167 insertions(+)

diff --git a/test/hilbert.go b/test/hilbert.go
new file mode 100644
index 0000000000..275b11997d
--- /dev/null
+++ b/test/hilbert.go
@@ -0,0 +1,167 @@
+// $G $D/$F.go && $L $F.$A && ./$A.out
+//
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+//
+// A little test program for rational arithmetics.
+// Computes a Hilbert matrix, its inverse, multiplies them
+// and verifies that the product is the identity matrix.
+
+package main
+
+import Big "bignum"
+import Fmt "fmt"
+
+
+func assert(p bool) {
+	if !p {
+		panic("assert failed");
+	}
+}
+
+
+var (
+	Zero = Big.Rat(0, 1);
+	One = Big.Rat(1, 1);
+)
+
+
+type Matrix struct {
+	n, m int;
+	a *[]*Big.Rational;
+}
+
+
+func (a *Matrix) at(i, j int) *Big.Rational {
+	assert(0 <= i && i < a.n && 0 <= j && j < a.m);
+	return a.a[i*a.m + j];
+}
+
+
+func (a *Matrix) set(i, j int, x *Big.Rational) {
+	assert(0 <= i && i < a.n && 0 <= j && j < a.m);
+	a.a[i*a.m + j] = x;
+}
+
+
+func NewMatrix(n, m int) *Matrix {
+	assert(0 <= n && 0 <= m);
+	a := new(Matrix);
+	a.n = n;
+	a.m = m;
+	a.a = new([]*Big.Rational, n*m);
+	return a;
+}
+
+
+func NewUnit(n int) *Matrix {
+	a := NewMatrix(n, n);
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			x := Zero;
+			if i == j {
+				x = One;
+			}
+			a.set(i, j, x);
+		}
+	}
+	return a;
+}
+
+
+func NewHilbert(n int) *Matrix {
+	a := NewMatrix(n, n);
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			x := Big.Rat(1, i + j + 1);
+			a.set(i, j, x);
+		}
+	}
+	return a;
+}
+
+
+func MakeRat(x *Big.Natural) *Big.Rational {
+	return Big.MakeRat(Big.MakeInt(false, x), Big.Nat(1));
+}
+
+
+func NewInverseHilbert(n int) *Matrix {
+	a := NewMatrix(n, n);
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			x0 := One;
+			if (i+j)&1 != 0 {
+				x0 = x0.Neg();
+			}
+			x1 := Big.Rat(i + j + 1, 1);
+			x2 := MakeRat(Big.Binomial(uint(n+i), uint(n-j-1)));
+			x3 := MakeRat(Big.Binomial(uint(n+j), uint(n-i-1)));
+			x4 := MakeRat(Big.Binomial(uint(i+j), uint(i)));
+			x4 = x4.Mul(x4);
+			a.set(i, j, x0.Mul(x1).Mul(x2).Mul(x3).Mul(x4));
+		}
+	}
+	return a;
+}
+
+
+func (a *Matrix) Mul(b *Matrix) *Matrix {
+	assert(a.m == b.n);
+	c := NewMatrix(a.n, b.m);
+	for i := 0; i < c.n; i++ {
+		for j := 0; j < c.m; j++ {
+			x := Zero;
+			for k := 0; k < a.m; k++ {
+				x = x.Add(a.at(i, k).Mul(b.at(k, j)));
+			}
+			c.set(i, j, x);
+		}
+	}
+	return c;
+}
+
+
+func (a *Matrix) Eql(b *Matrix) bool {
+	if a.n != b.n || a.m != b.m {
+		return false;
+	}
+	for i := 0; i < a.n; i++ {
+		for j := 0; j < a.m; j++ {
+			if a.at(i, j).Cmp(b.at(i,j)) != 0 {
+				return false;
+			}
+		}
+	}
+	return true;
+}
+
+
+func (a *Matrix) String() string {
+	s := "";
+	for i := 0; i < a.n; i++ {
+		for j := 0; j < a.m; j++ {
+			x := a.at(i, j);  // BUG 6g bug
+			s += Fmt.sprintf("\t%s", x);
+		}
+		s += "\n";
+	}
+	return s;
+}
+
+
+func main() {
+	n := 10;
+	a := NewHilbert(n);
+	b := NewInverseHilbert(n);
+	I := NewUnit(n);
+	ab := a.Mul(b);
+	if !ab.Eql(I) {
+		Fmt.println("a =", a);
+		Fmt.println("b =", b);
+		Fmt.println("a*b =", ab);
+		Fmt.println("I =", I);
+		panic("FAILED");
+	}
+}

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

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

元コミット内容

A recreational programming exercise:

Multiplication of a Hilbert matrix with its inverse using
Bignum.Rationals as a test case for rational arithmetic.

R=r
OCL=18706
CL=18706

変更の背景

このコミットは、Go言語の初期開発段階において、bignumパッケージで実装された有理数（Rational）演算の正確性と信頼性を検証するために導入されました。特に、浮動小数点数演算では精度問題が生じやすい数値計算において、正確な有理数演算が正しく機能するかを確認することが目的でした。

ヒルベルト行列は、その逆行列の要素が非常に大きな整数になるため、数値的に非常に不安定（悪条件）な行列として知られています。このような悪条件の行列に対する演算は、浮動小数点数では容易に丸め誤差が蓄積し、正確な結果を得ることが困難です。しかし、有理数演算を用いることで、無限の精度で計算が可能となり、理論的には正確な結果が得られるはずです。

この「レクリエーション的なプログラミング演習」という表現は、単なるテストケースとしてだけでなく、Go言語の新しいbignumパッケージの能力を示すデモンストレーションとしての側面も持っていたことを示唆しています。

前提知識の解説

ヒルベルト行列 (Hilbert Matrix)

ヒルベルト行列 $H$ は、要素 $H_{ij}$ が $1/(i+j-1)$ で与えられる正方行列です（1-indexedの場合）。例えば、3x3のヒルベルト行列は以下のようになります。

$$ H = \begin{pmatrix} 1 & 1/2 & 1/3 \ 1/2 & 1/3 & 1/4 \ 1/3 & 1/4 & 1/5 \end{pmatrix} $$

この行列は、その性質上、非常に悪条件（ill-conditioned）であることが知られています。これは、行列式が非常に小さく、逆行列の要素が非常に大きくなる傾向があるため、数値計算において小さな誤差が結果に大きく影響することを意味します。

ヒルベルト行列の逆行列 (Inverse Hilbert Matrix)

ヒルベルト行列の逆行列 $H^{-1}$ の要素 $(H^{-1})_{ij}$ は、以下の式で与えられます。

$$ (H^{-1})_{ij} = (-1)^{i+j} (i+j-1) \binom{n+i-1}{n-j} \binom{n+j-1}{n-i} \left(\binom{i+j-2}{i-1}\right)^2 $$

ここで、$\binom{n}{k}$ は二項係数（"n choose k"）を表します。この式からわかるように、逆行列の要素はすべて整数であり、行列のサイズ $n$ が大きくなると、その値は非常に急速に増大します。

有理数演算 (Rational Arithmetic)

有理数とは、整数 $p$ とゼロではない整数 $q$ を用いて $p/q$ の分数として表せる数のことです。有理数演算は、浮動小数点数演算とは異なり、計算中に丸め誤差が発生しません。これにより、理論的には無限の精度で計算を進めることができます。

Go言語のbignumパッケージは、任意精度の整数（Big.Int）や有理数（Big.Rational）を扱う機能を提供します。これにより、標準のfloat64型では表現できないような非常に大きな数や、正確な分数計算が必要な場合に利用されます。

行列の積 (Matrix Multiplication)

2つの行列 $A$ と $B$ の積 $C = AB$ は、行列 $A$ の列数と行列 $B$ の行数が等しい場合に定義されます。結果行列 $C$ の要素 $C_{ij}$ は、行列 $A$ の $i$ 行目と行列 $B$ の $j$ 列目の要素の積の和として計算されます。

$$ C_{ij} = \sum_{k=1}^{m} A_{ik} B_{kj} $$

ここで、$A$ は $n \times m$ 行列、$B$ は $m \times p$ 行列、$C$ は $n \times p$ 行列です。

単位行列 (Identity Matrix)

単位行列 $I$ は、主対角線上の要素がすべて1で、それ以外の要素がすべて0である正方行列です。任意の正方行列 $A$ に対して、$AI = IA = A$ が成り立ちます。行列とその逆行列の積は、常に単位行列になります。

$$ I = \begin{pmatrix} 1 & 0 & 0 \ 0 & 1 & 0 \ 0 & 0 & 1 \end{pmatrix} $$

技術的詳細

このコミットで追加されたtest/hilbert.goは、Go言語のbignumパッケージを利用して、ヒルベルト行列とその逆行列を生成し、それらの積が単位行列になることを検証するプログラムです。

プログラムは以下の主要な構造体と関数で構成されています。

  • Matrix構造体: 行列を表現するための構造体で、行数n、列数m、そしてBig.Rational型のポインタのスライスa（行列の要素を格納）を持ちます。
  • assert(p bool)関数: 条件pが偽の場合にパニックを発生させるシンプルなアサート関数です。テストの途中で予期せぬ状態になった場合にプログラムを停止させます。
  • ZeroとOne: Big.Rat(0, 1)とBig.Rat(1, 1)として定義された有理数のゼロとイチの定数です。
  • at(i, j)とset(i, j, x)メソッド: Matrix構造体の要素にアクセスしたり設定したりするためのヘルパーメソッドです。
  • NewMatrix(n, m int)関数: 指定されたサイズの新しい行列を初期化して返します。
  • NewUnit(n int)関数: 指定されたサイズの単位行列を生成して返します。
  • NewHilbert(n int)関数: 指定されたサイズのヒルベルト行列を生成して返します。要素はBig.Rat(1, i + j + 1)で計算されます。
  • MakeRat(x *Big.Natural)関数: Big.Natural型（符号なし任意精度整数）をBig.Rational型に変換するヘルパー関数です。分母を1として有理数を生成します。
  • NewInverseHilbert(n int)関数: 指定されたサイズのヒルベルト行列の逆行列を生成して返します。この関数は、前述の逆行列の公式を実装しており、Big.Binomial（二項係数）関数をbignumパッケージから利用しています。
    • x0: $(-1)^{i+j}$ の部分を計算します。i+jが奇数なら負の1、偶数なら正の1です。
    • x1: $(i+j+1)$ の部分を計算します。
    • x2: $\binom{n+i}{n-j-1}$ の部分を計算します。
    • x3: $\binom{n+j}{n-i-1}$ の部分を計算します。
    • x4: $\left(\binom{i+j}{i}\right)^2$ の部分を計算します。
    • これらすべての項をBig.Rationalの乗算メソッドMulで掛け合わせることで、逆行列の要素を計算しています。
  • Mul(b *Matrix)メソッド: 行列の積を計算します。標準的な行列乗算のアルゴリズム（$C_{ij} = \sum_{k} A_{ik} B_{kj}$）を実装しています。
  • Eql(b *Matrix)メソッド: 2つの行列が等しいかどうかを比較します。各要素がBig.Rationalの比較メソッドCmpで比較されます。
  • String()メソッド: 行列を文字列として整形して出力します。デバッグ目的で使用されます。コメントに// BUG 6g bugとあるのは、当時のGoコンパイラ（6g）のバグを示唆しており、特定の状況下でx := a.at(i, j)のような代入が正しく機能しない可能性があったことを示しています。
  • main()関数: プログラムのエントリポイントです。
    • n := 10; で行列のサイズを10に設定します。
    • NewHilbert(n)でヒルベルト行列aを生成します。
    • NewInverseHilbert(n)で逆ヒルベルト行列bを生成します。
    • NewUnit(n)で単位行列Iを生成します。
    • a.Mul(b)でaとbの積abを計算します。
    • ab.Eql(I)でabが単位行列Iと等しいか検証します。
    • もし等しくなければ、各行列の内容を出力し、panic("FAILED")でプログラムを終了させます。

このテストは、bignum.Rationalが正確な有理数演算を提供し、悪条件の行列計算においても期待通りの結果（積が単位行列になること）を導き出せることを確認するために設計されています。

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

このコミットでは、test/hilbert.goという新しいファイルが追加されました。このファイル全体が変更箇所であり、167行のコードが新規に挿入されています。

コアとなるコードの解説

コアとなるコードは、NewInverseHilbert関数とMulメソッド、そしてmain関数における検証ロジックです。

NewInverseHilbert(n int) *Matrix

この関数は、ヒルベルト行列の逆行列を正確に計算する部分であり、bignumパッケージのBig.Binomial（二項係数）やMul（乗算）、Neg（符号反転）といった有理数演算の機能がフル活用されています。特に、逆行列の要素が非常に大きな整数になるため、任意精度有理数型であるBig.Rationalが不可欠です。

func NewInverseHilbert(n int) *Matrix {
	a := NewMatrix(n, n);
	for i := 0; i < n; i++ {
		for j := 0; j < n; j++ {
			x0 := One;
			if (i+j)&1 != 0 { // (-1)^(i+j) の計算
				x0 = x0.Neg();
			}
			x1 := Big.Rat(i + j + 1, 1); // (i+j+1) の計算
			// 以下の x2, x3, x4 は二項係数を含む項の計算
			x2 := MakeRat(Big.Binomial(uint(n+i), uint(n-j-1)));
			x3 := MakeRat(Big.Binomial(uint(n+j), uint(n-i-1)));
			x4 := MakeRat(Big.Binomial(uint(i+j), uint(i)));
			x4 = x4.Mul(x4); // (二項係数)^2 の計算
			a.set(i, j, x0.Mul(x1).Mul(x2).Mul(x3).Mul(x4)); // 全ての項を乗算
		}
	}
	return a;
}

(a *Matrix) Mul(b *Matrix) *Matrix

このメソッドは、行列の乗算を実装しており、Big.Rational型の加算Addと乗算Mulメソッドを繰り返し使用します。これにより、中間結果もすべて有理数として正確に保持され、最終的な積の精度が保証されます。

func (a *Matrix) Mul(b *Matrix) *Matrix {
	assert(a.m == b.n); // 行列乗算の条件チェック
	c := NewMatrix(a.n, b.m); // 結果行列の初期化
	for i := 0; i < c.n; i++ {
		for j := 0; j < c.m; j++ {
			x := Zero; // 各要素の初期値はゼロ
			for k := 0; k < a.m; k++ {
				// A_ik * B_kj を計算し、x に加算
				x = x.Add(a.at(i, k).Mul(b.at(k, j)));
			}
			c.set(i, j, x); // 結果行列の要素を設定
		}
	}
	return c;
}

main()関数における検証

main関数では、生成されたヒルベルト行列とその逆行列の積が、理論的に単位行列と一致するかどうかをEqlメソッドで厳密に比較しています。この比較が成功することで、bignum.Rationalパッケージが提供する有理数演算が、複雑な数値計算においても正確であることを証明しています。

func main() {
	n := 10; // 行列のサイズ
	a := NewHilbert(n); // ヒルベルト行列
	b := NewInverseHilbert(n); // 逆ヒルベルト行列
	I := NewUnit(n); // 単位行列
	ab := a.Mul(b); // 積を計算
	if !ab.Eql(I) { // 単位行列との比較
		// 失敗した場合、詳細を出力してパニック
		Fmt.println("a =", a);
		Fmt.println("b =", b);
		Fmt.println("a*b =", ab);
		Fmt.println("I =", I);
		panic("FAILED");
	}
}

これらのコードは、Go言語のbignumパッケージが、数値的に不安定な問題に対しても高精度な計算能力を提供できることを示す、具体的なテストケースとして機能しています。

関連リンク

参考にした情報源リンク

  • Web search results for "Hilbert matrix inverse formula" (stata.com, proofwiki.org, mathworks.com, scipy.org) - ヒルベルト行列の逆行列の公式と、その数値的な不安定性に関する情報。
  • Go言語のソースコードリポジトリ（コミット履歴）
  • Go言語の公式ドキュメント（math/bigパッケージに関する情報）
  • 一般的な線形代数および数値解析の知識。