Subversion Repositories openrisc

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