1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
|
// 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
}
|
