github.com/remobjects/goldbaselibrary@v0.0.0-20230924164425-d458680a936b/Source/Gold/math/big/nat_elements.go (about) 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 implements unsigned multi-precision integers (natural 6 // numbers). They are the building blocks for the implementation 7 // of signed integers, rationals, and floating-point numbers. 8 // 9 // Caution: This implementation relies on the function "alias" 10 // which assumes that (nat) slice capacities are never 11 // changed (no 3-operand slice expressions). If that 12 // changes, alias needs to be updated for correctness. 13 14 package big 15 16 import ( 17 "encoding/binary" 18 "math/bits" 19 "math/rand" 20 "sync" 21 ) 22 23 // An unsigned integer x of the form 24 // 25 // x = x[n-1]*_B^(n-1) + x[n-2]*_B^(n-2) + ... + x[1]*_B + x[0] 26 // 27 // with 0 <= x[i] < _B and 0 <= i < n is stored in a slice of length n, 28 // with the digits x[i] as the slice elements. 29 // 30 // A number is normalized if the slice contains no leading 0 digits. 31 // During arithmetic operations, denormalized values may occur but are 32 // always normalized before returning the final result. The normalized 33 // representation of 0 is the empty or nil slice (length = 0). 34 // 35 type nat []Word 36 37 var ( 38 natOne = nat{1} 39 natTwo = nat{2} 40 natFive = nat{5} 41 natTen = nat{10} 42 ) 43 44 func (z nat) clear() { 45 for i := range z { 46 z[i] = 0 47 } 48 } 49 50 func (z nat) norm() nat { 51 i := len(z) 52 for i > 0 && z[i-1] == 0 { 53 i-- 54 } 55 return z[0:i] 56 } 57 58 func (z nat) make(n int) nat { 59 if n <= cap(z) { 60 return z[:n] // reuse z 61 } 62 if n == 1 { 63 // Most nats start small and stay that way; don't over-allocate. 64 return make(nat, 1) 65 } 66 // Choosing a good value for e has significant performance impact 67 // because it increases the chance that a value can be reused. 68 const e = 4 // extra capacity 69 return make(nat, n, n+e) 70 } 71 72 func (z nat) setWord(x Word) nat { 73 if x == 0 { 74 return z[:0] 75 } 76 z = z.make(1) 77 z[0] = x 78 return z 79 } 80 81 func (z nat) setUint64(x uint64) nat { 82 // single-word value 83 if w := Word(x); uint64(w) == x { 84 return z.setWord(w) 85 } 86 // 2-word value 87 z = z.make(2) 88 z[1] = Word(x >> 32) 89 z[0] = Word(x) 90 return z 91 } 92 93 func (z nat) set(x nat) nat { 94 z = z.make(len(x)) 95 copy(z, x) 96 return z 97 } 98 99 func (z nat) add(x, y nat) nat { 100 m := len(x) 101 n := len(y) 102 103 switch { 104 case m < n: 105 return z.add(y, x) 106 case m == 0: 107 // n == 0 because m >= n; result is 0 108 return z[:0] 109 case n == 0: 110 // result is x 111 return z.set(x) 112 } 113 // m > 0 114 115 z = z.make(m + 1) 116 c := addVV(z[0:n], x, y) 117 if m > n { 118 c = addVW(z[n:m], x[n:], c) 119 } 120 z[m] = c 121 122 return z.norm() 123 } 124 125 func (z nat) sub(x, y nat) nat { 126 m := len(x) 127 n := len(y) 128 129 switch { 130 case m < n: 131 panic("underflow") 132 case m == 0: 133 // n == 0 because m >= n; result is 0 134 return z[:0] 135 case n == 0: 136 // result is x 137 return z.set(x) 138 } 139 // m > 0 140 141 z = z.make(m) 142 c := subVV(z[0:n], x, y) 143 if m > n { 144 c = subVW(z[n:], x[n:], c) 145 } 146 if c != 0 { 147 panic("underflow") 148 } 149 150 return z.norm() 151 } 152 153 func (x nat) cmp(y nat) (r int) { 154 m := len(x) 155 n := len(y) 156 if m != n || m == 0 { 157 switch { 158 case m < n: 159 r = -1 160 case m > n: 161 r = 1 162 } 163 return 164 } 165 166 i := m - 1 167 for i > 0 && x[i] == y[i] { 168 i-- 169 } 170 171 switch { 172 case x[i] < y[i]: 173 r = -1 174 case x[i] > y[i]: 175 r = 1 176 } 177 return 178 } 179 180 func (z nat) mulAddWW(x nat, y, r Word) nat { 181 m := len(x) 182 if m == 0 || y == 0 { 183 return z.setWord(r) // result is r 184 } 185 // m > 0 186 187 z = z.make(m + 1) 188 z[m] = mulAddVWW(z[0:m], x, y, r) 189 190 return z.norm() 191 } 192 193 // basicMul multiplies x and y and leaves the result in z. 194 // The (non-normalized) result is placed in z[0 : len(x) + len(y)]. 195 func basicMul(z, x, y nat) { 196 z[0 : len(x)+len(y)].clear() // initialize z 197 for i, d := range y { 198 if d != 0 { 199 z[len(x)+i] = addMulVVW(z[i:i+len(x)], x, d) 200 } 201 } 202 } 203 204 // montgomery computes z mod m = x*y*2**(-n*_W) mod m, 205 // assuming k = -1/m mod 2**_W. 206 // z is used for storing the result which is returned; 207 // z must not alias x, y or m. 208 // See Gueron, "Efficient Software Implementations of Modular Exponentiation". 209 // https://eprint.iacr.org/2011/239.pdf 210 // In the terminology of that paper, this is an "Almost Montgomery Multiplication": 211 // x and y are required to satisfy 0 <= z < 2**(n*_W) and then the result 212 // z is guaranteed to satisfy 0 <= z < 2**(n*_W), but it may not be < m. 213 func (z nat) montgomery(x, y, m nat, k Word, n int) nat { 214 // This code assumes x, y, m are all the same length, n. 215 // (required by addMulVVW and the for loop). 216 // It also assumes that x, y are already reduced mod m, 217 // or else the result will not be properly reduced. 218 if len(x) != n || len(y) != n || len(m) != n { 219 panic("math/big: mismatched montgomery number lengths") 220 } 221 z = z.make(n * 2) 222 z.clear() 223 var c Word 224 for i := 0; i < n; i++ { 225 d := y[i] 226 c2 := addMulVVW(z[i:n+i], x, d) 227 t := z[i] * k 228 c3 := addMulVVW(z[i:n+i], m, t) 229 cx := c + c2 230 cy := cx + c3 231 z[n+i] = cy 232 if cx < c2 || cy < c3 { 233 c = 1 234 } else { 235 c = 0 236 } 237 } 238 if c != 0 { 239 subVV(z[:n], z[n:], m) 240 } else { 241 copy(z[:n], z[n:]) 242 } 243 return z[:n] 244 } 245 246 // Fast version of z[0:n+n>>1].add(z[0:n+n>>1], x[0:n]) w/o bounds checks. 247 // Factored out for readability - do not use outside karatsuba. 248 func karatsubaAdd(z, x nat, n int) { 249 if c := addVV(z[0:n], z, x); c != 0 { 250 addVW(z[n:n+n>>1], z[n:], c) 251 } 252 } 253 254 // Like karatsubaAdd, but does subtract. 255 func karatsubaSub(z, x nat, n int) { 256 if c := subVV(z[0:n], z, x); c != 0 { 257 subVW(z[n:n+n>>1], z[n:], c) 258 } 259 } 260 261 // Operands that are shorter than karatsubaThreshold are multiplied using 262 // "grade school" multiplication; for longer operands the Karatsuba algorithm 263 // is used. 264 var karatsubaThreshold = 40 // computed by calibrate_test.go 265 266 // karatsuba multiplies x and y and leaves the result in z. 267 // Both x and y must have the same length n and n must be a 268 // power of 2. The result vector z must have len(z) >= 6*n. 269 // The (non-normalized) result is placed in z[0 : 2*n]. 270 func karatsuba(z, x, y nat) { 271 n := len(y) 272 273 // Switch to basic multiplication if numbers are odd or small. 274 // (n is always even if karatsubaThreshold is even, but be 275 // conservative) 276 if n&1 != 0 || n < karatsubaThreshold || n < 2 { 277 basicMul(z, x, y) 278 return 279 } 280 // n&1 == 0 && n >= karatsubaThreshold && n >= 2 281 282 // Karatsuba multiplication is based on the observation that 283 // for two numbers x and y with: 284 // 285 // x = x1*b + x0 286 // y = y1*b + y0 287 // 288 // the product x*y can be obtained with 3 products z2, z1, z0 289 // instead of 4: 290 // 291 // x*y = x1*y1*b*b + (x1*y0 + x0*y1)*b + x0*y0 292 // = z2*b*b + z1*b + z0 293 // 294 // with: 295 // 296 // xd = x1 - x0 297 // yd = y0 - y1 298 // 299 // z1 = xd*yd + z2 + z0 300 // = (x1-x0)*(y0 - y1) + z2 + z0 301 // = x1*y0 - x1*y1 - x0*y0 + x0*y1 + z2 + z0 302 // = x1*y0 - z2 - z0 + x0*y1 + z2 + z0 303 // = x1*y0 + x0*y1 304 305 // split x, y into "digits" 306 n2 := n >> 1 // n2 >= 1 307 x1, x0 := x[n2:], x[0:n2] // x = x1*b + y0 308 y1, y0 := y[n2:], y[0:n2] // y = y1*b + y0 309 310 // z is used for the result and temporary storage: 311 // 312 // 6*n 5*n 4*n 3*n 2*n 1*n 0*n 313 // z = [z2 copy|z0 copy| xd*yd | yd:xd | x1*y1 | x0*y0 ] 314 // 315 // For each recursive call of karatsuba, an unused slice of 316 // z is passed in that has (at least) half the length of the 317 // caller's z. 318 319 // compute z0 and z2 with the result "in place" in z 320 karatsuba(z, x0, y0) // z0 = x0*y0 321 karatsuba(z[n:], x1, y1) // z2 = x1*y1 322 323 // compute xd (or the negative value if underflow occurs) 324 s := 1 // sign of product xd*yd 325 xd := z[2*n : 2*n+n2] 326 if subVV(xd, x1, x0) != 0 { // x1-x0 327 s = -s 328 subVV(xd, x0, x1) // x0-x1 329 } 330 331 // compute yd (or the negative value if underflow occurs) 332 yd := z[2*n+n2 : 3*n] 333 if subVV(yd, y0, y1) != 0 { // y0-y1 334 s = -s 335 subVV(yd, y1, y0) // y1-y0 336 } 337 338 // p = (x1-x0)*(y0-y1) == x1*y0 - x1*y1 - x0*y0 + x0*y1 for s > 0 339 // p = (x0-x1)*(y0-y1) == x0*y0 - x0*y1 - x1*y0 + x1*y1 for s < 0 340 p := z[n*3:] 341 karatsuba(p, xd, yd) 342 343 // save original z2:z0 344 // (ok to use upper half of z since we're done recursing) 345 r := z[n*4:] 346 copy(r, z[:n*2]) 347 348 // add up all partial products 349 // 350 // 2*n n 0 351 // z = [ z2 | z0 ] 352 // + [ z0 ] 353 // + [ z2 ] 354 // + [ p ] 355 // 356 karatsubaAdd(z[n2:], r, n) 357 karatsubaAdd(z[n2:], r[n:], n) 358 if s > 0 { 359 karatsubaAdd(z[n2:], p, n) 360 } else { 361 karatsubaSub(z[n2:], p, n) 362 } 363 } 364 365 // alias reports whether x and y share the same base array. 366 // Note: alias assumes that the capacity of underlying arrays 367 // is never changed for nat values; i.e. that there are 368 // no 3-operand slice expressions in this code (or worse, 369 // reflect-based operations to the same effect). 370 func alias(x, y nat) bool { 371 //return cap(x) > 0 && cap(y) > 0 && &x[0:cap(x)][cap(x)-1] == &y[0:cap(y)][cap(y)-1] 372 // elements change 373 return (x != nil) && (y != nil) && (x.fArray == y.fArray) 374 } 375 376 // addAt implements z += x<<(_W*i); z must be long enough. 377 // (we don't use nat.add because we need z to stay the same 378 // slice, and we don't need to normalize z after each addition) 379 func addAt(z, x nat, i int) { 380 if n := len(x); n > 0 { 381 if c := addVV(z[i:i+n], z[i:], x); c != 0 { 382 j := i + n 383 if j < len(z) { 384 addVW(z[j:], z[j:], c) 385 } 386 } 387 } 388 } 389 390 func max(x, y int) int { 391 if x > y { 392 return x 393 } 394 return y 395 } 396 397 // karatsubaLen computes an approximation to the maximum k <= n such that 398 // k = p<<i for a number p <= threshold and an i >= 0. Thus, the 399 // result is the largest number that can be divided repeatedly by 2 before 400 // becoming about the value of threshold. 401 func karatsubaLen(n, threshold int) int { 402 i := uint(0) 403 for n > threshold { 404 n >>= 1 405 i++ 406 } 407 return n << i 408 } 409 410 func (z nat) mul(x, y nat) nat { 411 m := len(x) 412 n := len(y) 413 414 switch { 415 case m < n: 416 return z.mul(y, x) 417 case m == 0 || n == 0: 418 return z[:0] 419 case n == 1: 420 return z.mulAddWW(x, y[0], 0) 421 } 422 // m >= n > 1 423 424 // determine if z can be reused 425 if alias(z, x) || alias(z, y) { 426 z = nil // z is an alias for x or y - cannot reuse 427 } 428 429 // use basic multiplication if the numbers are small 430 if n < karatsubaThreshold { 431 z = z.make(m + n) 432 basicMul(z, x, y) 433 return z.norm() 434 } 435 // m >= n && n >= karatsubaThreshold && n >= 2 436 437 // determine Karatsuba length k such that 438 // 439 // x = xh*b + x0 (0 <= x0 < b) 440 // y = yh*b + y0 (0 <= y0 < b) 441 // b = 1<<(_W*k) ("base" of digits xi, yi) 442 // 443 k := karatsubaLen(n, karatsubaThreshold) 444 // k <= n 445 446 // multiply x0 and y0 via Karatsuba 447 x0 := x[0:k] // x0 is not normalized 448 y0 := y[0:k] // y0 is not normalized 449 z = z.make(max(6*k, m+n)) // enough space for karatsuba of x0*y0 and full result of x*y 450 karatsuba(z, x0, y0) 451 z = z[0 : m+n] // z has final length but may be incomplete 452 z[2*k:].clear() // upper portion of z is garbage (and 2*k <= m+n since k <= n <= m) 453 454 // If xh != 0 or yh != 0, add the missing terms to z. For 455 // 456 // xh = xi*b^i + ... + x2*b^2 + x1*b (0 <= xi < b) 457 // yh = y1*b (0 <= y1 < b) 458 // 459 // the missing terms are 460 // 461 // x0*y1*b and xi*y0*b^i, xi*y1*b^(i+1) for i > 0 462 // 463 // since all the yi for i > 1 are 0 by choice of k: If any of them 464 // were > 0, then yh >= b^2 and thus y >= b^2. Then k' = k*2 would 465 // be a larger valid threshold contradicting the assumption about k. 466 // 467 if k < n || m != n { 468 var t nat 469 470 // add x0*y1*b 471 x0 := x0.norm() 472 y1 := y[k:] // y1 is normalized because y is 473 t = t.mul(x0, y1) // update t so we don't lose t's underlying array 474 addAt(z, t, k) 475 476 // add xi*y0<<i, xi*y1*b<<(i+k) 477 y0 := y0.norm() 478 for i := k; i < len(x); i += k { 479 xi := x[i:] 480 if len(xi) > k { 481 xi = xi[:k] 482 } 483 xi = xi.norm() 484 t = t.mul(xi, y0) 485 addAt(z, t, i) 486 t = t.mul(xi, y1) 487 addAt(z, t, i+k) 488 } 489 } 490 491 return z.norm() 492 } 493 494 // basicSqr sets z = x*x and is asymptotically faster than basicMul 495 // by about a factor of 2, but slower for small arguments due to overhead. 496 // Requirements: len(x) > 0, len(z) == 2*len(x) 497 // The (non-normalized) result is placed in z. 498 func basicSqr(z, x nat) { 499 n := len(x) 500 t := make(nat, 2*n) // temporary variable to hold the products 501 z[1], z[0] = mulWW(x[0], x[0]) // the initial square 502 for i := 1; i < n; i++ { 503 d := x[i] 504 // z collects the squares x[i] * x[i] 505 z[2*i+1], z[2*i] = mulWW(d, d) 506 // t collects the products x[i] * x[j] where j < i 507 t[2*i] = addMulVVW(t[i:2*i], x[0:i], d) 508 } 509 t[2*n-1] = shlVU(t[1:2*n-1], t[1:2*n-1], 1) // double the j < i products 510 addVV(z, z, t) // combine the result 511 } 512 513 // karatsubaSqr squares x and leaves the result in z. 514 // len(x) must be a power of 2 and len(z) >= 6*len(x). 515 // The (non-normalized) result is placed in z[0 : 2*len(x)]. 516 // 517 // The algorithm and the layout of z are the same as for karatsuba. 518 func karatsubaSqr(z, x nat) { 519 n := len(x) 520 521 if n&1 != 0 || n < karatsubaSqrThreshold || n < 2 { 522 basicSqr(z[:2*n], x) 523 return 524 } 525 526 n2 := n >> 1 527 x1, x0 := x[n2:], x[0:n2] 528 529 karatsubaSqr(z, x0) 530 karatsubaSqr(z[n:], x1) 531 532 // s = sign(xd*yd) == -1 for xd != 0; s == 1 for xd == 0 533 xd := z[2*n : 2*n+n2] 534 if subVV(xd, x1, x0) != 0 { 535 subVV(xd, x0, x1) 536 } 537 538 p := z[n*3:] 539 karatsubaSqr(p, xd) 540 541 r := z[n*4:] 542 copy(r, z[:n*2]) 543 544 karatsubaAdd(z[n2:], r, n) 545 karatsubaAdd(z[n2:], r[n:], n) 546 karatsubaSub(z[n2:], p, n) // s == -1 for p != 0; s == 1 for p == 0 547 } 548 549 // Operands that are shorter than basicSqrThreshold are squared using 550 // "grade school" multiplication; for operands longer than karatsubaSqrThreshold 551 // we use the Karatsuba algorithm optimized for x == y. 552 var basicSqrThreshold = 20 // computed by calibrate_test.go 553 var karatsubaSqrThreshold = 260 // computed by calibrate_test.go 554 555 // z = x*x 556 func (z nat) sqr(x nat) nat { 557 n := len(x) 558 switch { 559 case n == 0: 560 return z[:0] 561 case n == 1: 562 d := x[0] 563 z = z.make(2) 564 z[1], z[0] = mulWW(d, d) 565 return z.norm() 566 } 567 568 if alias(z, x) { 569 z = nil // z is an alias for x - cannot reuse 570 } 571 572 if n < basicSqrThreshold { 573 z = z.make(2 * n) 574 basicMul(z, x, x) 575 return z.norm() 576 } 577 if n < karatsubaSqrThreshold { 578 z = z.make(2 * n) 579 basicSqr(z, x) 580 return z.norm() 581 } 582 583 // Use Karatsuba multiplication optimized for x == y. 584 // The algorithm and layout of z are the same as for mul. 585 586 // z = (x1*b + x0)^2 = x1^2*b^2 + 2*x1*x0*b + x0^2 587 588 k := karatsubaLen(n, karatsubaSqrThreshold) 589 590 x0 := x[0:k] 591 z = z.make(max(6*k, 2*n)) 592 karatsubaSqr(z, x0) // z = x0^2 593 z = z[0 : 2*n] 594 z[2*k:].clear() 595 596 if k < n { 597 var t nat 598 x0 := x0.norm() 599 x1 := x[k:] 600 t = t.mul(x0, x1) 601 addAt(z, t, k) 602 addAt(z, t, k) // z = 2*x1*x0*b + x0^2 603 t = t.sqr(x1) 604 addAt(z, t, 2*k) // z = x1^2*b^2 + 2*x1*x0*b + x0^2 605 } 606 607 return z.norm() 608 } 609 610 // mulRange computes the product of all the unsigned integers in the 611 // range [a, b] inclusively. If a > b (empty range), the result is 1. 612 func (z nat) mulRange(a, b uint64) nat { 613 switch { 614 case a == 0: 615 // cut long ranges short (optimization) 616 return z.setUint64(0) 617 case a > b: 618 return z.setUint64(1) 619 case a == b: 620 return z.setUint64(a) 621 case a+1 == b: 622 return z.mul(nat(nil).setUint64(a), nat(nil).setUint64(b)) 623 } 624 m := (a + b) / 2 625 return z.mul(nat(nil).mulRange(a, m), nat(nil).mulRange(m+1, b)) 626 } 627 628 // q = (x-r)/y, with 0 <= r < y 629 func (z nat) divW(x nat, y Word) (q nat, r Word) { 630 m := len(x) 631 switch { 632 case y == 0: 633 panic("division by zero") 634 case y == 1: 635 q = z.set(x) // result is x 636 return 637 case m == 0: 638 q = z[:0] // result is 0 639 return 640 } 641 // m > 0 642 z = z.make(m) 643 r = divWVW(z, 0, x, y) 644 q = z.norm() 645 return 646 } 647 648 func (z nat) div(z2, u, v nat) (q, r nat) { 649 if len(v) == 0 { 650 panic("division by zero") 651 } 652 653 if u.cmp(v) < 0 { 654 q = z[:0] 655 r = z2.set(u) 656 return 657 } 658 659 if len(v) == 1 { 660 var r2 Word 661 q, r2 = z.divW(u, v[0]) 662 r = z2.setWord(r2) 663 return 664 } 665 666 q, r = z.divLarge(z2, u, v) 667 return 668 } 669 670 // getNat returns a *nat of len n. The contents may not be zero. 671 // The pool holds *nat to avoid allocation when converting to interface{}. 672 func getNat(n int) *nat { 673 var z *nat 674 if v := natPool.Get(); v != nil { 675 z = v.(*nat) 676 } 677 if z == nil { 678 z = new(nat) 679 } 680 *z = z.make(n) 681 return z 682 } 683 684 func putNat(x *nat) { 685 natPool.Put(x) 686 } 687 688 var natPool sync.Pool 689 690 // q = (uIn-r)/vIn, with 0 <= r < y 691 // Uses z as storage for q, and u as storage for r if possible. 692 // See Knuth, Volume 2, section 4.3.1, Algorithm D. 693 // Preconditions: 694 // len(vIn) >= 2 695 // len(uIn) >= len(vIn) 696 // u must not alias z 697 func (z nat) divLarge(u, uIn, vIn nat) (q, r nat) { 698 n := len(vIn) 699 m := len(uIn) - n 700 701 // D1. 702 shift := nlz(vIn[n-1]) 703 // do not modify vIn, it may be used by another goroutine simultaneously 704 vp := getNat(n) 705 v := *vp 706 shlVU(v, vIn, shift) 707 708 // u may safely alias uIn or vIn, the value of uIn is used to set u and vIn was already used 709 u = u.make(len(uIn) + 1) 710 u[len(uIn)] = shlVU(u[0:len(uIn)], uIn, shift) 711 712 // z may safely alias uIn or vIn, both values were used already 713 if alias(z, u) { 714 z = nil // z is an alias for u - cannot reuse 715 } 716 q = z.make(m + 1) 717 718 qhatvp := getNat(n + 1) 719 qhatv := *qhatvp 720 721 // D2. 722 vn1 := v[n-1] 723 for j := m; j >= 0; j-- { 724 // D3. 725 qhat := Word(_M) 726 if ujn := u[j+n]; ujn != vn1 { 727 var rhat Word 728 qhat, rhat = divWW(ujn, u[j+n-1], vn1) 729 730 // x1 | x2 = q̂v_{n-2} 731 vn2 := v[n-2] 732 x1, x2 := mulWW(qhat, vn2) 733 // test if q̂v_{n-2} > br̂ + u_{j+n-2} 734 ujn2 := u[j+n-2] 735 for greaterThan(x1, x2, rhat, ujn2) { 736 qhat-- 737 prevRhat := rhat 738 rhat += vn1 739 // v[n-1] >= 0, so this tests for overflow. 740 if rhat < prevRhat { 741 break 742 } 743 x1, x2 = mulWW(qhat, vn2) 744 } 745 } 746 747 // D4. 748 qhatv[n] = mulAddVWW(qhatv[0:n], v, qhat, 0) 749 750 c := subVV(u[j:j+len(qhatv)], u[j:], qhatv) 751 if c != 0 { 752 c := addVV(u[j:j+n], u[j:], v) 753 u[j+n] += c 754 qhat-- 755 } 756 757 q[j] = qhat 758 } 759 760 putNat(vp) 761 putNat(qhatvp) 762 763 q = q.norm() 764 shrVU(u, u, shift) 765 r = u.norm() 766 767 return q, r 768 } 769 770 // Length of x in bits. x must be normalized. 771 func (x nat) bitLen() int { 772 if i := len(x) - 1; i >= 0 { 773 return i*_W + bits.Len(uint(x[i])) 774 } 775 return 0 776 } 777 778 // trailingZeroBits returns the number of consecutive least significant zero 779 // bits of x. 780 func (x nat) trailingZeroBits() uint { 781 if len(x) == 0 { 782 return 0 783 } 784 var i uint 785 for x[i] == 0 { 786 i++ 787 } 788 // x[i] != 0 789 return i*_W + uint(bits.TrailingZeros(uint(x[i]))) 790 } 791 792 func same(x, y nat) bool { 793 return len(x) == len(y) && len(x) > 0 && &x[0] == &y[0] 794 } 795 796 // z = x << s 797 func (z nat) shl(x nat, s uint) nat { 798 if s == 0 { 799 if same(z, x) { 800 return z 801 } 802 if !alias(z, x) { 803 return z.set(x) 804 } 805 } 806 807 m := len(x) 808 if m == 0 { 809 return z[:0] 810 } 811 // m > 0 812 813 n := m + int(s/_W) 814 z = z.make(n + 1) 815 z[n] = shlVU(z[n-m:n], x, s%_W) 816 z[0 : n-m].clear() 817 818 return z.norm() 819 } 820 821 // z = x >> s 822 func (z nat) shr(x nat, s uint) nat { 823 if s == 0 { 824 if same(z, x) { 825 return z 826 } 827 if !alias(z, x) { 828 return z.set(x) 829 } 830 } 831 832 m := len(x) 833 n := m - int(s/_W) 834 if n <= 0 { 835 return z[:0] 836 } 837 // n > 0 838 839 z = z.make(n) 840 shrVU(z, x[m-n:], s%_W) 841 842 return z.norm() 843 } 844 845 func (z nat) setBit(x nat, i uint, b uint) nat { 846 j := int(i / _W) 847 m := Word(1) << (i % _W) 848 n := len(x) 849 switch b { 850 case 0: 851 z = z.make(n) 852 copy(z, x) 853 if j >= n { 854 // no need to grow 855 return z 856 } 857 z[j] &^= m 858 return z.norm() 859 case 1: 860 if j >= n { 861 z = z.make(j + 1) 862 z[n:].clear() 863 } else { 864 z = z.make(n) 865 } 866 copy(z, x) 867 z[j] |= m 868 // no need to normalize 869 return z 870 } 871 panic("set bit is not 0 or 1") 872 } 873 874 // bit returns the value of the i'th bit, with lsb == bit 0. 875 func (x nat) bit(i uint) uint { 876 j := i / _W 877 if j >= uint(len(x)) { 878 return 0 879 } 880 // 0 <= j < len(x) 881 return uint(x[j] >> (i % _W) & 1) 882 } 883 884 // sticky returns 1 if there's a 1 bit within the 885 // i least significant bits, otherwise it returns 0. 886 func (x nat) sticky(i uint) uint { 887 j := i / _W 888 if j >= uint(len(x)) { 889 if len(x) == 0 { 890 return 0 891 } 892 return 1 893 } 894 // 0 <= j < len(x) 895 for _, x := range x[:j] { 896 if x != 0 { 897 return 1 898 } 899 } 900 if x[j]<<(_W-i%_W) != 0 { 901 return 1 902 } 903 return 0 904 } 905 906 func (z nat) and(x, y nat) nat { 907 m := len(x) 908 n := len(y) 909 if m > n { 910 m = n 911 } 912 // m <= n 913 914 z = z.make(m) 915 for i := 0; i < m; i++ { 916 z[i] = x[i] & y[i] 917 } 918 919 return z.norm() 920 } 921 922 func (z nat) andNot(x, y nat) nat { 923 m := len(x) 924 n := len(y) 925 if n > m { 926 n = m 927 } 928 // m >= n 929 930 z = z.make(m) 931 for i := 0; i < n; i++ { 932 z[i] = x[i] &^ y[i] 933 } 934 copy(z[n:m], x[n:m]) 935 936 return z.norm() 937 } 938 939 func (z nat) or(x, y nat) nat { 940 m := len(x) 941 n := len(y) 942 s := x 943 if m < n { 944 n, m = m, n 945 s = y 946 } 947 // m >= n 948 949 z = z.make(m) 950 for i := 0; i < n; i++ { 951 z[i] = x[i] | y[i] 952 } 953 copy(z[n:m], s[n:m]) 954 955 return z.norm() 956 } 957 958 func (z nat) xor(x, y nat) nat { 959 m := len(x) 960 n := len(y) 961 s := x 962 if m < n { 963 n, m = m, n 964 s = y 965 } 966 // m >= n 967 968 z = z.make(m) 969 for i := 0; i < n; i++ { 970 z[i] = x[i] ^ y[i] 971 } 972 copy(z[n:m], s[n:m]) 973 974 return z.norm() 975 } 976 977 // greaterThan reports whether (x1<<_W + x2) > (y1<<_W + y2) 978 func greaterThan(x1, x2, y1, y2 Word) bool { 979 return x1 > y1 || x1 == y1 && x2 > y2 980 } 981 982 // modW returns x % d. 983 func (x nat) modW(d Word) (r Word) { 984 // TODO(agl): we don't actually need to store the q value. 985 var q nat 986 q = q.make(len(x)) 987 return divWVW(q, 0, x, d) 988 } 989 990 // random creates a random integer in [0..limit), using the space in z if 991 // possible. n is the bit length of limit. 992 func (z nat) random(rand *rand.Rand, limit nat, n int) nat { 993 if alias(z, limit) { 994 z = nil // z is an alias for limit - cannot reuse 995 } 996 z = z.make(len(limit)) 997 998 bitLengthOfMSW := uint(n % _W) 999 if bitLengthOfMSW == 0 { 1000 bitLengthOfMSW = _W 1001 } 1002 mask := Word((1 << bitLengthOfMSW) - 1) 1003 1004 for { 1005 switch _W { 1006 case 32: 1007 for i := range z { 1008 z[i] = Word(rand.Uint32()) 1009 } 1010 case 64: 1011 for i := range z { 1012 z[i] = Word(rand.Uint32()) | Word(rand.Uint32())<<32 1013 } 1014 default: 1015 panic("unknown word size") 1016 } 1017 z[len(limit)-1] &= mask 1018 if z.cmp(limit) < 0 { 1019 break 1020 } 1021 } 1022 1023 return z.norm() 1024 } 1025 1026 // If m != 0 (i.e., len(m) != 0), expNN sets z to x**y mod m; 1027 // otherwise it sets z to x**y. The result is the value of z. 1028 func (z nat) expNN(x, y, m nat) nat { 1029 if alias(z, x) || alias(z, y) { 1030 // We cannot allow in-place modification of x or y. 1031 z = nil 1032 } 1033 1034 // x**y mod 1 == 0 1035 if len(m) == 1 && m[0] == 1 { 1036 return z.setWord(0) 1037 } 1038 // m == 0 || m > 1 1039 1040 // x**0 == 1 1041 if len(y) == 0 { 1042 return z.setWord(1) 1043 } 1044 // y > 0 1045 1046 // x**1 mod m == x mod m 1047 if len(y) == 1 && y[0] == 1 && len(m) != 0 { 1048 _, z = nat(nil).div(z, x, m) 1049 return z 1050 } 1051 // y > 1 1052 1053 if len(m) != 0 { 1054 // We likely end up being as long as the modulus. 1055 z = z.make(len(m)) 1056 } 1057 z = z.set(x) 1058 1059 // If the base is non-trivial and the exponent is large, we use 1060 // 4-bit, windowed exponentiation. This involves precomputing 14 values 1061 // (x^2...x^15) but then reduces the number of multiply-reduces by a 1062 // third. Even for a 32-bit exponent, this reduces the number of 1063 // operations. Uses Montgomery method for odd moduli. 1064 if x.cmp(natOne) > 0 && len(y) > 1 && len(m) > 0 { 1065 if m[0]&1 == 1 { 1066 return z.expNNMontgomery(x, y, m) 1067 } 1068 return z.expNNWindowed(x, y, m) 1069 } 1070 1071 v := y[len(y)-1] // v > 0 because y is normalized and y > 0 1072 shift := nlz(v) + 1 1073 v <<= shift 1074 var q nat 1075 1076 const mask = 1 << (_W - 1) 1077 1078 // We walk through the bits of the exponent one by one. Each time we 1079 // see a bit, we square, thus doubling the power. If the bit is a one, 1080 // we also multiply by x, thus adding one to the power. 1081 1082 w := _W - int(shift) 1083 // zz and r are used to avoid allocating in mul and div as 1084 // otherwise the arguments would alias. 1085 var zz, r nat 1086 for j := 0; j < w; j++ { 1087 zz = zz.sqr(z) 1088 zz, z = z, zz 1089 1090 if v&mask != 0 { 1091 zz = zz.mul(z, x) 1092 zz, z = z, zz 1093 } 1094 1095 if len(m) != 0 { 1096 zz, r = zz.div(r, z, m) 1097 zz, r, q, z = q, z, zz, r 1098 } 1099 1100 v <<= 1 1101 } 1102 1103 for i := len(y) - 2; i >= 0; i-- { 1104 v = y[i] 1105 1106 for j := 0; j < _W; j++ { 1107 zz = zz.sqr(z) 1108 zz, z = z, zz 1109 1110 if v&mask != 0 { 1111 zz = zz.mul(z, x) 1112 zz, z = z, zz 1113 } 1114 1115 if len(m) != 0 { 1116 zz, r = zz.div(r, z, m) 1117 zz, r, q, z = q, z, zz, r 1118 } 1119 1120 v <<= 1 1121 } 1122 } 1123 1124 return z.norm() 1125 } 1126 1127 // expNNWindowed calculates x**y mod m using a fixed, 4-bit window. 1128 func (z nat) expNNWindowed(x, y, m nat) nat { 1129 // zz and r are used to avoid allocating in mul and div as otherwise 1130 // the arguments would alias. 1131 var zz, r nat 1132 1133 const n = 4 1134 // powers[i] contains x^i. 1135 var powers [1 << n]nat 1136 powers[0] = natOne 1137 powers[1] = x 1138 for i := 2; i < 1<<n; i += 2 { 1139 p2, p, p1 := &powers[i/2], &powers[i], &powers[i+1] 1140 *p = p.sqr(*p2) 1141 zz, r = zz.div(r, *p, m) 1142 *p, r = r, *p 1143 *p1 = p1.mul(*p, x) 1144 zz, r = zz.div(r, *p1, m) 1145 *p1, r = r, *p1 1146 } 1147 1148 z = z.setWord(1) 1149 1150 for i := len(y) - 1; i >= 0; i-- { 1151 yi := y[i] 1152 for j := 0; j < _W; j += n { 1153 if i != len(y)-1 || j != 0 { 1154 // Unrolled loop for significant performance 1155 // gain. Use go test -bench=".*" in crypto/rsa 1156 // to check performance before making changes. 1157 zz = zz.sqr(z) 1158 zz, z = z, zz 1159 zz, r = zz.div(r, z, m) 1160 z, r = r, z 1161 1162 zz = zz.sqr(z) 1163 zz, z = z, zz 1164 zz, r = zz.div(r, z, m) 1165 z, r = r, z 1166 1167 zz = zz.sqr(z) 1168 zz, z = z, zz 1169 zz, r = zz.div(r, z, m) 1170 z, r = r, z 1171 1172 zz = zz.sqr(z) 1173 zz, z = z, zz 1174 zz, r = zz.div(r, z, m) 1175 z, r = r, z 1176 } 1177 1178 zz = zz.mul(z, powers[yi>>(_W-n)]) 1179 zz, z = z, zz 1180 zz, r = zz.div(r, z, m) 1181 z, r = r, z 1182 1183 yi <<= n 1184 } 1185 } 1186 1187 return z.norm() 1188 } 1189 1190 // expNNMontgomery calculates x**y mod m using a fixed, 4-bit window. 1191 // Uses Montgomery representation. 1192 func (z nat) expNNMontgomery(x, y, m nat) nat { 1193 numWords := len(m) 1194 1195 // We want the lengths of x and m to be equal. 1196 // It is OK if x >= m as long as len(x) == len(m). 1197 if len(x) > numWords { 1198 _, x = nat(nil).div(nil, x, m) 1199 // Note: now len(x) <= numWords, not guaranteed ==. 1200 } 1201 if len(x) < numWords { 1202 rr := make(nat, numWords) 1203 copy(rr, x) 1204 x = rr 1205 } 1206 1207 // Ideally the precomputations would be performed outside, and reused 1208 // k0 = -m**-1 mod 2**_W. Algorithm from: Dumas, J.G. "On Newton–Raphson 1209 // Iteration for Multiplicative Inverses Modulo Prime Powers". 1210 k0 := 2 - m[0] 1211 t := m[0] - 1 1212 for i := 1; i < _W; i <<= 1 { 1213 t *= t 1214 k0 *= (t + 1) 1215 } 1216 k0 = -k0 1217 1218 // RR = 2**(2*_W*len(m)) mod m 1219 RR := nat(nil).setWord(1) 1220 zz := nat(nil).shl(RR, uint(2*numWords*_W)) 1221 _, RR = nat(nil).div(RR, zz, m) 1222 if len(RR) < numWords { 1223 zz = zz.make(numWords) 1224 copy(zz, RR) 1225 RR = zz 1226 } 1227 // one = 1, with equal length to that of m 1228 one := make(nat, numWords) 1229 one[0] = 1 1230 1231 const n = 4 1232 // powers[i] contains x^i 1233 var powers [1 << n]nat 1234 powers[0] = powers[0].montgomery(one, RR, m, k0, numWords) 1235 powers[1] = powers[1].montgomery(x, RR, m, k0, numWords) 1236 for i := 2; i < 1<<n; i++ { 1237 powers[i] = powers[i].montgomery(powers[i-1], powers[1], m, k0, numWords) 1238 } 1239 1240 // initialize z = 1 (Montgomery 1) 1241 z = z.make(numWords) 1242 copy(z, powers[0]) 1243 1244 zz = zz.make(numWords) 1245 1246 // same windowed exponent, but with Montgomery multiplications 1247 for i := len(y) - 1; i >= 0; i-- { 1248 yi := y[i] 1249 for j := 0; j < _W; j += n { 1250 if i != len(y)-1 || j != 0 { 1251 zz = zz.montgomery(z, z, m, k0, numWords) 1252 z = z.montgomery(zz, zz, m, k0, numWords) 1253 zz = zz.montgomery(z, z, m, k0, numWords) 1254 z = z.montgomery(zz, zz, m, k0, numWords) 1255 } 1256 zz = zz.montgomery(z, powers[yi>>(_W-n)], m, k0, numWords) 1257 z, zz = zz, z 1258 yi <<= n 1259 } 1260 } 1261 // convert to regular number 1262 zz = zz.montgomery(z, one, m, k0, numWords) 1263 1264 // One last reduction, just in case. 1265 // See golang.org/issue/13907. 1266 if zz.cmp(m) >= 0 { 1267 // Common case is m has high bit set; in that case, 1268 // since zz is the same length as m, there can be just 1269 // one multiple of m to remove. Just subtract. 1270 // We think that the subtract should be sufficient in general, 1271 // so do that unconditionally, but double-check, 1272 // in case our beliefs are wrong. 1273 // The div is not expected to be reached. 1274 zz = zz.sub(zz, m) 1275 if zz.cmp(m) >= 0 { 1276 _, zz = nat(nil).div(nil, zz, m) 1277 } 1278 } 1279 1280 return zz.norm() 1281 } 1282 1283 // bytes writes the value of z into buf using big-endian encoding. 1284 // len(buf) must be >= len(z)*_S. The value of z is encoded in the 1285 // slice buf[i:]. The number i of unused bytes at the beginning of 1286 // buf is returned as result. 1287 func (z nat) bytes(buf []byte) (i int) { 1288 i = len(buf) 1289 for _, d := range z { 1290 for j := 0; j < _S; j++ { 1291 i-- 1292 buf[i] = byte(d) 1293 d >>= 8 1294 } 1295 } 1296 1297 for i < len(buf) && buf[i] == 0 { 1298 i++ 1299 } 1300 1301 return 1302 } 1303 1304 // bigEndianWord returns the contents of buf interpreted as a big-endian encoded Word value. 1305 func bigEndianWord(buf []byte) Word { 1306 if _W == 64 { 1307 return Word(binary.BigEndian.Uint64(buf)) 1308 } 1309 return Word(binary.BigEndian.Uint32(buf)) 1310 } 1311 1312 // setBytes interprets buf as the bytes of a big-endian unsigned 1313 // integer, sets z to that value, and returns z. 1314 func (z nat) setBytes(buf []byte) nat { 1315 z = z.make((len(buf) + _S - 1) / _S) 1316 1317 i := len(buf) 1318 for k := 0; i >= _S; k++ { 1319 z[k] = bigEndianWord(buf[i-_S : i]) 1320 i -= _S 1321 } 1322 if i > 0 { 1323 var d Word 1324 for s := uint(0); i > 0; s += 8 { 1325 d |= Word(buf[i-1]) << s 1326 i-- 1327 } 1328 z[len(z)-1] = d 1329 } 1330 1331 return z.norm() 1332 } 1333 1334 // sqrt sets z = ⌊√x⌋ 1335 func (z nat) sqrt(x nat) nat { 1336 if x.cmp(natOne) <= 0 { 1337 return z.set(x) 1338 } 1339 if alias(z, x) { 1340 z = nil 1341 } 1342 1343 // Start with value known to be too large and repeat "z = ⌊(z + ⌊x/z⌋)/2⌋" until it stops getting smaller. 1344 // See Brent and Zimmermann, Modern Computer Arithmetic, Algorithm 1.13 (SqrtInt). 1345 // https://members.loria.fr/PZimmermann/mca/pub226.html 1346 // If x is one less than a perfect square, the sequence oscillates between the correct z and z+1; 1347 // otherwise it converges to the correct z and stays there. 1348 var z1, z2 nat 1349 z1 = z 1350 z1 = z1.setUint64(1) 1351 z1 = z1.shl(z1, uint(x.bitLen()+1)/2) // must be ≥ √x 1352 for n := 0; ; n++ { 1353 z2, _ = z2.div(nil, x, z1) 1354 z2 = z2.add(z2, z1) 1355 z2 = z2.shr(z2, 1) 1356 if z2.cmp(z1) >= 0 { 1357 // z1 is answer. 1358 // Figure out whether z1 or z2 is currently aliased to z by looking at loop count. 1359 if n&1 == 0 { 1360 return z1 1361 } 1362 return z.set(z1) 1363 } 1364 z1, z2 = z2, z1 1365 } 1366 }