blob: 24cb4d5587f6b3e42eba780b4341aedba1bf5154 [file] [log] [blame]
 // 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< 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 }