From 06c2c8945223b70111544103191c91ac0ecbf9e1 Mon Sep 17 00:00:00 2001
From: Robert Griesemer <gri@golang.org>
Date: Wed, 26 Aug 2009 09:46:12 -0700
Subject: [PATCH] added Newton-Raphson Division as an additional bignum
 testcase

R=rsc
DELTA=192  (192 added, 0 deleted, 0 changed)
OCL=33853
CL=33864
---
 src/pkg/bignum/nrdiv_test.go | 192 +++++++++++++++++++++++++++++++++++
 1 file changed, 192 insertions(+)
 create mode 100644 src/pkg/bignum/nrdiv_test.go

diff --git a/src/pkg/bignum/nrdiv_test.go b/src/pkg/bignum/nrdiv_test.go
new file mode 100644
index 00000000000..24cb4d5587f
--- /dev/null
+++ b/src/pkg/bignum/nrdiv_test.go
@@ -0,0 +1,192 @@
+// 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.
+
+// This file implements Newton-Raphson division and uses
+// it as an additional test case for bignum.
+//
+// Division of x/y is achieved by computing r = 1/y to
+// obtain the quotient q = x*r = x*(1/y) = x/y. The
+// reciprocal r is the solution for f(x) = 1/x - y and
+// the solution is approximated through iteration. The
+// iteration does not require division.
+
+package bignum
+
+import "testing"
+
+
+// An fpNat is a Natural scaled by a power of two
+// (an unsigned floating point representation). The
+// value of an fpNat x is x.m * 2^x.e .
+//
+type fpNat struct {
+	m Natural;
+	e int;
+}
+
+
+// sub computes x - y.
+func (x fpNat) sub(y fpNat) fpNat {
+	switch d := x.e - y.e; {
+	case d < 0:
+		return fpNat{x.m.Sub(y.m.Shl(uint(-d))), x.e};
+	case d > 0:
+		return fpNat{x.m.Shl(uint(d)).Sub(y.m), y.e};
+	}
+	return fpNat{x.m.Sub(y.m), x.e};
+}
+
+
+// mul2 computes x*2.
+func (x fpNat) mul2() fpNat {
+	return fpNat{x.m, x.e+1};
+}
+
+
+// mul computes x*y.
+func (x fpNat) mul(y fpNat) fpNat {
+	return fpNat{x.m.Mul(y.m), x.e + y.e};
+}
+
+
+// mant computes the (possibly truncated) Natural representation
+// of an fpNat x.
+//
+func (x fpNat) mant() Natural {
+	switch {
+	case x.e > 0: return x.m.Shl(uint(x.e));
+	case x.e < 0: return x.m.Shr(uint(-x.e));
+	}
+	return x.m;
+}
+
+
+// nrDivEst computes an estimate of the quotient q = x0/y0 and returns q.
+// q may be too small (usually by 1).
+//
+func nrDivEst(x0, y0 Natural) Natural {
+	if y0.IsZero() {
+		panic("division by zero");
+		return nil;
+	}
+	// y0 > 0
+
+	if y0.Cmp(Nat(1)) == 0 {
+		return x0;
+	}
+	// y0 > 1
+
+	switch d := x0.Cmp(y0); {
+	case d < 0:
+		return Nat(0);
+	case d == 0:
+		return Nat(1);
+	}
+	// x0 > y0 > 1
+
+	// Determine maximum result length.
+	maxLen := int(x0.Log2() - y0.Log2() + 1);
+
+	// In the following, each number x is represented
+	// as a mantissa x.m and an exponent x.e such that
+	// x = xm * 2^x.e.
+	x := fpNat{x0, 0};
+	y := fpNat{y0, 0};
+
+	// Determine a scale factor f = 2^e such that
+	// 0.5 <= y/f == y*(2^-e) < 1.0
+	// and scale y accordingly.
+	e := int(y.m.Log2()) + 1;
+	y.e -= e;
+
+	// t1
+	var c = 2.9142;
+	const n = 14;
+	t1 := fpNat{Nat(uint64(c*(1<<n))), -n};
+
+	// Compute initial value r0 for the reciprocal of y/f.
+	// r0 = t1 - 2*y
+	r := t1.sub(y.mul2());
+	two := fpNat{Nat(2), 0};
+
+	// Newton-Raphson iteration
+	p := Nat(0);
+	for i := 0; ; i++ {
+		// check if we are done
+		// TODO: Need to come up with a better test here
+		//       as it will reduce computation time significantly.
+		// q = x*r/f
+		q := x.mul(r);
+		q.e -= e;
+		res := q.mant();
+		if res.Cmp(p) == 0 {
+			return res;
+		}
+		p = res;
+
+		// r' = r*(2 - y*r)
+		r = r.mul(two.sub(y.mul(r)));
+
+		// reduce mantissa size
+		// TODO: Find smaller bound as it will reduce
+		//       computation time massively.
+		d := int(r.m.Log2()+1) - maxLen;
+		if d > 0 {
+			r = fpNat{r.m.Shr(uint(d)), r.e+d};
+		}
+	}
+
+	panic("unreachable");
+	return nil;
+}
+
+
+func nrdiv(x, y Natural) (q, r Natural) {
+	q = nrDivEst(x, y);
+	r = x.Sub(y.Mul(q));
+	// if r is too large, correct q and r
+	// (usually one iteration)
+	for r.Cmp(y) >= 0 {
+		q = q.Add(Nat(1));
+		r = r.Sub(y);
+	}
+	return;
+}
+
+
+func div(t *testing.T, x, y Natural) {
+	q, r := nrdiv(x, y);
+	qx, rx := x.DivMod(y);
+	if q.Cmp(qx) != 0 {
+		t.Errorf("x = %s, y = %s, got q = %s, want q = %s", x, y, q, qx);
+	}
+	if r.Cmp(rx) != 0 {
+		t.Errorf("x = %s, y = %s, got r = %s, want r = %s", x, y, r, rx);
+	}
+}
+
+
+func idiv(t *testing.T, x0, y0 uint64) {
+	div(t, Nat(x0), Nat(y0));
+}
+
+
+func TestNRDiv(t *testing.T) {
+	idiv(t, 17, 18);
+	idiv(t, 17, 17);
+	idiv(t, 17, 1);
+	idiv(t, 17, 16);
+	idiv(t, 17, 10);
+	idiv(t, 17, 9);
+	idiv(t, 17, 8);
+	idiv(t, 17, 5);
+	idiv(t, 17, 3);
+	idiv(t, 1025, 512);
+	idiv(t, 7489595, 2);
+	idiv(t, 5404679459, 78495);
+	idiv(t, 7484890589595, 7484890589594);
+	div(t, Fact(100), Fact(91));
+	div(t, Fact(1000), Fact(991));
+	//div(t, Fact(10000), Fact(9991));  // takes too long - disabled for now
+}