// $G $D/$F.go && $L $F.$A && ./$A.out

// 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 for rational arithmetics.
// Computes a Hilbert matrix, its inverse, multiplies them
// and verifies that the product is the identity matrix.

package main

import Big "exp/bignum"
import Fmt "fmt"


func assert(p bool) {
	if !p {
		panic("assert failed");
	}
}


var (
	Zero = Big.Rat(0, 1);
	One = Big.Rat(1, 1);
)


type Matrix struct {
	n, m int;
	a []*Big.Rational;
}


func (a *Matrix) at(i, j int) *Big.Rational {
	assert(0 <= i && i < a.n && 0 <= j && j < a.m);
	return a.a[i*a.m + j];
}


func (a *Matrix) set(i, j int, x *Big.Rational) {
	assert(0 <= i && i < a.n && 0 <= j && j < a.m);
	a.a[i*a.m + j] = x;
}


func NewMatrix(n, m int) *Matrix {
	assert(0 <= n && 0 <= m);
	a := new(Matrix);
	a.n = n;
	a.m = m;
	a.a = make([]*Big.Rational, 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 := Zero;
			if i == j {
				x = One;
			}
			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++ {
			x := Big.Rat(1, int64(i + j + 1));
			a.set(i, j, x);
		}
	}
	return a;
}


func MakeRat(x Big.Natural) *Big.Rational {
	return Big.MakeRat(Big.MakeInt(false, x), Big.Nat(1));
}


func NewInverseHilbert(n int) *Matrix {
	a := NewMatrix(n, n);
	for i := 0; i < n; i++ {
		for j := 0; j < n; j++ {
			x0 := One;
			if (i+j)&1 != 0 {
				x0 = x0.Neg();
			}
			x1 := Big.Rat(int64(i + j + 1), 1);
			x2 := MakeRat(Big.Binomial(uint(n+i), uint(n-j-1)));
			x3 := MakeRat(Big.Binomial(uint(n+j), uint(n-i-1)));
			x4 := MakeRat(Big.Binomial(uint(i+j), uint(i)));
			x4 = x4.Mul(x4);
			a.set(i, j, x0.Mul(x1).Mul(x2).Mul(x3).Mul(x4));
		}
	}
	return a;
}


func (a *Matrix) Mul(b *Matrix) *Matrix {
	assert(a.m == b.n);
	c := NewMatrix(a.n, b.m);
	for i := 0; i < c.n; i++ {
		for j := 0; j < c.m; j++ {
			x := Zero;
			for k := 0; k < a.m; k++ {
				x = x.Add(a.at(i, k).Mul(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 main() {
	n := 10;
	a := NewHilbert(n);
	b := NewInverseHilbert(n);
	I := NewUnit(n);
	ab := a.Mul(b);
	if !ab.Eql(I) {
		Fmt.Println("a =", a);
		Fmt.Println("b =", b);
		Fmt.Println("a*b =", ab);
		Fmt.Println("I =", I);
		panic("FAILED");
	}
}