diff options
Diffstat (limited to 'libgo/go/math/big/int.go')
-rw-r--r-- | libgo/go/math/big/int.go | 371 |
1 files changed, 239 insertions, 132 deletions
diff --git a/libgo/go/math/big/int.go b/libgo/go/math/big/int.go index 0eda9cd..47a288a 100644 --- a/libgo/go/math/big/int.go +++ b/libgo/go/math/big/int.go @@ -442,24 +442,28 @@ func (x *Int) BitLen() int { } // Exp sets z = x**y mod |m| (i.e. the sign of m is ignored), and returns z. -// If y <= 0, the result is 1 mod |m|; if m == nil or m == 0, z = x**y. +// If m == nil or m == 0, z = x**y unless y <= 0 then z = 1. // // Modular exponentation of inputs of a particular size is not a // cryptographically constant-time operation. func (z *Int) Exp(x, y, m *Int) *Int { // See Knuth, volume 2, section 4.6.3. - var yWords nat - if !y.neg { - yWords = y.abs + 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 + xWords = new(Int).ModInverse(x, m).abs } - // y >= 0 + yWords := y.abs var mWords nat if m != nil { mWords = m.abs // m.abs may be nil for m == 0 } - z.abs = z.abs.expNN(x.abs, yWords, mWords) + z.abs = z.abs.expNN(xWords, yWords, mWords) 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 @@ -485,162 +489,223 @@ func (z *Int) GCD(x, y, a, b *Int) *Int { } return z } - if x == nil && y == nil { - return z.lehmerGCD(a, b) - } - A := new(Int).Set(a) - B := new(Int).Set(b) - - X := new(Int) - lastX := new(Int).SetInt64(1) + return z.lehmerGCD(x, y, a, b) +} - q := new(Int) - temp := new(Int) +// 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]<<h | A.abs[n-2]>>(_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]<<h | B.abs[n-2]>>(_W-h) + case n == m+1: + a2 = B.abs[n-2] >> (_W - h) + default: + a2 = 0 + } - r := new(Int) - for len(B.abs) > 0 { - q, r = q.QuoRem(A, B, r) + // 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 +} - A, B, r = B, r, A +// 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, s, t *Int, u0, u1, v0, v1 Word, even bool) { + + t.abs = t.abs.setWord(u0) + s.abs = s.abs.setWord(v0) + t.neg = !even + s.neg = even + + t.Mul(A, t) + s.Mul(B, s) + + r.abs = r.abs.setWord(u1) + q.abs = q.abs.setWord(v1) + r.neg = even + q.neg = !even + + r.Mul(A, r) + q.Mul(B, q) + + A.Add(t, s) + B.Add(r, q) +} - temp.Set(X) - X.Mul(X, q) - X.Sub(lastX, X) - lastX.Set(temp) - } +// euclidUpdate performs a single step of the Euclidean GCD algorithm +// if extended is true, it also updates the cosequence Ua, Ub +func euclidUpdate(A, B, Ua, Ub, q, r, s, t *Int, extended bool) { + q, r = q.QuoRem(A, B, r) - if x != nil { - *x = *lastX - } + *A, *B, *r = *B, *r, *A - if y != nil { - // y = (z - a*x)/b - y.Mul(a, lastX) - y.Sub(A, y) - y.Div(y, b) + if extended { + // Ua, Ub = Ub, Ua - q*Ub + t.Set(Ub) + s.Mul(Ub, q) + Ub.Sub(Ua, s) + Ua.Set(t) } - - *z = *A - return z } // 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. -func (z *Int) lehmerGCD(a, b *Int) *Int { - // ensure a >= b - if a.abs.cmp(b.abs) < 0 { - a, b = b, a - } +// 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 { + var A, B, Ua, Ub *Int - // don't destroy incoming values of a and b - B := new(Int).Set(b) // must be set first in case b is an alias of z - A := z.Set(a) + A = new(Int).Set(a) + B = new(Int).Set(b) + + extended := x != nil || y != nil + + if extended { + // Ua (Ub) tracks how many times input a has been accumulated into A (B). + Ua = new(Int).SetInt64(1) + Ub = new(Int) + } // temp variables for multiprecision update - t := new(Int) + q := new(Int) r := new(Int) s := new(Int) - w := new(Int) + t := new(Int) + + // 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 { - // initialize the digits - var a1, a2, u0, u1, u2, v0, v1, 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] << h) | (A.abs[n-2] >> (_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] << h) | (B.abs[n-2] >> (_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 := a1 / a2 - a1, a2 = a2, a1-q*a2 - u0, u1, u2 = u1, u2, u1+q*u2 - v0, v1, v2 = v1, v2, v1+q*v2 - even = !even - } + // Attempt to calculate in single-precision using leading words of A and B. + u0, u1, v0, v1, even := lehmerSimulate(A, B) - // multiprecision step + // multiprecision Step if v0 != 0 { - // simulate the effect of the single precision steps using the cosequences + // 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, s, t, u0, u1, v0, v1, even) - t.abs = t.abs.setWord(u0) - s.abs = s.abs.setWord(v0) - t.neg = !even - s.neg = even - - t.Mul(A, t) - s.Mul(B, s) - - r.abs = r.abs.setWord(u1) - w.abs = w.abs.setWord(v1) - r.neg = even - w.neg = !even - - r.Mul(A, r) - w.Mul(B, w) - - A.Add(t, s) - B.Add(r, w) + if extended { + // Ua = u0*Ua + v0*Ub + // Ub = u1*Ua + v1*Ub + lehmerUpdate(Ua, Ub, q, r, s, t, u0, u1, v0, v1, even) + } } else { - // single-digit calculations failed to simluate any quotients - // do a standard Euclidean step - t.Rem(A, B) - A, B, t = B, t, A + // Single-digit calculations failed to simulate any quotients. + // Do a standard Euclidean step. + euclidUpdate(A, B, Ua, Ub, q, r, s, t, extended) } } if len(B.abs) > 0 { - // standard Euclidean algorithm base case for B a single Word + // extended Euclidean algorithm base case if B is a single Word if len(A.abs) > 1 { - // A is longer than a single Word - t.Rem(A, B) - A, B, t = B, t, A + // A is longer than a single Word, so one update is needed. + euclidUpdate(A, B, Ua, Ub, q, r, s, t, extended) } if len(B.abs) > 0 { - // A and B are both a single Word - a1, a2 := A.abs[0], B.abs[0] - for a2 != 0 { - a1, a2 = a2, a1%a2 + // 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 + } + + t.abs = t.abs.setWord(ua) + s.abs = s.abs.setWord(va) + t.neg = !even + s.neg = even + + t.Mul(Ua, t) + s.Mul(Ub, s) + + Ua.Add(t, s) + } else { + for bWord != 0 { + aWord, bWord = bWord, aWord%bWord + } } - A.abs[0] = a1 + A.abs[0] = aWord } } + + if x != nil { + *x = *Ua + } + + if y != nil { + // y = (z - a*x)/b + y.Mul(a, Ua) + y.Sub(A, y) + y.Div(y, b) + } + *z = *A + return z } @@ -659,20 +724,33 @@ func (z *Int) Rand(rnd *rand.Rand, n *Int) *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, the result is undefined. +// 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. 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 { - // GCD expects parameters a and b to be > 0. var g2 Int g = g2.Mod(g, n) } - var d Int - d.GCD(z, nil, g, n) - // x and y are such that g*x + n*y = d. Since g and n are - // relatively prime, d = 1. Taking that modulo n results in - // g*x = 1, therefore x is the inverse element. - if z.neg { - z.Add(z, n) + var d, x Int + 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 } @@ -745,6 +823,30 @@ func (z *Int) modSqrt3Mod4Prime(x, p *Int) *Int { return z } +// modSqrt5Mod8 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 { @@ -811,12 +913,17 @@ func (z *Int) ModSqrt(x, p *Int) *Int { x = new(Int).Mod(x, p) } - // Check whether p is 3 mod 4, and if so, use the faster algorithm. - if len(p.abs) > 0 && p.abs[0]%4 == 3 { + 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) } - // Otherwise, use Tonelli-Shanks. - return z.modSqrtTonelliShanks(x, p) } // Lsh sets z = x << n and returns z. |