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