aboutsummaryrefslogtreecommitdiff
path: root/src/secp256k1/sage
diff options
context:
space:
mode:
Diffstat (limited to 'src/secp256k1/sage')
-rw-r--r--src/secp256k1/sage/gen_exhaustive_groups.sage129
-rw-r--r--src/secp256k1/sage/group_prover.sage322
-rw-r--r--src/secp256k1/sage/secp256k1.sage306
-rw-r--r--src/secp256k1/sage/weierstrass_prover.sage264
4 files changed, 1021 insertions, 0 deletions
diff --git a/src/secp256k1/sage/gen_exhaustive_groups.sage b/src/secp256k1/sage/gen_exhaustive_groups.sage
new file mode 100644
index 0000000000..3c3c984811
--- /dev/null
+++ b/src/secp256k1/sage/gen_exhaustive_groups.sage
@@ -0,0 +1,129 @@
+# 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))
+
+orders_done = set()
+results = {}
+first = True
+for b in range(1, P):
+ # There are only 6 curves (up to isomorphism) of the form y^2=x^3+B. Stop once we have tried all.
+ if len(orders_done) == 6:
+ break
+
+ E = EllipticCurve(F, [0, b])
+ print("Analyzing curve y^2 = x^3 + %i" % b)
+ n = E.order()
+ # Skip curves with an order we've already tried
+ if n in orders_done:
+ print("- Isomorphic to earlier curve")
+ continue
+ orders_done.add(n)
+ # Skip curves isomorphic to the real secp256k1
+ if n.is_pseudoprime():
+ print(" - Isomorphic to secp256k1")
+ continue
+
+ print("- Finding subgroups")
+
+ # Find what prime subgroups exist
+ for f, _ in n.factor():
+ print("- Analyzing subgroup of order %i" % f)
+ # Skip subgroups of order >1000
+ if f < 4 or f > 1000:
+ print(" - Bad size")
+ continue
+
+ # Iterate over X coordinates until we find one that is on the curve, has order f,
+ # and for which curve isomorphism exists that maps it to X coordinate 1.
+ for x in range(1, P):
+ # Skip X coordinates not on the curve, and construct the full point otherwise.
+ if not E.is_x_coord(x):
+ continue
+ G = E.lift_x(F(x))
+
+ print(" - Analyzing (multiples of) point with X=%i" % x)
+
+ # Skip points whose order is not a multiple of f. Project the point to have
+ # order f otherwise.
+ if (G.order() % f):
+ print(" - Bad order")
+ continue
+ G = G * (G.order() // f)
+
+ # Find lambda for endomorphism. Skip if none can be found.
+ lam = None
+ for l in Integers(f)(1).nth_root(3, all=True):
+ if int(l)*G == E(BETA*G[0], G[1]):
+ lam = int(l)
+ break
+ if lam is None:
+ print(" - No endomorphism for this subgroup")
+ break
+
+ # Now look for an isomorphism of the curve that gives this point an X
+ # coordinate equal to 1.
+ # If (x,y) is on y^2 = x^3 + b, then (a^2*x, a^3*y) is on y^2 = x^3 + a^6*b.
+ # So look for m=a^2=1/x.
+ m = F(1)/G[0]
+ if not m.is_square():
+ print(" - No curve isomorphism maps it to a point with X=1")
+ continue
+ a = m.sqrt()
+ rb = a^6*b
+ RE = EllipticCurve(F, [0, rb])
+
+ # Use as generator twice the image of G under the above isormorphism.
+ # This means that generator*(1/2 mod f) will have X coordinate 1.
+ RG = RE(1, a^3*G[1]) * 2
+ # And even Y coordinate.
+ if int(RG[1]) % 2:
+ RG = -RG
+ assert(RG.order() == f)
+ assert(lam*RG == RE(BETA*RG[0], RG[1]))
+
+ # We have found curve RE:y^2=x^3+rb with generator RG of order f. Remember it
+ results[f] = {"b": rb, "G": RG, "lambda": lam}
+ print(" - Found solution")
+ break
+
+ print("")
+
+print("")
+print("")
+print("/* To be put in src/group_impl.h: */")
+first = True
+for f in sorted(results.keys()):
+ b = results[f]["b"]
+ G = results[f]["G"]
+ print("# %s EXHAUSTIVE_TEST_ORDER == %i" % ("if" if first else "elif", f))
+ first = False
+ print("static const secp256k1_ge secp256k1_ge_const_g = SECP256K1_GE_CONST(")
+ print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x," % tuple((int(G[0]) >> (32 * (7 - i))) & 0xffffffff for i in range(4)))
+ print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x," % tuple((int(G[0]) >> (32 * (7 - i))) & 0xffffffff for i in range(4, 8)))
+ print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x," % tuple((int(G[1]) >> (32 * (7 - i))) & 0xffffffff for i in range(4)))
+ print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x" % tuple((int(G[1]) >> (32 * (7 - i))) & 0xffffffff for i in range(4, 8)))
+ print(");")
+ print("static const secp256k1_fe secp256k1_fe_const_b = SECP256K1_FE_CONST(")
+ print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x," % tuple((int(b) >> (32 * (7 - i))) & 0xffffffff for i in range(4)))
+ print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x" % tuple((int(b) >> (32 * (7 - i))) & 0xffffffff for i in range(4, 8)))
+ print(");")
+print("# else")
+print("# error No known generator for the specified exhaustive test group order.")
+print("# endif")
+
+print("")
+print("")
+print("/* To be put in src/scalar_impl.h: */")
+first = True
+for f in sorted(results.keys()):
+ lam = results[f]["lambda"]
+ print("# %s EXHAUSTIVE_TEST_ORDER == %i" % ("if" if first else "elif", f))
+ first = False
+ print("# define EXHAUSTIVE_TEST_LAMBDA %i" % lam)
+print("# else")
+print("# error No known lambda for the specified exhaustive test group order.")
+print("# endif")
+print("")
diff --git a/src/secp256k1/sage/group_prover.sage b/src/secp256k1/sage/group_prover.sage
new file mode 100644
index 0000000000..8521f07999
--- /dev/null
+++ b/src/secp256k1/sage/group_prover.sage
@@ -0,0 +1,322 @@
+# This code supports verifying group implementations which have branches
+# or conditional statements (like cmovs), by allowing each execution path
+# to independently set assumptions on input or intermediary variables.
+#
+# The general approach is:
+# * A constraint is a tuple of two sets of symbolic expressions:
+# the first of which are required to evaluate to zero, the second of which
+# are required to evaluate to nonzero.
+# - A constraint is said to be conflicting if any of its nonzero expressions
+# is in the ideal with basis the zero expressions (in other words: when the
+# zero expressions imply that one of the nonzero expressions are zero).
+# * There is a list of laws that describe the intended behaviour, including
+# laws for addition and doubling. Each law is called with the symbolic point
+# coordinates as arguments, and returns:
+# - A constraint describing the assumptions under which it is applicable,
+# called "assumeLaw"
+# - A constraint describing the requirements of the law, called "require"
+# * Implementations are transliterated into functions that operate as well on
+# algebraic input points, and are called once per combination of branches
+# executed. Each execution returns:
+# - A constraint describing the assumptions this implementation requires
+# (such as Z1=1), called "assumeFormula"
+# - A constraint describing the assumptions this specific branch requires,
+# but which is by construction guaranteed to cover the entire space by
+# merging the results from all branches, called "assumeBranch"
+# - The result of the computation
+# * All combinations of laws with implementation branches are tried, and:
+# - If the combination of assumeLaw, assumeFormula, and assumeBranch results
+# in a conflict, it means this law does not apply to this branch, and it is
+# skipped.
+# - For others, we try to prove the require constraints hold, assuming the
+# information in assumeLaw + assumeFormula + assumeBranch, and if this does
+# not succeed, we fail.
+# + To prove an expression is zero, we check whether it belongs to the
+# ideal with the assumed zero expressions as basis. This test is exact.
+# + To prove an expression is nonzero, we check whether each of its
+# factors is contained in the set of nonzero assumptions' factors.
+# This test is not exact, so various combinations of original and
+# reduced expressions' factors are tried.
+# - If we succeed, we print out the assumptions from assumeFormula that
+# weren't implied by assumeLaw already. Those from assumeBranch are skipped,
+# 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
+
+class fastfrac:
+ """Fractions over rings."""
+
+ def __init__(self,R,top,bot=1):
+ """Construct a fractional, given a ring, a numerator, and denominator."""
+ self.R = R
+ if parent(top) == ZZ or parent(top) == R:
+ self.top = R(top)
+ self.bot = R(bot)
+ elif top.__class__ == fastfrac:
+ self.top = top.top
+ self.bot = top.bot * bot
+ else:
+ self.top = R(numerator(top))
+ self.bot = R(denominator(top)) * bot
+
+ def iszero(self,I):
+ """Return whether this fraction is zero given an ideal."""
+ return self.top in I and self.bot not in I
+
+ def reduce(self,assumeZero):
+ zero = self.R.ideal(map(numerator, assumeZero))
+ return fastfrac(self.R, zero.reduce(self.top)) / fastfrac(self.R, zero.reduce(self.bot))
+
+ def __add__(self,other):
+ """Add two fractions."""
+ if parent(other) == ZZ:
+ return fastfrac(self.R,self.top + self.bot * other,self.bot)
+ if other.__class__ == fastfrac:
+ return fastfrac(self.R,self.top * other.bot + self.bot * other.top,self.bot * other.bot)
+ return NotImplemented
+
+ def __sub__(self,other):
+ """Subtract two fractions."""
+ if parent(other) == ZZ:
+ return fastfrac(self.R,self.top - self.bot * other,self.bot)
+ if other.__class__ == fastfrac:
+ return fastfrac(self.R,self.top * other.bot - self.bot * other.top,self.bot * other.bot)
+ return NotImplemented
+
+ def __neg__(self):
+ """Return the negation of a fraction."""
+ return fastfrac(self.R,-self.top,self.bot)
+
+ def __mul__(self,other):
+ """Multiply two fractions."""
+ if parent(other) == ZZ:
+ return fastfrac(self.R,self.top * other,self.bot)
+ if other.__class__ == fastfrac:
+ return fastfrac(self.R,self.top * other.top,self.bot * other.bot)
+ return NotImplemented
+
+ def __rmul__(self,other):
+ """Multiply something else with a fraction."""
+ return self.__mul__(other)
+
+ def __div__(self,other):
+ """Divide two fractions."""
+ if parent(other) == ZZ:
+ return fastfrac(self.R,self.top,self.bot * other)
+ if other.__class__ == fastfrac:
+ return fastfrac(self.R,self.top * other.bot,self.bot * other.top)
+ return NotImplemented
+
+ def __pow__(self,other):
+ """Compute a power of a fraction."""
+ if parent(other) == ZZ:
+ if other < 0:
+ # Negative powers require flipping top and bottom
+ return fastfrac(self.R,self.bot ^ (-other),self.top ^ (-other))
+ else:
+ return fastfrac(self.R,self.top ^ other,self.bot ^ other)
+ return NotImplemented
+
+ def __str__(self):
+ return "fastfrac((" + str(self.top) + ") / (" + str(self.bot) + "))"
+ def __repr__(self):
+ return "%s" % self
+
+ def numerator(self):
+ return self.top
+
+class constraints:
+ """A set of constraints, consisting of zero and nonzero expressions.
+
+ Constraints can either be used to express knowledge or a requirement.
+
+ Both the fields zero and nonzero are maps from expressions to description
+ strings. The expressions that are the keys in zero are required to be zero,
+ and the expressions that are the keys in nonzero are required to be nonzero.
+
+ Note that (a != 0) and (b != 0) is the same as (a*b != 0), so all keys in
+ nonzero could be multiplied into a single key. This is often much less
+ efficient to work with though, so we keep them separate inside the
+ constraints. This allows higher-level code to do fast checks on the individual
+ nonzero elements, or combine them if needed for stronger checks.
+
+ We can't multiply the different zero elements, as it would suffice for one of
+ the factors to be zero, instead of all of them. Instead, the zero elements are
+ typically combined into an ideal first.
+ """
+
+ def __init__(self, **kwargs):
+ if 'zero' in kwargs:
+ self.zero = dict(kwargs['zero'])
+ else:
+ self.zero = dict()
+ if 'nonzero' in kwargs:
+ self.nonzero = dict(kwargs['nonzero'])
+ else:
+ self.nonzero = dict()
+
+ def negate(self):
+ return constraints(zero=self.nonzero, nonzero=self.zero)
+
+ def __add__(self, other):
+ zero = self.zero.copy()
+ zero.update(other.zero)
+ nonzero = self.nonzero.copy()
+ nonzero.update(other.nonzero)
+ return constraints(zero=zero, nonzero=nonzero)
+
+ def __str__(self):
+ return "constraints(zero=%s,nonzero=%s)" % (self.zero, self.nonzero)
+
+ def __repr__(self):
+ return "%s" % self
+
+
+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))
+ if 1 in zero:
+ return True
+ # First a cheap check whether any of the individual nonzero terms conflict on
+ # their own.
+ for nonzero in con.nonzero:
+ if nonzero.iszero(zero):
+ return True
+ # It can be the case that entries in the nonzero set do not individually
+ # conflict with the zero set, but their combination does. For example, knowing
+ # that either x or y is zero is equivalent to having x*y in the zero set.
+ # Having x or y individually in the nonzero set is not a conflict, but both
+ # simultaneously is, so that is the right thing to check for.
+ if reduce(lambda a,b: a * b, con.nonzero, fastfrac(R, 1)).iszero(zero):
+ return True
+ return False
+
+
+def get_nonzero_set(R, assume):
+ """Calculate a simple set of nonzero expressions"""
+ zero = R.ideal(map(numerator, assume.zero))
+ nonzero = set()
+ for nz in map(numerator, assume.nonzero):
+ for (f,n) in nz.factor():
+ nonzero.add(f)
+ rnz = zero.reduce(nz)
+ for (f,n) in rnz.factor():
+ nonzero.add(f)
+ return nonzero
+
+
+def prove_nonzero(R, exprs, assume):
+ """Check whether an expression is provably nonzero, given assumptions"""
+ zero = R.ideal(map(numerator, assume.zero))
+ nonzero = get_nonzero_set(R, assume)
+ expl = set()
+ ok = True
+ for expr in exprs:
+ if numerator(expr) in zero:
+ return (False, [exprs[expr]])
+ allexprs = reduce(lambda a,b: numerator(a)*numerator(b), exprs, 1)
+ for (f, n) in allexprs.factor():
+ if f not in nonzero:
+ ok = False
+ if ok:
+ return (True, None)
+ ok = True
+ for (f, n) in zero.reduce(numerator(allexprs)).factor():
+ if f not in nonzero:
+ ok = False
+ if ok:
+ return (True, None)
+ ok = True
+ for expr in exprs:
+ for (f,n) in numerator(expr).factor():
+ if f not in nonzero:
+ ok = False
+ if ok:
+ return (True, None)
+ ok = True
+ for expr in exprs:
+ for (f,n) in zero.reduce(numerator(expr)).factor():
+ if f not in nonzero:
+ expl.add(exprs[expr])
+ if expl:
+ return (False, list(expl))
+ else:
+ return (True, None)
+
+
+def prove_zero(R, exprs, assume):
+ """Check whether all of the passed expressions are provably zero, given assumptions"""
+ 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))
+ nonzero = prod(x for x in assume.nonzero)
+ expl = []
+ for expr in exprs:
+ if not expr.iszero(zero):
+ expl.append(exprs[expr])
+ if not expl:
+ return (True, None)
+ return (False, expl)
+
+
+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))
+ nonzero = get_nonzero_set(R, assume)
+ ret = set()
+ # Iterate over the extra zero expressions
+ for base in assumeExtra.zero:
+ if base not in zero:
+ add = []
+ for (f, n) in numerator(base).factor():
+ if f not in nonzero:
+ add += ["%s" % f]
+ if add:
+ ret.add((" * ".join(add)) + " = 0 [%s]" % assumeExtra.zero[base])
+ # Iterate over the extra nonzero expressions
+ for nz in assumeExtra.nonzero:
+ nzr = zeroextra.reduce(numerator(nz))
+ if nzr not in zeroextra:
+ for (f,n) in nzr.factor():
+ if zeroextra.reduce(f) not in nonzero:
+ ret.add("%s != 0" % zeroextra.reduce(f))
+ return ", ".join(x for x in ret)
+
+
+def check_symbolic(R, assumeLaw, assumeAssert, assumeBranch, require):
+ """Check a set of zero and nonzero requirements, given a set of zero and nonzero assumptions"""
+ assume = assumeLaw + assumeAssert + assumeBranch
+
+ if conflicts(R, assume):
+ # This formula does not apply
+ return None
+
+ describe = describe_extra(R, assumeLaw + assumeBranch, assumeAssert)
+
+ ok, msg = prove_zero(R, require.zero, assume)
+ if not ok:
+ return "FAIL, %s fails (assuming %s)" % (str(msg), describe)
+
+ res, expl = prove_nonzero(R, require.nonzero, assume)
+ if not res:
+ return "FAIL, %s fails (assuming %s)" % (str(expl), describe)
+
+ if describe != "":
+ return "OK (assuming %s)" % describe
+ else:
+ return "OK"
+
+
+def concrete_verify(c):
+ for k in c.zero:
+ if k != 0:
+ return (False, c.zero[k])
+ for k in c.nonzero:
+ if k == 0:
+ return (False, c.nonzero[k])
+ return (True, None)
diff --git a/src/secp256k1/sage/secp256k1.sage b/src/secp256k1/sage/secp256k1.sage
new file mode 100644
index 0000000000..a97e732f7f
--- /dev/null
+++ b/src/secp256k1/sage/secp256k1.sage
@@ -0,0 +1,306 @@
+# Test libsecp256k1' group operation implementations using prover.sage
+
+import sys
+
+load("group_prover.sage")
+load("weierstrass_prover.sage")
+
+def formula_secp256k1_gej_double_var(a):
+ """libsecp256k1's secp256k1_gej_double_var, used by various addition functions"""
+ rz = a.Z * a.Y
+ rz = rz * 2
+ t1 = a.X^2
+ t1 = t1 * 3
+ t2 = t1^2
+ t3 = a.Y^2
+ t3 = t3 * 2
+ t4 = t3^2
+ t4 = t4 * 2
+ t3 = t3 * a.X
+ rx = t3
+ rx = rx * 4
+ rx = -rx
+ rx = rx + t2
+ t2 = -t2
+ t3 = t3 * 6
+ t3 = t3 + t2
+ ry = t1 * t3
+ t2 = -t4
+ ry = ry + t2
+ return jacobianpoint(rx, ry, rz)
+
+def formula_secp256k1_gej_add_var(branch, a, b):
+ """libsecp256k1's secp256k1_gej_add_var"""
+ if branch == 0:
+ return (constraints(), constraints(nonzero={a.Infinity : 'a_infinite'}), b)
+ if branch == 1:
+ return (constraints(), constraints(zero={a.Infinity : 'a_finite'}, nonzero={b.Infinity : 'b_infinite'}), a)
+ z22 = b.Z^2
+ z12 = a.Z^2
+ u1 = a.X * z22
+ u2 = b.X * z12
+ s1 = a.Y * z22
+ s1 = s1 * b.Z
+ s2 = b.Y * z12
+ s2 = s2 * a.Z
+ h = -u1
+ h = h + u2
+ i = -s1
+ i = i + s2
+ if branch == 2:
+ r = formula_secp256k1_gej_double_var(a)
+ return (constraints(), constraints(zero={h : 'h=0', i : 'i=0', a.Infinity : 'a_finite', b.Infinity : 'b_finite'}), r)
+ if branch == 3:
+ return (constraints(), constraints(zero={h : 'h=0', a.Infinity : 'a_finite', b.Infinity : 'b_finite'}, nonzero={i : 'i!=0'}), point_at_infinity())
+ i2 = i^2
+ h2 = h^2
+ h3 = h2 * h
+ h = h * b.Z
+ rz = a.Z * h
+ t = u1 * h2
+ rx = t
+ rx = rx * 2
+ rx = rx + h3
+ rx = -rx
+ rx = rx + i2
+ ry = -rx
+ ry = ry + t
+ ry = ry * i
+ h3 = h3 * s1
+ h3 = -h3
+ ry = ry + h3
+ return (constraints(), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite'}, nonzero={h : 'h!=0'}), jacobianpoint(rx, ry, rz))
+
+def formula_secp256k1_gej_add_ge_var(branch, a, b):
+ """libsecp256k1's secp256k1_gej_add_ge_var, which assume bz==1"""
+ if branch == 0:
+ return (constraints(zero={b.Z - 1 : 'b.z=1'}), constraints(nonzero={a.Infinity : 'a_infinite'}), b)
+ if branch == 1:
+ return (constraints(zero={b.Z - 1 : 'b.z=1'}), constraints(zero={a.Infinity : 'a_finite'}, nonzero={b.Infinity : 'b_infinite'}), a)
+ z12 = a.Z^2
+ u1 = a.X
+ u2 = b.X * z12
+ s1 = a.Y
+ s2 = b.Y * z12
+ s2 = s2 * a.Z
+ h = -u1
+ h = h + u2
+ i = -s1
+ i = i + s2
+ if (branch == 2):
+ r = formula_secp256k1_gej_double_var(a)
+ return (constraints(zero={b.Z - 1 : 'b.z=1'}), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite', h : 'h=0', i : 'i=0'}), r)
+ if (branch == 3):
+ return (constraints(zero={b.Z - 1 : 'b.z=1'}), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite', h : 'h=0'}, nonzero={i : 'i!=0'}), point_at_infinity())
+ i2 = i^2
+ h2 = h^2
+ h3 = h * h2
+ rz = a.Z * h
+ t = u1 * h2
+ rx = t
+ rx = rx * 2
+ rx = rx + h3
+ rx = -rx
+ rx = rx + i2
+ ry = -rx
+ ry = ry + t
+ ry = ry * i
+ h3 = h3 * s1
+ h3 = -h3
+ ry = ry + h3
+ return (constraints(zero={b.Z - 1 : 'b.z=1'}), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite'}, nonzero={h : 'h!=0'}), jacobianpoint(rx, ry, rz))
+
+def formula_secp256k1_gej_add_zinv_var(branch, a, b):
+ """libsecp256k1's secp256k1_gej_add_zinv_var"""
+ bzinv = b.Z^(-1)
+ if branch == 0:
+ return (constraints(), constraints(nonzero={b.Infinity : 'b_infinite'}), a)
+ if branch == 1:
+ bzinv2 = bzinv^2
+ bzinv3 = bzinv2 * bzinv
+ rx = b.X * bzinv2
+ ry = b.Y * bzinv3
+ rz = 1
+ return (constraints(), constraints(zero={b.Infinity : 'b_finite'}, nonzero={a.Infinity : 'a_infinite'}), jacobianpoint(rx, ry, rz))
+ azz = a.Z * bzinv
+ z12 = azz^2
+ u1 = a.X
+ u2 = b.X * z12
+ s1 = a.Y
+ s2 = b.Y * z12
+ s2 = s2 * azz
+ h = -u1
+ h = h + u2
+ i = -s1
+ i = i + s2
+ if branch == 2:
+ r = formula_secp256k1_gej_double_var(a)
+ return (constraints(), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite', h : 'h=0', i : 'i=0'}), r)
+ if branch == 3:
+ return (constraints(), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite', h : 'h=0'}, nonzero={i : 'i!=0'}), point_at_infinity())
+ i2 = i^2
+ h2 = h^2
+ h3 = h * h2
+ rz = a.Z
+ rz = rz * h
+ t = u1 * h2
+ rx = t
+ rx = rx * 2
+ rx = rx + h3
+ rx = -rx
+ rx = rx + i2
+ ry = -rx
+ ry = ry + t
+ ry = ry * i
+ h3 = h3 * s1
+ h3 = -h3
+ ry = ry + h3
+ return (constraints(), constraints(zero={a.Infinity : 'a_finite', b.Infinity : 'b_finite'}, nonzero={h : 'h!=0'}), jacobianpoint(rx, ry, rz))
+
+def formula_secp256k1_gej_add_ge(branch, a, b):
+ """libsecp256k1's secp256k1_gej_add_ge"""
+ zeroes = {}
+ nonzeroes = {}
+ a_infinity = False
+ if (branch & 4) != 0:
+ nonzeroes.update({a.Infinity : 'a_infinite'})
+ a_infinity = True
+ else:
+ zeroes.update({a.Infinity : 'a_finite'})
+ zz = a.Z^2
+ u1 = a.X
+ u2 = b.X * zz
+ s1 = a.Y
+ s2 = b.Y * zz
+ s2 = s2 * a.Z
+ t = u1
+ t = t + u2
+ m = s1
+ m = m + s2
+ rr = t^2
+ m_alt = -u2
+ tt = u1 * m_alt
+ rr = rr + tt
+ degenerate = (branch & 3) == 3
+ if (branch & 1) != 0:
+ zeroes.update({m : 'm_zero'})
+ else:
+ nonzeroes.update({m : 'm_nonzero'})
+ if (branch & 2) != 0:
+ zeroes.update({rr : 'rr_zero'})
+ else:
+ nonzeroes.update({rr : 'rr_nonzero'})
+ rr_alt = s1
+ rr_alt = rr_alt * 2
+ m_alt = m_alt + u1
+ if not degenerate:
+ rr_alt = rr
+ m_alt = m
+ n = m_alt^2
+ q = n * t
+ n = n^2
+ if degenerate:
+ n = m
+ t = rr_alt^2
+ rz = a.Z * m_alt
+ infinity = False
+ if (branch & 8) != 0:
+ if not a_infinity:
+ infinity = True
+ zeroes.update({rz : 'r.z=0'})
+ else:
+ nonzeroes.update({rz : 'r.z!=0'})
+ rz = rz * 2
+ q = -q
+ t = t + q
+ rx = t
+ t = t * 2
+ t = t + q
+ t = t * rr_alt
+ t = t + n
+ ry = -t
+ rx = rx * 4
+ ry = ry * 4
+ if a_infinity:
+ rx = b.X
+ ry = b.Y
+ rz = 1
+ if infinity:
+ return (constraints(zero={b.Z - 1 : 'b.z=1', b.Infinity : 'b_finite'}), constraints(zero=zeroes, nonzero=nonzeroes), point_at_infinity())
+ return (constraints(zero={b.Z - 1 : 'b.z=1', b.Infinity : 'b_finite'}), constraints(zero=zeroes, nonzero=nonzeroes), jacobianpoint(rx, ry, rz))
+
+def formula_secp256k1_gej_add_ge_old(branch, a, b):
+ """libsecp256k1's old secp256k1_gej_add_ge, which fails when ay+by=0 but ax!=bx"""
+ a_infinity = (branch & 1) != 0
+ zero = {}
+ nonzero = {}
+ if a_infinity:
+ nonzero.update({a.Infinity : 'a_infinite'})
+ else:
+ zero.update({a.Infinity : 'a_finite'})
+ zz = a.Z^2
+ u1 = a.X
+ u2 = b.X * zz
+ s1 = a.Y
+ s2 = b.Y * zz
+ s2 = s2 * a.Z
+ z = a.Z
+ t = u1
+ t = t + u2
+ m = s1
+ m = m + s2
+ n = m^2
+ q = n * t
+ n = n^2
+ rr = t^2
+ t = u1 * u2
+ t = -t
+ rr = rr + t
+ t = rr^2
+ rz = m * z
+ infinity = False
+ if (branch & 2) != 0:
+ if not a_infinity:
+ infinity = True
+ else:
+ return (constraints(zero={b.Z - 1 : 'b.z=1', b.Infinity : 'b_finite'}), constraints(nonzero={z : 'conflict_a'}, zero={z : 'conflict_b'}), point_at_infinity())
+ zero.update({rz : 'r.z=0'})
+ else:
+ nonzero.update({rz : 'r.z!=0'})
+ rz = rz * (0 if a_infinity else 2)
+ rx = t
+ q = -q
+ rx = rx + q
+ q = q * 3
+ t = t * 2
+ t = t + q
+ t = t * rr
+ t = t + n
+ ry = -t
+ rx = rx * (0 if a_infinity else 4)
+ ry = ry * (0 if a_infinity else 4)
+ t = b.X
+ t = t * (1 if a_infinity else 0)
+ rx = rx + t
+ t = b.Y
+ t = t * (1 if a_infinity else 0)
+ ry = ry + t
+ t = (1 if a_infinity else 0)
+ rz = rz + t
+ if infinity:
+ return (constraints(zero={b.Z - 1 : 'b.z=1', b.Infinity : 'b_finite'}), constraints(zero=zero, nonzero=nonzero), point_at_infinity())
+ return (constraints(zero={b.Z - 1 : 'b.z=1', b.Infinity : 'b_finite'}), constraints(zero=zero, nonzero=nonzero), jacobianpoint(rx, ry, rz))
+
+if __name__ == "__main__":
+ check_symbolic_jacobian_weierstrass("secp256k1_gej_add_var", 0, 7, 5, formula_secp256k1_gej_add_var)
+ check_symbolic_jacobian_weierstrass("secp256k1_gej_add_ge_var", 0, 7, 5, formula_secp256k1_gej_add_ge_var)
+ check_symbolic_jacobian_weierstrass("secp256k1_gej_add_zinv_var", 0, 7, 5, formula_secp256k1_gej_add_zinv_var)
+ check_symbolic_jacobian_weierstrass("secp256k1_gej_add_ge", 0, 7, 16, formula_secp256k1_gej_add_ge)
+ check_symbolic_jacobian_weierstrass("secp256k1_gej_add_ge_old [should fail]", 0, 7, 4, formula_secp256k1_gej_add_ge_old)
+
+ if len(sys.argv) >= 2 and sys.argv[1] == "--exhaustive":
+ check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_var", 0, 7, 5, formula_secp256k1_gej_add_var, 43)
+ check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_ge_var", 0, 7, 5, formula_secp256k1_gej_add_ge_var, 43)
+ check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_zinv_var", 0, 7, 5, formula_secp256k1_gej_add_zinv_var, 43)
+ check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_ge", 0, 7, 16, formula_secp256k1_gej_add_ge, 43)
+ check_exhaustive_jacobian_weierstrass("secp256k1_gej_add_ge_old [should fail]", 0, 7, 4, formula_secp256k1_gej_add_ge_old, 43)
diff --git a/src/secp256k1/sage/weierstrass_prover.sage b/src/secp256k1/sage/weierstrass_prover.sage
new file mode 100644
index 0000000000..03ef2ec901
--- /dev/null
+++ b/src/secp256k1/sage/weierstrass_prover.sage
@@ -0,0 +1,264 @@
+# Prover implementation for Weierstrass curves of the form
+# y^2 = x^3 + A * x + B, specifically with a = 0 and b = 7, with group laws
+# operating on affine and Jacobian coordinates, including the point at infinity
+# represented by a 4th variable in coordinates.
+
+load("group_prover.sage")
+
+
+class affinepoint:
+ def __init__(self, x, y, infinity=0):
+ self.x = x
+ self.y = y
+ self.infinity = infinity
+ def __str__(self):
+ return "affinepoint(x=%s,y=%s,inf=%s)" % (self.x, self.y, self.infinity)
+
+
+class jacobianpoint:
+ def __init__(self, x, y, z, infinity=0):
+ self.X = x
+ self.Y = y
+ self.Z = z
+ self.Infinity = infinity
+ def __str__(self):
+ return "jacobianpoint(X=%s,Y=%s,Z=%s,inf=%s)" % (self.X, self.Y, self.Z, self.Infinity)
+
+
+def point_at_infinity():
+ return jacobianpoint(1, 1, 1, 1)
+
+
+def negate(p):
+ if p.__class__ == affinepoint:
+ return affinepoint(p.x, -p.y)
+ if p.__class__ == jacobianpoint:
+ return jacobianpoint(p.X, -p.Y, p.Z)
+ assert(False)
+
+
+def on_weierstrass_curve(A, B, p):
+ """Return a set of zero-expressions for an affine point to be on the curve"""
+ return constraints(zero={p.x^3 + A*p.x + B - p.y^2: 'on_curve'})
+
+
+def tangential_to_weierstrass_curve(A, B, p12, p3):
+ """Return a set of zero-expressions for ((x12,y12),(x3,y3)) to be a line that is tangential to the curve at (x12,y12)"""
+ return constraints(zero={
+ (p12.y - p3.y) * (p12.y * 2) - (p12.x^2 * 3 + A) * (p12.x - p3.x): 'tangential_to_curve'
+ })
+
+
+def colinear(p1, p2, p3):
+ """Return a set of zero-expressions for ((x1,y1),(x2,y2),(x3,y3)) to be collinear"""
+ return constraints(zero={
+ (p1.y - p2.y) * (p1.x - p3.x) - (p1.y - p3.y) * (p1.x - p2.x): 'colinear_1',
+ (p2.y - p3.y) * (p2.x - p1.x) - (p2.y - p1.y) * (p2.x - p3.x): 'colinear_2',
+ (p3.y - p1.y) * (p3.x - p2.x) - (p3.y - p2.y) * (p3.x - p1.x): 'colinear_3'
+ })
+
+
+def good_affine_point(p):
+ return constraints(nonzero={p.x : 'nonzero_x', p.y : 'nonzero_y'})
+
+
+def good_jacobian_point(p):
+ return constraints(nonzero={p.X : 'nonzero_X', p.Y : 'nonzero_Y', p.Z^6 : 'nonzero_Z'})
+
+
+def good_point(p):
+ return constraints(nonzero={p.Z^6 : 'nonzero_X'})
+
+
+def finite(p, *affine_fns):
+ con = good_point(p) + constraints(zero={p.Infinity : 'finite_point'})
+ if p.Z != 0:
+ return con + reduce(lambda a, b: a + b, (f(affinepoint(p.X / p.Z^2, p.Y / p.Z^3)) for f in affine_fns), con)
+ else:
+ return con
+
+def infinite(p):
+ return constraints(nonzero={p.Infinity : 'infinite_point'})
+
+
+def law_jacobian_weierstrass_add(A, B, pa, pb, pA, pB, pC):
+ """Check whether the passed set of coordinates is a valid Jacobian add, given assumptions"""
+ assumeLaw = (good_affine_point(pa) +
+ good_affine_point(pb) +
+ good_jacobian_point(pA) +
+ good_jacobian_point(pB) +
+ on_weierstrass_curve(A, B, pa) +
+ on_weierstrass_curve(A, B, pb) +
+ finite(pA) +
+ finite(pB) +
+ constraints(nonzero={pa.x - pb.x : 'different_x'}))
+ require = (finite(pC, lambda pc: on_weierstrass_curve(A, B, pc) +
+ colinear(pa, pb, negate(pc))))
+ return (assumeLaw, require)
+
+
+def law_jacobian_weierstrass_double(A, B, pa, pb, pA, pB, pC):
+ """Check whether the passed set of coordinates is a valid Jacobian doubling, given assumptions"""
+ assumeLaw = (good_affine_point(pa) +
+ good_affine_point(pb) +
+ good_jacobian_point(pA) +
+ good_jacobian_point(pB) +
+ on_weierstrass_curve(A, B, pa) +
+ on_weierstrass_curve(A, B, pb) +
+ finite(pA) +
+ finite(pB) +
+ constraints(zero={pa.x - pb.x : 'equal_x', pa.y - pb.y : 'equal_y'}))
+ require = (finite(pC, lambda pc: on_weierstrass_curve(A, B, pc) +
+ tangential_to_weierstrass_curve(A, B, pa, negate(pc))))
+ return (assumeLaw, require)
+
+
+def law_jacobian_weierstrass_add_opposites(A, B, pa, pb, pA, pB, pC):
+ assumeLaw = (good_affine_point(pa) +
+ good_affine_point(pb) +
+ good_jacobian_point(pA) +
+ good_jacobian_point(pB) +
+ on_weierstrass_curve(A, B, pa) +
+ on_weierstrass_curve(A, B, pb) +
+ finite(pA) +
+ finite(pB) +
+ constraints(zero={pa.x - pb.x : 'equal_x', pa.y + pb.y : 'opposite_y'}))
+ require = infinite(pC)
+ return (assumeLaw, require)
+
+
+def law_jacobian_weierstrass_add_infinite_a(A, B, pa, pb, pA, pB, pC):
+ assumeLaw = (good_affine_point(pa) +
+ good_affine_point(pb) +
+ good_jacobian_point(pA) +
+ good_jacobian_point(pB) +
+ on_weierstrass_curve(A, B, pb) +
+ infinite(pA) +
+ finite(pB))
+ require = finite(pC, lambda pc: constraints(zero={pc.x - pb.x : 'c.x=b.x', pc.y - pb.y : 'c.y=b.y'}))
+ return (assumeLaw, require)
+
+
+def law_jacobian_weierstrass_add_infinite_b(A, B, pa, pb, pA, pB, pC):
+ assumeLaw = (good_affine_point(pa) +
+ good_affine_point(pb) +
+ good_jacobian_point(pA) +
+ good_jacobian_point(pB) +
+ on_weierstrass_curve(A, B, pa) +
+ infinite(pB) +
+ finite(pA))
+ require = finite(pC, lambda pc: constraints(zero={pc.x - pa.x : 'c.x=a.x', pc.y - pa.y : 'c.y=a.y'}))
+ return (assumeLaw, require)
+
+
+def law_jacobian_weierstrass_add_infinite_ab(A, B, pa, pb, pA, pB, pC):
+ assumeLaw = (good_affine_point(pa) +
+ good_affine_point(pb) +
+ good_jacobian_point(pA) +
+ good_jacobian_point(pB) +
+ infinite(pA) +
+ infinite(pB))
+ require = infinite(pC)
+ return (assumeLaw, require)
+
+
+laws_jacobian_weierstrass = {
+ 'add': law_jacobian_weierstrass_add,
+ 'double': law_jacobian_weierstrass_double,
+ 'add_opposite': law_jacobian_weierstrass_add_opposites,
+ 'add_infinite_a': law_jacobian_weierstrass_add_infinite_a,
+ 'add_infinite_b': law_jacobian_weierstrass_add_infinite_b,
+ 'add_infinite_ab': law_jacobian_weierstrass_add_infinite_ab
+}
+
+
+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)
+ points = []
+ for x in xrange(0, p):
+ for y in xrange(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 pa in points:
+ for pb in points:
+ for ia in xrange(2):
+ for ib in xrange(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):
+ assumeAssert, assumeBranch, pC = formula(branch, pA, pB)
+ pC.X = F(pC.X)
+ pC.Y = F(pC.Y)
+ pC.Z = F(pC.Z)
+ pC.Infinity = F(pC.Infinity)
+ r, e = concrete_verify(assumeAssert + assumeBranch)
+ if r:
+ match = False
+ for key in laws_jacobian_weierstrass:
+ assumeLaw, require = laws_jacobian_weierstrass[key](A, B, pa, pb, pA, pB, pC)
+ 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)
+ 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
+
+
+def check_symbolic_function(R, assumeAssert, assumeBranch, f, A, B, pa, pb, pA, pB, pC):
+ assumeLaw, require = f(A, B, pa, pb, pA, pB, pC)
+ return check_symbolic(R, assumeLaw, assumeAssert, assumeBranch, require)
+
+def check_symbolic_jacobian_weierstrass(name, A, B, branches, formula):
+ """Verify an implementation of addition of Jacobian points on a Weierstrass curve symbolically"""
+ R.<ax,bx,ay,by,Az,Bz,Ai,Bi> = PolynomialRing(QQ,8,order='invlex')
+ lift = lambda x: fastfrac(R,x)
+ ax = lift(ax)
+ ay = lift(ay)
+ Az = lift(Az)
+ bx = lift(bx)
+ by = lift(by)
+ Bz = lift(Bz)
+ Ai = lift(Ai)
+ Bi = lift(Bi)
+
+ pa = affinepoint(ax, ay, Ai)
+ pb = affinepoint(bx, by, Bi)
+ pA = jacobianpoint(ax * Az^2, ay * Az^3, Az, Ai)
+ pB = jacobianpoint(bx * Bz^2, by * Bz^3, Bz, Bi)
+
+ res = {}
+
+ for key in laws_jacobian_weierstrass:
+ res[key] = []
+
+ print ("Formula " + name + ":")
+ count = 0
+ for branch in xrange(branches):
+ assumeFormula, assumeBranch, pC = formula(branch, pA, pB)
+ pC.X = lift(pC.X)
+ pC.Y = lift(pC.Y)
+ pC.Z = lift(pC.Z)
+ pC.Infinity = lift(pC.Infinity)
+
+ for key in laws_jacobian_weierstrass:
+ 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
+ val = res[key]
+ for x in val:
+ if x[0] is not None:
+ print " branch %i: %s" % (x[1], x[0])
+
+ print