Browse Source

Add sage verification script for the group laws

master
Pieter Wuille 7 years ago
parent
commit
03d4611c81
  1. 322
      sage/group_prover.sage
  2. 306
      sage/secp256k1.sage
  3. 264
      sage/weierstrass_prover.sage

322
sage/group_prover.sage

@ -0,0 +1,322 @@ @@ -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 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
# exectured. 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)

306
sage/secp256k1.sage

@ -0,0 +1,306 @@ @@ -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)

264
sage/weierstrass_prover.sage

@ -0,0 +1,264 @@ @@ -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
Loading…
Cancel
Save