mini_ecdsa.py

#Elliptic curve basics, tools for finding rational points, and ECDSA implementation.
#Brendan Cordy, 2015
#modified for Python 3.2.5

from fractions import Fraction
from math import ceil, sqrt
from random import SystemRandom, randrange
from hashlib import sha256
from time import clock

def pow_mod(x, y, z): #fast pow mod (need to fast check is point on the curve on not)
    "Calculate (x ** y) % z efficiently."
    number = 1
    while y:
        if y & 1:
            number = number * x % z
        y >>= 1
        x = x * x % z
    return number

#Affine Point (+Infinity) on an Elliptic Curve ---------------------------------------------------

class Point(object):
    #Construct a point with two given coordindates.
    def __init__(self, x, y):
        self.x, self.y = x, y
        self.inf = False

    #Construct the point at infinity.
    @classmethod
    def atInfinity(cls):
        P = cls(0, 0)
        P.inf = True
        return P

    #The secp256k1 generator.
    @classmethod
    def secp256k1(cls):
        return cls(55066263022277343669578718895168534326250603453777594175500187360389116729240,
                   32670510020758816978083085130507043184471273380659243275938904335757337482424)

    def __str__(self):
        if self.inf:
            return 'Inf'
        else:
            return '(' + str(self.x) + ',' + str(self.y) + ')'

    def __eq__(self,other):
        if self.inf:
            return other.inf
        elif other.inf:
            return self.inf
        else:
            return self.x == other.x and self.y == other.y

    def is_infinite(self):
        return self.inf

#Elliptic Curves over any Field ------------------------------------------------------------------

class Curve(object):
    #Set attributes of a general Weierstrass cubic y^2 = x^3 + ax^2 + bx + c over any field.
    def __init__(self, a, b, c, char, exp):
        self.a, self.b, self.c = a, b, c
        self.char, self.exp = char, exp
        print(self)

    def __str__(self):
        #Cases for 0, 1, -1, and general coefficients in the x^2 term.
        if self.a == 0:
            aTerm = ''
        elif self.a == 1:
            aTerm = ' + x^2'
        elif self.a == -1:
            aTerm = ' - x^2'
        elif self.a < 0:
            aTerm = " - " + str(-self.a) + 'x^2'
        else:
            aTerm = " + " + str(self.a) + 'x^2'
        #Cases for 0, 1, -1, and general coefficients in the x term.
        if self.b == 0:
            bTerm = ''
        elif self.b == 1:
            bTerm = ' + x'
        elif self.b == -1:
            bTerm = ' - x'
        elif self.b < 0:
            bTerm = " - " + str(-self.b) + 'x'
        else:
            bTerm = " + " + str(self.b) + 'x'
        #Cases for 0, 1, -1, and general coefficients in the constant term.
        if self.c == 0:
            cTerm = ''
        elif self.c < 0:
            cTerm = " - " + str(-self.c)
        else:
            cTerm = " + " + str(self.c)
        #Write out the nicely formatted Weierstrass equation.
        self.eq = 'y^2 = x^3' + aTerm + bTerm + cTerm
        #Print prettily.
        if self.char == 0:
            return self.eq + ' over Q'
        elif self.exp == 1:
            return self.eq + ' over ' + 'F_' + str(self.char)
        else:
            return self.eq + ' over ' + 'F_' + str(self.char) + '^' + str(self.exp)

    #Compute the discriminant.
    def discriminant(self):
        a, b, c = self.a, self.b, self.c
        return -4*a*a*a*c + a*a*b*b + 18*a*b*c - 4*b*b*b - 27*c*c

    #Compute the order of a point on the curve.
    def order(self, P):
        Q = P
        orderP = 1
        #Add P to Q repeatedly until obtaining the identity (point at infinity).
        while not Q.is_infinite():
            Q = self.add(P,Q)
            orderP += 1
        return orderP

    #List all multiples of a point on the curve.
    def generate(self, P):
        Q = P
        orbit = [str(Point.atInfinity())]
        #Repeatedly add P to Q, appending each (pretty printed) result.
        while not Q.is_infinite():
            orbit.append(str(Q))
            Q = self.add(P,Q)
        return orbit

    #Double a point on the curve.
    def double(self, P):
        return self.add(P,P)

    #Add P to itself k times.
    def mult(self, P, k):
        if P.is_infinite():
            return P
        elif k == 0:
            return Point.atInfinity()
        elif k < 0:
            return self.mult(self.invert(P), -k)
        else:
            #Convert k to a bitstring and use peasant multiplication to compute the product quickly.
            b = bin(k)[2:]
            return self.repeat_additions(P, b, 1)

    #Add efficiently by repeatedly doubling the given point, and adding the result to a running
    #total when, after the ith doubling, the ith digit in the bitstring b is a one.
    def repeat_additions(self, P, b, n):
        if b == '0':
            return Point.atInfinity()
        elif b == '1':
            return P
        elif b[-1] == '0':
            return self.repeat_additions(self.double(P), b[:-1], n+1)
        elif b[-1] == '1':
            return self.add(P, self.repeat_additions(self.double(P), b[:-1], n+1))

    #Returns a pretty printed list of points.
    def show_points(self):
        return [str(P) for P in self.get_points()]

