Source file src/pkg/big/hilbert_test.go

     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	}
The Go Programming Language

Source file src/pkg/big/hilbert_test.go
