github.com/goproxy0/go@v0.0.0-20171111080102-49cc0c489d2c/src/math/big/int.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 signed multi-precision integers. 6 7 package big 8 9 import ( 10 "fmt" 11 "io" 12 "math/rand" 13 "strings" 14 ) 15 16 // An Int represents a signed multi-precision integer. 17 // The zero value for an Int represents the value 0. 18 type Int struct { 19 neg bool // sign 20 abs nat // absolute value of the integer 21 } 22 23 var intOne = &Int{false, natOne} 24 25 // Sign returns: 26 // 27 // -1 if x < 0 28 // 0 if x == 0 29 // +1 if x > 0 30 // 31 func (x *Int) Sign() int { 32 if len(x.abs) == 0 { 33 return 0 34 } 35 if x.neg { 36 return -1 37 } 38 return 1 39 } 40 41 // SetInt64 sets z to x and returns z. 42 func (z *Int) SetInt64(x int64) *Int { 43 neg := false 44 if x < 0 { 45 neg = true 46 x = -x 47 } 48 z.abs = z.abs.setUint64(uint64(x)) 49 z.neg = neg 50 return z 51 } 52 53 // SetUint64 sets z to x and returns z. 54 func (z *Int) SetUint64(x uint64) *Int { 55 z.abs = z.abs.setUint64(x) 56 z.neg = false 57 return z 58 } 59 60 // NewInt allocates and returns a new Int set to x. 61 func NewInt(x int64) *Int { 62 return new(Int).SetInt64(x) 63 } 64 65 // Set sets z to x and returns z. 66 func (z *Int) Set(x *Int) *Int { 67 if z != x { 68 z.abs = z.abs.set(x.abs) 69 z.neg = x.neg 70 } 71 return z 72 } 73 74 // Bits provides raw (unchecked but fast) access to x by returning its 75 // absolute value as a little-endian Word slice. The result and x share 76 // the same underlying array. 77 // Bits is intended to support implementation of missing low-level Int 78 // functionality outside this package; it should be avoided otherwise. 79 func (x *Int) Bits() []Word { 80 return x.abs 81 } 82 83 // SetBits provides raw (unchecked but fast) access to z by setting its 84 // value to abs, interpreted as a little-endian Word slice, and returning 85 // z. The result and abs share the same underlying array. 86 // SetBits is intended to support implementation of missing low-level Int 87 // functionality outside this package; it should be avoided otherwise. 88 func (z *Int) SetBits(abs []Word) *Int { 89 z.abs = nat(abs).norm() 90 z.neg = false 91 return z 92 } 93 94 // Abs sets z to |x| (the absolute value of x) and returns z. 95 func (z *Int) Abs(x *Int) *Int { 96 z.Set(x) 97 z.neg = false 98 return z 99 } 100 101 // Neg sets z to -x and returns z. 102 func (z *Int) Neg(x *Int) *Int { 103 z.Set(x) 104 z.neg = len(z.abs) > 0 && !z.neg // 0 has no sign 105 return z 106 } 107 108 // Add sets z to the sum x+y and returns z. 109 func (z *Int) Add(x, y *Int) *Int { 110 neg := x.neg 111 if x.neg == y.neg { 112 // x + y == x + y 113 // (-x) + (-y) == -(x + y) 114 z.abs = z.abs.add(x.abs, y.abs) 115 } else { 116 // x + (-y) == x - y == -(y - x) 117 // (-x) + y == y - x == -(x - y) 118 if x.abs.cmp(y.abs) >= 0 { 119 z.abs = z.abs.sub(x.abs, y.abs) 120 } else { 121 neg = !neg 122 z.abs = z.abs.sub(y.abs, x.abs) 123 } 124 } 125 z.neg = len(z.abs) > 0 && neg // 0 has no sign 126 return z 127 } 128 129 // Sub sets z to the difference x-y and returns z. 130 func (z *Int) Sub(x, y *Int) *Int { 131 neg := x.neg 132 if x.neg != y.neg { 133 // x - (-y) == x + y 134 // (-x) - y == -(x + y) 135 z.abs = z.abs.add(x.abs, y.abs) 136 } else { 137 // x - y == x - y == -(y - x) 138 // (-x) - (-y) == y - x == -(x - y) 139 if x.abs.cmp(y.abs) >= 0 { 140 z.abs = z.abs.sub(x.abs, y.abs) 141 } else { 142 neg = !neg 143 z.abs = z.abs.sub(y.abs, x.abs) 144 } 145 } 146 z.neg = len(z.abs) > 0 && neg // 0 has no sign 147 return z 148 } 149 150 // Mul sets z to the product x*y and returns z. 151 func (z *Int) Mul(x, y *Int) *Int { 152 // x * y == x * y 153 // x * (-y) == -(x * y) 154 // (-x) * y == -(x * y) 155 // (-x) * (-y) == x * y 156 if x == y { 157 z.abs = z.abs.sqr(x.abs) 158 z.neg = false 159 return z 160 } 161 z.abs = z.abs.mul(x.abs, y.abs) 162 z.neg = len(z.abs) > 0 && x.neg != y.neg // 0 has no sign 163 return z 164 } 165 166 // MulRange sets z to the product of all integers 167 // in the range [a, b] inclusively and returns z. 168 // If a > b (empty range), the result is 1. 169 func (z *Int) MulRange(a, b int64) *Int { 170 switch { 171 case a > b: 172 return z.SetInt64(1) // empty range 173 case a <= 0 && b >= 0: 174 return z.SetInt64(0) // range includes 0 175 } 176 // a <= b && (b < 0 || a > 0) 177 178 neg := false 179 if a < 0 { 180 neg = (b-a)&1 == 0 181 a, b = -b, -a 182 } 183 184 z.abs = z.abs.mulRange(uint64(a), uint64(b)) 185 z.neg = neg 186 return z 187 } 188 189 // Binomial sets z to the binomial coefficient of (n, k) and returns z. 190 func (z *Int) Binomial(n, k int64) *Int { 191 // reduce the number of multiplications by reducing k 192 if n/2 < k && k <= n { 193 k = n - k // Binomial(n, k) == Binomial(n, n-k) 194 } 195 var a, b Int 196 a.MulRange(n-k+1, n) 197 b.MulRange(1, k) 198 return z.Quo(&a, &b) 199 } 200 201 // Quo sets z to the quotient x/y for y != 0 and returns z. 202 // If y == 0, a division-by-zero run-time panic occurs. 203 // Quo implements truncated division (like Go); see QuoRem for more details. 204 func (z *Int) Quo(x, y *Int) *Int { 205 z.abs, _ = z.abs.div(nil, x.abs, y.abs) 206 z.neg = len(z.abs) > 0 && x.neg != y.neg // 0 has no sign 207 return z 208 } 209 210 // Rem sets z to the remainder x%y for y != 0 and returns z. 211 // If y == 0, a division-by-zero run-time panic occurs. 212 // Rem implements truncated modulus (like Go); see QuoRem for more details. 213 func (z *Int) Rem(x, y *Int) *Int { 214 _, z.abs = nat(nil).div(z.abs, x.abs, y.abs) 215 z.neg = len(z.abs) > 0 && x.neg // 0 has no sign 216 return z 217 } 218 219 // QuoRem sets z to the quotient x/y and r to the remainder x%y 220 // and returns the pair (z, r) for y != 0. 221 // If y == 0, a division-by-zero run-time panic occurs. 222 // 223 // QuoRem implements T-division and modulus (like Go): 224 // 225 // q = x/y with the result truncated to zero 226 // r = x - y*q 227 // 228 // (See Daan Leijen, ``Division and Modulus for Computer Scientists''.) 229 // See DivMod for Euclidean division and modulus (unlike Go). 230 // 231 func (z *Int) QuoRem(x, y, r *Int) (*Int, *Int) { 232 z.abs, r.abs = z.abs.div(r.abs, x.abs, y.abs) 233 z.neg, r.neg = len(z.abs) > 0 && x.neg != y.neg, len(r.abs) > 0 && x.neg // 0 has no sign 234 return z, r 235 } 236 237 // Div sets z to the quotient x/y for y != 0 and returns z. 238 // If y == 0, a division-by-zero run-time panic occurs. 239 // Div implements Euclidean division (unlike Go); see DivMod for more details. 240 func (z *Int) Div(x, y *Int) *Int { 241 y_neg := y.neg // z may be an alias for y 242 var r Int 243 z.QuoRem(x, y, &r) 244 if r.neg { 245 if y_neg { 246 z.Add(z, intOne) 247 } else { 248 z.Sub(z, intOne) 249 } 250 } 251 return z 252 } 253 254 // Mod sets z to the modulus x%y for y != 0 and returns z. 255 // If y == 0, a division-by-zero run-time panic occurs. 256 // Mod implements Euclidean modulus (unlike Go); see DivMod for more details. 257 func (z *Int) Mod(x, y *Int) *Int { 258 y0 := y // save y 259 if z == y || alias(z.abs, y.abs) { 260 y0 = new(Int).Set(y) 261 } 262 var q Int 263 q.QuoRem(x, y, z) 264 if z.neg { 265 if y0.neg { 266 z.Sub(z, y0) 267 } else { 268 z.Add(z, y0) 269 } 270 } 271 return z 272 } 273 274 // DivMod sets z to the quotient x div y and m to the modulus x mod y 275 // and returns the pair (z, m) for y != 0. 276 // If y == 0, a division-by-zero run-time panic occurs. 277 // 278 // DivMod implements Euclidean division and modulus (unlike Go): 279 // 280 // q = x div y such that 281 // m = x - y*q with 0 <= m < |y| 282 // 283 // (See Raymond T. Boute, ``The Euclidean definition of the functions 284 // div and mod''. ACM Transactions on Programming Languages and 285 // Systems (TOPLAS), 14(2):127-144, New York, NY, USA, 4/1992. 286 // ACM press.) 287 // See QuoRem for T-division and modulus (like Go). 288 // 289 func (z *Int) DivMod(x, y, m *Int) (*Int, *Int) { 290 y0 := y // save y 291 if z == y || alias(z.abs, y.abs) { 292 y0 = new(Int).Set(y) 293 } 294 z.QuoRem(x, y, m) 295 if m.neg { 296 if y0.neg { 297 z.Add(z, intOne) 298 m.Sub(m, y0) 299 } else { 300 z.Sub(z, intOne) 301 m.Add(m, y0) 302 } 303 } 304 return z, m 305 } 306 307 // Cmp compares x and y and returns: 308 // 309 // -1 if x < y 310 // 0 if x == y 311 // +1 if x > y 312 // 313 func (x *Int) Cmp(y *Int) (r int) { 314 // x cmp y == x cmp y 315 // x cmp (-y) == x 316 // (-x) cmp y == y 317 // (-x) cmp (-y) == -(x cmp y) 318 switch { 319 case x.neg == y.neg: 320 r = x.abs.cmp(y.abs) 321 if x.neg { 322 r = -r 323 } 324 case x.neg: 325 r = -1 326 default: 327 r = 1 328 } 329 return 330 } 331 332 // CmpAbs compares the absolute values of x and y and returns: 333 // 334 // -1 if |x| < |y| 335 // 0 if |x| == |y| 336 // +1 if |x| > |y| 337 // 338 func (x *Int) CmpAbs(y *Int) int { 339 return x.abs.cmp(y.abs) 340 } 341 342 // low32 returns the least significant 32 bits of x. 343 func low32(x nat) uint32 { 344 if len(x) == 0 { 345 return 0 346 } 347 return uint32(x[0]) 348 } 349 350 // low64 returns the least significant 64 bits of x. 351 func low64(x nat) uint64 { 352 if len(x) == 0 { 353 return 0 354 } 355 v := uint64(x[0]) 356 if _W == 32 && len(x) > 1 { 357 return uint64(x[1])<<32 | v 358 } 359 return v 360 } 361 362 // Int64 returns the int64 representation of x. 363 // If x cannot be represented in an int64, the result is undefined. 364 func (x *Int) Int64() int64 { 365 v := int64(low64(x.abs)) 366 if x.neg { 367 v = -v 368 } 369 return v 370 } 371 372 // Uint64 returns the uint64 representation of x. 373 // If x cannot be represented in a uint64, the result is undefined. 374 func (x *Int) Uint64() uint64 { 375 return low64(x.abs) 376 } 377 378 // IsInt64 reports whether x can be represented as an int64. 379 func (x *Int) IsInt64() bool { 380 if len(x.abs) <= 64/_W { 381 w := int64(low64(x.abs)) 382 return w >= 0 || x.neg && w == -w 383 } 384 return false 385 } 386 387 // IsUint64 reports whether x can be represented as a uint64. 388 func (x *Int) IsUint64() bool { 389 return !x.neg && len(x.abs) <= 64/_W 390 } 391 392 // SetString sets z to the value of s, interpreted in the given base, 393 // and returns z and a boolean indicating success. The entire string 394 // (not just a prefix) must be valid for success. If SetString fails, 395 // the value of z is undefined but the returned value is nil. 396 // 397 // The base argument must be 0 or a value between 2 and MaxBase. If the base 398 // is 0, the string prefix determines the actual conversion base. A prefix of 399 // ``0x'' or ``0X'' selects base 16; the ``0'' prefix selects base 8, and a 400 // ``0b'' or ``0B'' prefix selects base 2. Otherwise the selected base is 10. 401 // 402 // For bases <= 36, lower and upper case letters are considered the same: 403 // The letters 'a' to 'z' and 'A' to 'Z' represent digit values 10 to 35. 404 // For bases > 36, the upper case letters 'A' to 'Z' represent the digit 405 // values 36 to 61. 406 // 407 func (z *Int) SetString(s string, base int) (*Int, bool) { 408 return z.setFromScanner(strings.NewReader(s), base) 409 } 410 411 // setFromScanner implements SetString given an io.BytesScanner. 412 // For documentation see comments of SetString. 413 func (z *Int) setFromScanner(r io.ByteScanner, base int) (*Int, bool) { 414 if _, _, err := z.scan(r, base); err != nil { 415 return nil, false 416 } 417 // entire content must have been consumed 418 if _, err := r.ReadByte(); err != io.EOF { 419 return nil, false 420 } 421 return z, true // err == io.EOF => scan consumed all content of r 422 } 423 424 // SetBytes interprets buf as the bytes of a big-endian unsigned 425 // integer, sets z to that value, and returns z. 426 func (z *Int) SetBytes(buf []byte) *Int { 427 z.abs = z.abs.setBytes(buf) 428 z.neg = false 429 return z 430 } 431 432 // Bytes returns the absolute value of x as a big-endian byte slice. 433 func (x *Int) Bytes() []byte { 434 buf := make([]byte, len(x.abs)*_S) 435 return buf[x.abs.bytes(buf):] 436 } 437 438 // BitLen returns the length of the absolute value of x in bits. 439 // The bit length of 0 is 0. 440 func (x *Int) BitLen() int { 441 return x.abs.bitLen() 442 } 443 444 // Exp sets z = x**y mod |m| (i.e. the sign of m is ignored), and returns z. 445 // If y <= 0, the result is 1 mod |m|; if m == nil or m == 0, z = x**y. 446 // 447 // Modular exponentation of inputs of a particular size is not a 448 // cryptographically constant-time operation. 449 func (z *Int) Exp(x, y, m *Int) *Int { 450 // See Knuth, volume 2, section 4.6.3. 451 var yWords nat 452 if !y.neg { 453 yWords = y.abs 454 } 455 // y >= 0 456 457 var mWords nat 458 if m != nil { 459 mWords = m.abs // m.abs may be nil for m == 0 460 } 461 462 z.abs = z.abs.expNN(x.abs, yWords, mWords) 463 z.neg = len(z.abs) > 0 && x.neg && len(yWords) > 0 && yWords[0]&1 == 1 // 0 has no sign 464 if z.neg && len(mWords) > 0 { 465 // make modulus result positive 466 z.abs = z.abs.sub(mWords, z.abs) // z == x**y mod |m| && 0 <= z < |m| 467 z.neg = false 468 } 469 470 return z 471 } 472 473 // GCD sets z to the greatest common divisor of a and b, which both must 474 // be > 0, and returns z. 475 // If x or y are not nil, GCD sets their value such that z = a*x + b*y. 476 // If either a or b is <= 0, GCD sets z = x = y = 0. 477 func (z *Int) GCD(x, y, a, b *Int) *Int { 478 if a.Sign() <= 0 || b.Sign() <= 0 { 479 z.SetInt64(0) 480 if x != nil { 481 x.SetInt64(0) 482 } 483 if y != nil { 484 y.SetInt64(0) 485 } 486 return z 487 } 488 if x == nil && y == nil { 489 return z.lehmerGCD(a, b) 490 } 491 492 A := new(Int).Set(a) 493 B := new(Int).Set(b) 494 495 X := new(Int) 496 lastX := new(Int).SetInt64(1) 497 498 q := new(Int) 499 temp := new(Int) 500 501 r := new(Int) 502 for len(B.abs) > 0 { 503 q, r = q.QuoRem(A, B, r) 504 505 A, B, r = B, r, A 506 507 temp.Set(X) 508 X.Mul(X, q) 509 X.Sub(lastX, X) 510 lastX.Set(temp) 511 } 512 513 if x != nil { 514 *x = *lastX 515 } 516 517 if y != nil { 518 // y = (z - a*x)/b 519 y.Mul(a, lastX) 520 y.Sub(A, y) 521 y.Div(y, b) 522 } 523 524 *z = *A 525 return z 526 } 527 528 // lehmerGCD sets z to the greatest common divisor of a and b, 529 // which both must be > 0, and returns z. 530 // See Knuth, The Art of Computer Programming, Vol. 2, Section 4.5.2, Algorithm L. 531 // This implementation uses the improved condition by Collins requiring only one 532 // quotient and avoiding the possibility of single Word overflow. 533 // See Jebelean, "Improving the multiprecision Euclidean algorithm", 534 // Design and Implementation of Symbolic Computation Systems, pp 45-58. 535 func (z *Int) lehmerGCD(a, b *Int) *Int { 536 537 // ensure a >= b 538 if a.abs.cmp(b.abs) < 0 { 539 a, b = b, a 540 } 541 542 // don't destroy incoming values of a and b 543 B := new(Int).Set(b) // must be set first in case b is an alias of z 544 A := z.Set(a) 545 546 // temp variables for multiprecision update 547 t := new(Int) 548 r := new(Int) 549 s := new(Int) 550 w := new(Int) 551 552 // loop invariant A >= B 553 for len(B.abs) > 1 { 554 555 // initialize the digits 556 var a1, a2, u0, u1, u2, v0, v1, v2 Word 557 558 m := len(B.abs) // m >= 2 559 n := len(A.abs) // n >= m >= 2 560 561 // extract the top Word of bits from A and B 562 h := nlz(A.abs[n-1]) 563 a1 = (A.abs[n-1] << h) | (A.abs[n-2] >> (_W - h)) 564 // B may have implicit zero words in the high bits if the lengths differ 565 switch { 566 case n == m: 567 a2 = (B.abs[n-1] << h) | (B.abs[n-2] >> (_W - h)) 568 case n == m+1: 569 a2 = (B.abs[n-2] >> (_W - h)) 570 default: 571 a2 = 0 572 } 573 574 // Since we are calculating with full words to avoid overflow, 575 // we use 'even' to track the sign of the cosequences. 576 // For even iterations: u0, v1 >= 0 && u1, v0 <= 0 577 // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0 578 // The first iteration starts with k=1 (odd). 579 even := false 580 // variables to track the cosequences 581 u0, u1, u2 = 0, 1, 0 582 v0, v1, v2 = 0, 0, 1 583 584 // calculate the quotient and cosequences using Collins' stopping condition 585 for a2 >= v2 && a1-a2 >= v1+v2 { 586 q := a1 / a2 587 a1, a2 = a2, a1-q*a2 588 u0, u1, u2 = u1, u2, u1+q*u2 589 v0, v1, v2 = v1, v2, v1+q*v2 590 even = !even 591 } 592 593 // multiprecision step 594 if v0 != 0 { 595 // simulate the effect of the single precision steps using the cosequences 596 // A = u0*A + v0*B 597 // B = u1*A + v1*B 598 599 t.abs = t.abs.setWord(u0) 600 s.abs = s.abs.setWord(v0) 601 t.neg = !even 602 s.neg = even 603 604 t.Mul(A, t) 605 s.Mul(B, s) 606 607 r.abs = r.abs.setWord(u1) 608 w.abs = w.abs.setWord(v1) 609 r.neg = even 610 w.neg = !even 611 612 r.Mul(A, r) 613 w.Mul(B, w) 614 615 A.Add(t, s) 616 B.Add(r, w) 617 618 } else { 619 // single-digit calculations failed to simluate any quotients 620 // do a standard Euclidean step 621 t.Rem(A, B) 622 A, B, t = B, t, A 623 } 624 } 625 626 if len(B.abs) > 0 { 627 // standard Euclidean algorithm base case for B a single Word 628 if len(A.abs) > 1 { 629 // A is longer than a single Word 630 t.Rem(A, B) 631 A, B, t = B, t, A 632 } 633 if len(B.abs) > 0 { 634 // A and B are both a single Word 635 a1, a2 := A.abs[0], B.abs[0] 636 for a2 != 0 { 637 a1, a2 = a2, a1%a2 638 } 639 A.abs[0] = a1 640 } 641 } 642 *z = *A 643 return z 644 } 645 646 // Rand sets z to a pseudo-random number in [0, n) and returns z. 647 // 648 // As this uses the math/rand package, it must not be used for 649 // security-sensitive work. Use crypto/rand.Int instead. 650 func (z *Int) Rand(rnd *rand.Rand, n *Int) *Int { 651 z.neg = false 652 if n.neg || len(n.abs) == 0 { 653 z.abs = nil 654 return z 655 } 656 z.abs = z.abs.random(rnd, n.abs, n.abs.bitLen()) 657 return z 658 } 659 660 // ModInverse sets z to the multiplicative inverse of g in the ring ℤ/nℤ 661 // and returns z. If g and n are not relatively prime, the result is undefined. 662 func (z *Int) ModInverse(g, n *Int) *Int { 663 if g.neg { 664 // GCD expects parameters a and b to be > 0. 665 var g2 Int 666 g = g2.Mod(g, n) 667 } 668 var d Int 669 d.GCD(z, nil, g, n) 670 // x and y are such that g*x + n*y = d. Since g and n are 671 // relatively prime, d = 1. Taking that modulo n results in 672 // g*x = 1, therefore x is the inverse element. 673 if z.neg { 674 z.Add(z, n) 675 } 676 return z 677 } 678 679 // Jacobi returns the Jacobi symbol (x/y), either +1, -1, or 0. 680 // The y argument must be an odd integer. 681 func Jacobi(x, y *Int) int { 682 if len(y.abs) == 0 || y.abs[0]&1 == 0 { 683 panic(fmt.Sprintf("big: invalid 2nd argument to Int.Jacobi: need odd integer but got %s", y)) 684 } 685 686 // We use the formulation described in chapter 2, section 2.4, 687 // "The Yacas Book of Algorithms": 688 // http://yacas.sourceforge.net/Algo.book.pdf 689 690 var a, b, c Int 691 a.Set(x) 692 b.Set(y) 693 j := 1 694 695 if b.neg { 696 if a.neg { 697 j = -1 698 } 699 b.neg = false 700 } 701 702 for { 703 if b.Cmp(intOne) == 0 { 704 return j 705 } 706 if len(a.abs) == 0 { 707 return 0 708 } 709 a.Mod(&a, &b) 710 if len(a.abs) == 0 { 711 return 0 712 } 713 // a > 0 714 715 // handle factors of 2 in 'a' 716 s := a.abs.trailingZeroBits() 717 if s&1 != 0 { 718 bmod8 := b.abs[0] & 7 719 if bmod8 == 3 || bmod8 == 5 { 720 j = -j 721 } 722 } 723 c.Rsh(&a, s) // a = 2^s*c 724 725 // swap numerator and denominator 726 if b.abs[0]&3 == 3 && c.abs[0]&3 == 3 { 727 j = -j 728 } 729 a.Set(&b) 730 b.Set(&c) 731 } 732 } 733 734 // modSqrt3Mod4 uses the identity 735 // (a^((p+1)/4))^2 mod p 736 // == u^(p+1) mod p 737 // == u^2 mod p 738 // to calculate the square root of any quadratic residue mod p quickly for 3 739 // mod 4 primes. 740 func (z *Int) modSqrt3Mod4Prime(x, p *Int) *Int { 741 e := new(Int).Add(p, intOne) // e = p + 1 742 e.Rsh(e, 2) // e = (p + 1) / 4 743 z.Exp(x, e, p) // z = x^e mod p 744 return z 745 } 746 747 // modSqrtTonelliShanks uses the Tonelli-Shanks algorithm to find the square 748 // root of a quadratic residue modulo any prime. 749 func (z *Int) modSqrtTonelliShanks(x, p *Int) *Int { 750 // Break p-1 into s*2^e such that s is odd. 751 var s Int 752 s.Sub(p, intOne) 753 e := s.abs.trailingZeroBits() 754 s.Rsh(&s, e) 755 756 // find some non-square n 757 var n Int 758 n.SetInt64(2) 759 for Jacobi(&n, p) != -1 { 760 n.Add(&n, intOne) 761 } 762 763 // Core of the Tonelli-Shanks algorithm. Follows the description in 764 // section 6 of "Square roots from 1; 24, 51, 10 to Dan Shanks" by Ezra 765 // Brown: 766 // https://www.maa.org/sites/default/files/pdf/upload_library/22/Polya/07468342.di020786.02p0470a.pdf 767 var y, b, g, t Int 768 y.Add(&s, intOne) 769 y.Rsh(&y, 1) 770 y.Exp(x, &y, p) // y = x^((s+1)/2) 771 b.Exp(x, &s, p) // b = x^s 772 g.Exp(&n, &s, p) // g = n^s 773 r := e 774 for { 775 // find the least m such that ord_p(b) = 2^m 776 var m uint 777 t.Set(&b) 778 for t.Cmp(intOne) != 0 { 779 t.Mul(&t, &t).Mod(&t, p) 780 m++ 781 } 782 783 if m == 0 { 784 return z.Set(&y) 785 } 786 787 t.SetInt64(0).SetBit(&t, int(r-m-1), 1).Exp(&g, &t, p) 788 // t = g^(2^(r-m-1)) mod p 789 g.Mul(&t, &t).Mod(&g, p) // g = g^(2^(r-m)) mod p 790 y.Mul(&y, &t).Mod(&y, p) 791 b.Mul(&b, &g).Mod(&b, p) 792 r = m 793 } 794 } 795 796 // ModSqrt sets z to a square root of x mod p if such a square root exists, and 797 // returns z. The modulus p must be an odd prime. If x is not a square mod p, 798 // ModSqrt leaves z unchanged and returns nil. This function panics if p is 799 // not an odd integer. 800 func (z *Int) ModSqrt(x, p *Int) *Int { 801 switch Jacobi(x, p) { 802 case -1: 803 return nil // x is not a square mod p 804 case 0: 805 return z.SetInt64(0) // sqrt(0) mod p = 0 806 case 1: 807 break 808 } 809 if x.neg || x.Cmp(p) >= 0 { // ensure 0 <= x < p 810 x = new(Int).Mod(x, p) 811 } 812 813 // Check whether p is 3 mod 4, and if so, use the faster algorithm. 814 if len(p.abs) > 0 && p.abs[0]%4 == 3 { 815 return z.modSqrt3Mod4Prime(x, p) 816 } 817 // Otherwise, use Tonelli-Shanks. 818 return z.modSqrtTonelliShanks(x, p) 819 } 820 821 // Lsh sets z = x << n and returns z. 822 func (z *Int) Lsh(x *Int, n uint) *Int { 823 z.abs = z.abs.shl(x.abs, n) 824 z.neg = x.neg 825 return z 826 } 827 828 // Rsh sets z = x >> n and returns z. 829 func (z *Int) Rsh(x *Int, n uint) *Int { 830 if x.neg { 831 // (-x) >> s == ^(x-1) >> s == ^((x-1) >> s) == -(((x-1) >> s) + 1) 832 t := z.abs.sub(x.abs, natOne) // no underflow because |x| > 0 833 t = t.shr(t, n) 834 z.abs = t.add(t, natOne) 835 z.neg = true // z cannot be zero if x is negative 836 return z 837 } 838 839 z.abs = z.abs.shr(x.abs, n) 840 z.neg = false 841 return z 842 } 843 844 // Bit returns the value of the i'th bit of x. That is, it 845 // returns (x>>i)&1. The bit index i must be >= 0. 846 func (x *Int) Bit(i int) uint { 847 if i == 0 { 848 // optimization for common case: odd/even test of x 849 if len(x.abs) > 0 { 850 return uint(x.abs[0] & 1) // bit 0 is same for -x 851 } 852 return 0 853 } 854 if i < 0 { 855 panic("negative bit index") 856 } 857 if x.neg { 858 t := nat(nil).sub(x.abs, natOne) 859 return t.bit(uint(i)) ^ 1 860 } 861 862 return x.abs.bit(uint(i)) 863 } 864 865 // SetBit sets z to x, with x's i'th bit set to b (0 or 1). 866 // That is, if b is 1 SetBit sets z = x | (1 << i); 867 // if b is 0 SetBit sets z = x &^ (1 << i). If b is not 0 or 1, 868 // SetBit will panic. 869 func (z *Int) SetBit(x *Int, i int, b uint) *Int { 870 if i < 0 { 871 panic("negative bit index") 872 } 873 if x.neg { 874 t := z.abs.sub(x.abs, natOne) 875 t = t.setBit(t, uint(i), b^1) 876 z.abs = t.add(t, natOne) 877 z.neg = len(z.abs) > 0 878 return z 879 } 880 z.abs = z.abs.setBit(x.abs, uint(i), b) 881 z.neg = false 882 return z 883 } 884 885 // And sets z = x & y and returns z. 886 func (z *Int) And(x, y *Int) *Int { 887 if x.neg == y.neg { 888 if x.neg { 889 // (-x) & (-y) == ^(x-1) & ^(y-1) == ^((x-1) | (y-1)) == -(((x-1) | (y-1)) + 1) 890 x1 := nat(nil).sub(x.abs, natOne) 891 y1 := nat(nil).sub(y.abs, natOne) 892 z.abs = z.abs.add(z.abs.or(x1, y1), natOne) 893 z.neg = true // z cannot be zero if x and y are negative 894 return z 895 } 896 897 // x & y == x & y 898 z.abs = z.abs.and(x.abs, y.abs) 899 z.neg = false 900 return z 901 } 902 903 // x.neg != y.neg 904 if x.neg { 905 x, y = y, x // & is symmetric 906 } 907 908 // x & (-y) == x & ^(y-1) == x &^ (y-1) 909 y1 := nat(nil).sub(y.abs, natOne) 910 z.abs = z.abs.andNot(x.abs, y1) 911 z.neg = false 912 return z 913 } 914 915 // AndNot sets z = x &^ y and returns z. 916 func (z *Int) AndNot(x, y *Int) *Int { 917 if x.neg == y.neg { 918 if x.neg { 919 // (-x) &^ (-y) == ^(x-1) &^ ^(y-1) == ^(x-1) & (y-1) == (y-1) &^ (x-1) 920 x1 := nat(nil).sub(x.abs, natOne) 921 y1 := nat(nil).sub(y.abs, natOne) 922 z.abs = z.abs.andNot(y1, x1) 923 z.neg = false 924 return z 925 } 926 927 // x &^ y == x &^ y 928 z.abs = z.abs.andNot(x.abs, y.abs) 929 z.neg = false 930 return z 931 } 932 933 if x.neg { 934 // (-x) &^ y == ^(x-1) &^ y == ^(x-1) & ^y == ^((x-1) | y) == -(((x-1) | y) + 1) 935 x1 := nat(nil).sub(x.abs, natOne) 936 z.abs = z.abs.add(z.abs.or(x1, y.abs), natOne) 937 z.neg = true // z cannot be zero if x is negative and y is positive 938 return z 939 } 940 941 // x &^ (-y) == x &^ ^(y-1) == x & (y-1) 942 y1 := nat(nil).sub(y.abs, natOne) 943 z.abs = z.abs.and(x.abs, y1) 944 z.neg = false 945 return z 946 } 947 948 // Or sets z = x | y and returns z. 949 func (z *Int) Or(x, y *Int) *Int { 950 if x.neg == y.neg { 951 if x.neg { 952 // (-x) | (-y) == ^(x-1) | ^(y-1) == ^((x-1) & (y-1)) == -(((x-1) & (y-1)) + 1) 953 x1 := nat(nil).sub(x.abs, natOne) 954 y1 := nat(nil).sub(y.abs, natOne) 955 z.abs = z.abs.add(z.abs.and(x1, y1), natOne) 956 z.neg = true // z cannot be zero if x and y are negative 957 return z 958 } 959 960 // x | y == x | y 961 z.abs = z.abs.or(x.abs, y.abs) 962 z.neg = false 963 return z 964 } 965 966 // x.neg != y.neg 967 if x.neg { 968 x, y = y, x // | is symmetric 969 } 970 971 // x | (-y) == x | ^(y-1) == ^((y-1) &^ x) == -(^((y-1) &^ x) + 1) 972 y1 := nat(nil).sub(y.abs, natOne) 973 z.abs = z.abs.add(z.abs.andNot(y1, x.abs), natOne) 974 z.neg = true // z cannot be zero if one of x or y is negative 975 return z 976 } 977 978 // Xor sets z = x ^ y and returns z. 979 func (z *Int) Xor(x, y *Int) *Int { 980 if x.neg == y.neg { 981 if x.neg { 982 // (-x) ^ (-y) == ^(x-1) ^ ^(y-1) == (x-1) ^ (y-1) 983 x1 := nat(nil).sub(x.abs, natOne) 984 y1 := nat(nil).sub(y.abs, natOne) 985 z.abs = z.abs.xor(x1, y1) 986 z.neg = false 987 return z 988 } 989 990 // x ^ y == x ^ y 991 z.abs = z.abs.xor(x.abs, y.abs) 992 z.neg = false 993 return z 994 } 995 996 // x.neg != y.neg 997 if x.neg { 998 x, y = y, x // ^ is symmetric 999 } 1000 1001 // x ^ (-y) == x ^ ^(y-1) == ^(x ^ (y-1)) == -((x ^ (y-1)) + 1) 1002 y1 := nat(nil).sub(y.abs, natOne) 1003 z.abs = z.abs.add(z.abs.xor(x.abs, y1), natOne) 1004 z.neg = true // z cannot be zero if only one of x or y is negative 1005 return z 1006 } 1007 1008 // Not sets z = ^x and returns z. 1009 func (z *Int) Not(x *Int) *Int { 1010 if x.neg { 1011 // ^(-x) == ^(^(x-1)) == x-1 1012 z.abs = z.abs.sub(x.abs, natOne) 1013 z.neg = false 1014 return z 1015 } 1016 1017 // ^x == -x-1 == -(x+1) 1018 z.abs = z.abs.add(x.abs, natOne) 1019 z.neg = true // z cannot be zero if x is positive 1020 return z 1021 } 1022 1023 // Sqrt sets z to ⌊√x⌋, the largest integer such that z² ≤ x, and returns z. 1024 // It panics if x is negative. 1025 func (z *Int) Sqrt(x *Int) *Int { 1026 if x.neg { 1027 panic("square root of negative number") 1028 } 1029 z.neg = false 1030 z.abs = z.abs.sqrt(x.abs) 1031 return z 1032 }