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