blob: 62e649be9fa5b51c755bd37eede1637ced570643 [file] [log] [blame]
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -07001// Copyright 2009 The Go Authors. All rights reserved.
2// Use of this source code is governed by a BSD-style
3// license that can be found in the LICENSE file.
4
5// This file contains operations on unsigned multi-precision integers.
6// These are the building blocks for the operations on signed integers
7// and rationals.
8
Robert Griesemer4fe7a382009-09-01 11:56:24 -07009// NOTE: PACKAGE UNDER CONSTRUCTION (use bignum for the time being)
10//
Robert Griesemer08a209f2009-08-26 12:55:54 -070011// This package implements multi-precision arithmetic (big numbers).
12// The following numeric types are supported:
13//
14// - Int signed integers
15//
16// All methods on Int take the result as the receiver; if it is one
17// of the operands it may be overwritten (and its memory reused).
18// To enable chaining of operations, the result is also returned.
19//
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070020package big
21
22// An unsigned integer x of the form
23//
24// x = x[n-1]*_B^(n-1) + x[n-2]*_B^(n-2) + ... + x[1]*_B + x[0]
25//
26// with 0 <= x[i] < _B and 0 <= i < n is stored in a slice of length n,
27// with the digits x[i] as the slice elements.
28//
29// A number is normalized if the slice contains no leading 0 digits.
30// During arithmetic operations, denormalized values may occur but are
31// always normalized before returning the final result. The normalized
32// representation of 0 is the empty or nil slice (length = 0).
33
Robert Griesemere5874222009-08-15 11:43:54 -070034// TODO(gri) - convert these routines into methods for type 'nat'
35// - decide if type 'nat' should be exported
36
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070037func normN(z []Word) []Word {
38 i := len(z);
39 for i > 0 && z[i-1] == 0 {
40 i--;
41 }
Russ Coxc62b3262009-10-06 11:42:55 -070042 z = z[0:i];
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070043 return z;
44}
45
46
Robert Griesemer88742ef2009-08-18 10:06:15 -070047func makeN(z []Word, m int, clear bool) []Word {
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070048 if len(z) > m {
Russ Coxc62b3262009-10-06 11:42:55 -070049 z = z[0:m]; // reuse z - has at least one extra word for a carry, if any
Robert Griesemer116b52d2009-08-18 11:48:47 -070050 if clear {
51 for i := range z {
52 z[i] = 0;
53 }
Robert Griesemer88742ef2009-08-18 10:06:15 -070054 }
Robert Griesemer116b52d2009-08-18 11:48:47 -070055 return z;
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070056 }
Robert Griesemer116b52d2009-08-18 11:48:47 -070057
Russ Coxc62b3262009-10-06 11:42:55 -070058 c := 4; // minimum capacity
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070059 if m > c {
60 c = m;
61 }
Russ Coxc62b3262009-10-06 11:42:55 -070062 return make([]Word, m, c+1); // +1: extra word for a carry, if any
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070063}
64
65
66func newN(z []Word, x uint64) []Word {
67 if x == 0 {
Robert Griesemer88742ef2009-08-18 10:06:15 -070068 return makeN(z, 0, false);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070069 }
70
71 // single-digit values
72 if x == uint64(Word(x)) {
Robert Griesemer88742ef2009-08-18 10:06:15 -070073 z = makeN(z, 1, false);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070074 z[0] = Word(x);
75 return z;
76 }
77
78 // compute number of words n required to represent x
79 n := 0;
80 for t := x; t > 0; t >>= _W {
81 n++;
82 }
83
84 // split x into n words
Robert Griesemer88742ef2009-08-18 10:06:15 -070085 z = makeN(z, n, false);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070086 for i := 0; i < n; i++ {
Russ Coxc62b3262009-10-06 11:42:55 -070087 z[i] = Word(x&_M);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070088 x >>= _W;
89 }
90
91 return z;
92}
93
94
95func setN(z, x []Word) []Word {
Robert Griesemer88742ef2009-08-18 10:06:15 -070096 z = makeN(z, len(x), false);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -070097 for i, d := range x {
98 z[i] = d;
99 }
100 return z;
101}
102
103
104func addNN(z, x, y []Word) []Word {
105 m := len(x);
106 n := len(y);
107
108 switch {
109 case m < n:
110 return addNN(z, y, x);
111 case m == 0:
112 // n == 0 because m >= n; result is 0
Robert Griesemer88742ef2009-08-18 10:06:15 -0700113 return makeN(z, 0, false);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700114 case n == 0:
115 // result is x
116 return setN(z, x);
117 }
Robert Griesemere5874222009-08-15 11:43:54 -0700118 // m > 0
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700119
Robert Griesemer88742ef2009-08-18 10:06:15 -0700120 z = makeN(z, m, false);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700121 c := addVV(&z[0], &x[0], &y[0], n);
122 if m > n {
123 c = addVW(&z[n], &x[n], c, m-n);
124 }
125 if c > 0 {
126 z = z[0 : m+1];
127 z[m] = c;
128 }
129
130 return z;
131}
132
133
134func subNN(z, x, y []Word) []Word {
135 m := len(x);
136 n := len(y);
137
138 switch {
139 case m < n:
140 panic("underflow");
141 case m == 0:
142 // n == 0 because m >= n; result is 0
Robert Griesemer88742ef2009-08-18 10:06:15 -0700143 return makeN(z, 0, false);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700144 case n == 0:
145 // result is x
146 return setN(z, x);
147 }
Robert Griesemere5874222009-08-15 11:43:54 -0700148 // m > 0
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700149
Robert Griesemer88742ef2009-08-18 10:06:15 -0700150 z = makeN(z, m, false);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700151 c := subVV(&z[0], &x[0], &y[0], n);
152 if m > n {
153 c = subVW(&z[n], &x[n], c, m-n);
154 }
155 if c != 0 {
156 panic("underflow");
157 }
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700158 z = normN(z);
Robert Griesemere5874222009-08-15 11:43:54 -0700159
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700160 return z;
161}
162
163
Robert Griesemer88742ef2009-08-18 10:06:15 -0700164func cmpNN(x, y []Word) (r int) {
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700165 m := len(x);
166 n := len(y);
167 if m != n || m == 0 {
Robert Griesemer88742ef2009-08-18 10:06:15 -0700168 switch {
Russ Coxc62b3262009-10-06 11:42:55 -0700169 case m < n:
170 r = -1;
171 case m > n:
172 r = 1;
Robert Griesemer88742ef2009-08-18 10:06:15 -0700173 }
174 return;
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700175 }
176
177 i := m-1;
178 for i > 0 && x[i] == y[i] {
179 i--;
180 }
181
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700182 switch {
Russ Coxc62b3262009-10-06 11:42:55 -0700183 case x[i] < y[i]:
184 r = -1;
185 case x[i] > y[i]:
186 r = 1;
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700187 }
Robert Griesemer88742ef2009-08-18 10:06:15 -0700188 return;
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700189}
190
191
Robert Griesemere5874222009-08-15 11:43:54 -0700192func mulAddNWW(z, x []Word, y, r Word) []Word {
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700193 m := len(x);
Robert Griesemere5874222009-08-15 11:43:54 -0700194 if m == 0 || y == 0 {
195 return newN(z, uint64(r)); // result is r
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700196 }
197 // m > 0
Robert Griesemere5874222009-08-15 11:43:54 -0700198
Robert Griesemer88742ef2009-08-18 10:06:15 -0700199 z = makeN(z, m, false);
Robert Griesemere5874222009-08-15 11:43:54 -0700200 c := mulAddVWW(&z[0], &x[0], y, r, m);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700201 if c > 0 {
202 z = z[0 : m+1];
203 z[m] = c;
204 }
Robert Griesemere5874222009-08-15 11:43:54 -0700205
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700206 return z;
207}
208
209
210func mulNN(z, x, y []Word) []Word {
Robert Griesemere5874222009-08-15 11:43:54 -0700211 m := len(x);
212 n := len(y);
213
214 switch {
215 case m < n:
Robert Griesemer88742ef2009-08-18 10:06:15 -0700216 return mulNN(z, y, x);
Robert Griesemere5874222009-08-15 11:43:54 -0700217 case m == 0 || n == 0:
Robert Griesemer88742ef2009-08-18 10:06:15 -0700218 return makeN(z, 0, false);
219 case n == 1:
220 return mulAddNWW(z, x, y[0], 0);
Robert Griesemere5874222009-08-15 11:43:54 -0700221 }
Robert Griesemer88742ef2009-08-18 10:06:15 -0700222 // m >= n && m > 1 && n > 1
Robert Griesemere5874222009-08-15 11:43:54 -0700223
Robert Griesemer88742ef2009-08-18 10:06:15 -0700224 z = makeN(z, m+n, true);
225 if &z[0] == &x[0] || &z[0] == &y[0] {
Russ Coxc62b3262009-10-06 11:42:55 -0700226 z = makeN(nil, m+n, true); // z is an alias for x or y - cannot reuse
Robert Griesemer88742ef2009-08-18 10:06:15 -0700227 }
228 for i := 0; i < n; i++ {
229 if f := y[i]; f != 0 {
230 z[m+i] = addMulVVW(&z[i], &x[0], f, m);
231 }
232 }
233 z = normN(z);
Robert Griesemere5874222009-08-15 11:43:54 -0700234
Russ Coxc62b3262009-10-06 11:42:55 -0700235 return z;
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700236}
237
238
239// q = (x-r)/y, with 0 <= r < y
240func divNW(z, x []Word, y Word) (q []Word, r Word) {
241 m := len(x);
242 switch {
243 case y == 0:
244 panic("division by zero");
245 case y == 1:
Russ Coxc62b3262009-10-06 11:42:55 -0700246 q = setN(z, x); // result is x
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700247 return;
248 case m == 0:
Russ Coxc62b3262009-10-06 11:42:55 -0700249 q = setN(z, nil); // result is 0
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700250 return;
251 }
252 // m > 0
Robert Griesemer88742ef2009-08-18 10:06:15 -0700253 z = makeN(z, m, false);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700254 r = divWVW(&z[0], 0, &x[0], y, m);
255 q = normN(z);
256 return;
257}
258
259
Robert Griesemer08a209f2009-08-26 12:55:54 -0700260// log2 computes the integer binary logarithm of x.
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700261// The result is the integer n for which 2^n <= x < 2^(n+1).
Robert Griesemer08a209f2009-08-26 12:55:54 -0700262// If x == 0, the result is -1.
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700263func log2(x Word) int {
264 n := 0;
265 for ; x > 0; x >>= 1 {
266 n++;
267 }
268 return n-1;
269}
270
271
Robert Griesemer08a209f2009-08-26 12:55:54 -0700272// log2N computes the integer binary logarithm of x.
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700273// The result is the integer n for which 2^n <= x < 2^(n+1).
Robert Griesemer08a209f2009-08-26 12:55:54 -0700274// If x == 0, the result is -1.
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700275func log2N(x []Word) int {
276 m := len(x);
277 if m > 0 {
278 return (m-1)*int(_W) + log2(x[m-1]);
279 }
280 return -1;
281}
282
283
284func hexValue(ch byte) int {
285 var d byte;
286 switch {
Russ Coxc62b3262009-10-06 11:42:55 -0700287 case '0' <= ch && ch <= '9':
288 d = ch-'0';
289 case 'a' <= ch && ch <= 'f':
290 d = ch-'a'+10;
291 case 'A' <= ch && ch <= 'F':
292 d = ch-'A'+10;
293 default:
294 return -1;
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700295 }
296 return int(d);
297}
298
299
300// scanN returns the natural number corresponding to the
301// longest possible prefix of s representing a natural number in a
302// given conversion base, the actual conversion base used, and the
303// prefix length. The syntax of natural numbers follows the syntax
304// of unsigned integer literals in Go.
305//
306// If the base argument is 0, the string prefix determines the actual
307// conversion base. A prefix of ``0x'' or ``0X'' selects base 16; the
308// ``0'' prefix selects base 8. Otherwise the selected base is 10.
309//
310func scanN(z []Word, s string, base int) ([]Word, int, int) {
311 // determine base if necessary
312 i, n := 0, len(s);
313 if base == 0 {
314 base = 10;
315 if n > 0 && s[0] == '0' {
316 if n > 1 && (s[1] == 'x' || s[1] == 'X') {
317 base, i = 16, 2;
318 } else {
319 base, i = 8, 1;
320 }
321 }
322 }
323 if base < 2 || 16 < base {
324 panic("illegal base");
325 }
326
327 // convert string
Robert Griesemer88742ef2009-08-18 10:06:15 -0700328 z = makeN(z, len(z), false);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700329 for ; i < n; i++ {
330 d := hexValue(s[i]);
331 if 0 <= d && d < base {
Robert Griesemere5874222009-08-15 11:43:54 -0700332 z = mulAddNWW(z, z, Word(base), Word(d));
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700333 } else {
334 break;
335 }
336 }
337
338 return z, base, i;
339}
340
341
342// string converts x to a string for a given base, with 2 <= base <= 16.
343// TODO(gri) in the style of the other routines, perhaps this should take
344// a []byte buffer and return it
345func stringN(x []Word, base int) string {
346 if base < 2 || 16 < base {
347 panic("illegal base");
348 }
349
350 if len(x) == 0 {
351 return "0";
352 }
353
354 // allocate buffer for conversion
Russ Coxc62b3262009-10-06 11:42:55 -0700355 i := (log2N(x)+1)/log2(Word(base)) + 1; // +1: round up
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700356 s := make([]byte, i);
357
358 // don't destroy x
359 q := setN(nil, x);
360
361 // convert
362 for len(q) > 0 {
363 i--;
364 var r Word;
365 q, r = divNW(q, q, 10);
366 s[i] = "0123456789abcdef"[r];
Russ Coxc62b3262009-10-06 11:42:55 -0700367 }
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700368
Russ Coxc62b3262009-10-06 11:42:55 -0700369 return string(s[i:len(s)]);
Robert Griesemerdb3bf9c2009-08-14 11:53:27 -0700370}