test/hilbert.go: convert to test case and benchmark for big.Rat

R=rsc
CC=golang-dev
https://golang.org/cl/1231044

3 files changed, 173 insertions, 169 deletions

diff --git a/src/pkg/big/hilbert_test.go b/src/pkg/big/hilbert_test.go
new file mode 100644
index 0000000000..66a21214d2
--- /dev/null
+++ b/src/pkg/big/hilbert_test.go
@@ -0,0 +1,173 @@
+// 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 and benchmark for rational arithmetics.
+// Computes a Hilbert matrix, its inverse, multiplies them
+// and verifies that the product is the identity matrix.
+
+package big
+
+import (
+	"fmt"
+	"testing"
+)
+
+
+type matrix struct {
+	n, m int
+	a    []*Rat
+}
+
+
+func (a *matrix) at(i, j int) *Rat {
+	if !(0 <= i && i < a.n && 0 <= j && j < a.m) {
+		panic("index out of range")
+	}
+	return a.a[i*a.m+j]
+}
+
+
+func (a *matrix) set(i, j int, x *Rat) {
+	if !(0 <= i && i < a.n && 0 <= j && j < a.m) {
+		panic("index out of range")
+	}
+	a.a[i*a.m+j] = x
+}
+
+
+func newMatrix(n, m int) *matrix {
+	if !(0 <= n && 0 <= m) {
+		panic("illegal matrix")
+	}
+	a := new(matrix)
+	a.n = n
+	a.m = m
+	a.a = make([]*Rat, 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 := NewRat(0, 1)
+			if i == j {
+				x.SetInt64(1)
+			}
+			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++ {
+			a.set(i, j, NewRat(1, int64(i+j+1)))
+		}
+	}
+	return a
+}
+
+
+func newInverseHilbert(n int) *matrix {
+	a := newMatrix(n, n)
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			x1 := new(Rat).SetInt64(int64(i + j + 1))
+			x2 := new(Rat).SetInt(new(Int).Binomial(int64(n+i), int64(n-j-1)))
+			x3 := new(Rat).SetInt(new(Int).Binomial(int64(n+j), int64(n-i-1)))
+			x4 := new(Rat).SetInt(new(Int).Binomial(int64(i+j), int64(i)))
+
+			x1.Mul(x1, x2)
+			x1.Mul(x1, x3)
+			x1.Mul(x1, x4)
+			x1.Mul(x1, x4)
+
+			if (i+j)&1 != 0 {
+				x1.Neg(x1)
+			}
+
+			a.set(i, j, x1)
+		}
+	}
+	return a
+}
+
+
+func (a *matrix) mul(b *matrix) *matrix {
+	if a.m != b.n {
+		panic("illegal matrix multiply")
+	}
+	c := newMatrix(a.n, b.m)
+	for i := 0; i < c.n; i++ {
+		for j := 0; j < c.m; j++ {
+			x := NewRat(0, 1)
+			for k := 0; k < a.m; k++ {
+				x.Add(x, new(Rat).Mul(a.at(i, k), 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++ {
+			s += fmt.Sprintf("\t%s", a.at(i, j))
+		}
+		s += "\n"
+	}
+	return s
+}
+
+
+func doHilbert(t *testing.T, n int) {
+	a := newHilbert(n)
+	b := newInverseHilbert(n)
+	I := newUnit(n)
+	ab := a.mul(b)
+	if !ab.eql(I) {
+		if t == nil {
+			panic("Hilbert failed")
+		}
+		t.Errorf("a   = %s\n", a)
+		t.Errorf("b   = %s\n", b)
+		t.Errorf("a*b = %s\n", ab)
+		t.Errorf("I   = %s\n", I)
+	}
+}
+
+
+func TestHilbert(t *testing.T) {
+	doHilbert(t, 10)
+}
+
+
+func BenchmarkHilbert(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		doHilbert(nil, 10)
+	}
+}
diff --git a/src/pkg/big/nat.go b/src/pkg/big/nat.go
index 19c3d88f73..dc066580a1 100755
--- a/src/pkg/big/nat.go
+++ b/src/pkg/big/nat.go
@@ -16,9 +16,6 @@
 // of the operands it may be overwritten (and its memory reused).
 // To enable chaining of operations, the result is also returned.
 //
-// If possible, one should use big over bignum as the latter is headed for
-// deprecation.
-//
 package big
 
 import "rand"
diff --git a/test/hilbert.go b/test/hilbert.go
deleted file mode 100644
index 07db353240..0000000000
--- a/test/hilbert.go
+++ /dev/null
@@ -1,166 +0,0 @@
-// $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 "exp/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 = make([]*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, int64(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(int64(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++ {
-			s += Fmt.Sprintf("\t%s", a.at(i, j));
-		}
-		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");
-	}
-}