123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160 |
- // 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)
- }
- }
|