mirror of
https://github.com/mozilla/gecko-dev.git
synced 2024-11-05 08:35:26 +00:00
a7f29f4b91
Obtained from https://github.com/amintos/PyECC.git (commit 1bfd3a41410c3650e57d58f7bd016bb0819af250)
382 lines
12 KiB
Python
382 lines
12 KiB
Python
|
|
# --- ELLIPTIC CURVE MATH ------------------------------------------------------
|
|
#
|
|
# curve definition: y^2 = x^3 - p*x - q
|
|
# over finite field: Z/nZ* (prime residue classes modulo a prime number n)
|
|
#
|
|
#
|
|
# COPYRIGHT (c) 2010 by Toni Mattis <solaris@live.de>
|
|
# ------------------------------------------------------------------------------
|
|
|
|
'''
|
|
Module for elliptic curve arithmetic over a prime field GF(n).
|
|
E(GF(n)) takes the form y**2 == x**3 - p*x - q (mod n) for a prime n.
|
|
|
|
0. Structures used by this module
|
|
|
|
PARAMETERS and SCALARS are non-negative (long) integers.
|
|
|
|
A POINT (x, y), usually denoted p1, p2, ...
|
|
is a pair of (long) integers where 0 <= x < n and 0 <= y < n
|
|
|
|
A POINT in PROJECTIVE COORDINATES, usually denoted jp1, jp2, ...
|
|
takes the form (X, Y, Z, Z**2, Z**3) where x = X / Z**2
|
|
and y = Y / z**3. This form is called Jacobian coordinates.
|
|
|
|
The NEUTRAL element "0" or "O" is represented by None
|
|
in both coordinate systems.
|
|
|
|
1. Basic Functions
|
|
|
|
euclid() Is the Extended Euclidean Algorithm.
|
|
inv() Computes the multiplicative inversion modulo n.
|
|
curve_q() Finds the curve parameter q (mod n)
|
|
when p and a point are given.
|
|
element() Tests whether a point (x, y) is on the curve.
|
|
|
|
2. Point transformations
|
|
|
|
to_projective() Converts a point (x, y) to projective coordinates.
|
|
from_projective() Converts a point from projective coordinates
|
|
to (x, y) using the transformation described above.
|
|
neg() Computes the inverse point -P in both coordinate
|
|
systems.
|
|
|
|
3. Slow point arithmetic
|
|
|
|
These algorithms make use of basic geometry and modular arithmetic
|
|
thus being suitable for small numbers and academic study.
|
|
|
|
add() Computes the sum of two (x, y)-points
|
|
mul() Perform scalar multiplication using "double & add"
|
|
|
|
4. Fast point arithmetic
|
|
|
|
These algorithms make use of projective coordinates, signed binary
|
|
expansion and a JSP-like approach (joint sparse form).
|
|
|
|
The following functions consume and return projective coordinates:
|
|
|
|
addf() Optimized point addition.
|
|
doublef() Optimized point doubling.
|
|
mulf() Highly optimized scalar multiplication.
|
|
muladdf() Highly optimized addition of two products.
|
|
|
|
The following functions use the optimized ones above but consume
|
|
and output (x, y)-coordinates for a more convenient usage:
|
|
|
|
mulp() Encapsulates mulf()
|
|
muladdp() Encapsulates muladdf()
|
|
|
|
For single additions add() is generally faster than an encapsulation of
|
|
addf() which would involve expensive coordinate transformations.
|
|
Hence there is no addp() and doublep().
|
|
'''
|
|
|
|
# BASIC MATH -------------------------------------------------------------------
|
|
|
|
def euclid(a, b):
|
|
'''Solve x*a + y*b = ggt(a, b) and return (x, y, ggt(a, b))'''
|
|
# Non-recursive approach hence suitable for large numbers
|
|
x = yy = 0
|
|
y = xx = 1
|
|
while b:
|
|
q = a // b
|
|
a, b = b, a % b
|
|
x, xx = xx - q * x, x
|
|
y, yy = yy - q * y, y
|
|
return xx, yy, a
|
|
|
|
def inv(a, n):
|
|
'''Perform inversion 1/a modulo n. a and n should be COPRIME.'''
|
|
# coprimality is not checked here in favour of performance
|
|
i = euclid(a, n)[0]
|
|
while i < 0:
|
|
i += n
|
|
return i
|
|
|
|
def curve_q(x, y, p, n):
|
|
'''Find curve parameter q mod n having point (x, y) and parameter p'''
|
|
return ((x * x - p) * x - y * y) % n
|
|
|
|
def element(point, p, q, n):
|
|
'''Test, whether the given point is on the curve (p, q, n)'''
|
|
if point:
|
|
x, y = point
|
|
return (x * x * x - p * x - q) % n == (y * y) % n
|
|
else:
|
|
return True
|
|
|
|
def to_projective(p):
|
|
'''Transform point p given as (x, y) to projective coordinates'''
|
|
if p:
|
|
return (p[0], p[1], 1, 1, 1)
|
|
else:
|
|
return None # Identity point (0)
|
|
|
|
def from_projective(jp, n):
|
|
'''Transform a point from projective coordinates to (x, y) mod n'''
|
|
if jp:
|
|
return (jp[0] * inv(jp[3], n)) % n, (jp[1] * inv(jp[4], n)) % n
|
|
else:
|
|
return None # Identity point (0)
|
|
|
|
def neg(p, n):
|
|
'''Compute the inverse point to p in any coordinate system'''
|
|
return (p[0], (n - p[1]) % n) + p[2:] if p else None
|
|
|
|
|
|
# POINT ADDITION ---------------------------------------------------------------
|
|
|
|
# addition of points in y**2 = x**3 - p*x - q over <Z/nZ*; +>
|
|
def add(p, q, n, p1, p2):
|
|
'''Add points p1 and p2 over curve (p, q, n)'''
|
|
if p1 and p2:
|
|
x1, y1 = p1
|
|
x2, y2 = p2
|
|
if (x1 - x2) % n:
|
|
s = ((y1 - y2) * inv(x1 - x2, n)) % n # slope
|
|
x = (s * s - x1 - x2) % n # intersection with curve
|
|
return (x, n - (y1 + s * (x - x1)) % n)
|
|
else:
|
|
if (y1 + y2) % n: # slope s calculated by derivation
|
|
s = ((3 * x1 * x1 - p) * inv(2 * y1, n)) % n
|
|
x = (s * s - 2 * x1) % n # intersection with curve
|
|
return (x, n - (y1 + s * (x - x1)) % n)
|
|
else:
|
|
return None
|
|
else: # either p1 is not none -> ret. p1, otherwiese p2, which may be
|
|
return p1 if p1 else p2 # none too.
|
|
|
|
|
|
# faster addition: redundancy in projective coordinates eliminates
|
|
# expensive inversions mod n.
|
|
def addf(p, q, n, jp1, jp2):
|
|
'''Add jp1 and jp2 in projective (jacobian) coordinates.'''
|
|
if jp1 and jp2:
|
|
|
|
x1, y1, z1, z1s, z1c = jp1
|
|
x2, y2, z2, z2s, z2c = jp2
|
|
|
|
s1 = (y1 * z2c) % n
|
|
s2 = (y2 * z1c) % n
|
|
|
|
u1 = (x1 * z2s) % n
|
|
u2 = (x2 * z1s) % n
|
|
|
|
if (u1 - u2) % n:
|
|
|
|
h = (u2 - u1) % n
|
|
r = (s2 - s1) % n
|
|
|
|
hs = (h * h) % n
|
|
hc = (hs * h) % n
|
|
|
|
x3 = (-hc - 2 * u1 * hs + r * r) % n
|
|
y3 = (-s1 * hc + r * (u1 * hs - x3)) % n
|
|
z3 = (z1 * z2 * h) % n
|
|
|
|
z3s = (z3 * z3) % n
|
|
z3c = (z3s * z3) % n
|
|
|
|
return (x3, y3, z3, z3s, z3c)
|
|
|
|
else:
|
|
if (s1 + s2) % n:
|
|
return doublef(p, q, n, jp1)
|
|
else:
|
|
return None
|
|
else:
|
|
return jp1 if jp1 else jp2
|
|
|
|
# explicit point doubling using redundant coordinates
|
|
def doublef(p, q, n, jp):
|
|
'''Double jp in projective (jacobian) coordinates'''
|
|
if not jp:
|
|
return None
|
|
x1, y1, z1, z1p2, z1p3 = jp
|
|
|
|
y1p2 = (y1 * y1) % n
|
|
a = (4 * x1 * y1p2) % n
|
|
b = (3 * x1 * x1 - p * z1p3 * z1) % n
|
|
x3 = (b * b - 2 * a) % n
|
|
y3 = (b * (a - x3) - 8 * y1p2 * y1p2) % n
|
|
z3 = (2 * y1 * z1) % n
|
|
z3p2 = (z3 * z3) % n
|
|
|
|
return x3, y3, z3, z3p2, (z3p2 * z3) % n
|
|
|
|
|
|
# SCALAR MULTIPLICATION --------------------------------------------------------
|
|
|
|
# scalar multiplication p1 * c = p1 + p1 + ... + p1 (c times) in O(log(n))
|
|
def mul(p, q, n, p1, c):
|
|
'''multiply point p1 by scalar c over curve (p, q, n)'''
|
|
res = None
|
|
while c > 0:
|
|
if c & 1:
|
|
res = add(p, q, n, res, p1)
|
|
c >>= 1 # c = c / 2
|
|
p1 = add(p, q, n, p1, p1) # p1 = p1 * 2
|
|
return res
|
|
|
|
|
|
# this method allows _signed_bin() to choose between 1 and -1. It will select
|
|
# the sign which leaves the higher number of zeroes in the binary
|
|
# representation (the higher GDB).
|
|
def _gbd(n):
|
|
'''Compute second greatest base-2 divisor'''
|
|
i = 1
|
|
if n <= 0: return 0
|
|
while not n % i:
|
|
i <<= 1
|
|
return i >> 2
|
|
|
|
|
|
# This method transforms n into a binary representation having signed bits.
|
|
# A signed binary expansion contains more zero-bits hence reducing the number
|
|
# of additions required by a multiplication algorithm.
|
|
#
|
|
# Example: 15 ( 0b1111 ) can be written as 16 - 1, resulting in (1,0,0,0,-1)
|
|
# and saving 2 additions. Subtraction can be performed as
|
|
# efficiently as addition.
|
|
def _signed_bin(n):
|
|
'''Transform n into an optimized signed binary representation'''
|
|
r = []
|
|
while n > 1:
|
|
if n & 1:
|
|
cp = _gbd(n + 1)
|
|
cn = _gbd(n - 1)
|
|
if cp > cn: # -1 leaves more zeroes -> subtract -1 (= +1)
|
|
r.append(-1)
|
|
n += 1
|
|
else: # +1 leaves more zeroes -> subtract +1 (= -1)
|
|
r.append(+1)
|
|
n -= 1
|
|
else:
|
|
r.append(0) # be glad about one more zero
|
|
n >>= 1
|
|
r.append(n)
|
|
return r[::-1]
|
|
|
|
|
|
# This multiplication algorithm combines signed binary expansion and
|
|
# fast addition using projective coordinates resulting in 5 to 10 times
|
|
# faster multiplication.
|
|
def mulf(p, q, n, jp1, c):
|
|
'''Multiply point jp1 by c in projective coordinates'''
|
|
sb = _signed_bin(c)
|
|
res = None
|
|
jp0 = neg(jp1, n) # additive inverse of jp1 to be used fot bit -1
|
|
for s in sb:
|
|
res = doublef(p, q, n, res)
|
|
if s:
|
|
res = addf(p, q, n, res, jp1) if s > 0 else \
|
|
addf(p, q, n, res, jp0)
|
|
return res
|
|
|
|
# Encapsulates mulf() in order to enable flat coordinates (x, y)
|
|
def mulp(p, q, n, p1, c):
|
|
'''Multiply point p by c using fast multiplication'''
|
|
return from_projective(mulf(p, q, n, to_projective(p1), c), n)
|
|
|
|
|
|
# Sum of two products using Shamir's trick and signed binary expansion
|
|
def muladdf(p, q, n, jp1, c1, jp2, c2):
|
|
'''Efficiently compute c1 * jp1 + c2 * jp2 in projective coordinates'''
|
|
s1 = _signed_bin(c1)
|
|
s2 = _signed_bin(c2)
|
|
diff = len(s2) - len(s1)
|
|
if diff > 0:
|
|
s1 = [0] * diff + s1
|
|
elif diff < 0:
|
|
s2 = [0] * -diff + s2
|
|
|
|
jp1p2 = addf(p, q, n, jp1, jp2)
|
|
jp1n2 = addf(p, q, n, jp1, neg(jp2, n))
|
|
|
|
precomp = ((None, jp2, neg(jp2, n)),
|
|
(jp1, jp1p2, jp1n2),
|
|
(neg(jp1, n), neg(jp1n2, n), neg(jp1p2, n)))
|
|
res = None
|
|
|
|
for i, j in zip(s1, s2):
|
|
res = doublef(p, q, n, res)
|
|
if i or j:
|
|
res = addf(p, q, n, res, precomp[i][j])
|
|
return res
|
|
|
|
# Encapsulate muladdf()
|
|
def muladdp(p, q, n, p1, c1, p2, c2):
|
|
'''Efficiently compute c1 * p1 + c2 * p2 in (x, y)-coordinates'''
|
|
return from_projective(muladdf(p, q, n,
|
|
to_projective(p1), c1,
|
|
to_projective(p2), c2), n)
|
|
|
|
# POINT COMPRESSION ------------------------------------------------------------
|
|
|
|
# Compute the square root modulo n
|
|
|
|
|
|
# Determine the sign-bit of a point allowing to reconstruct y-coordinates
|
|
# when x and the sign-bit are given:
|
|
def sign_bit(p1):
|
|
'''Return the signedness of a point p1'''
|
|
return p1[1] % 2 if p1 else 0
|
|
|
|
# Reconstruct the y-coordinate when curve parameters, x and the sign-bit of
|
|
# the y coordinate are given:
|
|
def y_from_x(x, p, q, n, sign):
|
|
'''Return the y coordinate over curve (p, q, n) for given (x, sign)'''
|
|
|
|
# optimized form of (x**3 - p*x - q) % n
|
|
a = (((x * x) % n - p) * x - q) % n
|
|
|
|
|
|
|
|
if __name__ == "__main__":
|
|
import rsa
|
|
import time
|
|
|
|
t = time.time()
|
|
n = rsa.get_prime(256/8, 20)
|
|
tp = time.time() - t
|
|
p = rsa.random.randint(1, n)
|
|
p1 = (rsa.random.randint(1, n), rsa.random.randint(1, n))
|
|
q = curve_q(p1[0], p1[1], p, n)
|
|
r1 = rsa.random.randint(1,n)
|
|
r2 = rsa.random.randint(1,n)
|
|
q1 = mulp(p, q, n, p1, r1)
|
|
q2 = mulp(p, q, n, p1, r2)
|
|
s1 = mulp(p, q, n, q1, r2)
|
|
s2 = mulp(p, q, n, q2, r1)
|
|
s1 == s2
|
|
tt = time.time() - t
|
|
|
|
def test(tcount, bits = 256):
|
|
n = rsa.get_prime(bits/8, 20)
|
|
p = rsa.random.randint(1, n)
|
|
p1 = (rsa.random.randint(1, n), rsa.random.randint(1, n))
|
|
q = curve_q(p1[0], p1[1], p, n)
|
|
p2 = mulp(p, q, n, p1, rsa.random.randint(1, n))
|
|
|
|
c1 = [rsa.random.randint(1, n) for i in xrange(tcount)]
|
|
c2 = [rsa.random.randint(1, n) for i in xrange(tcount)]
|
|
c = zip(c1, c2)
|
|
|
|
t = time.time()
|
|
for i, j in c:
|
|
from_projective(addf(p, q, n,
|
|
mulf(p, q, n, to_projective(p1), i),
|
|
mulf(p, q, n, to_projective(p2), j)), n)
|
|
t1 = time.time() - t
|
|
t = time.time()
|
|
for i, j in c:
|
|
muladdp(p, q, n, p1, i, p2, j)
|
|
t2 = time.time() - t
|
|
|
|
return tcount, t1, t2
|
|
|
|
|
|
|