#Elliptic Curves over Q --------------------------------------------------------------------------

class CurveOverQ(Curve):
    #Construct a Weierstrass cubic y^2 = x^3 + ax^2 + bx + c over Q.
    def __init__(self, a, b, c):
        Curve.__init__(self, a, b, c, 0, 0)

    def contains(self, P):
        if P.is_infinite():
            return True
        else:
            return P.y*P.y == P.x*P.x*P.x + self.a*P.x*P.x + self.b*P.x + self.c

    def get_points(self):
        #Start with the point at infinity.
        points = [Point.atInfinity()]
        #The only possible y values are divisors of the discriminant.
        for y in divisors(self.discriminant()):
            #Each possible y value yields a monic cubic polynomial in x, whose roots
            #must divide the constant term.
            const_term = self.c - y*y
            if const_term != 0:
                for x in divisors(const_term):
                    P = Point(x,y)
                    if 0 == x*x*x + self.a*x*x + self.b*x + const_term and self.has_finite_order(P):
                        points.append(P)
            #If the constant term is zero, factor out x and look for rational roots
            #of the resulting quadratic polynomial. Any such roots must divide b.
            elif self.b != 0:
                for x in divisors(self.b):
                    P = Point(x,y)
                    if 0 == x*x*x+self.a*x*x+self.b*x+const_term and self.has_finite_order(P):
                        points.append(P)
            #If the constant term and b are both zero, factor out x^2 and look for rational
            #roots of the resulting linear polynomial. Any such roots must divide a.
            elif self.a != 0:
                for x in divisors(self.a):
                    P = Point(x,y)
                    if 0 == x*x*x+self.a*x*x+self.b*x+const_term and self.has_finite_order(P):
                        points.append(P)
            #If the constant term, b, and a are all zero, we have 0 = x^3 + c - y^2 with
            #const_term = c - y^2 = 0, so (0,y) is a point on the curve.
            else:
                points.append(Point(0,y))
        #Ensure that there are no duplicates in our list of points.
        unique_points = []
        for P in points:
            addP = True
            for Q in unique_points:
                if P == Q:
                    addP = False
            if addP:
                unique_points.append(P)
        return unique_points

    def invert(self, P):
        if P.is_infinite():
            return P
        else:
            return Point(P.x, -P.y)

    def add(self, P_1, P_2):
        #Compute the differences in the coordinates.
        y_diff = P_2.y - P_1.y
        x_diff = P_2.x - P_1.x
        #Cases involving the point at infinity.
        if P_1.is_infinite():
            return P_2
        elif P_2.is_infinite():
            return P_1
        #Case for adding an affine point to its inverse.
        elif x_diff == 0 and y_diff != 0:
            return Point.atInfinity()
        #Case for adding an affine point to itself.
        elif x_diff == 0 and y_diff == 0:
            #If the point is on the x-axis, there's a vertical tangent there (assuming
            #the curve in nonsingular) so we obtain the point at infinity.
            if P_1.y == 0:
                return Point.atInfinity()
            #Otherwise the result is an affine point on the curve we can arrive at by
            #following the tangent line, whose slope is given below.
            else:
                ld = Fraction(3*P_1.x*P_1.x + 2*self.a*P_1.x + self.b, 2*P_1.y)
        #Case for adding two distinct affine points, where we compute the slope of
        #the secant line through the two points.
        else:
            ld = Fraction(y_diff, x_diff)
        #Use the slope of the tangent line or secant line to compute the result.
        nu = P_1.y - ld*P_1.x
        x = ld*ld - self.a  -P_1.x - P_2.x
        y = -ld*x - nu
        return Point(x,y)

    #Use the Nagell-Lutz Theorem and Mazur's Theorem to potentially save time.
    def order(self, P):
        Q = P
        orderP = 1
        #Add P to Q repeatedly until obtaining the point at infinity.
        while not Q.is_infinite():
            Q = self.add(P,Q)
            orderP += 1
            #If we ever obtain non integer coordinates, the point has infinite order.
            if Q.x != int(Q.x) or Q.y != int(Q.y):
                return -1
            #Moreover, all finite order points have order at most 12.
            if orderP > 12:
                return -1
        return orderP

    def has_finite_order(self, P):
            return not self.order(P) == -1

    def torsion_group(self):
        highest_order = 1
        #Find the rational point with the highest order.
        for P in self.get_points():
            if self.order(P) > highest_order:
                highest_order = self.order(P)
        #If this point generates the entire torsion group, the torsion group is cyclic.
        if highest_order == len(self.get_points()):
            print('Z/' + str(highest_order) + 'Z')
        #If not, by Mazur's Theorem the torsion group must be a direct product of Z/2Z
        #with the cyclic group generated by the highest order point.
        else:
            print('Z/2Z x ' + 'Z/' + str(highest_order) + 'Z')
        print(C.show_points())

