...
Run Format

Source file src/math/big/hilbert_test.go

Documentation: math/big

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

View as plain text