OLD | NEW |
1 // Copyright 2009 The Go Authors. All rights reserved. | 1 // Copyright 2009 The Go Authors. All rights reserved. |
2 // Use of this source code is governed by a BSD-style | 2 // Use of this source code is governed by a BSD-style |
3 // license that can be found in the LICENSE file. | 3 // license that can be found in the LICENSE file. |
4 | 4 |
5 // This file contains operations on unsigned multi-precision integers. | 5 // This file contains operations on unsigned multi-precision integers. |
6 // These are the building blocks for the operations on signed integers | 6 // These are the building blocks for the operations on signed integers |
7 // and rationals. | 7 // and rationals. |
8 | 8 |
9 // This package implements multi-precision arithmetic (big numbers). | 9 // This package implements multi-precision arithmetic (big numbers). |
10 // The following numeric types are supported: | 10 // The following numeric types are supported: |
(...skipping 20 matching lines...) Expand all Loading... |
31 // | 31 // |
32 // A number is normalized if the slice contains no leading 0 digits. | 32 // A number is normalized if the slice contains no leading 0 digits. |
33 // During arithmetic operations, denormalized values may occur but are | 33 // During arithmetic operations, denormalized values may occur but are |
34 // always normalized before returning the final result. The normalized | 34 // always normalized before returning the final result. The normalized |
35 // representation of 0 is the empty or nil slice (length = 0). | 35 // representation of 0 is the empty or nil slice (length = 0). |
36 | 36 |
37 // TODO(gri) - convert these routines into methods for type 'nat' | 37 // TODO(gri) - convert these routines into methods for type 'nat' |
38 // - decide if type 'nat' should be exported | 38 // - decide if type 'nat' should be exported |
39 | 39 |
40 func normN(z []Word) []Word { | 40 func normN(z []Word) []Word { |
41 » i := len(z); | 41 » i := len(z) |
42 for i > 0 && z[i-1] == 0 { | 42 for i > 0 && z[i-1] == 0 { |
43 i-- | 43 i-- |
44 } | 44 } |
45 » z = z[0:i]; | 45 » z = z[0:i] |
46 » return z; | 46 » return z |
47 } | 47 } |
48 | 48 |
49 | 49 |
50 func makeN(z []Word, m int, clear bool) []Word { | 50 func makeN(z []Word, m int, clear bool) []Word { |
51 if len(z) > m { | 51 if len(z) > m { |
52 » » z = z[0:m];» // reuse z - has at least one extra word for a c
arry, if any | 52 » » z = z[0:m] // reuse z - has at least one extra word for a carry,
if any |
53 if clear { | 53 if clear { |
54 for i := range z { | 54 for i := range z { |
55 z[i] = 0 | 55 z[i] = 0 |
56 } | 56 } |
57 } | 57 } |
58 » » return z; | 58 » » return z |
59 } | 59 } |
60 | 60 |
61 » c := 4;»// minimum capacity | 61 » c := 4 // minimum capacity |
62 if m > c { | 62 if m > c { |
63 c = m | 63 c = m |
64 } | 64 } |
65 » return make([]Word, m, c+1);» // +1: extra word for a carry, if any | 65 » return make([]Word, m, c+1) // +1: extra word for a carry, if any |
66 } | 66 } |
67 | 67 |
68 | 68 |
69 func newN(z []Word, x uint64) []Word { | 69 func newN(z []Word, x uint64) []Word { |
70 if x == 0 { | 70 if x == 0 { |
71 return makeN(z, 0, false) | 71 return makeN(z, 0, false) |
72 } | 72 } |
73 | 73 |
74 // single-digit values | 74 // single-digit values |
75 if x == uint64(Word(x)) { | 75 if x == uint64(Word(x)) { |
76 » » z = makeN(z, 1, false); | 76 » » z = makeN(z, 1, false) |
77 » » z[0] = Word(x); | 77 » » z[0] = Word(x) |
78 » » return z; | 78 » » return z |
79 } | 79 } |
80 | 80 |
81 // compute number of words n required to represent x | 81 // compute number of words n required to represent x |
82 » n := 0; | 82 » n := 0 |
83 for t := x; t > 0; t >>= _W { | 83 for t := x; t > 0; t >>= _W { |
84 n++ | 84 n++ |
85 } | 85 } |
86 | 86 |
87 // split x into n words | 87 // split x into n words |
88 » z = makeN(z, n, false); | 88 » z = makeN(z, n, false) |
89 for i := 0; i < n; i++ { | 89 for i := 0; i < n; i++ { |
90 » » z[i] = Word(x & _M); | 90 » » z[i] = Word(x & _M) |
91 » » x >>= _W; | 91 » » x >>= _W |
92 } | 92 } |
93 | 93 |
94 » return z; | 94 » return z |
95 } | 95 } |
96 | 96 |
97 | 97 |
98 func setN(z, x []Word) []Word { | 98 func setN(z, x []Word) []Word { |
99 » z = makeN(z, len(x), false); | 99 » z = makeN(z, len(x), false) |
100 for i, d := range x { | 100 for i, d := range x { |
101 z[i] = d | 101 z[i] = d |
102 } | 102 } |
103 » return z; | 103 » return z |
104 } | 104 } |
105 | 105 |
106 | 106 |
107 func addNN(z, x, y []Word) []Word { | 107 func addNN(z, x, y []Word) []Word { |
108 » m := len(x); | 108 » m := len(x) |
109 » n := len(y); | 109 » n := len(y) |
110 | 110 |
111 switch { | 111 switch { |
112 case m < n: | 112 case m < n: |
113 return addNN(z, y, x) | 113 return addNN(z, y, x) |
114 case m == 0: | 114 case m == 0: |
115 // n == 0 because m >= n; result is 0 | 115 // n == 0 because m >= n; result is 0 |
116 return makeN(z, 0, false) | 116 return makeN(z, 0, false) |
117 case n == 0: | 117 case n == 0: |
118 // result is x | 118 // result is x |
119 return setN(z, x) | 119 return setN(z, x) |
120 } | 120 } |
121 // m > 0 | 121 // m > 0 |
122 | 122 |
123 » z = makeN(z, m, false); | 123 » z = makeN(z, m, false) |
124 » c := addVV(&z[0], &x[0], &y[0], n); | 124 » c := addVV(&z[0], &x[0], &y[0], n) |
125 if m > n { | 125 if m > n { |
126 c = addVW(&z[n], &x[n], c, m-n) | 126 c = addVW(&z[n], &x[n], c, m-n) |
127 } | 127 } |
128 if c > 0 { | 128 if c > 0 { |
129 » » z = z[0 : m+1]; | 129 » » z = z[0 : m+1] |
130 » » z[m] = c; | 130 » » z[m] = c |
131 } | 131 } |
132 | 132 |
133 » return z; | 133 » return z |
134 } | 134 } |
135 | 135 |
136 | 136 |
137 func subNN(z, x, y []Word) []Word { | 137 func subNN(z, x, y []Word) []Word { |
138 » m := len(x); | 138 » m := len(x) |
139 » n := len(y); | 139 » n := len(y) |
140 | 140 |
141 switch { | 141 switch { |
142 case m < n: | 142 case m < n: |
143 panic("underflow") | 143 panic("underflow") |
144 case m == 0: | 144 case m == 0: |
145 // n == 0 because m >= n; result is 0 | 145 // n == 0 because m >= n; result is 0 |
146 return makeN(z, 0, false) | 146 return makeN(z, 0, false) |
147 case n == 0: | 147 case n == 0: |
148 // result is x | 148 // result is x |
149 return setN(z, x) | 149 return setN(z, x) |
150 } | 150 } |
151 // m > 0 | 151 // m > 0 |
152 | 152 |
153 » z = makeN(z, m, false); | 153 » z = makeN(z, m, false) |
154 » c := subVV(&z[0], &x[0], &y[0], n); | 154 » c := subVV(&z[0], &x[0], &y[0], n) |
155 if m > n { | 155 if m > n { |
156 c = subVW(&z[n], &x[n], c, m-n) | 156 c = subVW(&z[n], &x[n], c, m-n) |
157 } | 157 } |
158 if c != 0 { | 158 if c != 0 { |
159 panic("underflow") | 159 panic("underflow") |
160 } | 160 } |
161 » z = normN(z); | 161 » z = normN(z) |
162 | 162 |
163 » return z; | 163 » return z |
164 } | 164 } |
165 | 165 |
166 | 166 |
167 func cmpNN(x, y []Word) (r int) { | 167 func cmpNN(x, y []Word) (r int) { |
168 » m := len(x); | 168 » m := len(x) |
169 » n := len(y); | 169 » n := len(y) |
170 if m != n || m == 0 { | 170 if m != n || m == 0 { |
171 switch { | 171 switch { |
172 case m < n: | 172 case m < n: |
173 r = -1 | 173 r = -1 |
174 case m > n: | 174 case m > n: |
175 r = 1 | 175 r = 1 |
176 } | 176 } |
177 » » return; | 177 » » return |
178 } | 178 } |
179 | 179 |
180 » i := m - 1; | 180 » i := m - 1 |
181 for i > 0 && x[i] == y[i] { | 181 for i > 0 && x[i] == y[i] { |
182 i-- | 182 i-- |
183 } | 183 } |
184 | 184 |
185 switch { | 185 switch { |
186 case x[i] < y[i]: | 186 case x[i] < y[i]: |
187 r = -1 | 187 r = -1 |
188 case x[i] > y[i]: | 188 case x[i] > y[i]: |
189 r = 1 | 189 r = 1 |
190 } | 190 } |
191 » return; | 191 » return |
192 } | 192 } |
193 | 193 |
194 | 194 |
195 func mulAddNWW(z, x []Word, y, r Word) []Word { | 195 func mulAddNWW(z, x []Word, y, r Word) []Word { |
196 » m := len(x); | 196 » m := len(x) |
197 if m == 0 || y == 0 { | 197 if m == 0 || y == 0 { |
198 » » return newN(z, uint64(r))» // result is r | 198 » » return newN(z, uint64(r)) // result is r |
199 } | 199 } |
200 // m > 0 | 200 // m > 0 |
201 | 201 |
202 » z = makeN(z, m, false); | 202 » z = makeN(z, m, false) |
203 » c := mulAddVWW(&z[0], &x[0], y, r, m); | 203 » c := mulAddVWW(&z[0], &x[0], y, r, m) |
204 if c > 0 { | 204 if c > 0 { |
205 » » z = z[0 : m+1]; | 205 » » z = z[0 : m+1] |
206 » » z[m] = c; | 206 » » z[m] = c |
207 } | 207 } |
208 | 208 |
209 » return z; | 209 » return z |
210 } | 210 } |
211 | 211 |
212 | 212 |
213 func mulNN(z, x, y []Word) []Word { | 213 func mulNN(z, x, y []Word) []Word { |
214 » m := len(x); | 214 » m := len(x) |
215 » n := len(y); | 215 » n := len(y) |
216 | 216 |
217 switch { | 217 switch { |
218 case m < n: | 218 case m < n: |
219 return mulNN(z, y, x) | 219 return mulNN(z, y, x) |
220 case m == 0 || n == 0: | 220 case m == 0 || n == 0: |
221 return makeN(z, 0, false) | 221 return makeN(z, 0, false) |
222 case n == 1: | 222 case n == 1: |
223 return mulAddNWW(z, x, y[0], 0) | 223 return mulAddNWW(z, x, y[0], 0) |
224 } | 224 } |
225 // m >= n && m > 1 && n > 1 | 225 // m >= n && m > 1 && n > 1 |
226 | 226 |
227 » z = makeN(z, m+n, true); | 227 » z = makeN(z, m+n, true) |
228 if &z[0] == &x[0] || &z[0] == &y[0] { | 228 if &z[0] == &x[0] || &z[0] == &y[0] { |
229 » » z = makeN(nil, m+n, true)» // z is an alias for x or y - ca
nnot reuse | 229 » » z = makeN(nil, m+n, true) // z is an alias for x or y - cannot r
euse |
230 } | 230 } |
231 for i := 0; i < n; i++ { | 231 for i := 0; i < n; i++ { |
232 if f := y[i]; f != 0 { | 232 if f := y[i]; f != 0 { |
233 z[m+i] = addMulVVW(&z[i], &x[0], f, m) | 233 z[m+i] = addMulVVW(&z[i], &x[0], f, m) |
234 } | 234 } |
235 } | 235 } |
236 » z = normN(z); | 236 » z = normN(z) |
237 | 237 |
238 » return z; | 238 » return z |
239 } | 239 } |
240 | 240 |
241 | 241 |
242 // q = (x-r)/y, with 0 <= r < y | 242 // q = (x-r)/y, with 0 <= r < y |
243 func divNW(z, x []Word, y Word) (q []Word, r Word) { | 243 func divNW(z, x []Word, y Word) (q []Word, r Word) { |
244 » m := len(x); | 244 » m := len(x) |
245 switch { | 245 switch { |
246 case y == 0: | 246 case y == 0: |
247 panic("division by zero") | 247 panic("division by zero") |
248 case y == 1: | 248 case y == 1: |
249 » » q = setN(z, x);»// result is x | 249 » » q = setN(z, x) // result is x |
250 » » return; | 250 » » return |
251 case m == 0: | 251 case m == 0: |
252 » » q = setN(z, nil);» // result is 0 | 252 » » q = setN(z, nil) // result is 0 |
253 » » return; | 253 » » return |
254 } | 254 } |
255 // m > 0 | 255 // m > 0 |
256 » z = makeN(z, m, false); | 256 » z = makeN(z, m, false) |
257 » r = divWVW(&z[0], 0, &x[0], y, m); | 257 » r = divWVW(&z[0], 0, &x[0], y, m) |
258 » q = normN(z); | 258 » q = normN(z) |
259 » return; | 259 » return |
260 } | 260 } |
261 | 261 |
262 | 262 |
263 func divNN(z, z2, u, v []Word) (q, r []Word) { | 263 func divNN(z, z2, u, v []Word) (q, r []Word) { |
264 if len(v) == 0 { | 264 if len(v) == 0 { |
265 panic("Divide by zero undefined") | 265 panic("Divide by zero undefined") |
266 } | 266 } |
267 | 267 |
268 if cmpNN(u, v) < 0 { | 268 if cmpNN(u, v) < 0 { |
269 » » q = makeN(z, 0, false); | 269 » » q = makeN(z, 0, false) |
270 » » r = setN(z2, u); | 270 » » r = setN(z2, u) |
271 » » return; | 271 » » return |
272 } | 272 } |
273 | 273 |
274 if len(v) == 1 { | 274 if len(v) == 1 { |
275 » » var rprime Word; | 275 » » var rprime Word |
276 » » q, rprime = divNW(z, u, v[0]); | 276 » » q, rprime = divNW(z, u, v[0]) |
277 if rprime > 0 { | 277 if rprime > 0 { |
278 » » » r = makeN(z2, 1, false); | 278 » » » r = makeN(z2, 1, false) |
279 » » » r[0] = rprime; | 279 » » » r[0] = rprime |
280 } else { | 280 } else { |
281 r = makeN(z2, 0, false) | 281 r = makeN(z2, 0, false) |
282 } | 282 } |
283 » » return; | 283 » » return |
284 } | 284 } |
285 | 285 |
286 » q, r = divLargeNN(z, z2, u, v); | 286 » q, r = divLargeNN(z, z2, u, v) |
287 » return; | 287 » return |
288 } | 288 } |
289 | 289 |
290 | 290 |
291 // q = (uIn-r)/v, with 0 <= r < y | 291 // q = (uIn-r)/v, with 0 <= r < y |
292 // See Knuth, Volume 2, section 4.3.1, Algorithm D. | 292 // See Knuth, Volume 2, section 4.3.1, Algorithm D. |
293 // Preconditions: | 293 // Preconditions: |
294 // len(v) >= 2 | 294 // len(v) >= 2 |
295 // len(uIn) >= len(v) | 295 // len(uIn) >= len(v) |
296 func divLargeNN(z, z2, uIn, v []Word) (q, r []Word) { | 296 func divLargeNN(z, z2, uIn, v []Word) (q, r []Word) { |
297 » n := len(v); | 297 » n := len(v) |
298 » m := len(uIn) - len(v); | 298 » m := len(uIn) - len(v) |
299 | 299 |
300 » u := makeN(z2, len(uIn)+1, false); | 300 » u := makeN(z2, len(uIn)+1, false) |
301 » qhatv := make([]Word, len(v)+1); | 301 » qhatv := make([]Word, len(v)+1) |
302 » q = makeN(z, m+1, false); | 302 » q = makeN(z, m+1, false) |
303 | 303 |
304 // D1. | 304 // D1. |
305 » shift := leadingZeroBits(v[n-1]); | 305 » shift := leadingZeroBits(v[n-1]) |
306 » shiftLeft(v, v, shift); | 306 » shiftLeft(v, v, shift) |
307 » shiftLeft(u, uIn, shift); | 307 » shiftLeft(u, uIn, shift) |
308 » u[len(uIn)] = uIn[len(uIn)-1] >> (_W - uint(shift)); | 308 » u[len(uIn)] = uIn[len(uIn)-1] >> (_W - uint(shift)) |
309 | 309 |
310 // D2. | 310 // D2. |
311 for j := m; j >= 0; j-- { | 311 for j := m; j >= 0; j-- { |
312 // D3. | 312 // D3. |
313 » » var qhat Word; | 313 » » var qhat Word |
314 if u[j+n] == v[n-1] { | 314 if u[j+n] == v[n-1] { |
315 qhat = _B - 1 | 315 qhat = _B - 1 |
316 } else { | 316 } else { |
317 » » » var rhat Word; | 317 » » » var rhat Word |
318 » » » qhat, rhat = divWW_g(u[j+n], u[j+n-1], v[n-1]); | 318 » » » qhat, rhat = divWW_g(u[j+n], u[j+n-1], v[n-1]) |
319 | 319 |
320 // x1 | x2 = q̂v_{n-2} | 320 // x1 | x2 = q̂v_{n-2} |
321 » » » x1, x2 := mulWW_g(qhat, v[n-2]); | 321 » » » x1, x2 := mulWW_g(qhat, v[n-2]) |
322 // test if q̂v_{n-2} > br̂ + u_{j+n-2} | 322 // test if q̂v_{n-2} > br̂ + u_{j+n-2} |
323 for greaterThan(x1, x2, rhat, u[j+n-2]) { | 323 for greaterThan(x1, x2, rhat, u[j+n-2]) { |
324 » » » » qhat--; | 324 » » » » qhat-- |
325 » » » » prevRhat := rhat; | 325 » » » » prevRhat := rhat |
326 » » » » rhat += v[n-1]; | 326 » » » » rhat += v[n-1] |
327 // v[n-1] >= 0, so this tests for overflow. | 327 // v[n-1] >= 0, so this tests for overflow. |
328 if rhat < prevRhat { | 328 if rhat < prevRhat { |
329 break | 329 break |
330 } | 330 } |
331 » » » » x1, x2 = mulWW_g(qhat, v[n-2]); | 331 » » » » x1, x2 = mulWW_g(qhat, v[n-2]) |
332 } | 332 } |
333 } | 333 } |
334 | 334 |
335 // D4. | 335 // D4. |
336 » » qhatv[len(v)] = mulAddVWW(&qhatv[0], &v[0], qhat, 0, len(v)); | 336 » » qhatv[len(v)] = mulAddVWW(&qhatv[0], &v[0], qhat, 0, len(v)) |
337 | 337 |
338 » » c := subVV(&u[j], &u[j], &qhatv[0], len(qhatv)); | 338 » » c := subVV(&u[j], &u[j], &qhatv[0], len(qhatv)) |
339 if c != 0 { | 339 if c != 0 { |
340 » » » c := addVV(&u[j], &u[j], &v[0], len(v)); | 340 » » » c := addVV(&u[j], &u[j], &v[0], len(v)) |
341 » » » u[j+len(v)] += c; | 341 » » » u[j+len(v)] += c |
342 » » » qhat--; | 342 » » » qhat-- |
343 } | 343 } |
344 | 344 |
345 » » q[j] = qhat; | 345 » » q[j] = qhat |
346 } | 346 } |
347 | 347 |
348 » q = normN(q); | 348 » q = normN(q) |
349 » shiftRight(u, u, shift); | 349 » shiftRight(u, u, shift) |
350 » shiftRight(v, v, shift); | 350 » shiftRight(v, v, shift) |
351 » r = normN(u); | 351 » r = normN(u) |
352 | 352 |
353 » return q, r; | 353 » return q, r |
354 } | 354 } |
355 | 355 |
356 | 356 |
357 // log2 computes the integer binary logarithm of x. | 357 // log2 computes the integer binary logarithm of x. |
358 // The result is the integer n for which 2^n <= x < 2^(n+1). | 358 // The result is the integer n for which 2^n <= x < 2^(n+1). |
359 // If x == 0, the result is -1. | 359 // If x == 0, the result is -1. |
360 func log2(x Word) int { | 360 func log2(x Word) int { |
361 » n := 0; | 361 » n := 0 |
362 for ; x > 0; x >>= 1 { | 362 for ; x > 0; x >>= 1 { |
363 n++ | 363 n++ |
364 } | 364 } |
365 » return n - 1; | 365 » return n - 1 |
366 } | 366 } |
367 | 367 |
368 | 368 |
369 // log2N computes the integer binary logarithm of x. | 369 // log2N computes the integer binary logarithm of x. |
370 // The result is the integer n for which 2^n <= x < 2^(n+1). | 370 // The result is the integer n for which 2^n <= x < 2^(n+1). |
371 // If x == 0, the result is -1. | 371 // If x == 0, the result is -1. |
372 func log2N(x []Word) int { | 372 func log2N(x []Word) int { |
373 » m := len(x); | 373 » m := len(x) |
374 if m > 0 { | 374 if m > 0 { |
375 return (m-1)*_W + log2(x[m-1]) | 375 return (m-1)*_W + log2(x[m-1]) |
376 } | 376 } |
377 » return -1; | 377 » return -1 |
378 } | 378 } |
379 | 379 |
380 | 380 |
381 func hexValue(ch byte) int { | 381 func hexValue(ch byte) int { |
382 » var d byte; | 382 » var d byte |
383 switch { | 383 switch { |
384 case '0' <= ch && ch <= '9': | 384 case '0' <= ch && ch <= '9': |
385 d = ch - '0' | 385 d = ch - '0' |
386 case 'a' <= ch && ch <= 'f': | 386 case 'a' <= ch && ch <= 'f': |
387 d = ch - 'a' + 10 | 387 d = ch - 'a' + 10 |
388 case 'A' <= ch && ch <= 'F': | 388 case 'A' <= ch && ch <= 'F': |
389 d = ch - 'A' + 10 | 389 d = ch - 'A' + 10 |
390 default: | 390 default: |
391 return -1 | 391 return -1 |
392 } | 392 } |
393 » return int(d); | 393 » return int(d) |
394 } | 394 } |
395 | 395 |
396 | 396 |
397 // scanN returns the natural number corresponding to the | 397 // scanN returns the natural number corresponding to the |
398 // longest possible prefix of s representing a natural number in a | 398 // longest possible prefix of s representing a natural number in a |
399 // given conversion base, the actual conversion base used, and the | 399 // given conversion base, the actual conversion base used, and the |
400 // prefix length. The syntax of natural numbers follows the syntax | 400 // prefix length. The syntax of natural numbers follows the syntax |
401 // of unsigned integer literals in Go. | 401 // of unsigned integer literals in Go. |
402 // | 402 // |
403 // If the base argument is 0, the string prefix determines the actual | 403 // If the base argument is 0, the string prefix determines the actual |
404 // conversion base. A prefix of ``0x'' or ``0X'' selects base 16; the | 404 // conversion base. A prefix of ``0x'' or ``0X'' selects base 16; the |
405 // ``0'' prefix selects base 8. Otherwise the selected base is 10. | 405 // ``0'' prefix selects base 8. Otherwise the selected base is 10. |
406 // | 406 // |
407 func scanN(z []Word, s string, base int) ([]Word, int, int) { | 407 func scanN(z []Word, s string, base int) ([]Word, int, int) { |
408 // determine base if necessary | 408 // determine base if necessary |
409 » i, n := 0, len(s); | 409 » i, n := 0, len(s) |
410 if base == 0 { | 410 if base == 0 { |
411 » » base = 10; | 411 » » base = 10 |
412 if n > 0 && s[0] == '0' { | 412 if n > 0 && s[0] == '0' { |
413 if n > 1 && (s[1] == 'x' || s[1] == 'X') { | 413 if n > 1 && (s[1] == 'x' || s[1] == 'X') { |
414 if n == 2 { | 414 if n == 2 { |
415 // Reject a string which is just '0x' as
nonsense. | 415 // Reject a string which is just '0x' as
nonsense. |
416 return nil, 0, 0 | 416 return nil, 0, 0 |
417 } | 417 } |
418 » » » » base, i = 16, 2; | 418 » » » » base, i = 16, 2 |
419 } else { | 419 } else { |
420 base, i = 8, 1 | 420 base, i = 8, 1 |
421 } | 421 } |
422 } | 422 } |
423 } | 423 } |
424 if base < 2 || 16 < base { | 424 if base < 2 || 16 < base { |
425 panic("illegal base") | 425 panic("illegal base") |
426 } | 426 } |
427 | 427 |
428 // convert string | 428 // convert string |
429 » z = makeN(z, len(z), false); | 429 » z = makeN(z, len(z), false) |
430 for ; i < n; i++ { | 430 for ; i < n; i++ { |
431 » » d := hexValue(s[i]); | 431 » » d := hexValue(s[i]) |
432 if 0 <= d && d < base { | 432 if 0 <= d && d < base { |
433 z = mulAddNWW(z, z, Word(base), Word(d)) | 433 z = mulAddNWW(z, z, Word(base), Word(d)) |
434 } else { | 434 } else { |
435 break | 435 break |
436 } | 436 } |
437 } | 437 } |
438 | 438 |
439 » return z, base, i; | 439 » return z, base, i |
440 } | 440 } |
441 | 441 |
442 | 442 |
443 // string converts x to a string for a given base, with 2 <= base <= 16. | 443 // string converts x to a string for a given base, with 2 <= base <= 16. |
444 // TODO(gri) in the style of the other routines, perhaps this should take | 444 // TODO(gri) in the style of the other routines, perhaps this should take |
445 // a []byte buffer and return it | 445 // a []byte buffer and return it |
446 func stringN(x []Word, base int) string { | 446 func stringN(x []Word, base int) string { |
447 if base < 2 || 16 < base { | 447 if base < 2 || 16 < base { |
448 panic("illegal base") | 448 panic("illegal base") |
449 } | 449 } |
450 | 450 |
451 if len(x) == 0 { | 451 if len(x) == 0 { |
452 return "0" | 452 return "0" |
453 } | 453 } |
454 | 454 |
455 // allocate buffer for conversion | 455 // allocate buffer for conversion |
456 » i := (log2N(x)+1)/log2(Word(base)) + 1;»// +1: round up | 456 » i := (log2N(x)+1)/log2(Word(base)) + 1 // +1: round up |
457 » s := make([]byte, i); | 457 » s := make([]byte, i) |
458 | 458 |
459 // don't destroy x | 459 // don't destroy x |
460 » q := setN(nil, x); | 460 » q := setN(nil, x) |
461 | 461 |
462 // convert | 462 // convert |
463 for len(q) > 0 { | 463 for len(q) > 0 { |
464 » » i--; | 464 » » i-- |
465 » » var r Word; | 465 » » var r Word |
466 » » q, r = divNW(q, q, Word(base)); | 466 » » q, r = divNW(q, q, Word(base)) |
467 » » s[i] = "0123456789abcdef"[r]; | 467 » » s[i] = "0123456789abcdef"[r] |
468 } | 468 } |
469 | 469 |
470 » return string(s[i:]); | 470 » return string(s[i:]) |
471 } | 471 } |
472 | 472 |
473 | 473 |
474 // leadingZeroBits returns the number of leading zero bits in x. | 474 // leadingZeroBits returns the number of leading zero bits in x. |
475 func leadingZeroBits(x Word) int { | 475 func leadingZeroBits(x Word) int { |
476 » c := 0; | 476 » c := 0 |
477 if x < 1<<(_W/2) { | 477 if x < 1<<(_W/2) { |
478 » » x <<= _W / 2; | 478 » » x <<= _W / 2 |
479 » » c = _W / 2; | 479 » » c = _W / 2 |
480 } | 480 } |
481 | 481 |
482 for i := 0; x != 0; i++ { | 482 for i := 0; x != 0; i++ { |
483 if x&(1<<(_W-1)) != 0 { | 483 if x&(1<<(_W-1)) != 0 { |
484 return i + c | 484 return i + c |
485 } | 485 } |
486 » » x <<= 1; | 486 » » x <<= 1 |
487 } | 487 } |
488 | 488 |
489 » return _W; | 489 » return _W |
490 } | 490 } |
491 | 491 |
492 const deBruijn32 = 0x077CB531 | 492 const deBruijn32 = 0x077CB531 |
493 | 493 |
494 var deBruijn32Lookup = []byte{ | 494 var deBruijn32Lookup = []byte{ |
495 0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8, | 495 0, 1, 28, 2, 29, 14, 24, 3, 30, 22, 20, 15, 25, 17, 4, 8, |
496 31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9, | 496 31, 27, 13, 23, 21, 19, 16, 7, 26, 12, 18, 6, 11, 5, 10, 9, |
497 } | 497 } |
498 | 498 |
499 const deBruijn64 = 0x03f79d71b4ca8b09 | 499 const deBruijn64 = 0x03f79d71b4ca8b09 |
(...skipping 19 matching lines...) Expand all Loading... |
519 // substring ended up at the top of the word. | 519 // substring ended up at the top of the word. |
520 switch _W { | 520 switch _W { |
521 case 32: | 521 case 32: |
522 return int(deBruijn32Lookup[((x&-x)*deBruijn32)>>27]) | 522 return int(deBruijn32Lookup[((x&-x)*deBruijn32)>>27]) |
523 case 64: | 523 case 64: |
524 return int(deBruijn64Lookup[((x&-x)*(deBruijn64&_M))>>58]) | 524 return int(deBruijn64Lookup[((x&-x)*(deBruijn64&_M))>>58]) |
525 default: | 525 default: |
526 panic("Unknown word size") | 526 panic("Unknown word size") |
527 } | 527 } |
528 | 528 |
529 » return 0; | 529 » return 0 |
530 } | 530 } |
531 | 531 |
532 | 532 |
533 func shiftLeft(dst, src []Word, n int) { | 533 func shiftLeft(dst, src []Word, n int) { |
534 if len(src) == 0 { | 534 if len(src) == 0 { |
535 return | 535 return |
536 } | 536 } |
537 | 537 |
538 » ñ := _W - uint(n); | 538 » ñ := _W - uint(n) |
539 for i := len(src) - 1; i >= 1; i-- { | 539 for i := len(src) - 1; i >= 1; i-- { |
540 » » dst[i] = src[i] << uint(n); | 540 » » dst[i] = src[i] << uint(n) |
541 » » dst[i] |= src[i-1] >> ñ; | 541 » » dst[i] |= src[i-1] >> ñ |
542 } | 542 } |
543 » dst[0] = src[0] << uint(n); | 543 » dst[0] = src[0] << uint(n) |
544 } | 544 } |
545 | 545 |
546 | 546 |
547 func shiftRight(dst, src []Word, n int) { | 547 func shiftRight(dst, src []Word, n int) { |
548 if len(src) == 0 { | 548 if len(src) == 0 { |
549 return | 549 return |
550 } | 550 } |
551 | 551 |
552 » ñ := _W - uint(n); | 552 » ñ := _W - uint(n) |
553 for i := 0; i < len(src)-1; i++ { | 553 for i := 0; i < len(src)-1; i++ { |
554 » » dst[i] = src[i] >> uint(n); | 554 » » dst[i] = src[i] >> uint(n) |
555 » » dst[i] |= src[i+1] << ñ; | 555 » » dst[i] |= src[i+1] << ñ |
556 } | 556 } |
557 » dst[len(src)-1] = src[len(src)-1] >> uint(n); | 557 » dst[len(src)-1] = src[len(src)-1] >> uint(n) |
558 } | 558 } |
559 | 559 |
560 | 560 |
561 // greaterThan returns true iff (x1<<_W + x2) > (y1<<_W + y2) | 561 // greaterThan returns true iff (x1<<_W + x2) > (y1<<_W + y2) |
562 func greaterThan(x1, x2, y1, y2 Word) bool» { return x1 > y1 || x1 == y1 &&
x2 > y2 } | 562 func greaterThan(x1, x2, y1, y2 Word) bool { return x1 > y1 || x1 == y1 && x2 >
y2 } |
563 | 563 |
564 | 564 |
565 // modNW returns x % d. | 565 // modNW returns x % d. |
566 func modNW(x []Word, d Word) (r Word) { | 566 func modNW(x []Word, d Word) (r Word) { |
567 // TODO(agl): we don't actually need to store the q value. | 567 // TODO(agl): we don't actually need to store the q value. |
568 » q := makeN(nil, len(x), false); | 568 » q := makeN(nil, len(x), false) |
569 » return divWVW(&q[0], 0, &x[0], d, len(x)); | 569 » return divWVW(&q[0], 0, &x[0], d, len(x)) |
570 } | 570 } |
571 | 571 |
572 | 572 |
573 // powersOfTwoDecompose finds q and k such that q * 1<<k = n and q is odd. | 573 // powersOfTwoDecompose finds q and k such that q * 1<<k = n and q is odd. |
574 func powersOfTwoDecompose(n []Word) (q []Word, k Word) { | 574 func powersOfTwoDecompose(n []Word) (q []Word, k Word) { |
575 if len(n) == 0 { | 575 if len(n) == 0 { |
576 return n, 0 | 576 return n, 0 |
577 } | 577 } |
578 | 578 |
579 » zeroWords := 0; | 579 » zeroWords := 0 |
580 for n[zeroWords] == 0 { | 580 for n[zeroWords] == 0 { |
581 zeroWords++ | 581 zeroWords++ |
582 } | 582 } |
583 // One of the words must be non-zero by invariant, therefore | 583 // One of the words must be non-zero by invariant, therefore |
584 // zeroWords < len(n). | 584 // zeroWords < len(n). |
585 » x := trailingZeroBits(n[zeroWords]); | 585 » x := trailingZeroBits(n[zeroWords]) |
586 | 586 |
587 » q = makeN(nil, len(n)-zeroWords, false); | 587 » q = makeN(nil, len(n)-zeroWords, false) |
588 » shiftRight(q, n[zeroWords:], x); | 588 » shiftRight(q, n[zeroWords:], x) |
589 | 589 |
590 » k = Word(_W*zeroWords + x); | 590 » k = Word(_W*zeroWords + x) |
591 » return; | 591 » return |
592 } | 592 } |
593 | 593 |
594 | 594 |
595 // randomN creates a random integer in [0..limit), using the space in z if | 595 // randomN creates a random integer in [0..limit), using the space in z if |
596 // possible. n is the bit length of limit. | 596 // possible. n is the bit length of limit. |
597 func randomN(z []Word, rand *rand.Rand, limit []Word, n int) []Word { | 597 func randomN(z []Word, rand *rand.Rand, limit []Word, n int) []Word { |
598 » bitLengthOfMSW := uint(n % _W); | 598 » bitLengthOfMSW := uint(n % _W) |
599 » mask := Word((1 << bitLengthOfMSW) - 1); | 599 » mask := Word((1 << bitLengthOfMSW) - 1) |
600 » z = makeN(z, len(limit), false); | 600 » z = makeN(z, len(limit), false) |
601 | 601 |
602 for { | 602 for { |
603 for i := range z { | 603 for i := range z { |
604 switch _W { | 604 switch _W { |
605 case 32: | 605 case 32: |
606 z[i] = Word(rand.Uint32()) | 606 z[i] = Word(rand.Uint32()) |
607 case 64: | 607 case 64: |
608 z[i] = Word(rand.Uint32()) | Word(rand.Uint32())
<<32 | 608 z[i] = Word(rand.Uint32()) | Word(rand.Uint32())
<<32 |
609 } | 609 } |
610 } | 610 } |
611 | 611 |
612 » » z[len(limit)-1] &= mask; | 612 » » z[len(limit)-1] &= mask |
613 | 613 |
614 if cmpNN(z, limit) < 0 { | 614 if cmpNN(z, limit) < 0 { |
615 break | 615 break |
616 } | 616 } |
617 } | 617 } |
618 | 618 |
619 » return z; | 619 » return z |
620 } | 620 } |
621 | 621 |
622 | 622 |
623 // If m != nil, expNNN calculates x**y mod m. Otherwise it calculates x**y. It | 623 // If m != nil, expNNN calculates x**y mod m. Otherwise it calculates x**y. It |
624 // reuses the storage of z if possible. | 624 // reuses the storage of z if possible. |
625 func expNNN(z, x, y, m []Word) []Word { | 625 func expNNN(z, x, y, m []Word) []Word { |
626 if len(y) == 0 { | 626 if len(y) == 0 { |
627 » » z = makeN(z, 1, false); | 627 » » z = makeN(z, 1, false) |
628 » » z[0] = 1; | 628 » » z[0] = 1 |
629 » » return z; | 629 » » return z |
630 } | 630 } |
631 | 631 |
632 if m != nil { | 632 if m != nil { |
633 // We likely end up being as long as the modulus. | 633 // We likely end up being as long as the modulus. |
634 z = makeN(z, len(m), false) | 634 z = makeN(z, len(m), false) |
635 } | 635 } |
636 » z = setN(z, x); | 636 » z = setN(z, x) |
637 » v := y[len(y)-1]; | 637 » v := y[len(y)-1] |
638 // It's invalid for the most significant word to be zero, therefore we | 638 // It's invalid for the most significant word to be zero, therefore we |
639 // will find a one bit. | 639 // will find a one bit. |
640 » shift := leadingZeros(v) + 1; | 640 » shift := leadingZeros(v) + 1 |
641 » v <<= shift; | 641 » v <<= shift |
642 » var q []Word; | 642 » var q []Word |
643 | 643 |
644 » const mask = 1 << (_W - 1); | 644 » const mask = 1 << (_W - 1) |
645 | 645 |
646 // We walk through the bits of the exponent one by one. Each time we | 646 // We walk through the bits of the exponent one by one. Each time we |
647 // see a bit, we square, thus doubling the power. If the bit is a one, | 647 // see a bit, we square, thus doubling the power. If the bit is a one, |
648 // we also multiply by x, thus adding one to the power. | 648 // we also multiply by x, thus adding one to the power. |
649 | 649 |
650 » w := _W - int(shift); | 650 » w := _W - int(shift) |
651 for j := 0; j < w; j++ { | 651 for j := 0; j < w; j++ { |
652 » » z = mulNN(z, z, z); | 652 » » z = mulNN(z, z, z) |
653 | 653 |
654 if v&mask != 0 { | 654 if v&mask != 0 { |
655 z = mulNN(z, z, x) | 655 z = mulNN(z, z, x) |
656 } | 656 } |
657 | 657 |
658 if m != nil { | 658 if m != nil { |
659 q, z = divNN(q, z, z, m) | 659 q, z = divNN(q, z, z, m) |
660 } | 660 } |
661 | 661 |
662 » » v <<= 1; | 662 » » v <<= 1 |
663 } | 663 } |
664 | 664 |
665 for i := len(y) - 2; i >= 0; i-- { | 665 for i := len(y) - 2; i >= 0; i-- { |
666 » » v = y[i]; | 666 » » v = y[i] |
667 | 667 |
668 for j := 0; j < _W; j++ { | 668 for j := 0; j < _W; j++ { |
669 » » » z = mulNN(z, z, z); | 669 » » » z = mulNN(z, z, z) |
670 | 670 |
671 if v&mask != 0 { | 671 if v&mask != 0 { |
672 z = mulNN(z, z, x) | 672 z = mulNN(z, z, x) |
673 } | 673 } |
674 | 674 |
675 if m != nil { | 675 if m != nil { |
676 q, z = divNN(q, z, z, m) | 676 q, z = divNN(q, z, z, m) |
677 } | 677 } |
678 | 678 |
679 » » » v <<= 1; | 679 » » » v <<= 1 |
680 } | 680 } |
681 } | 681 } |
682 | 682 |
683 » return z; | 683 » return z |
684 } | 684 } |
685 | 685 |
686 | 686 |
687 // lenN returns the bit length of z. | 687 // lenN returns the bit length of z. |
688 func lenN(z []Word) int { | 688 func lenN(z []Word) int { |
689 if len(z) == 0 { | 689 if len(z) == 0 { |
690 return 0 | 690 return 0 |
691 } | 691 } |
692 | 692 |
693 » return (len(z)-1)*_W + (_W - leadingZeroBits(z[len(z)-1])); | 693 » return (len(z)-1)*_W + (_W - leadingZeroBits(z[len(z)-1])) |
694 } | 694 } |
695 | 695 |
696 | 696 |
697 const ( | 697 const ( |
698 » primesProduct32»= 0xC0CFD797;» » // Π {p ∈ primes, 2 < p <= 29} | 698 » primesProduct32 = 0xC0CFD797 // Π {p ∈ primes, 2 < p <= 29} |
699 » primesProduct64»= 0xE221F97C30E94E1D;» // Π {p ∈ primes, 2 < p <= 53} | 699 » primesProduct64 = 0xE221F97C30E94E1D // Π {p ∈ primes, 2 < p <= 53} |
700 ) | 700 ) |
701 | 701 |
702 var bigOne = []Word{1} | 702 var bigOne = []Word{1} |
703 var bigTwo = []Word{2} | 703 var bigTwo = []Word{2} |
704 | 704 |
705 // ProbablyPrime performs n Miller-Rabin tests to check whether n is prime. | 705 // ProbablyPrime performs n Miller-Rabin tests to check whether n is prime. |
706 // If it returns true, n is prime with probability 1 - 1/4^n. | 706 // If it returns true, n is prime with probability 1 - 1/4^n. |
707 // If it returns false, n is not prime. | 707 // If it returns false, n is not prime. |
708 func probablyPrime(n []Word, reps int) bool { | 708 func probablyPrime(n []Word, reps int) bool { |
709 if len(n) == 0 { | 709 if len(n) == 0 { |
710 return false | 710 return false |
711 } | 711 } |
712 | 712 |
713 if len(n) == 1 { | 713 if len(n) == 1 { |
714 if n[0]%2 == 0 { | 714 if n[0]%2 == 0 { |
715 return n[0] == 2 | 715 return n[0] == 2 |
716 } | 716 } |
717 | 717 |
718 // We have to exclude these cases because we reject all | 718 // We have to exclude these cases because we reject all |
719 // multiples of these numbers below. | 719 // multiples of these numbers below. |
720 if n[0] == 3 || n[0] == 5 || n[0] == 7 || n[0] == 11 || | 720 if n[0] == 3 || n[0] == 5 || n[0] == 7 || n[0] == 11 || |
721 n[0] == 13 || n[0] == 17 || n[0] == 19 || n[0] == 23 || | 721 n[0] == 13 || n[0] == 17 || n[0] == 19 || n[0] == 23 || |
722 n[0] == 29 || n[0] == 31 || n[0] == 37 || n[0] == 41 || | 722 n[0] == 29 || n[0] == 31 || n[0] == 37 || n[0] == 41 || |
723 n[0] == 43 || n[0] == 47 || n[0] == 53 { | 723 n[0] == 43 || n[0] == 47 || n[0] == 53 { |
724 return true | 724 return true |
725 } | 725 } |
726 } | 726 } |
727 | 727 |
728 » var r Word; | 728 » var r Word |
729 switch _W { | 729 switch _W { |
730 case 32: | 730 case 32: |
731 r = modNW(n, primesProduct32) | 731 r = modNW(n, primesProduct32) |
732 case 64: | 732 case 64: |
733 r = modNW(n, primesProduct64&_M) | 733 r = modNW(n, primesProduct64&_M) |
734 default: | 734 default: |
735 panic("Unknown word size") | 735 panic("Unknown word size") |
736 } | 736 } |
737 | 737 |
738 if r%3 == 0 || r%5 == 0 || r%7 == 0 || r%11 == 0 || | 738 if r%3 == 0 || r%5 == 0 || r%7 == 0 || r%11 == 0 || |
739 r%13 == 0 || r%17 == 0 || r%19 == 0 || r%23 == 0 || r%29 == 0 { | 739 r%13 == 0 || r%17 == 0 || r%19 == 0 || r%23 == 0 || r%29 == 0 { |
740 return false | 740 return false |
741 } | 741 } |
742 | 742 |
743 if _W == 64 && (r%31 == 0 || r%37 == 0 || r%41 == 0 || | 743 if _W == 64 && (r%31 == 0 || r%37 == 0 || r%41 == 0 || |
744 r%43 == 0 || r%47 == 0 || r%53 == 0) { | 744 r%43 == 0 || r%47 == 0 || r%53 == 0) { |
745 return false | 745 return false |
746 } | 746 } |
747 | 747 |
748 » nm1 := subNN(nil, n, bigOne); | 748 » nm1 := subNN(nil, n, bigOne) |
749 // 1<<k * q = nm1; | 749 // 1<<k * q = nm1; |
750 » q, k := powersOfTwoDecompose(nm1); | 750 » q, k := powersOfTwoDecompose(nm1) |
751 | 751 |
752 » nm3 := subNN(nil, nm1, bigTwo); | 752 » nm3 := subNN(nil, nm1, bigTwo) |
753 » rand := rand.New(rand.NewSource(int64(n[0]))); | 753 » rand := rand.New(rand.NewSource(int64(n[0]))) |
754 | 754 |
755 » var x, y, quotient []Word; | 755 » var x, y, quotient []Word |
756 » nm3Len := lenN(nm3); | 756 » nm3Len := lenN(nm3) |
757 | 757 |
758 NextRandom: | 758 NextRandom: |
759 for i := 0; i < reps; i++ { | 759 for i := 0; i < reps; i++ { |
760 » » x = randomN(x, rand, nm3, nm3Len); | 760 » » x = randomN(x, rand, nm3, nm3Len) |
761 » » addNN(x, x, bigTwo); | 761 » » addNN(x, x, bigTwo) |
762 » » y = expNNN(y, x, q, n); | 762 » » y = expNNN(y, x, q, n) |
763 if cmpNN(y, bigOne) == 0 || cmpNN(y, nm1) == 0 { | 763 if cmpNN(y, bigOne) == 0 || cmpNN(y, nm1) == 0 { |
764 continue | 764 continue |
765 } | 765 } |
766 for j := Word(1); j < k; j++ { | 766 for j := Word(1); j < k; j++ { |
767 » » » y = mulNN(y, y, y); | 767 » » » y = mulNN(y, y, y) |
768 » » » quotient, y = divNN(quotient, y, y, n); | 768 » » » quotient, y = divNN(quotient, y, y, n) |
769 if cmpNN(y, nm1) == 0 { | 769 if cmpNN(y, nm1) == 0 { |
770 continue NextRandom | 770 continue NextRandom |
771 } | 771 } |
772 if cmpNN(y, bigOne) == 0 { | 772 if cmpNN(y, bigOne) == 0 { |
773 return false | 773 return false |
774 } | 774 } |
775 } | 775 } |
776 » » return false; | 776 » » return false |
777 } | 777 } |
778 | 778 |
779 » return true; | 779 » return true |
780 } | 780 } |
OLD | NEW |