blob: 1a84341b3c0b63bfa7040c5631c2af9c8670578a [file] [log] [blame]
 // 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) } }