| |
| # --- 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 |
| |
| |
| |
