github.com/c12o16h1/go/src@v0.0.0-20200114212001-5a151c0f00ed/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 == y: 327 // nothing to do 328 case x.neg == y.neg: 329 r = x.abs.cmp(y.abs) 330 if x.neg { 331 r = -r 332 } 333 case x.neg: 334 r = -1 335 default: 336 r = 1 337 } 338 return 339 } 340 341 // CmpAbs compares the absolute values of x and y and returns: 342 // 343 // -1 if |x| < |y| 344 // 0 if |x| == |y| 345 // +1 if |x| > |y| 346 // 347 func (x *Int) CmpAbs(y *Int) int { 348 return x.abs.cmp(y.abs) 349 } 350 351 // low32 returns the least significant 32 bits of x. 352 func low32(x nat) uint32 { 353 if len(x) == 0 { 354 return 0 355 } 356 return uint32(x[0]) 357 } 358 359 // low64 returns the least significant 64 bits of x. 360 func low64(x nat) uint64 { 361 if len(x) == 0 { 362 return 0 363 } 364 v := uint64(x[0]) 365 if _W == 32 && len(x) > 1 { 366 return uint64(x[1])<<32 | v 367 } 368 return v 369 } 370 371 // Int64 returns the int64 representation of x. 372 // If x cannot be represented in an int64, the result is undefined. 373 func (x *Int) Int64() int64 { 374 v := int64(low64(x.abs)) 375 if x.neg { 376 v = -v 377 } 378 return v 379 } 380 381 // Uint64 returns the uint64 representation of x. 382 // If x cannot be represented in a uint64, the result is undefined. 383 func (x *Int) Uint64() uint64 { 384 return low64(x.abs) 385 } 386 387 // IsInt64 reports whether x can be represented as an int64. 388 func (x *Int) IsInt64() bool { 389 if len(x.abs) <= 64/_W { 390 w := int64(low64(x.abs)) 391 return w >= 0 || x.neg && w == -w 392 } 393 return false 394 } 395 396 // IsUint64 reports whether x can be represented as a uint64. 397 func (x *Int) IsUint64() bool { 398 return !x.neg && len(x.abs) <= 64/_W 399 } 400 401 // SetString sets z to the value of s, interpreted in the given base, 402 // and returns z and a boolean indicating success. The entire string 403 // (not just a prefix) must be valid for success. If SetString fails, 404 // the value of z is undefined but the returned value is nil. 405 // 406 // The base argument must be 0 or a value between 2 and MaxBase. 407 // For base 0, the number prefix determines the actual base: A prefix of 408 // ``0b'' or ``0B'' selects base 2, ``0'', ``0o'' or ``0O'' selects base 8, 409 // and ``0x'' or ``0X'' selects base 16. Otherwise, the selected base is 10 410 // and no prefix is accepted. 411 // 412 // For bases <= 36, lower and upper case letters are considered the same: 413 // The letters 'a' to 'z' and 'A' to 'Z' represent digit values 10 to 35. 414 // For bases > 36, the upper case letters 'A' to 'Z' represent the digit 415 // values 36 to 61. 416 // 417 // For base 0, an underscore character ``_'' may appear between a base 418 // prefix and an adjacent digit, and between successive digits; such 419 // underscores do not change the value of the number. 420 // Incorrect placement of underscores is reported as an error if there 421 // are no other errors. If base != 0, underscores are not recognized 422 // and act like any other character that is not a valid digit. 423 // 424 func (z *Int) SetString(s string, base int) (*Int, bool) { 425 return z.setFromScanner(strings.NewReader(s), base) 426 } 427 428 // setFromScanner implements SetString given an io.BytesScanner. 429 // For documentation see comments of SetString. 430 func (z *Int) setFromScanner(r io.ByteScanner, base int) (*Int, bool) { 431 if _, _, err := z.scan(r, base); err != nil { 432 return nil, false 433 } 434 // entire content must have been consumed 435 if _, err := r.ReadByte(); err != io.EOF { 436 return nil, false 437 } 438 return z, true // err == io.EOF => scan consumed all content of r 439 } 440 441 // SetBytes interprets buf as the bytes of a big-endian unsigned 442 // integer, sets z to that value, and returns z. 443 func (z *Int) SetBytes(buf []byte) *Int { 444 z.abs = z.abs.setBytes(buf) 445 z.neg = false 446 return z 447 } 448 449 // Bytes returns the absolute value of x as a big-endian byte slice. 450 func (x *Int) Bytes() []byte { 451 buf := make([]byte, len(x.abs)*_S) 452 return buf[x.abs.bytes(buf):] 453 } 454 455 // BitLen returns the length of the absolute value of x in bits. 456 // The bit length of 0 is 0. 457 func (x *Int) BitLen() int { 458 return x.abs.bitLen() 459 } 460 461 // TrailingZeroBits returns the number of consecutive least significant zero 462 // bits of |x|. 463 func (x *Int) TrailingZeroBits() uint { 464 return x.abs.trailingZeroBits() 465 } 466 467 // Exp sets z = x**y mod |m| (i.e. the sign of m is ignored), and returns z. 468 // If m == nil or m == 0, z = x**y unless y <= 0 then z = 1. If m > 0, y < 0, 469 // and x and n are not relatively prime, z is unchanged and nil is returned. 470 // 471 // Modular exponentiation of inputs of a particular size is not a 472 // cryptographically constant-time operation. 473 func (z *Int) Exp(x, y, m *Int) *Int { 474 // See Knuth, volume 2, section 4.6.3. 475 xWords := x.abs 476 if y.neg { 477 if m == nil || len(m.abs) == 0 { 478 return z.SetInt64(1) 479 } 480 // for y < 0: x**y mod m == (x**(-1))**|y| mod m 481 inverse := new(Int).ModInverse(x, m) 482 if inverse == nil { 483 return nil 484 } 485 xWords = inverse.abs 486 } 487 yWords := y.abs 488 489 var mWords nat 490 if m != nil { 491 mWords = m.abs // m.abs may be nil for m == 0 492 } 493 494 z.abs = z.abs.expNN(xWords, yWords, mWords) 495 z.neg = len(z.abs) > 0 && x.neg && len(yWords) > 0 && yWords[0]&1 == 1 // 0 has no sign 496 if z.neg && len(mWords) > 0 { 497 // make modulus result positive 498 z.abs = z.abs.sub(mWords, z.abs) // z == x**y mod |m| && 0 <= z < |m| 499 z.neg = false 500 } 501 502 return z 503 } 504 505 // GCD sets z to the greatest common divisor of a and b and returns z. 506 // If x or y are not nil, GCD sets their value such that z = a*x + b*y. 507 // Regardless of the signs of a and b, z is always >= 0. 508 // If a == b == 0, GCD sets z = x = y = 0. 509 // If a == 0 and b != 0, GCD sets z = |b|, x = 0, y = sign(b) * 1. 510 // If a != 0 and b == 0, GCD sets z = |a|, x = sign(a) * 1, y = 0. 511 func (z *Int) GCD(x, y, a, b *Int) *Int { 512 if len(a.abs) == 0 || len(b.abs) == 0 { 513 lenA, lenB, negA, negB := len(a.abs), len(b.abs), a.neg, b.neg 514 if lenA == 0 { 515 z.Set(b) 516 } else { 517 z.Set(a) 518 } 519 z.neg = false 520 if x != nil { 521 if lenA == 0 { 522 x.SetUint64(0) 523 } else { 524 x.SetUint64(1) 525 x.neg = negA 526 } 527 } 528 if y != nil { 529 if lenB == 0 { 530 y.SetUint64(0) 531 } else { 532 y.SetUint64(1) 533 y.neg = negB 534 } 535 } 536 return z 537 } 538 539 return z.lehmerGCD(x, y, a, b) 540 } 541 542 // lehmerSimulate attempts to simulate several Euclidean update steps 543 // using the leading digits of A and B. It returns u0, u1, v0, v1 544 // such that A and B can be updated as: 545 // A = u0*A + v0*B 546 // B = u1*A + v1*B 547 // Requirements: A >= B and len(B.abs) >= 2 548 // Since we are calculating with full words to avoid overflow, 549 // we use 'even' to track the sign of the cosequences. 550 // For even iterations: u0, v1 >= 0 && u1, v0 <= 0 551 // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0 552 func lehmerSimulate(A, B *Int) (u0, u1, v0, v1 Word, even bool) { 553 // initialize the digits 554 var a1, a2, u2, v2 Word 555 556 m := len(B.abs) // m >= 2 557 n := len(A.abs) // n >= m >= 2 558 559 // extract the top Word of bits from A and B 560 h := nlz(A.abs[n-1]) 561 a1 = A.abs[n-1]<<h | A.abs[n-2]>>(_W-h) 562 // B may have implicit zero words in the high bits if the lengths differ 563 switch { 564 case n == m: 565 a2 = B.abs[n-1]<<h | B.abs[n-2]>>(_W-h) 566 case n == m+1: 567 a2 = B.abs[n-2] >> (_W - h) 568 default: 569 a2 = 0 570 } 571 572 // Since we are calculating with full words to avoid overflow, 573 // we use 'even' to track the sign of the cosequences. 574 // For even iterations: u0, v1 >= 0 && u1, v0 <= 0 575 // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0 576 // The first iteration starts with k=1 (odd). 577 even = false 578 // variables to track the cosequences 579 u0, u1, u2 = 0, 1, 0 580 v0, v1, v2 = 0, 0, 1 581 582 // Calculate the quotient and cosequences using Collins' stopping condition. 583 // Note that overflow of a Word is not possible when computing the remainder 584 // sequence and cosequences since the cosequence size is bounded by the input size. 585 // See section 4.2 of Jebelean for details. 586 for a2 >= v2 && a1-a2 >= v1+v2 { 587 q, r := a1/a2, a1%a2 588 a1, a2 = a2, r 589 u0, u1, u2 = u1, u2, u1+q*u2 590 v0, v1, v2 = v1, v2, v1+q*v2 591 even = !even 592 } 593 return 594 } 595 596 // lehmerUpdate updates the inputs A and B such that: 597 // A = u0*A + v0*B 598 // B = u1*A + v1*B 599 // where the signs of u0, u1, v0, v1 are given by even 600 // For even == true: u0, v1 >= 0 && u1, v0 <= 0 601 // For even == false: u0, v1 <= 0 && u1, v0 >= 0 602 // q, r, s, t are temporary variables to avoid allocations in the multiplication 603 func lehmerUpdate(A, B, q, r, s, t *Int, u0, u1, v0, v1 Word, even bool) { 604 605 t.abs = t.abs.setWord(u0) 606 s.abs = s.abs.setWord(v0) 607 t.neg = !even 608 s.neg = even 609 610 t.Mul(A, t) 611 s.Mul(B, s) 612 613 r.abs = r.abs.setWord(u1) 614 q.abs = q.abs.setWord(v1) 615 r.neg = even 616 q.neg = !even 617 618 r.Mul(A, r) 619 q.Mul(B, q) 620 621 A.Add(t, s) 622 B.Add(r, q) 623 } 624 625 // euclidUpdate performs a single step of the Euclidean GCD algorithm 626 // if extended is true, it also updates the cosequence Ua, Ub 627 func euclidUpdate(A, B, Ua, Ub, q, r, s, t *Int, extended bool) { 628 q, r = q.QuoRem(A, B, r) 629 630 *A, *B, *r = *B, *r, *A 631 632 if extended { 633 // Ua, Ub = Ub, Ua - q*Ub 634 t.Set(Ub) 635 s.Mul(Ub, q) 636 Ub.Sub(Ua, s) 637 Ua.Set(t) 638 } 639 } 640 641 // lehmerGCD sets z to the greatest common divisor of a and b, 642 // which both must be != 0, and returns z. 643 // If x or y are not nil, their values are set such that z = a*x + b*y. 644 // See Knuth, The Art of Computer Programming, Vol. 2, Section 4.5.2, Algorithm L. 645 // This implementation uses the improved condition by Collins requiring only one 646 // quotient and avoiding the possibility of single Word overflow. 647 // See Jebelean, "Improving the multiprecision Euclidean algorithm", 648 // Design and Implementation of Symbolic Computation Systems, pp 45-58. 649 // The cosequences are updated according to Algorithm 10.45 from 650 // Cohen et al. "Handbook of Elliptic and Hyperelliptic Curve Cryptography" pp 192. 651 func (z *Int) lehmerGCD(x, y, a, b *Int) *Int { 652 var A, B, Ua, Ub *Int 653 654 A = new(Int).Abs(a) 655 B = new(Int).Abs(b) 656 657 extended := x != nil || y != nil 658 659 if extended { 660 // Ua (Ub) tracks how many times input a has been accumulated into A (B). 661 Ua = new(Int).SetInt64(1) 662 Ub = new(Int) 663 } 664 665 // temp variables for multiprecision update 666 q := new(Int) 667 r := new(Int) 668 s := new(Int) 669 t := new(Int) 670 671 // ensure A >= B 672 if A.abs.cmp(B.abs) < 0 { 673 A, B = B, A 674 Ub, Ua = Ua, Ub 675 } 676 677 // loop invariant A >= B 678 for len(B.abs) > 1 { 679 // Attempt to calculate in single-precision using leading words of A and B. 680 u0, u1, v0, v1, even := lehmerSimulate(A, B) 681 682 // multiprecision Step 683 if v0 != 0 { 684 // Simulate the effect of the single-precision steps using the cosequences. 685 // A = u0*A + v0*B 686 // B = u1*A + v1*B 687 lehmerUpdate(A, B, q, r, s, t, u0, u1, v0, v1, even) 688 689 if extended { 690 // Ua = u0*Ua + v0*Ub 691 // Ub = u1*Ua + v1*Ub 692 lehmerUpdate(Ua, Ub, q, r, s, t, u0, u1, v0, v1, even) 693 } 694 695 } else { 696 // Single-digit calculations failed to simulate any quotients. 697 // Do a standard Euclidean step. 698 euclidUpdate(A, B, Ua, Ub, q, r, s, t, extended) 699 } 700 } 701 702 if len(B.abs) > 0 { 703 // extended Euclidean algorithm base case if B is a single Word 704 if len(A.abs) > 1 { 705 // A is longer than a single Word, so one update is needed. 706 euclidUpdate(A, B, Ua, Ub, q, r, s, t, extended) 707 } 708 if len(B.abs) > 0 { 709 // A and B are both a single Word. 710 aWord, bWord := A.abs[0], B.abs[0] 711 if extended { 712 var ua, ub, va, vb Word 713 ua, ub = 1, 0 714 va, vb = 0, 1 715 even := true 716 for bWord != 0 { 717 q, r := aWord/bWord, aWord%bWord 718 aWord, bWord = bWord, r 719 ua, ub = ub, ua+q*ub 720 va, vb = vb, va+q*vb 721 even = !even 722 } 723 724 t.abs = t.abs.setWord(ua) 725 s.abs = s.abs.setWord(va) 726 t.neg = !even 727 s.neg = even 728 729 t.Mul(Ua, t) 730 s.Mul(Ub, s) 731 732 Ua.Add(t, s) 733 } else { 734 for bWord != 0 { 735 aWord, bWord = bWord, aWord%bWord 736 } 737 } 738 A.abs[0] = aWord 739 } 740 } 741 negA := a.neg 742 if y != nil { 743 // avoid aliasing b needed in the division below 744 if y == b { 745 B.Set(b) 746 } else { 747 B = b 748 } 749 // y = (z - a*x)/b 750 y.Mul(a, Ua) // y can safely alias a 751 if negA { 752 y.neg = !y.neg 753 } 754 y.Sub(A, y) 755 y.Div(y, B) 756 } 757 758 if x != nil { 759 *x = *Ua 760 if negA { 761 x.neg = !x.neg 762 } 763 } 764 765 *z = *A 766 767 return z 768 } 769 770 // Rand sets z to a pseudo-random number in [0, n) and returns z. 771 // 772 // As this uses the math/rand package, it must not be used for 773 // security-sensitive work. Use crypto/rand.Int instead. 774 func (z *Int) Rand(rnd *rand.Rand, n *Int) *Int { 775 z.neg = false 776 if n.neg || len(n.abs) == 0 { 777 z.abs = nil 778 return z 779 } 780 z.abs = z.abs.random(rnd, n.abs, n.abs.bitLen()) 781 return z 782 } 783 784 // ModInverse sets z to the multiplicative inverse of g in the ring ℤ/nℤ 785 // and returns z. If g and n are not relatively prime, g has no multiplicative 786 // inverse in the ring ℤ/nℤ. In this case, z is unchanged and the return value 787 // is nil. 788 func (z *Int) ModInverse(g, n *Int) *Int { 789 // GCD expects parameters a and b to be > 0. 790 if n.neg { 791 var n2 Int 792 n = n2.Neg(n) 793 } 794 if g.neg { 795 var g2 Int 796 g = g2.Mod(g, n) 797 } 798 var d, x Int 799 d.GCD(&x, nil, g, n) 800 801 // if and only if d==1, g and n are relatively prime 802 if d.Cmp(intOne) != 0 { 803 return nil 804 } 805 806 // x and y are such that g*x + n*y = 1, therefore x is the inverse element, 807 // but it may be negative, so convert to the range 0 <= z < |n| 808 if x.neg { 809 z.Add(&x, n) 810 } else { 811 z.Set(&x) 812 } 813 return z 814 } 815 816 // Jacobi returns the Jacobi symbol (x/y), either +1, -1, or 0. 817 // The y argument must be an odd integer. 818 func Jacobi(x, y *Int) int { 819 if len(y.abs) == 0 || y.abs[0]&1 == 0 { 820 panic(fmt.Sprintf("big: invalid 2nd argument to Int.Jacobi: need odd integer but got %s", y)) 821 } 822 823 // We use the formulation described in chapter 2, section 2.4, 824 // "The Yacas Book of Algorithms": 825 // http://yacas.sourceforge.net/Algo.book.pdf 826 827 var a, b, c Int 828 a.Set(x) 829 b.Set(y) 830 j := 1 831 832 if b.neg { 833 if a.neg { 834 j = -1 835 } 836 b.neg = false 837 } 838 839 for { 840 if b.Cmp(intOne) == 0 { 841 return j 842 } 843 if len(a.abs) == 0 { 844 return 0 845 } 846 a.Mod(&a, &b) 847 if len(a.abs) == 0 { 848 return 0 849 } 850 // a > 0 851 852 // handle factors of 2 in 'a' 853 s := a.abs.trailingZeroBits() 854 if s&1 != 0 { 855 bmod8 := b.abs[0] & 7 856 if bmod8 == 3 || bmod8 == 5 { 857 j = -j 858 } 859 } 860 c.Rsh(&a, s) // a = 2^s*c 861 862 // swap numerator and denominator 863 if b.abs[0]&3 == 3 && c.abs[0]&3 == 3 { 864 j = -j 865 } 866 a.Set(&b) 867 b.Set(&c) 868 } 869 } 870 871 // modSqrt3Mod4 uses the identity 872 // (a^((p+1)/4))^2 mod p 873 // == u^(p+1) mod p 874 // == u^2 mod p 875 // to calculate the square root of any quadratic residue mod p quickly for 3 876 // mod 4 primes. 877 func (z *Int) modSqrt3Mod4Prime(x, p *Int) *Int { 878 e := new(Int).Add(p, intOne) // e = p + 1 879 e.Rsh(e, 2) // e = (p + 1) / 4 880 z.Exp(x, e, p) // z = x^e mod p 881 return z 882 } 883 884 // modSqrt5Mod8 uses Atkin's observation that 2 is not a square mod p 885 // alpha == (2*a)^((p-5)/8) mod p 886 // beta == 2*a*alpha^2 mod p is a square root of -1 887 // b == a*alpha*(beta-1) mod p is a square root of a 888 // to calculate the square root of any quadratic residue mod p quickly for 5 889 // mod 8 primes. 890 func (z *Int) modSqrt5Mod8Prime(x, p *Int) *Int { 891 // p == 5 mod 8 implies p = e*8 + 5 892 // e is the quotient and 5 the remainder on division by 8 893 e := new(Int).Rsh(p, 3) // e = (p - 5) / 8 894 tx := new(Int).Lsh(x, 1) // tx = 2*x 895 alpha := new(Int).Exp(tx, e, p) 896 beta := new(Int).Mul(alpha, alpha) 897 beta.Mod(beta, p) 898 beta.Mul(beta, tx) 899 beta.Mod(beta, p) 900 beta.Sub(beta, intOne) 901 beta.Mul(beta, x) 902 beta.Mod(beta, p) 903 beta.Mul(beta, alpha) 904 z.Mod(beta, p) 905 return z 906 } 907 908 // modSqrtTonelliShanks uses the Tonelli-Shanks algorithm to find the square 909 // root of a quadratic residue modulo any prime. 910 func (z *Int) modSqrtTonelliShanks(x, p *Int) *Int { 911 // Break p-1 into s*2^e such that s is odd. 912 var s Int 913 s.Sub(p, intOne) 914 e := s.abs.trailingZeroBits() 915 s.Rsh(&s, e) 916 917 // find some non-square n 918 var n Int 919 n.SetInt64(2) 920 for Jacobi(&n, p) != -1 { 921 n.Add(&n, intOne) 922 } 923 924 // Core of the Tonelli-Shanks algorithm. Follows the description in 925 // section 6 of "Square roots from 1; 24, 51, 10 to Dan Shanks" by Ezra 926 // Brown: 927 // https://www.maa.org/sites/default/files/pdf/upload_library/22/Polya/07468342.di020786.02p0470a.pdf 928 var y, b, g, t Int 929 y.Add(&s, intOne) 930 y.Rsh(&y, 1) 931 y.Exp(x, &y, p) // y = x^((s+1)/2) 932 b.Exp(x, &s, p) // b = x^s 933 g.Exp(&n, &s, p) // g = n^s 934 r := e 935 for { 936 // find the least m such that ord_p(b) = 2^m 937 var m uint 938 t.Set(&b) 939 for t.Cmp(intOne) != 0 { 940 t.Mul(&t, &t).Mod(&t, p) 941 m++ 942 } 943 944 if m == 0 { 945 return z.Set(&y) 946 } 947 948 t.SetInt64(0).SetBit(&t, int(r-m-1), 1).Exp(&g, &t, p) 949 // t = g^(2^(r-m-1)) mod p 950 g.Mul(&t, &t).Mod(&g, p) // g = g^(2^(r-m)) mod p 951 y.Mul(&y, &t).Mod(&y, p) 952 b.Mul(&b, &g).Mod(&b, p) 953 r = m 954 } 955 } 956 957 // ModSqrt sets z to a square root of x mod p if such a square root exists, and 958 // returns z. The modulus p must be an odd prime. If x is not a square mod p, 959 // ModSqrt leaves z unchanged and returns nil. This function panics if p is 960 // not an odd integer. 961 func (z *Int) ModSqrt(x, p *Int) *Int { 962 switch Jacobi(x, p) { 963 case -1: 964 return nil // x is not a square mod p 965 case 0: 966 return z.SetInt64(0) // sqrt(0) mod p = 0 967 case 1: 968 break 969 } 970 if x.neg || x.Cmp(p) >= 0 { // ensure 0 <= x < p 971 x = new(Int).Mod(x, p) 972 } 973 974 switch { 975 case p.abs[0]%4 == 3: 976 // Check whether p is 3 mod 4, and if so, use the faster algorithm. 977 return z.modSqrt3Mod4Prime(x, p) 978 case p.abs[0]%8 == 5: 979 // Check whether p is 5 mod 8, use Atkin's algorithm. 980 return z.modSqrt5Mod8Prime(x, p) 981 default: 982 // Otherwise, use Tonelli-Shanks. 983 return z.modSqrtTonelliShanks(x, p) 984 } 985 } 986 987 // Lsh sets z = x << n and returns z. 988 func (z *Int) Lsh(x *Int, n uint) *Int { 989 z.abs = z.abs.shl(x.abs, n) 990 z.neg = x.neg 991 return z 992 } 993 994 // Rsh sets z = x >> n and returns z. 995 func (z *Int) Rsh(x *Int, n uint) *Int { 996 if x.neg { 997 // (-x) >> s == ^(x-1) >> s == ^((x-1) >> s) == -(((x-1) >> s) + 1) 998 t := z.abs.sub(x.abs, natOne) // no underflow because |x| > 0 999 t = t.shr(t, n) 1000 z.abs = t.add(t, natOne) 1001 z.neg = true // z cannot be zero if x is negative 1002 return z 1003 } 1004 1005 z.abs = z.abs.shr(x.abs, n) 1006 z.neg = false 1007 return z 1008 } 1009 1010 // Bit returns the value of the i'th bit of x. That is, it 1011 // returns (x>>i)&1. The bit index i must be >= 0. 1012 func (x *Int) Bit(i int) uint { 1013 if i == 0 { 1014 // optimization for common case: odd/even test of x 1015 if len(x.abs) > 0 { 1016 return uint(x.abs[0] & 1) // bit 0 is same for -x 1017 } 1018 return 0 1019 } 1020 if i < 0 { 1021 panic("negative bit index") 1022 } 1023 if x.neg { 1024 t := nat(nil).sub(x.abs, natOne) 1025 return t.bit(uint(i)) ^ 1 1026 } 1027 1028 return x.abs.bit(uint(i)) 1029 } 1030 1031 // SetBit sets z to x, with x's i'th bit set to b (0 or 1). 1032 // That is, if b is 1 SetBit sets z = x | (1 << i); 1033 // if b is 0 SetBit sets z = x &^ (1 << i). If b is not 0 or 1, 1034 // SetBit will panic. 1035 func (z *Int) SetBit(x *Int, i int, b uint) *Int { 1036 if i < 0 { 1037 panic("negative bit index") 1038 } 1039 if x.neg { 1040 t := z.abs.sub(x.abs, natOne) 1041 t = t.setBit(t, uint(i), b^1) 1042 z.abs = t.add(t, natOne) 1043 z.neg = len(z.abs) > 0 1044 return z 1045 } 1046 z.abs = z.abs.setBit(x.abs, uint(i), b) 1047 z.neg = false 1048 return z 1049 } 1050 1051 // And sets z = x & y and returns z. 1052 func (z *Int) And(x, y *Int) *Int { 1053 if x.neg == y.neg { 1054 if x.neg { 1055 // (-x) & (-y) == ^(x-1) & ^(y-1) == ^((x-1) | (y-1)) == -(((x-1) | (y-1)) + 1) 1056 x1 := nat(nil).sub(x.abs, natOne) 1057 y1 := nat(nil).sub(y.abs, natOne) 1058 z.abs = z.abs.add(z.abs.or(x1, y1), natOne) 1059 z.neg = true // z cannot be zero if x and y are negative 1060 return z 1061 } 1062 1063 // x & y == x & y 1064 z.abs = z.abs.and(x.abs, y.abs) 1065 z.neg = false 1066 return z 1067 } 1068 1069 // x.neg != y.neg 1070 if x.neg { 1071 x, y = y, x // & is symmetric 1072 } 1073 1074 // x & (-y) == x & ^(y-1) == x &^ (y-1) 1075 y1 := nat(nil).sub(y.abs, natOne) 1076 z.abs = z.abs.andNot(x.abs, y1) 1077 z.neg = false 1078 return z 1079 } 1080 1081 // AndNot sets z = x &^ y and returns z. 1082 func (z *Int) AndNot(x, y *Int) *Int { 1083 if x.neg == y.neg { 1084 if x.neg { 1085 // (-x) &^ (-y) == ^(x-1) &^ ^(y-1) == ^(x-1) & (y-1) == (y-1) &^ (x-1) 1086 x1 := nat(nil).sub(x.abs, natOne) 1087 y1 := nat(nil).sub(y.abs, natOne) 1088 z.abs = z.abs.andNot(y1, x1) 1089 z.neg = false 1090 return z 1091 } 1092 1093 // x &^ y == x &^ y 1094 z.abs = z.abs.andNot(x.abs, y.abs) 1095 z.neg = false 1096 return z 1097 } 1098 1099 if x.neg { 1100 // (-x) &^ y == ^(x-1) &^ y == ^(x-1) & ^y == ^((x-1) | y) == -(((x-1) | y) + 1) 1101 x1 := nat(nil).sub(x.abs, natOne) 1102 z.abs = z.abs.add(z.abs.or(x1, y.abs), natOne) 1103 z.neg = true // z cannot be zero if x is negative and y is positive 1104 return z 1105 } 1106 1107 // x &^ (-y) == x &^ ^(y-1) == x & (y-1) 1108 y1 := nat(nil).sub(y.abs, natOne) 1109 z.abs = z.abs.and(x.abs, y1) 1110 z.neg = false 1111 return z 1112 } 1113 1114 // Or sets z = x | y and returns z. 1115 func (z *Int) Or(x, y *Int) *Int { 1116 if x.neg == y.neg { 1117 if x.neg { 1118 // (-x) | (-y) == ^(x-1) | ^(y-1) == ^((x-1) & (y-1)) == -(((x-1) & (y-1)) + 1) 1119 x1 := nat(nil).sub(x.abs, natOne) 1120 y1 := nat(nil).sub(y.abs, natOne) 1121 z.abs = z.abs.add(z.abs.and(x1, y1), natOne) 1122 z.neg = true // z cannot be zero if x and y are negative 1123 return z 1124 } 1125 1126 // x | y == x | y 1127 z.abs = z.abs.or(x.abs, y.abs) 1128 z.neg = false 1129 return z 1130 } 1131 1132 // x.neg != y.neg 1133 if x.neg { 1134 x, y = y, x // | is symmetric 1135 } 1136 1137 // x | (-y) == x | ^(y-1) == ^((y-1) &^ x) == -(^((y-1) &^ x) + 1) 1138 y1 := nat(nil).sub(y.abs, natOne) 1139 z.abs = z.abs.add(z.abs.andNot(y1, x.abs), natOne) 1140 z.neg = true // z cannot be zero if one of x or y is negative 1141 return z 1142 } 1143 1144 // Xor sets z = x ^ y and returns z. 1145 func (z *Int) Xor(x, y *Int) *Int { 1146 if x.neg == y.neg { 1147 if x.neg { 1148 // (-x) ^ (-y) == ^(x-1) ^ ^(y-1) == (x-1) ^ (y-1) 1149 x1 := nat(nil).sub(x.abs, natOne) 1150 y1 := nat(nil).sub(y.abs, natOne) 1151 z.abs = z.abs.xor(x1, y1) 1152 z.neg = false 1153 return z 1154 } 1155 1156 // x ^ y == x ^ y 1157 z.abs = z.abs.xor(x.abs, y.abs) 1158 z.neg = false 1159 return z 1160 } 1161 1162 // x.neg != y.neg 1163 if x.neg { 1164 x, y = y, x // ^ is symmetric 1165 } 1166 1167 // x ^ (-y) == x ^ ^(y-1) == ^(x ^ (y-1)) == -((x ^ (y-1)) + 1) 1168 y1 := nat(nil).sub(y.abs, natOne) 1169 z.abs = z.abs.add(z.abs.xor(x.abs, y1), natOne) 1170 z.neg = true // z cannot be zero if only one of x or y is negative 1171 return z 1172 } 1173 1174 // Not sets z = ^x and returns z. 1175 func (z *Int) Not(x *Int) *Int { 1176 if x.neg { 1177 // ^(-x) == ^(^(x-1)) == x-1 1178 z.abs = z.abs.sub(x.abs, natOne) 1179 z.neg = false 1180 return z 1181 } 1182 1183 // ^x == -x-1 == -(x+1) 1184 z.abs = z.abs.add(x.abs, natOne) 1185 z.neg = true // z cannot be zero if x is positive 1186 return z 1187 } 1188 1189 // Sqrt sets z to ⌊√x⌋, the largest integer such that z² ≤ x, and returns z. 1190 // It panics if x is negative. 1191 func (z *Int) Sqrt(x *Int) *Int { 1192 if x.neg { 1193 panic("square root of negative number") 1194 } 1195 z.neg = false 1196 z.abs = z.abs.sqrt(x.abs) 1197 return z 1198 }