libgo/go/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
  17 type matrix struct {
  18         n, m int
  19         a    []*Rat
  20 }
  21
  22
  23 func (a *matrix) at(i, j int) *Rat {
  24         if !(0 <= i && i < a.n && 0 <= j && j < a.m) {
  25                 panic("index out of range")
  26         }
  27         return a.a[i*a.m+j]
  28 }
  29
  30
  31 func (a *matrix) set(i, j int, x *Rat) {
  32         if !(0 <= i && i < a.n && 0 <= j && j < a.m) {
  33                 panic("index out of range")
  34         }
  35         a.a[i*a.m+j] = x
  36 }
  37
  38
  39 func newMatrix(n, m int) *matrix {
  40         if !(0 <= n && 0 <= m) {
  41                 panic("illegal matrix")
  42         }
  43         a := new(matrix)
  44         a.n = n
  45         a.m = m
  46         a.a = make([]*Rat, n*m)
  47         return a
  48 }
  49
  50
  51 func newUnit(n int) *matrix {
  52         a := newMatrix(n, n)
  53         for i := 0; i < n; i++ {
  54                 for j := 0; j < n; j++ {
  55                         x := NewRat(0, 1)
  56                         if i == j {
  57                                 x.SetInt64(1)
  58                         }
  59                         a.set(i, j, x)
  60                 }
  61         }
  62         return a
  63 }
  64
  65
  66 func newHilbert(n int) *matrix {
  67         a := newMatrix(n, n)
  68         for i := 0; i < n; i++ {
  69                 for j := 0; j < n; j++ {
  70                         a.set(i, j, NewRat(1, int64(i+j+1)))
  71                 }
  72         }
  73         return a
  74 }
  75
  76
  77 func newInverseHilbert(n int) *matrix {
  78         a := newMatrix(n, n)
  79         for i := 0; i < n; i++ {
  80                 for j := 0; j < n; j++ {
  81                         x1 := new(Rat).SetInt64(int64(i + j + 1))
  82                         x2 := new(Rat).SetInt(new(Int).Binomial(int64(n+i), int64(n-j-1)))
  83                         x3 := new(Rat).SetInt(new(Int).Binomial(int64(n+j), int64(n-i-1)))
  84                         x4 := new(Rat).SetInt(new(Int).Binomial(int64(i+j), int64(i)))
  85
  86                         x1.Mul(x1, x2)
  87                         x1.Mul(x1, x3)
  88                         x1.Mul(x1, x4)
  89                         x1.Mul(x1, x4)
  90
  91                         if (i+j)&1 != 0 {
  92                                 x1.Neg(x1)
  93                         }
  94
  95                         a.set(i, j, x1)
  96                 }
  97         }
  98         return a
  99 }
 100
 101
 102 func (a *matrix) mul(b *matrix) *matrix {
 103         if a.m != b.n {
 104                 panic("illegal matrix multiply")
 105         }
 106         c := newMatrix(a.n, b.m)
 107         for i := 0; i < c.n; i++ {
 108                 for j := 0; j < c.m; j++ {
 109                         x := NewRat(0, 1)
 110                         for k := 0; k < a.m; k++ {
 111                                 x.Add(x, new(Rat).Mul(a.at(i, k), b.at(k, j)))
 112                         }
 113                         c.set(i, j, x)
 114                 }
 115         }
 116         return c
 117 }
 118
 119
 120 func (a *matrix) eql(b *matrix) bool {
 121         if a.n != b.n || a.m != b.m {
 122                 return false
 123         }
 124         for i := 0; i < a.n; i++ {
 125                 for j := 0; j < a.m; j++ {
 126                         if a.at(i, j).Cmp(b.at(i, j)) != 0 {
 127                                 return false
 128                         }
 129                 }
 130         }
 131         return true
 132 }
 133
 134
 135 func (a *matrix) String() string {
 136         s := ""
 137         for i := 0; i < a.n; i++ {
 138                 for j := 0; j < a.m; j++ {
 139                         s += fmt.Sprintf("\t%s", a.at(i, j))
 140                 }
 141                 s += "\n"
 142         }
 143         return s
 144 }
 145
 146
 147 func doHilbert(t *testing.T, n int) {
 148         a := newHilbert(n)
 149         b := newInverseHilbert(n)
 150         I := newUnit(n)
 151         ab := a.mul(b)
 152         if !ab.eql(I) {
 153                 if t == nil {
 154                         panic("Hilbert failed")
 155                 }
 156                 t.Errorf("a   = %s\n", a)
 157                 t.Errorf("b   = %s\n", b)
 158                 t.Errorf("a*b = %s\n", ab)
 159                 t.Errorf("I   = %s\n", I)
 160         }
 161 }
 162
 163
 164 func TestHilbert(t *testing.T) {
 165         doHilbert(t, 10)
 166 }
 167
 168
 169 func BenchmarkHilbert(b *testing.B) {
 170         for i := 0; i < b.N; i++ {
 171                 doHilbert(nil, 10)
 172         }
 173 }