#Elliptic Curves over Prime Order Fields ---------------------------------------------------------

class CurveOverFp(Curve):
    #Construct a Weierstrass cubic y^2 = x^3 + ax^2 + bx + c over Fp.
    def __init__(self, a, b, c, p):
        Curve.__init__(self, a, b, c, p, 1)

    #The secp256k1 curve.
    @classmethod
    def secp256k1(cls):
        return cls(0, 0, 7, 2**256-2**32-2**9-2**8-2**7-2**6-2**4-1)

    def contains(self, P, pow_mod=0):
        if P.is_infinite():
            return True
        elif(pow_mod==0):		#	check:	(y^2) mod p == (x^3 + ax^2 + bx + c) mod p
            return (P.y*P.y) % self.char == (P.x*P.x*P.x + self.a*P.x*P.x + self.b*P.x + self.c) % self.char	# true/false
        else:					#  check this faster, using pow_mod
            return (pow_mod(P.y, 2, self.char) == (pow_mod(P.x, 3, self.char) + pow_mod(self.a*P.x, 2, self.char) + (self.b * P.x) + self.c)%self.char); #true/false

    def get_points(self):
        #Start with the point at infinity.
        points = [Point.atInfinity()]

        #Just brute force the rest.
        for x in range(self.char):
                for y in range(self.char):
                    P = Point(x,y)
                    if (y*y) % self.char == (x*x*x + self.a*x*x + self.b*x + self.c) % self.char:
                        points.append(P)
        return points

    def invert(self, P):	#additive inversion
        if P.is_infinite():
            return P
        else:
            return Point(P.x, -P.y % self.char)	#-Q = (x_Q, -y_Q%p)

    def add(self, P_1, P_2):
        #Adding points over Fp and can be done in exactly the same way as adding over Q,
        #but with of the all arithmetic now happening in Fp.
        y_diff = (P_2.y - P_1.y) % self.char
        x_diff = (P_2.x - P_1.x) % self.char
        if P_1.is_infinite():
            return P_2
        elif P_2.is_infinite():
            return P_1
        elif x_diff == 0 and y_diff != 0:
            return Point.atInfinity()
        elif x_diff == 0 and y_diff == 0:
            if P_1.y == 0:
                return Point.atInfinity()
            else:
                ld = ((3*P_1.x*P_1.x + 2*self.a*P_1.x + self.b) * mult_inv(2*P_1.y, self.char)) % self.char
        else:
            ld = (y_diff * mult_inv(x_diff, self.char)) % self.char
        nu = (P_1.y - ld*P_1.x) % self.char
        x = (ld*ld - self.a - P_1.x - P_2.x) % self.char
        y = (-ld*x - nu) % self.char
        return Point(x,y)

    def subtract(self, point, another): #C.subtract(from, minus)
        if(another.is_infinite()):
            return point;
        else:
            return self.add(point, self.invert(another));	# (P - Q) = P + (-Q)
        
    def divide_point(self, point, k, n):
