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