Home | History | Annotate | Download | only in 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