#        print("k", k, "n", n, "mult_inv(k, n)", mult_inv(k, n));
        return self.mult(point, mult_inv(k, n))	# (P / k) = P * (k^-1 mod n)

    def getY(self, X, Y_parity_bit=0):						# for each X there is two Y on curve, odd and even.
        a = (pow_mod(X, 3, self.char) + 7) % self.char		# a = ( ( ( (x^3) mod n ) ax^2 + bx + c ) mod n )
        Y = pow_mod(a, (self.char+1)//4, self.char)			# y = a^{(n+1)//4} mod n
        if Y % 2 != Y_parity_bit:
            Y = -Y % self.char								# invert Y
#        return Y;
        return Y%self.char;

    def get_point_by_X(self, X, Y_parity_bit=0):			# each point can be encoded as X and parity_bit for one from two Y's.
#        return Point(X, self.getY(X, Y_parity_bit));		#return this.
        return Point(X%self.char, self.getY(X, Y_parity_bit));		#return this.

#Elliptic Curves over Prime Power Order Fields ---------------------------------------------------

class CurveOverFq(Curve):
    #Construct a Weierstrass cubic y^2 = x^3 + ax^2 + bx + c over Fp^n.
    def __init__(self, a, b, c, p, n, irred_poly):
        self.irred_poly = irred_poly
        Curve.__init__(self, a, b, c, p, n)

    #TODO: Implement it!

#Number Theoretic Functions ----------------------------------------------------------------------

def divisors(n):
    divs = [0]
    for i in range(1, abs(n) + 1):
        if n % i == 0:
            divs.append(i)
            divs.append(-i)
    return divs

#Extended Euclidean algorithm.
def euclid(sml, big):
    #When the smaller value is zero, it's done, gcd = b = 0*sml + 1*big.
    if sml == 0:
        return (big, 0, 1)
    else:
        #Repeat with sml and the remainder, big%sml.
        g, y, x = euclid(big % sml, sml)
        #Backtrack through the calculation, rewriting the gcd as we go. From the values just
        #returned above, we have gcd = y*(big%sml) + x*sml, and rewriting big%sml we obtain
        #gcd = y*(big - (big//sml)*sml) + x*sml = (x - (big//sml)*y)*sml + y*big.
        return (g, x - (big//sml)*y, y)

#    Compute the multiplicative modular inverse mod n of a with 0 < a < n.
def mult_inv(a, n):	#    return b, for which: a * b ===  1 (mod n), or return tuple with gcd
    g, x, y = euclid(a, n)
    #If gcd(a,n) is not one, then a has no multiplicative inverse.
    if g != 1:
        #raise ValueError('multiplicative inverse does not exist for value a mod n =', a, "mod", n);
        #return str('multiplicative inverse does not exist for value a mod n = ')+str(a)+str(" mod ")+str(n);
        return ("gcd!=1", g);    # if value have no modular inverse - don't show error and just return tuple, with gcd.
    #If gcd(a,n) = 1, and gcd(a,n) = x*a + y*n, x is the multiplicative inverse of a.
    else:
        return x % n

#ECDSA functions ---------------------------------------------------------------------------------

#Use sha256 to hash a message, and return the hash value as an interger.
def hash(message):
    return int(sha256(message).hexdigest(), 16)

#Hash the message and return integer whose binary representation is the the L leftmost bits
#of the hash value, where L is the bit length of n.
def hash_and_truncate(message, n):
    h = hash(message)
    b = bin(h)[2:len(bin(n))]
    return int(b, 2)

#Generate a keypair using the point P of order n on the given curve. The private key is a
#positive integer d smaller than n, and the public key is Q = dP.
def generate_keypair(curve, P, n):
    sysrand = SystemRandom()
    d = sysrand.randrange(1, n)
    Q = curve.mult(P, d)
    print("Priv key: d = " + str(d))
    print("Publ key: Q = " + str(Q))
    return (d, Q)

#Create a digital signature for the string message using a given curve with a distinguished
#point P which generates a prime order subgroup of size n.
def sign(message, curve, P, n, keypair):
    #Extract the private and public keys, and compute z by hashing the message.
    d, Q = keypair
    z = hash_and_truncate(message, n)
    #Choose a randomly selected secret point kP then compute r and s.
    r, s = 0, 0
    while r == 0 or s == 0:
        k = 4
        R = curve.mult(P, k)
        r = R.x % n
        s = (mult_inv(k, n) * (z + r*d)) % n
    print('ECDSA sig: (Q, r, s) = (' + str(Q) + ', ' + str(r) + ', ' + str(s) + ')')
    return (Q, r, s)

#Verify the string message is authentic, given an ECDSA signature generated using a curve with
#a distinguished point P that generates a prime order subgroup of size n.
def verify(message, curve, P, n, sig):
    Q, r, s = sig
    #Confirm that Q is on the curve.
    if Q.is_infinite() or not curve.contains(Q):
        return False
    #Confirm that Q has order that divides n.
    if not curve.mult(Q,n).is_infinite():
        return False
    #Confirm that r and s are at least in the acceptable range.
    if r > n or s > n:
        return False
    #Compute z in the same manner used in the signing procedure,
    #and verify the message is authentic.
    z = hash_and_truncate(message, n)
    w = mult_inv(s, n) % n
    u_1, u_2 = z * w % n, r * w % n
    C_1, C_2 = curve.mult(P, u_1), curve.mult(Q, u_2)
    C = curve.add(C_1, C_2)
    return r % n == C.x % n

#Key Cracking Functions --------------------------------------------------------------------------

#Find d for which Q = dP by simply trying all possibilities
def crack_brute_force(curve, P, n, Q, start=0, show_each=0):	#it can be parallelized
    start_time = clock()
    print("start: ", start, "up to:", n);
    for d in range(start, n):
        if(show_each!=0 and (d%show_each==0)): print("crack_brute_force: last incorrect d = ", d); #this number can be used as "start" value to continue.
        if curve.mult(P,d) == Q:
            end_time = clock()
            print("Priv key: d = " + str(d))
            print("Time: " + str(round(end_time - start_time, 3)) + " secs")
            break

#Find d for which Q = dP using the baby-step giant-step algortihm.
def crack_baby_giant(curve, P, n, Q, m=0):	#m can be specified, 0 by default sqrt(n)
    start_time = clock()
    if m==0: m = int(ceil(sqrt(n))); #for m = sqrt(n), Baby-Step-Giant-Step working for O(sqrt(n)) multiplies, and eating sqrt(n) memory.
    #else - O(n/m) multiplies, and eating m memory for storing all points for baby-steps.
    #Build a hash table with all bP with 0 < b < m using a dictionary. The dicitonary value
    #stores b so that it can be quickly recovered after a matching giant step hash.
    baby_table = {}
    for b in range(m):
        bP = curve.mult(P,b)
        baby_table[str(bP)] = b
    #Check if Q - gmP is in the hash table for all 0 < g < m. If we get such a matching hash,
    #we have Q - gmP = bP, so extract b from the dictionary, then Q = (b + gm)P.
    for g in range((int(n/m)+1)):
        R = curve.add(Q, curve.invert(curve.mult(P, g*m)))
        if str(R) in baby_table.keys():
            b = baby_table[str(R)]
            end_time = clock()
            print("Priv key: d = " + str((b + g*m) % n))
            print("Time: " + str(round(end_time - start_time, 3)) + " secs")
            break

#Find d for which Q = dP using Pollard's rho algorithm.
def crack_rho(curve, P, n, Q, bits):
    start_time = clock()
    while True:
        R_list = []
        #Compute 2^bits randomly selected linear combinations of P and Q, storing them as triples
        #of the form (aP + bQ, a, b) in R_list.
        for i in range(2**bits):
            a, b = randrange(0,n), randrange(0,n)
            R_list.append((curve.add(curve.mult(P,a), curve.mult(Q,b)), a, b))
        #Compute a new random linear combination of P and Q to start the cycle-finding.
        aT, bT = randrange(0,n), randrange(0,n)
        aH, bH = aT, bT
        T = curve.add(curve.mult(P,aT), curve.mult(Q,bT))
        H = curve.add(curve.mult(P,aH), curve.mult(Q,bH))
        while True:
            #Advance the tortoise one step, by adding a point in R_list determined by the last b
            #bits in the binary explansion of the x coordinate of the current position.
            j = int(bin(T.x)[len(bin(T.x)) - bits : len(bin(T.x))], 2)
            T, aT, bT = curve.add(T, R_list[j][0]), (aT + R_list[j][1]) % n, (bT + R_list[j][2]) % n
            #Advance the hare two steps, again by adding points in R_list determined by the last
            #b bits in the binary explansion of the x coordinate of the current position.
            for i in range(2):
                j = int(bin(H.x)[len(bin(H.x)) - bits : len(bin(H.x))], 2)
                H, aH, bH = curve.add(H, R_list[j][0]), (aH + R_list[j][1]) % n, (bH + R_list[j][2]) % n
            #If the tortoise and hare arrive at the same point, a cycle has been found.
            if(T == H):
                break
        #It is possible that the tortoise and hare arrive at exactly the same linear combination.
        if bH == bT:
            end_time = clock()
            print("Rho failed with identical linear combinations")
            print(str(end_time - start_time) + " secs")
        else:
            end_time = clock()
            print("\nPollard-rho, results:");
            modular_inverse = mult_inv((bH - bT) % n, n);        #Some values have no modular inverse, and will be returned gcd.

            if isinstance(modular_inverse, int):
                modular_inverse =  modular_inverse % n
                priv = ((aT - aH) * modular_inverse) % n;
                pub = curve.mult(P, priv);
                if(pub.__eq__(Q)):
                    print("Priv key: d = " + str(priv) + str(", pub: ")+str(pub));
                    print("Time: " + str(round(end_time - start_time, 3)) + " secs")
                    break;
                else:
                    print("pub not equal Q", "pub", pub, "Q", Q);
                    continue;
            elif ( (isinstance(modular_inverse, tuple)) and (modular_inverse[0] == "gcd!=1") ) :
                print("invalid modinv... gcd = ", modular_inverse[1], "Try again...");
                continue;

#Find d by exploiting two messages signed with the same value of k.
def crack_from_ECDSA_repeat_k(curve, P, n, m1, sig1, m2, sig2):
    Q1, r1, s1 = sig1
    Q2, r2, s2 = sig2
    #Check that the two messages were signed with the same k. This check may pass even if m1
    #and m2 were signed with distinct k, in which case the value of d computed below will be
    #wrong, but this is a very unlikely scenario.
    if not r1 == r2:
        print("Messages signed with distinct k")
    else:
        z1 = hash_and_truncate(m1, n)
        z2 = hash_and_truncate(m2, n)
        k = (z1 - z2) * mult_inv((s1 - s2) % n, n) % n
        d = mult_inv(r1, n) * ((s1 * k) % n - z1) % n
        print("Priv key: d = " + str(d))

# Run tests_mini_ecdsa.py to test the functions in this mini_ecdsa.py
# Run tests_ECC.py to test Elliptic-Curve Encryption functions from ECC.py