161 lines
2.9 KiB
Go
161 lines
2.9 KiB
Go
// 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)
|
|
}
|
|
}
|