--- /dev/null
+// 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)
+ }
+}
+++ /dev/null
-// $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");
- }
-}