diff options
Diffstat (limited to 'src/secp256k1/sage')
| -rw-r--r-- | src/secp256k1/sage/gen_exhaustive_groups.sage | 129 | ||||
| -rw-r--r-- | src/secp256k1/sage/group_prover.sage | 322 | ||||
| -rw-r--r-- | src/secp256k1/sage/secp256k1.sage | 306 | ||||
| -rw-r--r-- | src/secp256k1/sage/weierstrass_prover.sage | 264 |
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 |