Skip to content

G
golang
Commit 38b2d10b authored by Robert Griesemer's avatar Robert Griesemer
Browse files
Options

test/hilbert.go: convert to test case and benchmark for big.Rat 

R=rsc
CC=golang-dev
https://golang.org/cl/1231044
parent 88b308a2
Hide whitespace changes
Inline Side-by-side
Showing with 173 additions and 169 deletions
src/pkg/big/hilbert_test.go 0 → 100644
View file @ 38b2d10b
// 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)
}
}
src/pkg/big/nat.go
View file @ 38b2d10b
......@@ -16,9 +16,6 @@
// of the operands it may be overwritten (and its memory reused).
// To enable chaining of operations, the result is also returned.
//
// If possible, one should use big over bignum as the latter is headed for
// deprecation.
//
package big
import "rand"
......
test/hilbert.go deleted 100644 → 0
View file @ 88b308a2
// $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");
}
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment