Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 1 | // 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 Griesemer | 4fe7a38 | 2009-09-01 11:56:24 -0700 | [diff] [blame] | 9 | // NOTE: PACKAGE UNDER CONSTRUCTION (use bignum for the time being) |
| 10 | // |
Robert Griesemer | 08a209f | 2009-08-26 12:55:54 -0700 | [diff] [blame] | 11 | // 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 Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 20 | package 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 Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 34 | // TODO(gri) - convert these routines into methods for type 'nat' |
| 35 | // - decide if type 'nat' should be exported |
| 36 | |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 37 | func normN(z []Word) []Word { |
| 38 | i := len(z); |
| 39 | for i > 0 && z[i-1] == 0 { |
| 40 | i--; |
| 41 | } |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 42 | z = z[0:i]; |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 43 | return z; |
| 44 | } |
| 45 | |
| 46 | |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 47 | func makeN(z []Word, m int, clear bool) []Word { |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 48 | if len(z) > m { |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 49 | z = z[0:m]; // reuse z - has at least one extra word for a carry, if any |
Robert Griesemer | 116b52d | 2009-08-18 11:48:47 -0700 | [diff] [blame] | 50 | if clear { |
| 51 | for i := range z { |
| 52 | z[i] = 0; |
| 53 | } |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 54 | } |
Robert Griesemer | 116b52d | 2009-08-18 11:48:47 -0700 | [diff] [blame] | 55 | return z; |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 56 | } |
Robert Griesemer | 116b52d | 2009-08-18 11:48:47 -0700 | [diff] [blame] | 57 | |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 58 | c := 4; // minimum capacity |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 59 | if m > c { |
| 60 | c = m; |
| 61 | } |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 62 | return make([]Word, m, c+1); // +1: extra word for a carry, if any |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 63 | } |
| 64 | |
| 65 | |
| 66 | func newN(z []Word, x uint64) []Word { |
| 67 | if x == 0 { |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 68 | return makeN(z, 0, false); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 69 | } |
| 70 | |
| 71 | // single-digit values |
| 72 | if x == uint64(Word(x)) { |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 73 | z = makeN(z, 1, false); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 74 | 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 Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 85 | z = makeN(z, n, false); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 86 | for i := 0; i < n; i++ { |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 87 | z[i] = Word(x&_M); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 88 | x >>= _W; |
| 89 | } |
| 90 | |
| 91 | return z; |
| 92 | } |
| 93 | |
| 94 | |
| 95 | func setN(z, x []Word) []Word { |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 96 | z = makeN(z, len(x), false); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 97 | for i, d := range x { |
| 98 | z[i] = d; |
| 99 | } |
| 100 | return z; |
| 101 | } |
| 102 | |
| 103 | |
| 104 | func 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 Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 113 | return makeN(z, 0, false); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 114 | case n == 0: |
| 115 | // result is x |
| 116 | return setN(z, x); |
| 117 | } |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 118 | // m > 0 |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 119 | |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 120 | z = makeN(z, m, false); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 121 | 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 | |
| 134 | func 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 Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 143 | return makeN(z, 0, false); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 144 | case n == 0: |
| 145 | // result is x |
| 146 | return setN(z, x); |
| 147 | } |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 148 | // m > 0 |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 149 | |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 150 | z = makeN(z, m, false); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 151 | 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 Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 158 | z = normN(z); |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 159 | |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 160 | return z; |
| 161 | } |
| 162 | |
| 163 | |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 164 | func cmpNN(x, y []Word) (r int) { |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 165 | m := len(x); |
| 166 | n := len(y); |
| 167 | if m != n || m == 0 { |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 168 | switch { |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 169 | case m < n: |
| 170 | r = -1; |
| 171 | case m > n: |
| 172 | r = 1; |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 173 | } |
| 174 | return; |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 175 | } |
| 176 | |
| 177 | i := m-1; |
| 178 | for i > 0 && x[i] == y[i] { |
| 179 | i--; |
| 180 | } |
| 181 | |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 182 | switch { |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 183 | case x[i] < y[i]: |
| 184 | r = -1; |
| 185 | case x[i] > y[i]: |
| 186 | r = 1; |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 187 | } |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 188 | return; |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 189 | } |
| 190 | |
| 191 | |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 192 | func mulAddNWW(z, x []Word, y, r Word) []Word { |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 193 | m := len(x); |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 194 | if m == 0 || y == 0 { |
| 195 | return newN(z, uint64(r)); // result is r |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 196 | } |
| 197 | // m > 0 |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 198 | |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 199 | z = makeN(z, m, false); |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 200 | c := mulAddVWW(&z[0], &x[0], y, r, m); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 201 | if c > 0 { |
| 202 | z = z[0 : m+1]; |
| 203 | z[m] = c; |
| 204 | } |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 205 | |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 206 | return z; |
| 207 | } |
| 208 | |
| 209 | |
| 210 | func mulNN(z, x, y []Word) []Word { |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 211 | m := len(x); |
| 212 | n := len(y); |
| 213 | |
| 214 | switch { |
| 215 | case m < n: |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 216 | return mulNN(z, y, x); |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 217 | case m == 0 || n == 0: |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 218 | return makeN(z, 0, false); |
| 219 | case n == 1: |
| 220 | return mulAddNWW(z, x, y[0], 0); |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 221 | } |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 222 | // m >= n && m > 1 && n > 1 |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 223 | |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 224 | z = makeN(z, m+n, true); |
| 225 | if &z[0] == &x[0] || &z[0] == &y[0] { |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 226 | z = makeN(nil, m+n, true); // z is an alias for x or y - cannot reuse |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 227 | } |
| 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 Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 234 | |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 235 | return z; |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 236 | } |
| 237 | |
| 238 | |
| 239 | // q = (x-r)/y, with 0 <= r < y |
| 240 | func 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 Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 246 | q = setN(z, x); // result is x |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 247 | return; |
| 248 | case m == 0: |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 249 | q = setN(z, nil); // result is 0 |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 250 | return; |
| 251 | } |
| 252 | // m > 0 |
Robert Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 253 | z = makeN(z, m, false); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 254 | r = divWVW(&z[0], 0, &x[0], y, m); |
| 255 | q = normN(z); |
| 256 | return; |
| 257 | } |
| 258 | |
| 259 | |
Robert Griesemer | 08a209f | 2009-08-26 12:55:54 -0700 | [diff] [blame] | 260 | // log2 computes the integer binary logarithm of x. |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 261 | // The result is the integer n for which 2^n <= x < 2^(n+1). |
Robert Griesemer | 08a209f | 2009-08-26 12:55:54 -0700 | [diff] [blame] | 262 | // If x == 0, the result is -1. |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 263 | func log2(x Word) int { |
| 264 | n := 0; |
| 265 | for ; x > 0; x >>= 1 { |
| 266 | n++; |
| 267 | } |
| 268 | return n-1; |
| 269 | } |
| 270 | |
| 271 | |
Robert Griesemer | 08a209f | 2009-08-26 12:55:54 -0700 | [diff] [blame] | 272 | // log2N computes the integer binary logarithm of x. |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 273 | // The result is the integer n for which 2^n <= x < 2^(n+1). |
Robert Griesemer | 08a209f | 2009-08-26 12:55:54 -0700 | [diff] [blame] | 274 | // If x == 0, the result is -1. |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 275 | func 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 | |
| 284 | func hexValue(ch byte) int { |
| 285 | var d byte; |
| 286 | switch { |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 287 | 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 Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 295 | } |
| 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 | // |
| 310 | func 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 Griesemer | 88742ef | 2009-08-18 10:06:15 -0700 | [diff] [blame] | 328 | z = makeN(z, len(z), false); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 329 | for ; i < n; i++ { |
| 330 | d := hexValue(s[i]); |
| 331 | if 0 <= d && d < base { |
Robert Griesemer | e587422 | 2009-08-15 11:43:54 -0700 | [diff] [blame] | 332 | z = mulAddNWW(z, z, Word(base), Word(d)); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 333 | } 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 |
| 345 | func 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 Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 355 | i := (log2N(x)+1)/log2(Word(base)) + 1; // +1: round up |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 356 | 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 Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 367 | } |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 368 | |
Russ Cox | c62b326 | 2009-10-06 11:42:55 -0700 | [diff] [blame^] | 369 | return string(s[i:len(s)]); |
Robert Griesemer | db3bf9c | 2009-08-14 11:53:27 -0700 | [diff] [blame] | 370 | } |