diff --git a/src/math/big/int.go b/src/math/big/int.go index 1e99a964665484..661813a2fd426f 100644 --- a/src/math/big/int.go +++ b/src/math/big/int.go @@ -14,1442 +14,4 @@ import ( "sync" ) -// An Int represents a signed multi-precision integer. -// The zero value for an Int represents the value 0. -// -// Operations always take pointer arguments (*Int) rather -// than Int values, and each unique Int value requires -// its own unique *Int pointer. To "copy" an Int value, -// an existing (or newly allocated) Int must be set to -// a new value using the [Int.Set] method; shallow copies -// of Ints are not supported and may lead to errors. -// -// Note that methods may leak the Int's value through timing side-channels. -// Because of this and because of the scope and complexity of the -// implementation, Int is not well-suited to implement cryptographic operations. -// The standard library avoids exposing non-trivial Int methods to -// attacker-controlled inputs and the determination of whether a bug in math/big -// is considered a security vulnerability might depend on the impact on the -// standard library. -type Int struct { - neg bool // sign - abs nat // absolute value of the integer -} - -var intOne = &Int{false, natOne} - -// Sign returns: -// - -1 if x < 0; -// - 0 if x == 0; -// - +1 if x > 0. -func (x *Int) Sign() int { - // This function is used in cryptographic operations. It must not leak - // anything but the Int's sign and bit size through side-channels. Any - // changes must be reviewed by a security expert. - if len(x.abs) == 0 { - return 0 - } - if x.neg { - return -1 - } - return 1 -} - -// SetInt64 sets z to x and returns z. -func (z *Int) SetInt64(x int64) *Int { - neg := false - if x < 0 { - neg = true - x = -x - } - z.abs = z.abs.setUint64(uint64(x)) - z.neg = neg - return z -} - -// SetUint64 sets z to x and returns z. -func (z *Int) SetUint64(x uint64) *Int { - z.abs = z.abs.setUint64(x) - z.neg = false - return z -} - -// NewInt allocates and returns a new [Int] set to x. -func NewInt(x int64) *Int { - // This code is arranged to be inlineable and produce - // zero allocations when inlined. See issue 29951. - u := uint64(x) - if x < 0 { - u = -u - } - var abs []Word - if x == 0 { - } else if _W == 32 && u>>32 != 0 { - abs = []Word{Word(u), Word(u >> 32)} - } else { - abs = []Word{Word(u)} - } - return &Int{neg: x < 0, abs: abs} -} - -// Set sets z to x and returns z. -func (z *Int) Set(x *Int) *Int { - if z != x { - z.abs = z.abs.set(x.abs) - z.neg = x.neg - } - return z -} - -// Bits provides raw (unchecked but fast) access to x by returning its -// absolute value as a little-endian [Word] slice. The result and x share -// the same underlying array. -// Bits is intended to support implementation of missing low-level [Int] -// functionality outside this package; it should be avoided otherwise. -func (x *Int) Bits() []Word { - // This function is used in cryptographic operations. It must not leak - // anything but the Int's sign and bit size through side-channels. Any - // changes must be reviewed by a security expert. - return x.abs -} - -// SetBits provides raw (unchecked but fast) access to z by setting its -// value to abs, interpreted as a little-endian [Word] slice, and returning -// z. The result and abs share the same underlying array. -// SetBits is intended to support implementation of missing low-level [Int] -// functionality outside this package; it should be avoided otherwise. -func (z *Int) SetBits(abs []Word) *Int { - z.abs = nat(abs).norm() - z.neg = false - return z -} - -// Abs sets z to |x| (the absolute value of x) and returns z. -func (z *Int) Abs(x *Int) *Int { - z.Set(x) - z.neg = false - return z -} - -// Neg sets z to -x and returns z. -func (z *Int) Neg(x *Int) *Int { - z.Set(x) - z.neg = len(z.abs) > 0 && !z.neg // 0 has no sign - return z -} - -// Add sets z to the sum x+y and returns z. -func (z *Int) Add(x, y *Int) *Int { - neg := x.neg - if x.neg == y.neg { - // x + y == x + y - // (-x) + (-y) == -(x + y) - z.abs = z.abs.add(x.abs, y.abs) - } else { - // x + (-y) == x - y == -(y - x) - // (-x) + y == y - x == -(x - y) - if x.abs.cmp(y.abs) >= 0 { - z.abs = z.abs.sub(x.abs, y.abs) - } else { - neg = !neg - z.abs = z.abs.sub(y.abs, x.abs) - } - } - z.neg = len(z.abs) > 0 && neg // 0 has no sign - return z -} - -// Sub sets z to the difference x-y and returns z. -func (z *Int) Sub(x, y *Int) *Int { - neg := x.neg - if x.neg != y.neg { - // x - (-y) == x + y - // (-x) - y == -(x + y) - z.abs = z.abs.add(x.abs, y.abs) - } else { - // x - y == x - y == -(y - x) - // (-x) - (-y) == y - x == -(x - y) - if x.abs.cmp(y.abs) >= 0 { - z.abs = z.abs.sub(x.abs, y.abs) - } else { - neg = !neg - z.abs = z.abs.sub(y.abs, x.abs) - } - } - z.neg = len(z.abs) > 0 && neg // 0 has no sign - return z -} - -// Mul sets z to the product x*y and returns z. -func (z *Int) Mul(x, y *Int) *Int { - z.mul(nil, x, y) - return z -} - -// mul is like Mul but takes an explicit stack to use, for internal use. -// It does not return a *Int because doing so makes the stack-allocated Ints -// used in natmul.go escape to the heap (even though the result is unused). -func (z *Int) mul(stk *stack, x, y *Int) { - // x * y == x * y - // x * (-y) == -(x * y) - // (-x) * y == -(x * y) - // (-x) * (-y) == x * y - if x == y { - z.abs = z.abs.sqr(stk, x.abs) - z.neg = false - return - } - z.abs = z.abs.mul(stk, x.abs, y.abs) - z.neg = len(z.abs) > 0 && x.neg != y.neg // 0 has no sign -} - -// MulRange sets z to the product of all integers -// in the range [a, b] inclusively and returns z. -// If a > b (empty range), the result is 1. -func (z *Int) MulRange(a, b int64) *Int { - switch { - case a > b: - return z.SetInt64(1) // empty range - case a <= 0 && b >= 0: - return z.SetInt64(0) // range includes 0 - } - // a <= b && (b < 0 || a > 0) - - neg := false - if a < 0 { - neg = (b-a)&1 == 0 - a, b = -b, -a - } - - z.abs = z.abs.mulRange(nil, uint64(a), uint64(b)) - z.neg = neg - return z -} - -// Binomial sets z to the binomial coefficient C(n, k) and returns z. -func (z *Int) Binomial(n, k int64) *Int { - if k > n || k < 0 { - return z.SetInt64(0) - } - // reduce the number of multiplications by reducing k - if k > n-k { - k = n - k // C(n, k) == C(n, n-k) - } - // C(n, k) == n * (n-1) * ... * (n-k+1) / k * (k-1) * ... * 1 - // == n * (n-1) * ... * (n-k+1) / 1 * (1+1) * ... * k - // - // Using the multiplicative formula produces smaller values - // at each step, requiring fewer allocations and computations: - // - // z = 1 - // for i := 0; i < k; i = i+1 { - // z *= n-i - // z /= i+1 - // } - // - // finally to avoid computing i+1 twice per loop: - // - // z = 1 - // i := 0 - // for i < k { - // z *= n-i - // i++ - // z /= i - // } - var N, K, i, t Int - N.SetInt64(n) - K.SetInt64(k) - z.Set(intOne) - for i.Cmp(&K) < 0 { - z.Mul(z, t.Sub(&N, &i)) - i.Add(&i, intOne) - z.Quo(z, &i) - } - return z -} - -// Quo sets z to the quotient x/y for y != 0 and returns z. -// If y == 0, a division-by-zero run-time panic occurs. -// Quo implements truncated division (like Go); see [Int.QuoRem] for more details. -func (z *Int) Quo(x, y *Int) *Int { - z.abs, _ = z.abs.div(nil, nil, x.abs, y.abs) - z.neg = len(z.abs) > 0 && x.neg != y.neg // 0 has no sign - return z -} - -// Rem sets z to the remainder x%y for y != 0 and returns z. -// If y == 0, a division-by-zero run-time panic occurs. -// Rem implements truncated modulus (like Go); see [Int.QuoRem] for more details. -func (z *Int) Rem(x, y *Int) *Int { - _, z.abs = nat(nil).div(nil, z.abs, x.abs, y.abs) - z.neg = len(z.abs) > 0 && x.neg // 0 has no sign - return z -} - -// QuoRem sets z to the quotient x/y and r to the remainder x%y -// and returns the pair (z, r) for y != 0. -// If y == 0, a division-by-zero run-time panic occurs. -// -// QuoRem implements T-division and modulus (like Go): -// -// q = x/y with the result truncated to zero -// r = x - y*q -// -// (See Daan Leijen, “Division and Modulus for Computer Scientists”.) -// See [Int.DivMod] for Euclidean division and modulus (unlike Go). -func (z *Int) QuoRem(x, y, r *Int) (*Int, *Int) { - z.abs, r.abs = z.abs.div(nil, r.abs, x.abs, y.abs) - z.neg, r.neg = len(z.abs) > 0 && x.neg != y.neg, len(r.abs) > 0 && x.neg // 0 has no sign - return z, r -} - -// Div sets z to the quotient x/y for y != 0 and returns z. -// If y == 0, a division-by-zero run-time panic occurs. -// Div implements Euclidean division (unlike Go); see [Int.DivMod] for more details. -func (z *Int) Div(x, y *Int) *Int { - y_neg := y.neg // z may be an alias for y - var r Int - z.QuoRem(x, y, &r) - if r.neg { - if y_neg { - z.Add(z, intOne) - } else { - z.Sub(z, intOne) - } - } - return z -} - -// Mod sets z to the modulus x%y for y != 0 and returns z. -// If y == 0, a division-by-zero run-time panic occurs. -// Mod implements Euclidean modulus (unlike Go); see [Int.DivMod] for more details. -func (z *Int) Mod(x, y *Int) *Int { - y0 := y // save y - if z == y || alias(z.abs, y.abs) { - y0 = new(Int).Set(y) - } - var q Int - q.QuoRem(x, y, z) - if z.neg { - if y0.neg { - z.Sub(z, y0) - } else { - z.Add(z, y0) - } - } - return z -} - -// DivMod sets z to the quotient x div y and m to the modulus x mod y -// and returns the pair (z, m) for y != 0. -// If y == 0, a division-by-zero run-time panic occurs. -// -// DivMod implements Euclidean division and modulus (unlike Go): -// -// q = x div y such that -// m = x - y*q with 0 <= m < |y| -// -// (See Raymond T. Boute, “The Euclidean definition of the functions -// div and mod”. ACM Transactions on Programming Languages and -// Systems (TOPLAS), 14(2):127-144, New York, NY, USA, 4/1992. -// ACM press.) -// See [Int.QuoRem] for T-division and modulus (like Go). -func (z *Int) DivMod(x, y, m *Int) (*Int, *Int) { - y0 := y // save y - if z == y || alias(z.abs, y.abs) { - y0 = new(Int).Set(y) - } - z.QuoRem(x, y, m) - if m.neg { - if y0.neg { - z.Add(z, intOne) - m.Sub(m, y0) - } else { - z.Sub(z, intOne) - m.Add(m, y0) - } - } - return z, m -} - -// Rounding modes that determine how the integer quotient is adjusted in an integer division. -// See Daan Leijen, “Division and Modulus for Computer Scientists”, for details. -const ( - Trunc = ToZero // T-division (same as Go division) - Floor = ToNegativeInf // F-division - Round = ToNearestEven // R-division - Ceil = ToPositiveInf // C-division -) - -// Divide computes the integer quotient q and remainder r such that -// -// q = f(x/y) -// r = x - y*q -// -// where f is described by the rounding mode, -// which must be one of [Trunc], [Floor], [Round] or [Ceil]. -// Divide sets z to q if z != nil, updates r if r != nil, -// and returns the pair (z, r) if y != 0. -// If y == 0, a division-by-zero run-time panic occurs. -func (z *Int) Divide(x, y, r *Int, mode RoundingMode) (*Int, *Int) { - // TODO: optimize the code where z or r is nil - var z_abs nat - if z != nil { - z_abs = z.abs - } - var r_neg bool - var r_abs nat - if r != nil { - r_abs = r.abs - } - y_abs := y.abs // save y - if z == y || alias(z_abs, y.abs) { - y_abs = nat(nil).set(y.abs) - } - neg := x.neg != y.neg - z_abs, r_abs = z_abs.div(nil, r_abs, x.abs, y.abs) - if len(r_abs) > 0 { - switch mode { - case Trunc: - r_neg = x.neg - case Floor: - r_neg = y.neg - if neg { - z_abs = z_abs.add(z_abs, natOne) - r_abs = r_abs.sub(y_abs, r_abs) - } - case Ceil: - r_neg = !y.neg - if !neg { - z_abs = z_abs.add(z_abs, natOne) - r_abs = r_abs.sub(y_abs, r_abs) - } - case Round: - switch nat(nil).mul(nil, r_abs, natTwo).cmp(y_abs) { - case -1: - r_neg = x.neg - case 0: - even := len(z_abs) == 0 || z_abs[0]&1 == 0 - if even { - r_neg = x.neg - break - } - fallthrough - case 1: - r_neg = !x.neg - z_abs = z_abs.add(z_abs, natOne) - r_abs = r_abs.sub(y_abs, r_abs) - } - default: - panic("unsupported rounding mode") - } - } - if z != nil { - z.abs = z_abs - z.neg = neg && len(z_abs) > 0 // 0 has no sign - } - if r != nil { - r.abs = r_abs - r.neg = r_neg - } - return z, r -} - -// Cmp compares x and y and returns: -// - -1 if x < y; -// - 0 if x == y; -// - +1 if x > y. -func (x *Int) Cmp(y *Int) (r int) { - // x cmp y == x cmp y - // x cmp (-y) == x - // (-x) cmp y == y - // (-x) cmp (-y) == -(x cmp y) - switch { - case x == y: - // nothing to do - case x.neg == y.neg: - r = x.abs.cmp(y.abs) - if x.neg { - r = -r - } - case x.neg: - r = -1 - default: - r = 1 - } - return -} - -// CmpAbs compares the absolute values of x and y and returns: -// - -1 if |x| < |y|; -// - 0 if |x| == |y|; -// - +1 if |x| > |y|. -func (x *Int) CmpAbs(y *Int) int { - return x.abs.cmp(y.abs) -} - -// low32 returns the least significant 32 bits of x. -func low32(x nat) uint32 { - if len(x) == 0 { - return 0 - } - return uint32(x[0]) -} - -// low64 returns the least significant 64 bits of x. -func low64(x nat) uint64 { - if len(x) == 0 { - return 0 - } - v := uint64(x[0]) - if _W == 32 && len(x) > 1 { - return uint64(x[1])<<32 | v - } - return v -} - -// Int64 returns the int64 representation of x. -// If x cannot be represented in an int64, the result is undefined. -func (x *Int) Int64() int64 { - v := int64(low64(x.abs)) - if x.neg { - v = -v - } - return v -} - -// Uint64 returns the uint64 representation of x. -// If x cannot be represented in a uint64, the result is undefined. -func (x *Int) Uint64() uint64 { - return low64(x.abs) -} - -// IsInt64 reports whether x can be represented as an int64. -func (x *Int) IsInt64() bool { - if len(x.abs) <= 64/_W { - w := int64(low64(x.abs)) - return w >= 0 || x.neg && w == -w - } - return false -} - -// IsUint64 reports whether x can be represented as a uint64. -func (x *Int) IsUint64() bool { - return !x.neg && len(x.abs) <= 64/_W -} - -// Float64 returns the float64 value nearest x, -// and an indication of any rounding that occurred. -func (x *Int) Float64() (float64, Accuracy) { - n := x.abs.bitLen() // NB: still uses slow crypto impl! - if n == 0 { - return 0.0, Exact - } - - // Fast path: no more than 53 significant bits. - if n <= 53 || n < 64 && n-int(x.abs.trailingZeroBits()) <= 53 { - f := float64(low64(x.abs)) - if x.neg { - f = -f - } - return f, Exact - } - - return new(Float).SetInt(x).Float64() -} - -// SetString sets z to the value of s, interpreted in the given base, -// and returns z and a boolean indicating success. The entire string -// (not just a prefix) must be valid for success. If SetString fails, -// the value of z is undefined but the returned value is nil. -// -// The base argument must be 0 or a value between 2 and [MaxBase]. -// For base 0, the number prefix determines the actual base: A prefix of -// “0b” or “0B” selects base 2, “0”, “0o” or “0O” selects base 8, -// and “0x” or “0X” selects base 16. Otherwise, the selected base is 10 -// and no prefix is accepted. -// -// For bases <= 36, lower and upper case letters are considered the same: -// The letters 'a' to 'z' and 'A' to 'Z' represent digit values 10 to 35. -// For bases > 36, the upper case letters 'A' to 'Z' represent the digit -// values 36 to 61. -// -// For base 0, an underscore character “_” may appear between a base -// prefix and an adjacent digit, and between successive digits; such -// underscores do not change the value of the number. -// Incorrect placement of underscores is reported as an error if there -// are no other errors. If base != 0, underscores are not recognized -// and act like any other character that is not a valid digit. -func (z *Int) SetString(s string, base int) (*Int, bool) { - return z.setFromScanner(strings.NewReader(s), base) -} - -// setFromScanner implements SetString given an io.ByteScanner. -// For documentation see comments of SetString. -func (z *Int) setFromScanner(r io.ByteScanner, base int) (*Int, bool) { - if _, _, err := z.scan(r, base); err != nil { - return nil, false - } - // entire content must have been consumed - if _, err := r.ReadByte(); err != io.EOF { - return nil, false - } - return z, true // err == io.EOF => scan consumed all content of r -} - -// SetBytes interprets buf as the bytes of a big-endian unsigned -// integer, sets z to that value, and returns z. -func (z *Int) SetBytes(buf []byte) *Int { - z.abs = z.abs.setBytes(buf) - z.neg = false - return z -} - -// Bytes returns the absolute value of x as a big-endian byte slice. -// -// To use a fixed length slice, or a preallocated one, use [Int.FillBytes]. -func (x *Int) Bytes() []byte { - // This function is used in cryptographic operations. It must not leak - // anything but the Int's sign and bit size through side-channels. Any - // changes must be reviewed by a security expert. - buf := make([]byte, len(x.abs)*_S) - return buf[x.abs.bytes(buf):] -} - -// FillBytes sets buf to the absolute value of x, storing it as a zero-extended -// big-endian byte slice, and returns buf. -// -// If the absolute value of x doesn't fit in buf, FillBytes will panic. -func (x *Int) FillBytes(buf []byte) []byte { - // Clear whole buffer. - clear(buf) - x.abs.bytes(buf) - return buf -} - -// BitLen returns the length of the absolute value of x in bits. -// The bit length of 0 is 0. -func (x *Int) BitLen() int { - // This function is used in cryptographic operations. It must not leak - // anything but the Int's sign and bit size through side-channels. Any - // changes must be reviewed by a security expert. - return x.abs.bitLen() -} - -// TrailingZeroBits returns the number of consecutive least significant zero -// bits of |x|. -func (x *Int) TrailingZeroBits() uint { - return x.abs.trailingZeroBits() -} - -// Exp sets z = x**y mod |m| (i.e. the sign of m is ignored), and returns z. -// If m == nil or m == 0, z = x**y unless y <= 0 then z = 1. If m != 0, y < 0, -// and x and m are not relatively prime, z is unchanged and nil is returned. -// -// Modular exponentiation of inputs of a particular size is not a -// cryptographically constant-time operation. -func (z *Int) Exp(x, y, m *Int) *Int { - return z.exp(x, y, m, false) -} - -func (z *Int) expSlow(x, y, m *Int) *Int { - return z.exp(x, y, m, true) -} - -func (z *Int) exp(x, y, m *Int, slow bool) *Int { - // See Knuth, volume 2, section 4.6.3. - xWords := x.abs - if y.neg { - if m == nil || len(m.abs) == 0 { - return z.SetInt64(1) - } - // for y < 0: x**y mod m == (x**(-1))**|y| mod m - inverse := new(Int).ModInverse(x, m) - if inverse == nil { - return nil - } - xWords = inverse.abs - } - yWords := y.abs - - var mWords nat - if m != nil { - if z == m || alias(z.abs, m.abs) { - m = new(Int).Set(m) - } - mWords = m.abs // m.abs may be nil for m == 0 - } - - z.abs = z.abs.expNN(nil, xWords, yWords, mWords, slow) - z.neg = len(z.abs) > 0 && x.neg && len(yWords) > 0 && yWords[0]&1 == 1 // 0 has no sign - if z.neg && len(mWords) > 0 { - // make modulus result positive - z.abs = z.abs.sub(mWords, z.abs) // z == x**y mod |m| && 0 <= z < |m| - z.neg = false - } - - return z -} - -// GCD sets z to the greatest common divisor of a and b and returns z. -// If x or y are not nil, GCD sets their value such that z = a*x + b*y. -// -// a and b may be positive, zero or negative. (Before Go 1.14 both had -// to be > 0.) Regardless of the signs of a and b, z is always >= 0. -// -// If a == b == 0, GCD sets z = x = y = 0. -// -// If a == 0 and b != 0, GCD sets z = |b|, x = 0, y = sign(b) * 1. -// -// If a != 0 and b == 0, GCD sets z = |a|, x = sign(a) * 1, y = 0. -func (z *Int) GCD(x, y, a, b *Int) *Int { - if len(a.abs) == 0 || len(b.abs) == 0 { - lenA, lenB, negA, negB := len(a.abs), len(b.abs), a.neg, b.neg - if lenA == 0 { - z.Set(b) - } else { - z.Set(a) - } - z.neg = false - if x != nil { - if lenA == 0 { - x.SetUint64(0) - } else { - x.SetUint64(1) - x.neg = negA - } - } - if y != nil { - if lenB == 0 { - y.SetUint64(0) - } else { - y.SetUint64(1) - y.neg = negB - } - } - return z - } - - return z.lehmerGCD(x, y, a, b) -} - -// lehmerSimulate attempts to simulate several Euclidean update steps -// using the leading digits of A and B. It returns u0, u1, v0, v1 -// such that A and B can be updated as: -// -// A = u0*A + v0*B -// B = u1*A + v1*B -// -// Requirements: A >= B and len(B.abs) >= 2 -// Since we are calculating with full words to avoid overflow, -// we use 'even' to track the sign of the cosequences. -// For even iterations: u0, v1 >= 0 && u1, v0 <= 0 -// For odd iterations: u0, v1 <= 0 && u1, v0 >= 0 -func lehmerSimulate(A, B *Int) (u0, u1, v0, v1 Word, even bool) { - // initialize the digits - var a1, a2, u2, v2 Word - - m := len(B.abs) // m >= 2 - n := len(A.abs) // n >= m >= 2 - - // extract the top Word of bits from A and B - h := nlz(A.abs[n-1]) - a1 = A.abs[n-1]<>(_W-h) - // B may have implicit zero words in the high bits if the lengths differ - switch { - case n == m: - a2 = B.abs[n-1]<>(_W-h) - case n == m+1: - a2 = B.abs[n-2] >> (_W - h) - default: - a2 = 0 - } - - // Since we are calculating with full words to avoid overflow, - // we use 'even' to track the sign of the cosequences. - // For even iterations: u0, v1 >= 0 && u1, v0 <= 0 - // For odd iterations: u0, v1 <= 0 && u1, v0 >= 0 - // The first iteration starts with k=1 (odd). - even = false - // variables to track the cosequences - u0, u1, u2 = 0, 1, 0 - v0, v1, v2 = 0, 0, 1 - - // Calculate the quotient and cosequences using Collins' stopping condition. - // Note that overflow of a Word is not possible when computing the remainder - // sequence and cosequences since the cosequence size is bounded by the input size. - // See section 4.2 of Jebelean for details. - for a2 >= v2 && a1-a2 >= v1+v2 { - q, r := a1/a2, a1%a2 - a1, a2 = a2, r - u0, u1, u2 = u1, u2, u1+q*u2 - v0, v1, v2 = v1, v2, v1+q*v2 - even = !even - } - return -} - -// lehmerUpdate updates the inputs A and B such that: -// -// A = u0*A + v0*B -// B = u1*A + v1*B -// -// where the signs of u0, u1, v0, v1 are given by even -// For even == true: u0, v1 >= 0 && u1, v0 <= 0 -// For even == false: u0, v1 <= 0 && u1, v0 >= 0 -// q, r, s, t are temporary variables to avoid allocations in the multiplication. -func lehmerUpdate(A, B, q, r *Int, u0, u1, v0, v1 Word, even bool) { - mulW(q, B, even, v0) - mulW(r, A, even, u1) - mulW(A, A, !even, u0) - mulW(B, B, !even, v1) - A.Add(A, q) - B.Add(B, r) -} - -// mulW sets z = x * (-?)w -// where the minus sign is present when neg is true. -func mulW(z, x *Int, neg bool, w Word) { - z.abs = z.abs.mulAddWW(x.abs, w, 0) - z.neg = x.neg != neg -} - -// euclidUpdate performs a single step of the Euclidean GCD algorithm -// if extended is true, it also updates the cosequence Ua, Ub. -// q and r are used as temporaries; the initial values are ignored. -func euclidUpdate(A, B, Ua, Ub, q, r *Int, extended bool) (nA, nB, nr, nUa, nUb *Int) { - q.QuoRem(A, B, r) - - if extended { - // Ua, Ub = Ub, Ua-q*Ub - q.Mul(q, Ub) - Ua, Ub = Ub, Ua - Ub.Sub(Ub, q) - } - - return B, r, A, Ua, Ub -} - -// sixIntPool is used to reduce allocation of limbs used in temporary integers -// used to calculate lehmerGCD. -type sixIntPool struct { - pool sync.Pool -} - -func (t *sixIntPool) put(data *[6]Int) { - data[0].SetInt64(0) - data[1].SetInt64(0) - data[2].SetInt64(0) - data[3].SetInt64(0) - data[4].SetInt64(0) - data[5].SetInt64(0) - t.pool.Put(data) -} - -func (t *sixIntPool) get() *[6]Int { - return t.pool.Get().(*[6]Int) -} - -var sixIntP = sixIntPool{ - sync.Pool{ - New: func() any { - return &[6]Int{} - }, - }, -} - -// lehmerGCD sets z to the greatest common divisor of a and b, -// which both must be != 0, and returns z. -// If x or y are not nil, their values are set such that z = a*x + b*y. -// See Knuth, The Art of Computer Programming, Vol. 2, Section 4.5.2, Algorithm L. -// This implementation uses the improved condition by Collins requiring only one -// quotient and avoiding the possibility of single Word overflow. -// See Jebelean, "Improving the multiprecision Euclidean algorithm", -// Design and Implementation of Symbolic Computation Systems, pp 45-58. -// The cosequences are updated according to Algorithm 10.45 from -// Cohen et al. "Handbook of Elliptic and Hyperelliptic Curve Cryptography" pp 192. -func (z *Int) lehmerGCD(x, y, a, b *Int) *Int { - // recycle limbs to reduce allocations. - data := sixIntP.get() - defer sixIntP.put(data) - - var A, B, Ua, Ub *Int = &data[0], &data[1], &data[2], &data[3] - - A.Abs(a) - B.Abs(b) - - extended := x != nil || y != nil - - if extended { - // Ua (Ub) tracks how many times input a has been accumulated into A (B). - Ua.SetInt64(1) - } - - // temp variables for multiprecision update - q := &data[4] - r := &data[5] - - // ensure A >= B - if A.abs.cmp(B.abs) < 0 { - A, B = B, A - Ub, Ua = Ua, Ub - } - - // loop invariant A >= B - for len(B.abs) > 1 { - // Attempt to calculate in single-precision using leading words of A and B. - u0, u1, v0, v1, even := lehmerSimulate(A, B) - - // multiprecision Step - if v0 != 0 { - // Simulate the effect of the single-precision steps using the cosequences. - // A = u0*A + v0*B - // B = u1*A + v1*B - lehmerUpdate(A, B, q, r, u0, u1, v0, v1, even) - - if extended { - // Ua = u0*Ua + v0*Ub - // Ub = u1*Ua + v1*Ub - lehmerUpdate(Ua, Ub, q, r, u0, u1, v0, v1, even) - } - - } else { - // Single-digit calculations failed to simulate any quotients. - // Do a standard Euclidean step. - A, B, r, Ua, Ub = euclidUpdate(A, B, Ua, Ub, q, r, extended) - } - } - - if len(B.abs) > 0 { - // extended Euclidean algorithm base case if B is a single Word - if len(A.abs) > 1 { - // A is longer than a single Word, so one update is needed. - A, B, r, Ua, Ub = euclidUpdate(A, B, Ua, Ub, q, r, extended) - } - if len(B.abs) > 0 { - // A and B are both a single Word. - aWord, bWord := A.abs[0], B.abs[0] - if extended { - var ua, ub, va, vb Word - ua, ub = 1, 0 - va, vb = 0, 1 - even := true - for bWord != 0 { - q, r := aWord/bWord, aWord%bWord - aWord, bWord = bWord, r - ua, ub = ub, ua+q*ub - va, vb = vb, va+q*vb - even = !even - } - - mulW(Ua, Ua, !even, ua) - mulW(Ub, Ub, even, va) - Ua.Add(Ua, Ub) - } else { - for bWord != 0 { - aWord, bWord = bWord, aWord%bWord - } - } - A.abs[0] = aWord - } - } - negA := a.neg - if y != nil { - // avoid aliasing b needed in the division below - if y == b { - B.Set(b) - } else { - B = b - } - // y = (z - a*x)/b - y.Mul(a, Ua) // y can safely alias a - if negA { - y.neg = !y.neg - } - y.Sub(A, y) - y.Div(y, B) - } - - if x != nil { - x.Set(Ua) - if negA { - x.neg = !x.neg - } - } - - z.Set(A) - - return z -} - -// Rand sets z to a pseudo-random number in [0, n) and returns z. -// -// As this uses the [math/rand] package, it must not be used for -// security-sensitive work. Use [crypto/rand.Int] instead. -func (z *Int) Rand(rnd *rand.Rand, n *Int) *Int { - // z.neg is not modified before the if check, because z and n might alias. - if n.neg || len(n.abs) == 0 { - z.neg = false - z.abs = nil - return z - } - z.neg = false - z.abs = z.abs.random(rnd, n.abs, n.abs.bitLen()) - return z -} - -// twoIntPool is used to reduce allocation of limbs used in temporary integers -// used to calculate ModInverse. -type twoIntPool struct { - pool sync.Pool -} - -func (t *twoIntPool) put(data *[2]Int) { - data[0].SetInt64(0) - data[1].SetInt64(0) - t.pool.Put(data) -} - -func (t *twoIntPool) get() *[2]Int { - return t.pool.Get().(*[2]Int) -} - -var twoIntP = twoIntPool{ - sync.Pool{ - New: func() any { - return &[2]Int{} - }, - }, -} - -// ModInverse sets z to the multiplicative inverse of g in the ring ℤ/nℤ -// and returns z. If g and n are not relatively prime, g has no multiplicative -// inverse in the ring ℤ/nℤ. In this case, z is unchanged and the return value -// is nil. If n == 0, a division-by-zero run-time panic occurs. -func (z *Int) ModInverse(g, n *Int) *Int { - // GCD expects parameters a and b to be > 0. - if n.neg { - var n2 Int - n = n2.Neg(n) - } - if g.neg { - var g2 Int - g = g2.Mod(g, n) - } - - // recycle limbs to reduce allocations. - data := twoIntP.get() - defer twoIntP.put(data) - - var d, x *Int = &data[0], &data[1] - d.GCD(x, nil, g, n) - - // if and only if d==1, g and n are relatively prime - if d.Cmp(intOne) != 0 { - return nil - } - - // x and y are such that g*x + n*y = 1, therefore x is the inverse element, - // but it may be negative, so convert to the range 0 <= z < |n| - if x.neg { - z.Add(x, n) - } else { - z.Set(x) - } - - return z -} - -func (z nat) modInverse(g, n nat) nat { - // TODO(rsc): ModInverse should be implemented in terms of this function. - return (&Int{abs: z}).ModInverse(&Int{abs: g}, &Int{abs: n}).abs -} - -// Jacobi returns the Jacobi symbol (x/y), either +1, -1, or 0. -// The y argument must be an odd integer. -func Jacobi(x, y *Int) int { - if len(y.abs) == 0 || y.abs[0]&1 == 0 { - panic(fmt.Sprintf("big: invalid 2nd argument to Int.Jacobi: need odd integer but got %s", y.String())) - } - - // We use the formulation described in chapter 2, section 2.4, - // "The Yacas Book of Algorithms": - // http://yacas.sourceforge.net/Algo.book.pdf - - var a, b, c Int - a.Set(x) - b.Set(y) - j := 1 - - if b.neg { - if a.neg { - j = -1 - } - b.neg = false - } - - for { - if b.Cmp(intOne) == 0 { - return j - } - if len(a.abs) == 0 { - return 0 - } - a.Mod(&a, &b) - if len(a.abs) == 0 { - return 0 - } - // a > 0 - - // handle factors of 2 in 'a' - s := a.abs.trailingZeroBits() - if s&1 != 0 { - bmod8 := b.abs[0] & 7 - if bmod8 == 3 || bmod8 == 5 { - j = -j - } - } - c.Rsh(&a, s) // a = 2^s*c - - // swap numerator and denominator - if b.abs[0]&3 == 3 && c.abs[0]&3 == 3 { - j = -j - } - a.Set(&b) - b.Set(&c) - } -} - -// modSqrt3Mod4 uses the identity -// -// (a^((p+1)/4))^2 mod p -// == u^(p+1) mod p -// == u^2 mod p -// -// to calculate the square root of any quadratic residue mod p quickly for 3 -// mod 4 primes. -func (z *Int) modSqrt3Mod4Prime(x, p *Int) *Int { - e := new(Int).Add(p, intOne) // e = p + 1 - e.Rsh(e, 2) // e = (p + 1) / 4 - z.Exp(x, e, p) // z = x^e mod p - return z -} - -// modSqrt5Mod8Prime uses Atkin's observation that 2 is not a square mod p -// -// alpha == (2*a)^((p-5)/8) mod p -// beta == 2*a*alpha^2 mod p is a square root of -1 -// b == a*alpha*(beta-1) mod p is a square root of a -// -// to calculate the square root of any quadratic residue mod p quickly for 5 -// mod 8 primes. -func (z *Int) modSqrt5Mod8Prime(x, p *Int) *Int { - // p == 5 mod 8 implies p = e*8 + 5 - // e is the quotient and 5 the remainder on division by 8 - e := new(Int).Rsh(p, 3) // e = (p - 5) / 8 - tx := new(Int).Lsh(x, 1) // tx = 2*x - alpha := new(Int).Exp(tx, e, p) - beta := new(Int).Mul(alpha, alpha) - beta.Mod(beta, p) - beta.Mul(beta, tx) - beta.Mod(beta, p) - beta.Sub(beta, intOne) - beta.Mul(beta, x) - beta.Mod(beta, p) - beta.Mul(beta, alpha) - z.Mod(beta, p) - return z -} - -// modSqrtTonelliShanks uses the Tonelli-Shanks algorithm to find the square -// root of a quadratic residue modulo any prime. -func (z *Int) modSqrtTonelliShanks(x, p *Int) *Int { - // Break p-1 into s*2^e such that s is odd. - var s Int - s.Sub(p, intOne) - e := s.abs.trailingZeroBits() - s.Rsh(&s, e) - - // find some non-square n - var n Int - n.SetInt64(2) - for Jacobi(&n, p) != -1 { - n.Add(&n, intOne) - } - - // Core of the Tonelli-Shanks algorithm. Follows the description in - // section 6 of "Square roots from 1; 24, 51, 10 to Dan Shanks" by Ezra - // Brown: - // https://www.maa.org/sites/default/files/pdf/upload_library/22/Polya/07468342.di020786.02p0470a.pdf - var y, b, g, t Int - y.Add(&s, intOne) - y.Rsh(&y, 1) - y.Exp(x, &y, p) // y = x^((s+1)/2) - b.Exp(x, &s, p) // b = x^s - g.Exp(&n, &s, p) // g = n^s - r := e - for { - // find the least m such that ord_p(b) = 2^m - var m uint - t.Set(&b) - for t.Cmp(intOne) != 0 { - t.Mul(&t, &t).Mod(&t, p) - m++ - } - - if m == 0 { - return z.Set(&y) - } - - t.SetInt64(0).SetBit(&t, int(r-m-1), 1).Exp(&g, &t, p) - // t = g^(2^(r-m-1)) mod p - g.Mul(&t, &t).Mod(&g, p) // g = g^(2^(r-m)) mod p - y.Mul(&y, &t).Mod(&y, p) - b.Mul(&b, &g).Mod(&b, p) - r = m - } -} - -// ModSqrt sets z to a square root of x mod p if such a square root exists, and -// returns z. The modulus p must be an odd prime. If x is not a square mod p, -// ModSqrt leaves z unchanged and returns nil. This function panics if p is -// not an odd integer, its behavior is undefined if p is odd but not prime. -func (z *Int) ModSqrt(x, p *Int) *Int { - switch Jacobi(x, p) { - case -1: - return nil // x is not a square mod p - case 0: - return z.SetInt64(0) // sqrt(0) mod p = 0 - case 1: - break - } - if x.neg || x.Cmp(p) >= 0 { // ensure 0 <= x < p - x = new(Int).Mod(x, p) - } - - switch { - case p.abs[0]%4 == 3: - // Check whether p is 3 mod 4, and if so, use the faster algorithm. - return z.modSqrt3Mod4Prime(x, p) - case p.abs[0]%8 == 5: - // Check whether p is 5 mod 8, use Atkin's algorithm. - return z.modSqrt5Mod8Prime(x, p) - default: - // Otherwise, use Tonelli-Shanks. - return z.modSqrtTonelliShanks(x, p) - } -} - -// Lsh sets z = x << n and returns z. -func (z *Int) Lsh(x *Int, n uint) *Int { - z.abs = z.abs.lsh(x.abs, n) - z.neg = x.neg - return z -} - -// Rsh sets z = x >> n and returns z. -func (z *Int) Rsh(x *Int, n uint) *Int { - if x.neg { - // (-x) >> s == ^(x-1) >> s == ^((x-1) >> s) == -(((x-1) >> s) + 1) - t := z.abs.sub(x.abs, natOne) // no underflow because |x| > 0 - t = t.rsh(t, n) - z.abs = t.add(t, natOne) - z.neg = true // z cannot be zero if x is negative - return z - } - - z.abs = z.abs.rsh(x.abs, n) - z.neg = false - return z -} - -// Bit returns the value of the i'th bit of x. That is, it -// returns (x>>i)&1. The bit index i must be >= 0. -func (x *Int) Bit(i int) uint { - if i == 0 { - // optimization for common case: odd/even test of x - if len(x.abs) > 0 { - return uint(x.abs[0] & 1) // bit 0 is same for -x - } - return 0 - } - if i < 0 { - panic("negative bit index") - } - if x.neg { - t := nat(nil).sub(x.abs, natOne) - return t.bit(uint(i)) ^ 1 - } - - return x.abs.bit(uint(i)) -} - -// SetBit sets z to x, with x's i'th bit set to b (0 or 1). -// That is, -// - if b is 1, SetBit sets z = x | (1 << i); -// - if b is 0, SetBit sets z = x &^ (1 << i); -// - if b is not 0 or 1, SetBit will panic. -func (z *Int) SetBit(x *Int, i int, b uint) *Int { - if i < 0 { - panic("negative bit index") - } - if x.neg { - t := z.abs.sub(x.abs, natOne) - t = t.setBit(t, uint(i), b^1) - z.abs = t.add(t, natOne) - z.neg = len(z.abs) > 0 - return z - } - z.abs = z.abs.setBit(x.abs, uint(i), b) - z.neg = false - return z -} - -// And sets z = x & y and returns z. -func (z *Int) And(x, y *Int) *Int { - if x.neg == y.neg { - if x.neg { - // (-x) & (-y) == ^(x-1) & ^(y-1) == ^((x-1) | (y-1)) == -(((x-1) | (y-1)) + 1) - x1 := nat(nil).sub(x.abs, natOne) - y1 := nat(nil).sub(y.abs, natOne) - z.abs = z.abs.add(z.abs.or(x1, y1), natOne) - z.neg = true // z cannot be zero if x and y are negative - return z - } - - // x & y == x & y - z.abs = z.abs.and(x.abs, y.abs) - z.neg = false - return z - } - - // x.neg != y.neg - if x.neg { - x, y = y, x // & is symmetric - } - - // x & (-y) == x & ^(y-1) == x &^ (y-1) - y1 := nat(nil).sub(y.abs, natOne) - z.abs = z.abs.andNot(x.abs, y1) - z.neg = false - return z -} - -// AndNot sets z = x &^ y and returns z. -func (z *Int) AndNot(x, y *Int) *Int { - if x.neg == y.neg { - if x.neg { - // (-x) &^ (-y) == ^(x-1) &^ ^(y-1) == ^(x-1) & (y-1) == (y-1) &^ (x-1) - x1 := nat(nil).sub(x.abs, natOne) - y1 := nat(nil).sub(y.abs, natOne) - z.abs = z.abs.andNot(y1, x1) - z.neg = false - return z - } - - // x &^ y == x &^ y - z.abs = z.abs.andNot(x.abs, y.abs) - z.neg = false - return z - } - - if x.neg { - // (-x) &^ y == ^(x-1) &^ y == ^(x-1) & ^y == ^((x-1) | y) == -(((x-1) | y) + 1) - x1 := nat(nil).sub(x.abs, natOne) - z.abs = z.abs.add(z.abs.or(x1, y.abs), natOne) - z.neg = true // z cannot be zero if x is negative and y is positive - return z - } - - // x &^ (-y) == x &^ ^(y-1) == x & (y-1) - y1 := nat(nil).sub(y.abs, natOne) - z.abs = z.abs.and(x.abs, y1) - z.neg = false - return z -} - -// Or sets z = x | y and returns z. -func (z *Int) Or(x, y *Int) *Int { - if x.neg == y.neg { - if x.neg { - // (-x) | (-y) == ^(x-1) | ^(y-1) == ^((x-1) & (y-1)) == -(((x-1) & (y-1)) + 1) - x1 := nat(nil).sub(x.abs, natOne) - y1 := nat(nil).sub(y.abs, natOne) - z.abs = z.abs.add(z.abs.and(x1, y1), natOne) - z.neg = true // z cannot be zero if x and y are negative - return z - } - - // x | y == x | y - z.abs = z.abs.or(x.abs, y.abs) - z.neg = false - return z - } - - // x.neg != y.neg - if x.neg { - x, y = y, x // | is symmetric - } - - // x | (-y) == x | ^(y-1) == ^((y-1) &^ x) == -(^((y-1) &^ x) + 1) - y1 := nat(nil).sub(y.abs, natOne) - z.abs = z.abs.add(z.abs.andNot(y1, x.abs), natOne) - z.neg = true // z cannot be zero if one of x or y is negative - return z -} - -// Xor sets z = x ^ y and returns z. -func (z *Int) Xor(x, y *Int) *Int { - if x.neg == y.neg { - if x.neg { - // (-x) ^ (-y) == ^(x-1) ^ ^(y-1) == (x-1) ^ (y-1) - x1 := nat(nil).sub(x.abs, natOne) - y1 := nat(nil).sub(y.abs, natOne) - z.abs = z.abs.xor(x1, y1) - z.neg = false - return z - } - - // x ^ y == x ^ y - z.abs = z.abs.xor(x.abs, y.abs) - z.neg = false - return z - } - - // x.neg != y.neg - if x.neg { - x, y = y, x // ^ is symmetric - } - - // x ^ (-y) == x ^ ^(y-1) == ^(x ^ (y-1)) == -((x ^ (y-1)) + 1) - y1 := nat(nil).sub(y.abs, natOne) - z.abs = z.abs.add(z.abs.xor(x.abs, y1), natOne) - z.neg = true // z cannot be zero if only one of x or y is negative - return z -} - -// Not sets z = ^x and returns z. -func (z *Int) Not(x *Int) *Int { - if x.neg { - // ^(-x) == ^(^(x-1)) == x-1 - z.abs = z.abs.sub(x.abs, natOne) - z.neg = false - return z - } - - // ^x == -x-1 == -(x+1) - z.abs = z.abs.add(x.abs, natOne) - z.neg = true // z cannot be zero if x is positive - return z -} - -// Sqrt sets z to ⌊√x⌋, the largest integer such that z² ≤ x, and returns z. -// It panics if x is negative. -func (z *Int) Sqrt(x *Int) *Int { - if x.neg { - panic("square root of negative number") - } - z.neg = false - z.abs = z.abs.sqrt(nil, x.abs) - return z -} +// Note: Handle Int.Divide remainder aliasing when r aliases y (#80882)