diff --git a/src/pkg/big/hilbert_test.go b/src/pkg/big/hilbert_test.go new file mode 100644 index 00000000000..66a21214d22 --- /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 19c3d88f737..dc066580a19 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 07db353240c..00000000000 --- 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"); - } -}