2009-08-26 10:46:12 -06:00
|
|
|
// 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 {
|
2009-10-06 15:55:39 -06:00
|
|
|
m Natural;
|
|
|
|
e int;
|
2009-08-26 10:46:12 -06:00
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
// sub computes x - y.
|
|
|
|
func (x fpNat) sub(y fpNat) fpNat {
|
|
|
|
switch d := x.e - y.e; {
|
|
|
|
case d < 0:
|
2009-11-09 13:07:39 -07:00
|
|
|
return fpNat{x.m.Sub(y.m.Shl(uint(-d))), x.e}
|
2009-08-26 10:46:12 -06:00
|
|
|
case d > 0:
|
2009-11-09 13:07:39 -07:00
|
|
|
return fpNat{x.m.Shl(uint(d)).Sub(y.m), y.e}
|
2009-08-26 10:46:12 -06:00
|
|
|
}
|
|
|
|
return fpNat{x.m.Sub(y.m), x.e};
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
// mul2 computes x*2.
|
2009-11-06 15:24:38 -07:00
|
|
|
func (x fpNat) mul2() fpNat { return fpNat{x.m, x.e + 1} }
|
2009-08-26 10:46:12 -06:00
|
|
|
|
|
|
|
|
|
|
|
// mul computes x*y.
|
2009-11-06 15:24:38 -07:00
|
|
|
func (x fpNat) mul(y fpNat) fpNat { return fpNat{x.m.Mul(y.m), x.e + y.e} }
|
2009-08-26 10:46:12 -06:00
|
|
|
|
|
|
|
|
|
|
|
// mant computes the (possibly truncated) Natural representation
|
|
|
|
// of an fpNat x.
|
|
|
|
//
|
|
|
|
func (x fpNat) mant() Natural {
|
|
|
|
switch {
|
2009-10-06 15:55:39 -06:00
|
|
|
case x.e > 0:
|
2009-11-09 13:07:39 -07:00
|
|
|
return x.m.Shl(uint(x.e))
|
2009-10-06 15:55:39 -06:00
|
|
|
case x.e < 0:
|
2009-11-09 13:07:39 -07:00
|
|
|
return x.m.Shr(uint(-x.e))
|
2009-08-26 10:46:12 -06:00
|
|
|
}
|
|
|
|
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 {
|
2009-11-09 13:07:39 -07:00
|
|
|
return x0
|
2009-08-26 10:46:12 -06:00
|
|
|
}
|
|
|
|
// y0 > 1
|
|
|
|
|
|
|
|
switch d := x0.Cmp(y0); {
|
|
|
|
case d < 0:
|
2009-11-09 13:07:39 -07:00
|
|
|
return Nat(0)
|
2009-08-26 10:46:12 -06:00
|
|
|
case d == 0:
|
2009-11-09 13:07:39 -07:00
|
|
|
return Nat(1)
|
2009-08-26 10:46:12 -06:00
|
|
|
}
|
|
|
|
// 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.
|
2009-11-09 22:09:34 -07:00
|
|
|
e := int(y.m.Log2()) + 1;
|
2009-08-26 10:46:12 -06:00
|
|
|
y.e -= e;
|
|
|
|
|
|
|
|
// t1
|
|
|
|
var c = 2.9142;
|
|
|
|
const n = 14;
|
2009-11-09 22:09:34 -07:00
|
|
|
t1 := fpNat{Nat(uint64(c * (1 << n))), -n};
|
2009-08-26 10:46:12 -06:00
|
|
|
|
|
|
|
// 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 {
|
2009-11-09 13:07:39 -07:00
|
|
|
return res
|
2009-08-26 10:46:12 -06:00
|
|
|
}
|
|
|
|
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.
|
2009-11-09 22:09:34 -07:00
|
|
|
d := int(r.m.Log2()+1) - maxLen;
|
2009-08-26 10:46:12 -06:00
|
|
|
if d > 0 {
|
2009-11-09 13:07:39 -07:00
|
|
|
r = fpNat{r.m.Shr(uint(d)), r.e + d}
|
2009-08-26 10:46:12 -06:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
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 {
|
2009-11-09 13:07:39 -07:00
|
|
|
t.Errorf("x = %s, y = %s, got q = %s, want q = %s", x, y, q, qx)
|
2009-08-26 10:46:12 -06:00
|
|
|
}
|
|
|
|
if r.Cmp(rx) != 0 {
|
2009-11-09 13:07:39 -07:00
|
|
|
t.Errorf("x = %s, y = %s, got r = %s, want r = %s", x, y, r, rx)
|
2009-08-26 10:46:12 -06:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
|
2009-11-06 15:24:38 -07:00
|
|
|
func idiv(t *testing.T, x0, y0 uint64) { div(t, Nat(x0), Nat(y0)) }
|
2009-08-26 10:46:12 -06:00
|
|
|
|
|
|
|
|
|
|
|
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));
|
2009-10-27 23:47:54 -06:00
|
|
|
//div(t, Fact(10000), Fact(9991)); // takes too long - disabled for now
|
2009-08-26 10:46:12 -06:00
|
|
|
}
|