hilbert_test.go

Documentation: math/big

		 1  // Copyright 2009 The Go Authors. All rights reserved.
		 2  // Use of this source code is governed by a BSD-style
		 3  // license that can be found in the LICENSE file.
		 4  
		 5  // A little test program and benchmark for rational arithmetics.
		 6  // Computes a Hilbert matrix, its inverse, multiplies them
		 7  // and verifies that the product is the identity matrix.
		 8  
		 9  package big
		10  
		11  import (
		12  	"fmt"
		13  	"testing"
		14  )
		15  
		16  type matrix struct {
		17  	n, m int
		18  	a		[]*Rat
		19  }
		20  
		21  func (a *matrix) at(i, j int) *Rat {
		22  	if !(0 <= i && i < a.n && 0 <= j && j < a.m) {
		23  		panic("index out of range")
		24  	}
		25  	return a.a[i*a.m+j]
		26  }
		27  
		28  func (a *matrix) set(i, j int, x *Rat) {
		29  	if !(0 <= i && i < a.n && 0 <= j && j < a.m) {
		30  		panic("index out of range")
		31  	}
		32  	a.a[i*a.m+j] = x
		33  }
		34  
		35  func newMatrix(n, m int) *matrix {
		36  	if !(0 <= n && 0 <= m) {
		37  		panic("illegal matrix")
		38  	}
		39  	a := new(matrix)
		40  	a.n = n
		41  	a.m = m
		42  	a.a = make([]*Rat, n*m)
		43  	return a
		44  }
		45  
		46  func newUnit(n int) *matrix {
		47  	a := newMatrix(n, n)
		48  	for i := 0; i < n; i++ {
		49  		for j := 0; j < n; j++ {
		50  			x := NewRat(0, 1)
		51  			if i == j {
		52  				x.SetInt64(1)
		53  			}
		54  			a.set(i, j, x)
		55  		}
		56  	}
		57  	return a
		58  }
		59  
		60  func newHilbert(n int) *matrix {
		61  	a := newMatrix(n, n)
		62  	for i := 0; i < n; i++ {
		63  		for j := 0; j < n; j++ {
		64  			a.set(i, j, NewRat(1, int64(i+j+1)))
		65  		}
		66  	}
		67  	return a
		68  }
		69  
		70  func newInverseHilbert(n int) *matrix {
		71  	a := newMatrix(n, n)
		72  	for i := 0; i < n; i++ {
		73  		for j := 0; j < n; j++ {
		74  			x1 := new(Rat).SetInt64(int64(i + j + 1))
		75  			x2 := new(Rat).SetInt(new(Int).Binomial(int64(n+i), int64(n-j-1)))
		76  			x3 := new(Rat).SetInt(new(Int).Binomial(int64(n+j), int64(n-i-1)))
		77  			x4 := new(Rat).SetInt(new(Int).Binomial(int64(i+j), int64(i)))
		78  
		79  			x1.Mul(x1, x2)
		80  			x1.Mul(x1, x3)
		81  			x1.Mul(x1, x4)
		82  			x1.Mul(x1, x4)
		83  
		84  			if (i+j)&1 != 0 {
		85  				x1.Neg(x1)
		86  			}
		87  
		88  			a.set(i, j, x1)
		89  		}
		90  	}
		91  	return a
		92  }
		93  
		94  func (a *matrix) mul(b *matrix) *matrix {
		95  	if a.m != b.n {
		96  		panic("illegal matrix multiply")
		97  	}
		98  	c := newMatrix(a.n, b.m)
		99  	for i := 0; i < c.n; i++ {
	 100  		for j := 0; j < c.m; j++ {
	 101  			x := NewRat(0, 1)
	 102  			for k := 0; k < a.m; k++ {
	 103  				x.Add(x, new(Rat).Mul(a.at(i, k), b.at(k, j)))
	 104  			}
	 105  			c.set(i, j, x)
	 106  		}
	 107  	}
	 108  	return c
	 109  }
	 110  
	 111  func (a *matrix) eql(b *matrix) bool {
	 112  	if a.n != b.n || a.m != b.m {
	 113  		return false
	 114  	}
	 115  	for i := 0; i < a.n; i++ {
	 116  		for j := 0; j < a.m; j++ {
	 117  			if a.at(i, j).Cmp(b.at(i, j)) != 0 {
	 118  				return false
	 119  			}
	 120  		}
	 121  	}
	 122  	return true
	 123  }
	 124  
	 125  func (a *matrix) String() string {
	 126  	s := ""
	 127  	for i := 0; i < a.n; i++ {
	 128  		for j := 0; j < a.m; j++ {
	 129  			s += fmt.Sprintf("\t%s", a.at(i, j))
	 130  		}
	 131  		s += "\n"
	 132  	}
	 133  	return s
	 134  }
	 135  
	 136  func doHilbert(t *testing.T, n int) {
	 137  	a := newHilbert(n)
	 138  	b := newInverseHilbert(n)
	 139  	I := newUnit(n)
	 140  	ab := a.mul(b)
	 141  	if !ab.eql(I) {
	 142  		if t == nil {
	 143  			panic("Hilbert failed")
	 144  		}
	 145  		t.Errorf("a	 = %s\n", a)
	 146  		t.Errorf("b	 = %s\n", b)
	 147  		t.Errorf("a*b = %s\n", ab)
	 148  		t.Errorf("I	 = %s\n", I)
	 149  	}
	 150  }
	 151  
	 152  func TestHilbert(t *testing.T) {
	 153  	doHilbert(t, 10)
	 154  }
	 155  
	 156  func BenchmarkHilbert(b *testing.B) {
	 157  	for i := 0; i < b.N; i++ {
	 158  		doHilbert(nil, 10)
	 159  	}
	 160  }
	 161
View as plain text