Update libsecp256k1 subtree to latest upstream master

6 files changed, 181 insertions, 31 deletions

diff --git a/src/secp256k1/sage/gen_exhaustive_groups.sage b/src/secp256k1/sage/gen_exhaustive_groups.sage
index 3c3c984811..01d15dcdea 100644
--- a/src/secp256k1/sage/gen_exhaustive_groups.sage
+++ b/src/secp256k1/sage/gen_exhaustive_groups.sage
@@ -1,9 +1,4 @@
-# Define field size and field
-P = 2^256 - 2^32 - 977
-F = GF(P)
-BETA = F(0x7ae96a2b657c07106e64479eac3434e99cf0497512f58995c1396c28719501ee)
-
-assert(BETA != F(1) and BETA^3 == F(1))
+load("secp256k1_params.sage")
 
 orders_done = set()
 results = {}
diff --git a/src/secp256k1/sage/gen_split_lambda_constants.sage b/src/secp256k1/sage/gen_split_lambda_constants.sage
new file mode 100644
index 0000000000..7d4359e0f6
--- /dev/null
+++ b/src/secp256k1/sage/gen_split_lambda_constants.sage
@@ -0,0 +1,114 @@
+""" Generates the constants used in secp256k1_scalar_split_lambda.
+
+See the comments for secp256k1_scalar_split_lambda in src/scalar_impl.h for detailed explanations.
+"""
+
+load("secp256k1_params.sage")
+
+def inf_norm(v):
+    """Returns the infinity norm of a vector."""
+    return max(map(abs, v))
+
+def gauss_reduction(i1, i2):
+    v1, v2 = i1.copy(), i2.copy()
+    while True:
+        if inf_norm(v2) < inf_norm(v1):
+            v1, v2 = v2, v1
+        # This is essentially
+        #    m = round((v1[0]*v2[0] + v1[1]*v2[1]) / (inf_norm(v1)**2))
+        # (rounding to the nearest integer) without relying on floating point arithmetic.
+        m = ((v1[0]*v2[0] + v1[1]*v2[1]) + (inf_norm(v1)**2) // 2) // (inf_norm(v1)**2)
+        if m == 0:
+            return v1, v2
+        v2[0] -= m*v1[0]
+        v2[1] -= m*v1[1]
+
+def find_split_constants_gauss():
+    """Find constants for secp256k1_scalar_split_lamdba using gauss reduction."""
+    (v11, v12), (v21, v22) = gauss_reduction([0, N], [1, int(LAMBDA)])
+
+    # We use related vectors in secp256k1_scalar_split_lambda.
+    A1, B1 = -v21, -v11
+    A2, B2 = v22, -v21
+
+    return A1, B1, A2, B2
+
+def find_split_constants_explicit_tof():
+    """Find constants for secp256k1_scalar_split_lamdba using the trace of Frobenius.
+
+    See Benjamin Smith: "Easy scalar decompositions for efficient scalar multiplication on
+    elliptic curves and genus 2 Jacobians" (https://eprint.iacr.org/2013/672), Example 2
+    """
+    assert P % 3 == 1 # The paper says P % 3 == 2 but that appears to be a mistake, see [10].
+    assert C.j_invariant() == 0
+
+    t = C.trace_of_frobenius()
+
+    c = Integer(sqrt((4*P - t**2)/3))
+    A1 = Integer((t - c)/2 - 1)
+    B1 = c
+
+    A2 = Integer((t + c)/2 - 1)
+    B2 = Integer(1 - (t - c)/2)
+
+    # We use a negated b values in secp256k1_scalar_split_lambda.
+    B1, B2 = -B1, -B2
+
+    return A1, B1, A2, B2
+
+A1, B1, A2, B2 = find_split_constants_explicit_tof()
+
+# For extra fun, use an independent method to recompute the constants.
+assert (A1, B1, A2, B2) == find_split_constants_gauss()
+
+# PHI : Z[l] -> Z_n where phi(a + b*l) == a + b*lambda mod n.
+def PHI(a,b):
+    return Z(a + LAMBDA*b)
+
+# Check that (A1, B1) and (A2, B2) are in the kernel of PHI.
+assert PHI(A1, B1) == Z(0)
+assert PHI(A2, B2) == Z(0)
+
+# Check that the parallelogram generated by (A1, A2) and (B1, B2)
+# is a fundamental domain by containing exactly N points.
+# Since the LHS is the determinant and N != 0, this also checks that
+# (A1, A2) and (B1, B2) are linearly independent. By the previous
+# assertions, (A1, A2) and (B1, B2) are a basis of the kernel.
+assert A1*B2 - B1*A2 == N
+
+# Check that their components are short enough.
+assert (A1 + A2)/2 < sqrt(N)
+assert B1 < sqrt(N)
+assert B2 < sqrt(N)
+
+G1 = round((2**384)*B2/N)
+G2 = round((2**384)*(-B1)/N)
+
+def rnddiv2(v):
+    if v & 1:
+        v += 1
+    return v >> 1
+
+def scalar_lambda_split(k):
+    """Equivalent to secp256k1_scalar_lambda_split()."""
+    c1 = rnddiv2((k * G1) >> 383)
+    c2 = rnddiv2((k * G2) >> 383)
+    c1 = (c1 * -B1) % N
+    c2 = (c2 * -B2) % N
+    r2 = (c1 + c2) % N
+    r1 = (k + r2 * -LAMBDA) % N
+    return (r1, r2)
+
+# The result of scalar_lambda_split can depend on the representation of k (mod n).
+SPECIAL = (2**383) // G2 + 1
+assert scalar_lambda_split(SPECIAL) != scalar_lambda_split(SPECIAL + N)
+
+print('  A1     =', hex(A1))
+print(' -B1     =', hex(-B1))
+print('  A2     =', hex(A2))
+print(' -B2     =', hex(-B2))
+print('         =', hex(Z(-B2)))
+print(' -LAMBDA =', hex(-LAMBDA))
+
+print('  G1     =', hex(G1))
+print('  G2     =', hex(G2))
diff --git a/src/secp256k1/sage/group_prover.sage b/src/secp256k1/sage/group_prover.sage
index 8521f07999..b200bfeae3 100644
--- a/src/secp256k1/sage/group_prover.sage
+++ b/src/secp256k1/sage/group_prover.sage
@@ -42,7 +42,7 @@
 #     as we assume that all constraints in it are complementary with each other.
 #
 # Based on the sage verification scripts used in the Explicit-Formulas Database
-# by Tanja Lange and others, see http://hyperelliptic.org/EFD
+# by Tanja Lange and others, see https://hyperelliptic.org/EFD
 
 class fastfrac:
   """Fractions over rings."""
@@ -65,7 +65,7 @@ class fastfrac:
     return self.top in I and self.bot not in I
 
   def reduce(self,assumeZero):
-    zero = self.R.ideal(map(numerator, assumeZero))
+    zero = self.R.ideal(list(map(numerator, assumeZero)))
     return fastfrac(self.R, zero.reduce(self.top)) / fastfrac(self.R, zero.reduce(self.bot))
 
   def __add__(self,other):
@@ -100,7 +100,7 @@ class fastfrac:
     """Multiply something else with a fraction."""
     return self.__mul__(other)
 
-  def __div__(self,other):
+  def __truediv__(self,other):
     """Divide two fractions."""
     if parent(other) == ZZ:
       return fastfrac(self.R,self.top,self.bot * other)
@@ -108,6 +108,11 @@ class fastfrac:
       return fastfrac(self.R,self.top * other.bot,self.bot * other.top)
     return NotImplemented
 
+  # Compatibility wrapper for Sage versions based on Python 2
+  def __div__(self,other):
+     """Divide two fractions."""
+     return self.__truediv__(other)
+
   def __pow__(self,other):
     """Compute a power of a fraction."""
     if parent(other) == ZZ:
@@ -175,7 +180,7 @@ class constraints:
 
 def conflicts(R, con):
   """Check whether any of the passed non-zero assumptions is implied by the zero assumptions"""
-  zero = R.ideal(map(numerator, con.zero))
+  zero = R.ideal(list(map(numerator, con.zero)))
   if 1 in zero:
     return True
   # First a cheap check whether any of the individual nonzero terms conflict on
@@ -195,7 +200,7 @@ def conflicts(R, con):
 
 def get_nonzero_set(R, assume):
   """Calculate a simple set of nonzero expressions"""
-  zero = R.ideal(map(numerator, assume.zero))
+  zero = R.ideal(list(map(numerator, assume.zero)))
   nonzero = set()
   for nz in map(numerator, assume.nonzero):
     for (f,n) in nz.factor():
@@ -208,7 +213,7 @@ def get_nonzero_set(R, assume):
 
 def prove_nonzero(R, exprs, assume):
   """Check whether an expression is provably nonzero, given assumptions"""
-  zero = R.ideal(map(numerator, assume.zero))
+  zero = R.ideal(list(map(numerator, assume.zero)))
   nonzero = get_nonzero_set(R, assume)
   expl = set()
   ok = True
@@ -250,7 +255,7 @@ def prove_zero(R, exprs, assume):
   r, e = prove_nonzero(R, dict(map(lambda x: (fastfrac(R, x.bot, 1), exprs[x]), exprs)), assume)
   if not r:
     return (False, map(lambda x: "Possibly zero denominator: %s" % x, e))
-  zero = R.ideal(map(numerator, assume.zero))
+  zero = R.ideal(list(map(numerator, assume.zero)))
   nonzero = prod(x for x in assume.nonzero)
   expl = []
   for expr in exprs:
@@ -265,8 +270,8 @@ def describe_extra(R, assume, assumeExtra):
   """Describe what assumptions are added, given existing assumptions"""
   zerox = assume.zero.copy()
   zerox.update(assumeExtra.zero)
-  zero = R.ideal(map(numerator, assume.zero))
-  zeroextra = R.ideal(map(numerator, zerox))
+  zero = R.ideal(list(map(numerator, assume.zero)))
+  zeroextra = R.ideal(list(map(numerator, zerox)))
   nonzero = get_nonzero_set(R, assume)
   ret = set()
   # Iterate over the extra zero expressions
diff --git a/src/secp256k1/sage/secp256k1.sage b/src/secp256k1/sage/prove_group_implementations.sage
index a97e732f7f..a97e732f7f 100644
--- a/src/secp256k1/sage/secp256k1.sage
+++ b/src/secp256k1/sage/prove_group_implementations.sage
diff --git a/src/secp256k1/sage/secp256k1_params.sage b/src/secp256k1/sage/secp256k1_params.sage
new file mode 100644
index 0000000000..4e000726ed
--- /dev/null
+++ b/src/secp256k1/sage/secp256k1_params.sage
@@ -0,0 +1,36 @@
+"""Prime order of finite field underlying secp256k1 (2^256 - 2^32 - 977)"""
+P = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F
+
+"""Finite field underlying secp256k1"""
+F = FiniteField(P)
+
+"""Elliptic curve secp256k1: y^2 = x^3 + 7"""
+C = EllipticCurve([F(0), F(7)])
+
+"""Base point of secp256k1"""
+G = C.lift_x(0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798)
+
+"""Prime order of secp256k1"""
+N = C.order()
+
+"""Finite field of scalars of secp256k1"""
+Z = FiniteField(N)
+
+""" Beta value of secp256k1 non-trivial endomorphism: lambda * (x, y) = (beta * x, y)"""
+BETA = F(2)^((P-1)/3)
+
+""" Lambda value of secp256k1 non-trivial endomorphism: lambda * (x, y) = (beta * x, y)"""
+LAMBDA = Z(3)^((N-1)/3)
+
+assert is_prime(P)
+assert is_prime(N)
+
+assert BETA != F(1)
+assert BETA^3 == F(1)
+assert BETA^2 + BETA + 1 == 0
+
+assert LAMBDA != Z(1)
+assert LAMBDA^3 == Z(1)
+assert LAMBDA^2 + LAMBDA + 1 == 0
+
+assert Integer(LAMBDA)*G == C(BETA*G[0], G[1])
diff --git a/src/secp256k1/sage/weierstrass_prover.sage b/src/secp256k1/sage/weierstrass_prover.sage
index 03ef2ec901..b770c6dafe 100644
--- a/src/secp256k1/sage/weierstrass_prover.sage
+++ b/src/secp256k1/sage/weierstrass_prover.sage
@@ -175,24 +175,24 @@ laws_jacobian_weierstrass = {
 def check_exhaustive_jacobian_weierstrass(name, A, B, branches, formula, p):
   """Verify an implementation of addition of Jacobian points on a Weierstrass curve, by executing and validating the result for every possible addition in a prime field"""
   F = Integers(p)
-  print "Formula %s on Z%i:" % (name, p)
+  print("Formula %s on Z%i:" % (name, p))
   points = []
-  for x in xrange(0, p):
-    for y in xrange(0, p):
+  for x in range(0, p):
+    for y in range(0, p):
       point = affinepoint(F(x), F(y))
       r, e = concrete_verify(on_weierstrass_curve(A, B, point))
       if r:
         points.append(point)
 
-  for za in xrange(1, p):
-    for zb in xrange(1, p):
+  for za in range(1, p):
+    for zb in range(1, p):
       for pa in points:
         for pb in points:
-          for ia in xrange(2):
-            for ib in xrange(2):
+          for ia in range(2):
+            for ib in range(2):
               pA = jacobianpoint(pa.x * F(za)^2, pa.y * F(za)^3, F(za), ia)
               pB = jacobianpoint(pb.x * F(zb)^2, pb.y * F(zb)^3, F(zb), ib)
-              for branch in xrange(0, branches):
+              for branch in range(0, branches):
                 assumeAssert, assumeBranch, pC = formula(branch, pA, pB)
                 pC.X = F(pC.X)
                 pC.Y = F(pC.Y)
@@ -206,13 +206,13 @@ def check_exhaustive_jacobian_weierstrass(name, A, B, branches, formula, p):
                     r, e = concrete_verify(assumeLaw)
                     if r:
                       if match:
-                        print "  multiple branches for (%s,%s,%s,%s) + (%s,%s,%s,%s)" % (pA.X, pA.Y, pA.Z, pA.Infinity, pB.X, pB.Y, pB.Z, pB.Infinity)
+                        print("  multiple branches for (%s,%s,%s,%s) + (%s,%s,%s,%s)" % (pA.X, pA.Y, pA.Z, pA.Infinity, pB.X, pB.Y, pB.Z, pB.Infinity))
                       else:
                         match = True
                       r, e = concrete_verify(require)
                       if not r:
-                        print "  failure in branch %i for (%s,%s,%s,%s) + (%s,%s,%s,%s) = (%s,%s,%s,%s): %s" % (branch, pA.X, pA.Y, pA.Z, pA.Infinity, pB.X, pB.Y, pB.Z, pB.Infinity, pC.X, pC.Y, pC.Z, pC.Infinity, e)
-  print
+                        print("  failure in branch %i for (%s,%s,%s,%s) + (%s,%s,%s,%s) = (%s,%s,%s,%s): %s" % (branch, pA.X, pA.Y, pA.Z, pA.Infinity, pB.X, pB.Y, pB.Z, pB.Infinity, pC.X, pC.Y, pC.Z, pC.Infinity, e))
+  print()
 
 
 def check_symbolic_function(R, assumeAssert, assumeBranch, f, A, B, pa, pb, pA, pB, pC):
@@ -242,9 +242,9 @@ def check_symbolic_jacobian_weierstrass(name, A, B, branches, formula):
   for key in laws_jacobian_weierstrass:
     res[key] = []
 
-  print ("Formula " + name + ":")
+  print("Formula " + name + ":")
   count = 0
-  for branch in xrange(branches):
+  for branch in range(branches):
     assumeFormula, assumeBranch, pC = formula(branch, pA, pB)
     pC.X = lift(pC.X)
     pC.Y = lift(pC.Y)
@@ -255,10 +255,10 @@ def check_symbolic_jacobian_weierstrass(name, A, B, branches, formula):
       res[key].append((check_symbolic_function(R, assumeFormula, assumeBranch, laws_jacobian_weierstrass[key], A, B, pa, pb, pA, pB, pC), branch))
 
   for key in res:
-    print "  %s:" % key
+    print("  %s:" % key)
     val = res[key]
     for x in val:
       if x[0] is not None:
-        print "    branch %i: %s" % (x[1], x[0])
+        print("    branch %i: %s" % (x[1], x[0]))
 
-  print
+  print()