github.com/hikaru7719/go@v0.0.0-20181025140707-c8b2ac68906a/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 m == nil or m == 0, z = x**y unless y <= 0 then z = 1. 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 xWords := x.abs 452 if y.neg { 453 if m == nil || len(m.abs) == 0 { 454 return z.SetInt64(1) 455 } 456 // for y < 0: x**y mod m == (x**(-1))**|y| mod m 457 xWords = new(Int).ModInverse(x, m).abs 458 } 459 yWords := y.abs 460 461 var mWords nat 462 if m != nil { 463 mWords = m.abs // m.abs may be nil for m == 0 464 } 465 466 z.abs = z.abs.expNN(xWords, yWords, mWords) 467 z.neg = len(z.abs) > 0 && x.neg && len(yWords) > 0 && yWords[0]&1 == 1 // 0 has no sign 468 if z.neg && len(mWords) > 0 { 469 // make modulus result positive 470 z.abs = z.abs.sub(mWords, z.abs) // z == x**y mod |m| && 0 <= z < |m| 471 z.neg = false 472 } 473 474 return z 475 } 476 477 // GCD sets z to the greatest common divisor of a and b, which both must 478 // be > 0, and returns z. 479 // If x or y are not nil, GCD sets their value such that z = a*x + b*y. 480 // If either a or b is <= 0, GCD sets z = x = y = 0. 481 func (z *Int) GCD(x, y, a, b *Int) *Int { 482 if a.Sign() <= 0 || b.Sign() <= 0 { 483 z.SetInt64(0) 484 if x != nil { 485 x.SetInt64(0) 486 } 487 if y != nil { 488 y.SetInt64(0) 489 } 490 return z 491 } 492 493 return z.lehmerGCD(x, y, a, b) 494 } 495 496 // lehmerSimulate attempts to simulate several Euclidean update steps 497 // using the leading digits of A and B. It returns u0, u1, v0, v1 498 // such that A and B can be updated as: 499 // A = u0*A + v0*B 500 // B = u1*A + v1*B 501 // Requirements: A >= B and len(B.abs) >= 2 502 // Since we are calculating with full words to avoid overflow, 503 // we use 'even' to track the sign of the cosequences. 504 // For even iterations: u0, v1 >= 0 && u1, v0 <= 0 505 // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0 506 func lehmerSimulate(A, B *Int) (u0, u1, v0, v1 Word, even bool) { 507 // initialize the digits 508 var a1, a2, u2, v2 Word 509 510 m := len(B.abs) // m >= 2 511 n := len(A.abs) // n >= m >= 2 512 513 // extract the top Word of bits from A and B 514 h := nlz(A.abs[n-1]) 515 a1 = A.abs[n-1]<<h | A.abs[n-2]>>(_W-h) 516 // B may have implicit zero words in the high bits if the lengths differ 517 switch { 518 case n == m: 519 a2 = B.abs[n-1]<<h | B.abs[n-2]>>(_W-h) 520 case n == m+1: 521 a2 = B.abs[n-2] >> (_W - h) 522 default: 523 a2 = 0 524 } 525 526 // Since we are calculating with full words to avoid overflow, 527 // we use 'even' to track the sign of the cosequences. 528 // For even iterations: u0, v1 >= 0 && u1, v0 <= 0 529 // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0 530 // The first iteration starts with k=1 (odd). 531 even = false 532 // variables to track the cosequences 533 u0, u1, u2 = 0, 1, 0 534 v0, v1, v2 = 0, 0, 1 535 536 // Calculate the quotient and cosequences using Collins' stopping condition. 537 // Note that overflow of a Word is not possible when computing the remainder 538 // sequence and cosequences since the cosequence size is bounded by the input size. 539 // See section 4.2 of Jebelean for details. 540 for a2 >= v2 && a1-a2 >= v1+v2 { 541 q, r := a1/a2, a1%a2 542 a1, a2 = a2, r 543 u0, u1, u2 = u1, u2, u1+q*u2 544 v0, v1, v2 = v1, v2, v1+q*v2 545 even = !even 546 } 547 return 548 } 549 550 // lehmerUpdate updates the inputs A and B such that: 551 // A = u0*A + v0*B 552 // B = u1*A + v1*B 553 // where the signs of u0, u1, v0, v1 are given by even 554 // For even == true: u0, v1 >= 0 && u1, v0 <= 0 555 // For even == false: u0, v1 <= 0 && u1, v0 >= 0 556 // q, r, s, t are temporary variables to avoid allocations in the multiplication 557 func lehmerUpdate(A, B, q, r, s, t *Int, u0, u1, v0, v1 Word, even bool) { 558 559 t.abs = t.abs.setWord(u0) 560 s.abs = s.abs.setWord(v0) 561 t.neg = !even 562 s.neg = even 563 564 t.Mul(A, t) 565 s.Mul(B, s) 566 567 r.abs = r.abs.setWord(u1) 568 q.abs = q.abs.setWord(v1) 569 r.neg = even 570 q.neg = !even 571 572 r.Mul(A, r) 573 q.Mul(B, q) 574 575 A.Add(t, s) 576 B.Add(r, q) 577 } 578 579 // euclidUpdate performs a single step of the Euclidean GCD algorithm 580 // if extended is true, it also updates the cosequence Ua, Ub 581 func euclidUpdate(A, B, Ua, Ub, q, r, s, t *Int, extended bool) { 582 q, r = q.QuoRem(A, B, r) 583 584 *A, *B, *r = *B, *r, *A 585 586 if extended { 587 // Ua, Ub = Ub, Ua - q*Ub 588 t.Set(Ub) 589 s.Mul(Ub, q) 590 Ub.Sub(Ua, s) 591 Ua.Set(t) 592 } 593 } 594 595 // lehmerGCD sets z to the greatest common divisor of a and b, 596 // which both must be > 0, and returns z. 597 // If x or y are not nil, their values are set such that z = a*x + b*y. 598 // See Knuth, The Art of Computer Programming, Vol. 2, Section 4.5.2, Algorithm L. 599 // This implementation uses the improved condition by Collins requiring only one 600 // quotient and avoiding the possibility of single Word overflow. 601 // See Jebelean, "Improving the multiprecision Euclidean algorithm", 602 // Design and Implementation of Symbolic Computation Systems, pp 45-58. 603 // The cosequences are updated according to Algorithm 10.45 from 604 // Cohen et al. "Handbook of Elliptic and Hyperelliptic Curve Cryptography" pp 192. 605 func (z *Int) lehmerGCD(x, y, a, b *Int) *Int { 606 var A, B, Ua, Ub *Int 607 608 A = new(Int).Set(a) 609 B = new(Int).Set(b) 610 611 extended := x != nil || y != nil 612 613 if extended { 614 // Ua (Ub) tracks how many times input a has been accumulated into A (B). 615 Ua = new(Int).SetInt64(1) 616 Ub = new(Int) 617 } 618 619 // temp variables for multiprecision update 620 q := new(Int) 621 r := new(Int) 622 s := new(Int) 623 t := new(Int) 624 625 // ensure A >= B 626 if A.abs.cmp(B.abs) < 0 { 627 A, B = B, A 628 Ub, Ua = Ua, Ub 629 } 630 631 // loop invariant A >= B 632 for len(B.abs) > 1 { 633 // Attempt to calculate in single-precision using leading words of A and B. 634 u0, u1, v0, v1, even := lehmerSimulate(A, B) 635 636 // multiprecision Step 637 if v0 != 0 { 638 // Simulate the effect of the single-precision steps using the cosequences. 639 // A = u0*A + v0*B 640 // B = u1*A + v1*B 641 lehmerUpdate(A, B, q, r, s, t, u0, u1, v0, v1, even) 642 643 if extended { 644 // Ua = u0*Ua + v0*Ub 645 // Ub = u1*Ua + v1*Ub 646 lehmerUpdate(Ua, Ub, q, r, s, t, u0, u1, v0, v1, even) 647 } 648 649 } else { 650 // Single-digit calculations failed to simulate any quotients. 651 // Do a standard Euclidean step. 652 euclidUpdate(A, B, Ua, Ub, q, r, s, t, extended) 653 } 654 } 655 656 if len(B.abs) > 0 { 657 // extended Euclidean algorithm base case if B is a single Word 658 if len(A.abs) > 1 { 659 // A is longer than a single Word, so one update is needed. 660 euclidUpdate(A, B, Ua, Ub, q, r, s, t, extended) 661 } 662 if len(B.abs) > 0 { 663 // A and B are both a single Word. 664 aWord, bWord := A.abs[0], B.abs[0] 665 if extended { 666 var ua, ub, va, vb Word 667 ua, ub = 1, 0 668 va, vb = 0, 1 669 even := true 670 for bWord != 0 { 671 q, r := aWord/bWord, aWord%bWord 672 aWord, bWord = bWord, r 673 ua, ub = ub, ua+q*ub 674 va, vb = vb, va+q*vb 675 even = !even 676 } 677 678 t.abs = t.abs.setWord(ua) 679 s.abs = s.abs.setWord(va) 680 t.neg = !even 681 s.neg = even 682 683 t.Mul(Ua, t) 684 s.Mul(Ub, s) 685 686 Ua.Add(t, s) 687 } else { 688 for bWord != 0 { 689 aWord, bWord = bWord, aWord%bWord 690 } 691 } 692 A.abs[0] = aWord 693 } 694 } 695 696 if x != nil { 697 *x = *Ua 698 } 699 700 if y != nil { 701 // y = (z - a*x)/b 702 y.Mul(a, Ua) 703 y.Sub(A, y) 704 y.Div(y, b) 705 } 706 707 *z = *A 708 709 return z 710 } 711 712 // Rand sets z to a pseudo-random number in [0, n) and returns z. 713 // 714 // As this uses the math/rand package, it must not be used for 715 // security-sensitive work. Use crypto/rand.Int instead. 716 func (z *Int) Rand(rnd *rand.Rand, n *Int) *Int { 717 z.neg = false 718 if n.neg || len(n.abs) == 0 { 719 z.abs = nil 720 return z 721 } 722 z.abs = z.abs.random(rnd, n.abs, n.abs.bitLen()) 723 return z 724 } 725 726 // ModInverse sets z to the multiplicative inverse of g in the ring ℤ/nℤ 727 // and returns z. If g and n are not relatively prime, g has no multiplicative 728 // inverse in the ring ℤ/nℤ. In this case, z is unchanged and the return value 729 // is nil. 730 func (z *Int) ModInverse(g, n *Int) *Int { 731 // GCD expects parameters a and b to be > 0. 732 if n.neg { 733 var n2 Int 734 n = n2.Neg(n) 735 } 736 if g.neg { 737 var g2 Int 738 g = g2.Mod(g, n) 739 } 740 var d, x Int 741 d.GCD(&x, nil, g, n) 742 743 // if and only if d==1, g and n are relatively prime 744 if d.Cmp(intOne) != 0 { 745 return nil 746 } 747 748 // x and y are such that g*x + n*y = 1, therefore x is the inverse element, 749 // but it may be negative, so convert to the range 0 <= z < |n| 750 if x.neg { 751 z.Add(&x, n) 752 } else { 753 z.Set(&x) 754 } 755 return z 756 } 757 758 // Jacobi returns the Jacobi symbol (x/y), either +1, -1, or 0. 759 // The y argument must be an odd integer. 760 func Jacobi(x, y *Int) int { 761 if len(y.abs) == 0 || y.abs[0]&1 == 0 { 762 panic(fmt.Sprintf("big: invalid 2nd argument to Int.Jacobi: need odd integer but got %s", y)) 763 } 764 765 // We use the formulation described in chapter 2, section 2.4, 766 // "The Yacas Book of Algorithms": 767 // http://yacas.sourceforge.net/Algo.book.pdf 768 769 var a, b, c Int 770 a.Set(x) 771 b.Set(y) 772 j := 1 773 774 if b.neg { 775 if a.neg { 776 j = -1 777 } 778 b.neg = false 779 } 780 781 for { 782 if b.Cmp(intOne) == 0 { 783 return j 784 } 785 if len(a.abs) == 0 { 786 return 0 787 } 788 a.Mod(&a, &b) 789 if len(a.abs) == 0 { 790 return 0 791 } 792 // a > 0 793 794 // handle factors of 2 in 'a' 795 s := a.abs.trailingZeroBits() 796 if s&1 != 0 { 797 bmod8 := b.abs[0] & 7 798 if bmod8 == 3 || bmod8 == 5 { 799 j = -j 800 } 801 } 802 c.Rsh(&a, s) // a = 2^s*c 803 804 // swap numerator and denominator 805 if b.abs[0]&3 == 3 && c.abs[0]&3 == 3 { 806 j = -j 807 } 808 a.Set(&b) 809 b.Set(&c) 810 } 811 } 812 813 // modSqrt3Mod4 uses the identity 814 // (a^((p+1)/4))^2 mod p 815 // == u^(p+1) mod p 816 // == u^2 mod p 817 // to calculate the square root of any quadratic residue mod p quickly for 3 818 // mod 4 primes. 819 func (z *Int) modSqrt3Mod4Prime(x, p *Int) *Int { 820 e := new(Int).Add(p, intOne) // e = p + 1 821 e.Rsh(e, 2) // e = (p + 1) / 4 822 z.Exp(x, e, p) // z = x^e mod p 823 return z 824 } 825 826 // modSqrt5Mod8 uses Atkin's observation that 2 is not a square mod p 827 // alpha == (2*a)^((p-5)/8) mod p 828 // beta == 2*a*alpha^2 mod p is a square root of -1 829 // b == a*alpha*(beta-1) mod p is a square root of a 830 // to calculate the square root of any quadratic residue mod p quickly for 5 831 // mod 8 primes. 832 func (z *Int) modSqrt5Mod8Prime(x, p *Int) *Int { 833 // p == 5 mod 8 implies p = e*8 + 5 834 // e is the quotient and 5 the remainder on division by 8 835 e := new(Int).Rsh(p, 3) // e = (p - 5) / 8 836 tx := new(Int).Lsh(x, 1) // tx = 2*x 837 alpha := new(Int).Exp(tx, e, p) 838 beta := new(Int).Mul(alpha, alpha) 839 beta.Mod(beta, p) 840 beta.Mul(beta, tx) 841 beta.Mod(beta, p) 842 beta.Sub(beta, intOne) 843 beta.Mul(beta, x) 844 beta.Mod(beta, p) 845 beta.Mul(beta, alpha) 846 z.Mod(beta, p) 847 return z 848 } 849 850 // modSqrtTonelliShanks uses the Tonelli-Shanks algorithm to find the square 851 // root of a quadratic residue modulo any prime. 852 func (z *Int) modSqrtTonelliShanks(x, p *Int) *Int { 853 // Break p-1 into s*2^e such that s is odd. 854 var s Int 855 s.Sub(p, intOne) 856 e := s.abs.trailingZeroBits() 857 s.Rsh(&s, e) 858 859 // find some non-square n 860 var n Int 861 n.SetInt64(2) 862 for Jacobi(&n, p) != -1 { 863 n.Add(&n, intOne) 864 } 865 866 // Core of the Tonelli-Shanks algorithm. Follows the description in 867 // section 6 of "Square roots from 1; 24, 51, 10 to Dan Shanks" by Ezra 868 // Brown: 869 // https://www.maa.org/sites/default/files/pdf/upload_library/22/Polya/07468342.di020786.02p0470a.pdf 870 var y, b, g, t Int 871 y.Add(&s, intOne) 872 y.Rsh(&y, 1) 873 y.Exp(x, &y, p) // y = x^((s+1)/2) 874 b.Exp(x, &s, p) // b = x^s 875 g.Exp(&n, &s, p) // g = n^s 876 r := e 877 for { 878 // find the least m such that ord_p(b) = 2^m 879 var m uint 880 t.Set(&b) 881 for t.Cmp(intOne) != 0 { 882 t.Mul(&t, &t).Mod(&t, p) 883 m++ 884 } 885 886 if m == 0 { 887 return z.Set(&y) 888 } 889 890 t.SetInt64(0).SetBit(&t, int(r-m-1), 1).Exp(&g, &t, p) 891 // t = g^(2^(r-m-1)) mod p 892 g.Mul(&t, &t).Mod(&g, p) // g = g^(2^(r-m)) mod p 893 y.Mul(&y, &t).Mod(&y, p) 894 b.Mul(&b, &g).Mod(&b, p) 895 r = m 896 } 897 } 898 899 // ModSqrt sets z to a square root of x mod p if such a square root exists, and 900 // returns z. The modulus p must be an odd prime. If x is not a square mod p, 901 // ModSqrt leaves z unchanged and returns nil. This function panics if p is 902 // not an odd integer. 903 func (z *Int) ModSqrt(x, p *Int) *Int { 904 switch Jacobi(x, p) { 905 case -1: 906 return nil // x is not a square mod p 907 case 0: 908 return z.SetInt64(0) // sqrt(0) mod p = 0 909 case 1: 910 break 911 } 912 if x.neg || x.Cmp(p) >= 0 { // ensure 0 <= x < p 913 x = new(Int).Mod(x, p) 914 } 915 916 switch { 917 case p.abs[0]%4 == 3: 918 // Check whether p is 3 mod 4, and if so, use the faster algorithm. 919 return z.modSqrt3Mod4Prime(x, p) 920 case p.abs[0]%8 == 5: 921 // Check whether p is 5 mod 8, use Atkin's algorithm. 922 return z.modSqrt5Mod8Prime(x, p) 923 default: 924 // Otherwise, use Tonelli-Shanks. 925 return z.modSqrtTonelliShanks(x, p) 926 } 927 } 928 929 // Lsh sets z = x << n and returns z. 930 func (z *Int) Lsh(x *Int, n uint) *Int { 931 z.abs = z.abs.shl(x.abs, n) 932 z.neg = x.neg 933 return z 934 } 935 936 // Rsh sets z = x >> n and returns z. 937 func (z *Int) Rsh(x *Int, n uint) *Int { 938 if x.neg { 939 // (-x) >> s == ^(x-1) >> s == ^((x-1) >> s) == -(((x-1) >> s) + 1) 940 t := z.abs.sub(x.abs, natOne) // no underflow because |x| > 0 941 t = t.shr(t, n) 942 z.abs = t.add(t, natOne) 943 z.neg = true // z cannot be zero if x is negative 944 return z 945 } 946 947 z.abs = z.abs.shr(x.abs, n) 948 z.neg = false 949 return z 950 } 951 952 // Bit returns the value of the i'th bit of x. That is, it 953 // returns (x>>i)&1. The bit index i must be >= 0. 954 func (x *Int) Bit(i int) uint { 955 if i == 0 { 956 // optimization for common case: odd/even test of x 957 if len(x.abs) > 0 { 958 return uint(x.abs[0] & 1) // bit 0 is same for -x 959 } 960 return 0 961 } 962 if i < 0 { 963 panic("negative bit index") 964 } 965 if x.neg { 966 t := nat(nil).sub(x.abs, natOne) 967 return t.bit(uint(i)) ^ 1 968 } 969 970 return x.abs.bit(uint(i)) 971 } 972 973 // SetBit sets z to x, with x's i'th bit set to b (0 or 1). 974 // That is, if b is 1 SetBit sets z = x | (1 << i); 975 // if b is 0 SetBit sets z = x &^ (1 << i). If b is not 0 or 1, 976 // SetBit will panic. 977 func (z *Int) SetBit(x *Int, i int, b uint) *Int { 978 if i < 0 { 979 panic("negative bit index") 980 } 981 if x.neg { 982 t := z.abs.sub(x.abs, natOne) 983 t = t.setBit(t, uint(i), b^1) 984 z.abs = t.add(t, natOne) 985 z.neg = len(z.abs) > 0 986 return z 987 } 988 z.abs = z.abs.setBit(x.abs, uint(i), b) 989 z.neg = false 990 return z 991 } 992 993 // And sets z = x & y and returns z. 994 func (z *Int) And(x, y *Int) *Int { 995 if x.neg == y.neg { 996 if x.neg { 997 // (-x) & (-y) == ^(x-1) & ^(y-1) == ^((x-1) | (y-1)) == -(((x-1) | (y-1)) + 1) 998 x1 := nat(nil).sub(x.abs, natOne) 999 y1 := nat(nil).sub(y.abs, natOne) 1000 z.abs = z.abs.add(z.abs.or(x1, y1), natOne) 1001 z.neg = true // z cannot be zero if x and y are negative 1002 return z 1003 } 1004 1005 // x & y == x & y 1006 z.abs = z.abs.and(x.abs, y.abs) 1007 z.neg = false 1008 return z 1009 } 1010 1011 // x.neg != y.neg 1012 if x.neg { 1013 x, y = y, x // & is symmetric 1014 } 1015 1016 // x & (-y) == x & ^(y-1) == x &^ (y-1) 1017 y1 := nat(nil).sub(y.abs, natOne) 1018 z.abs = z.abs.andNot(x.abs, y1) 1019 z.neg = false 1020 return z 1021 } 1022 1023 // AndNot sets z = x &^ y and returns z. 1024 func (z *Int) AndNot(x, y *Int) *Int { 1025 if x.neg == y.neg { 1026 if x.neg { 1027 // (-x) &^ (-y) == ^(x-1) &^ ^(y-1) == ^(x-1) & (y-1) == (y-1) &^ (x-1) 1028 x1 := nat(nil).sub(x.abs, natOne) 1029 y1 := nat(nil).sub(y.abs, natOne) 1030 z.abs = z.abs.andNot(y1, x1) 1031 z.neg = false 1032 return z 1033 } 1034 1035 // x &^ y == x &^ y 1036 z.abs = z.abs.andNot(x.abs, y.abs) 1037 z.neg = false 1038 return z 1039 } 1040 1041 if x.neg { 1042 // (-x) &^ y == ^(x-1) &^ y == ^(x-1) & ^y == ^((x-1) | y) == -(((x-1) | y) + 1) 1043 x1 := nat(nil).sub(x.abs, natOne) 1044 z.abs = z.abs.add(z.abs.or(x1, y.abs), natOne) 1045 z.neg = true // z cannot be zero if x is negative and y is positive 1046 return z 1047 } 1048 1049 // x &^ (-y) == x &^ ^(y-1) == x & (y-1) 1050 y1 := nat(nil).sub(y.abs, natOne) 1051 z.abs = z.abs.and(x.abs, y1) 1052 z.neg = false 1053 return z 1054 } 1055 1056 // Or sets z = x | y and returns z. 1057 func (z *Int) Or(x, y *Int) *Int { 1058 if x.neg == y.neg { 1059 if x.neg { 1060 // (-x) | (-y) == ^(x-1) | ^(y-1) == ^((x-1) & (y-1)) == -(((x-1) & (y-1)) + 1) 1061 x1 := nat(nil).sub(x.abs, natOne) 1062 y1 := nat(nil).sub(y.abs, natOne) 1063 z.abs = z.abs.add(z.abs.and(x1, y1), natOne) 1064 z.neg = true // z cannot be zero if x and y are negative 1065 return z 1066 } 1067 1068 // x | y == x | y 1069 z.abs = z.abs.or(x.abs, y.abs) 1070 z.neg = false 1071 return z 1072 } 1073 1074 // x.neg != y.neg 1075 if x.neg { 1076 x, y = y, x // | is symmetric 1077 } 1078 1079 // x | (-y) == x | ^(y-1) == ^((y-1) &^ x) == -(^((y-1) &^ x) + 1) 1080 y1 := nat(nil).sub(y.abs, natOne) 1081 z.abs = z.abs.add(z.abs.andNot(y1, x.abs), natOne) 1082 z.neg = true // z cannot be zero if one of x or y is negative 1083 return z 1084 } 1085 1086 // Xor sets z = x ^ y and returns z. 1087 func (z *Int) Xor(x, y *Int) *Int { 1088 if x.neg == y.neg { 1089 if x.neg { 1090 // (-x) ^ (-y) == ^(x-1) ^ ^(y-1) == (x-1) ^ (y-1) 1091 x1 := nat(nil).sub(x.abs, natOne) 1092 y1 := nat(nil).sub(y.abs, natOne) 1093 z.abs = z.abs.xor(x1, y1) 1094 z.neg = false 1095 return z 1096 } 1097 1098 // x ^ y == x ^ y 1099 z.abs = z.abs.xor(x.abs, y.abs) 1100 z.neg = false 1101 return z 1102 } 1103 1104 // x.neg != y.neg 1105 if x.neg { 1106 x, y = y, x // ^ is symmetric 1107 } 1108 1109 // x ^ (-y) == x ^ ^(y-1) == ^(x ^ (y-1)) == -((x ^ (y-1)) + 1) 1110 y1 := nat(nil).sub(y.abs, natOne) 1111 z.abs = z.abs.add(z.abs.xor(x.abs, y1), natOne) 1112 z.neg = true // z cannot be zero if only one of x or y is negative 1113 return z 1114 } 1115 1116 // Not sets z = ^x and returns z. 1117 func (z *Int) Not(x *Int) *Int { 1118 if x.neg { 1119 // ^(-x) == ^(^(x-1)) == x-1 1120 z.abs = z.abs.sub(x.abs, natOne) 1121 z.neg = false 1122 return z 1123 } 1124 1125 // ^x == -x-1 == -(x+1) 1126 z.abs = z.abs.add(x.abs, natOne) 1127 z.neg = true // z cannot be zero if x is positive 1128 return z 1129 } 1130 1131 // Sqrt sets z to ⌊√x⌋, the largest integer such that z² ≤ x, and returns z. 1132 // It panics if x is negative. 1133 func (z *Int) Sqrt(x *Int) *Int { 1134 if x.neg { 1135 panic("square root of negative number") 1136 } 1137 z.neg = false 1138 z.abs = z.abs.sqrt(x.abs) 1139 return z 1140 }