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