mirror of
https://github.com/bitcoin/bitcoin.git
synced 2024-11-20 02:25:40 +01:00
Merge pull request #302
03d4611 Add sage verification script for the group laws (Pieter Wuille)This commit is contained in:
Pieter Wuille
commit
3026daa095
322
sage/group_prover.sage
Normal file
322
sage/group_prover.sage
Normal file
@ -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
Normal file
306
sage/secp256k1.sage
Normal file
@ -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
Normal file
264
sage/weierstrass_prover.sage
Normal file
@ -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…
Reference in New Issue
Block a user