python implementation of ecdsa calculations, demonstrating how to recover a private key from two signatures with identical 'r', and demonstrating how to find the public key from a signature and message, or from two signatures.

ecdsa_demo.py

from __future__ import print_function, division
"""
By Willem Hengeveld <[email protected]>

ecdsa implementation in python

demonstrating several 'unconventional' calculations,
like finding a public key from a signature,
and finding a private key from 2 signatures with identical 'r'
"""
# (gcd,c,d) = GCD(a, b)  ===> a*c+b*d==gcd
def GCD(a, b):
    prevx, x = 1, 0
    prevy, y = 0, 1
    while b:
        q = a//b
        x, prevx = prevx - q*x, x
        y, prevy = prevy - q*y, y
        a, b = b, a%b
    return a, prevx, prevy


def modinv(x, m):
    (gcd, c, d) = GCD(x,m)
    return c

def samefield(a,b): 
    """
    determine if a uses the same field 
    """
    if a.field!=b.field: 
        raise Exception("field mismatch")
    return True

class FiniteField:
    """
    FiniteField implements a value modulus a number.
    """
    class Value:
        """
        represent a value in the FiniteField

        this class forwards all operations to the FiniteField class
        """
        def __init__(self, field, value):
            self.field = field
            self.value = int(value)

        # Value * int
        def __add__(self, rhs): return self.field.add(self, self.field.value(rhs))
        def __sub__(self, rhs): return self.field.sub(self, self.field.value(rhs))
        def __mul__(self, rhs): return self.field.mul(self, self.field.value(rhs))
        def __div__(self, rhs): return self.field.div(self, self.field.value(rhs))
        def __truediv__(self, rhs): return self.__div__(rhs)
        def __floordiv__(self, rhs): return self.__div__(rhs)
        def __pow__(self, rhs): return self.field.pow(self, rhs)

        # int * Value
        def __radd__(self, lhs): return self.field.add(self.field.value(lhs), self)
        def __rsub__(self, lhs): return self.field.sub(self.field.value(lhs), self)
        def __rmul__(self, lhs): return self.field.mul(self.field.value(lhs), self)
        def __rdiv__(self, lhs): return self.field.div(self.field.value(lhs), self)
        def __rtruediv__(self, lhs): return self.__rdiv__(lhs)
        def __rfloordiv__(self, lhs): return self.__rdiv__(lhs)
        def __rpow__(self, lhs): return self.field.pow(self.field.value(lhs), self)

        def __eq__(self, rhs): return self.field.eq(self, self.field.value(rhs))
        def __ne__(self, rhs): return not (self==rhs)

        def __str__(self): return "0x%x" % self.value
        def __neg__(self): return self.field.neg(self)


        def sqrt(self, flag): return self.field.sqrt(self, flag)
        def inverse(self):  return self.field.inverse(self)
        def __nonzero__(self): return self.field.nonzero(self)
        def __bool__(self): return self.__nonzero__() != 0
        def __int__(self): return self.field.intvalue(self)

    def __init__(self, p):
        self.p = p

    """
    several basic operators
    """
    def add(self, lhs, rhs): return samefield(lhs, rhs) and self.value((lhs.value + rhs.value) % self.p)
    def sub(self, lhs, rhs): return samefield(lhs, rhs) and self.value((lhs.value - rhs.value) % self.p)
    def mul(self, lhs, rhs): return samefield(lhs, rhs) and self.value((lhs.value * rhs.value) % self.p)
    def div(self, lhs, rhs): return samefield(lhs, rhs) and self.value((lhs.value * rhs.inverse()) % self.p)
    def pow(self, lhs, rhs): return self.value(pow(lhs.value, int(rhs), self.p))
    def eq(self, lhs, rhs): return (lhs.value-rhs.value) % self.p == 0
    def neg(self, val): return self.value(self.p-val.value)

    def sqrt(self, val, flag):
        """
        calculate the square root modulus p
        """
        if not val:
            return val
        sw = self.p % 8
        if sw==3 or sw==7:
            res = val**((self.p+1)//4)
        elif sw==5:
            x = val**((self.p+1)//4)
            if x==1:
                res = val**((self.p+3)//8)
            else:
                res = (4*val)**((self.p-5)//8)*2*val
        else:
            raise Exception("modsqrt non supported for (p%8)==1")
        if res.value%2==flag:
            return res
        else:
            return -res

    def inverse(self, value):
        """
        calculate the multiplicative inverse
        """
        return modinv(value.value, self.p)

    def nonzero(self, x):
        return 1 if not (x.value % self.p)==0 else 0

    def value(self, x):
        """
        converts an integer or FinitField.Value to a value of this FiniteField.
        """
        return x if isinstance(x, FiniteField.Value) and x.field==self else FiniteField.Value(self, x)

    def zero(self):
        """
        returns the additive identity value

        meaning:  a + 0 = a
        """
        return FiniteField.Value(self, 0)
    def one(self):
        """
        returns the multiplicative identity value

        meaning a * 1 = a
        """
        return FiniteField.Value(self, 1)

    def intvalue(self, x):
        return x.value % self.p


class EllipticCurve:
    """
    EllipticCurve implements a point on a elliptic curve
    """
    class Point:
        """
        represent a value in the EllipticCurve

        this class forwards all operations to the EllipticCurve class
        """
        def __init__(self, curve, x, y):
            self.curve = curve
            self.x = x
            self.y = y
        # Point + Point
        def __add__(self, rhs): return self.curve.add(self, rhs)
        def __sub__(self, rhs): return self.curve.sub(self, rhs)

        # Point * int   or Point * Value
        def __mul__(self, rhs): return self.curve.mul(self, rhs)
        def __div__(self, rhs): return self.curve.div(self, rhs)
        def __truediv__(self, rhs): return self.__div__(rhs)
        def __floordiv__(self, rhs): return self.__div__(rhs)

        def __eq__(self, rhs): return self.curve.eq(self, rhs)
        def __ne__(self, rhs): return not (self==rhs)

        def __str__(self): return "(%s,%s)" % (self.x, self.y)
        def __neg__(self): return self.curve.neg(self)

        def __nonzero__(self): return self.curve.nonzero(self)
        def __bool__(self): return self.__nonzero__() != 0
        def isoncurve(self):
            return self.curve.isoncurve(self)

    def __init__(self, field, a, b):
        self.field = field
        self.a = field.value(a)
        self.b = field.value(b)

    def add(self, p, q):
        """
        perform elliptic curve addition
        """
        if not p: return q
        if not q: return p

        # calculate the slope of the intersection line
        if p==q:
            if p.y==0:
                return self.zero()
            l = (3* p.x**2 + self.a) // (2* p.y)
        elif p.x==q.x:
            return self.zero()
        else:
            l = (p.y-q.y)//(p.x-q.x)

        # calculate the intersection point
        x = l**2 - ( p.x + q.x )
        y = l*(p.x-x)-p.y
        return self.point(x,y)

    # subtraction is :  a - b  =  a + -b
    def sub(self, lhs, rhs): return lhs + -rhs

    # scalar multiplication is implemented like repeated addition
    def mul(self, pt, scalar): 
        scalar = int(scalar)
        accumulator = self.zero()
        shifter = pt
        while scalar != 0:
            bit = scalar % 2
            if bit:
                accumulator += shifter
            shifter += shifter
            scalar //= 2

        return accumulator

    def div(self, pt, scalar):
        """
        scalar division:  P / a = P * (1/a)

        scalar is assumed to be of type FiniteField(grouporder)
        """
        return pt * (1//scalar)

    def eq(self, lhs, rhs): return lhs.x==rhs.x and lhs.y==rhs.y
    def neg(self, pt):
        return self.point(pt.x, -pt.y)
    def nonzero(self, pt):
        return 1 if pt.x or pt.y else 0
    def zero(self):
        """
        Return the additive identity point ( aka '0' )

        P + 0 = P
        """
        return self.point(self.field.zero(), self.field.zero())

    def point(self, x, y):
        """
        construct a point from 2 values
        """
        return EllipticCurve.Point(self, self.field.value(x), self.field.value(y))

    def isoncurve(self, p):
        """
        verifies if a point is on the curve
        """
        return not p or (p.y**2 == p.x**3 + self.a*p.x + self.b)

    def decompress(self, x, flag):
        """
        calculate the y coordinate given only the x value.
        there are 2 possible solutions, use 'flag' to select.
        """
        x = self.field.value(x)
        ysquare = x**3 + self.a*x+self.b

        return self.point(x, ysquare.sqrt(flag))




class ECDSA:
    """
    Digital Signature Algorithm using Elliptic Curves
    """
    def __init__(self, ec, G, n):
        self.ec = ec
        self.G = G
        self.GFn = FiniteField(n)

    def calcpub(self, privkey):
        """
        calculate the public key for private key x

        return G*x
        """
        return self.G * self.GFn.value(privkey)

    def sign(self, message, privkey, secret):
        """
        sign the message using private key and sign secret

        for signsecret k, message m, privatekey x
        return (G*k,  (m+x*r)/k)
        """
        m = self.GFn.value(message)
        x = self.GFn.value(privkey)
        k = self.GFn.value(secret)

        R = self.G * k

        r = self.GFn.value(R.x)
        s = (m + x*r) // k

        return (r, s)

     
    def verify(self, message, pubkey, rnum, snum):
        """
        Verify the signature

        for message m, pubkey Y, signature (r,s)

        r = xcoord(R)

        verify that :  G*m+Y*r=R*s

        this is true because: { Y=G*x, and R=G*k, s=(m+x*r)/k }
         
        G*m+G*x*r = G*k*(m+x*r)/k  ->
        G*(m+x*r) = G*(m+x*r)

        several ways to do the verification:
            r == xcoord[ G*(m/s) + Y*(r/s) ]  <<< the standard way
            R * s == G*m + Y*r
            r == xcoord[ (G*m + Y*r)/s) ]
 
        """
        m = self.GFn.value(message)
        r = self.GFn.value(rnum)
        s = self.GFn.value(snum)

        R = self.G * (m//s) + pubkey * (r//s)

        # alternative methods of verifying
        #RORG = self.ec.decompress(r, 0)
        #RR = self.G * m + pubkey * r
        #print("#1: %s .. %s"  % (RR, RORG*s))
        #print("#2: %s .. %s"  % (RR*(1//s), r))
        #print("#3: %s .. %s"  % (R, r))

        return R.x == r

    def findpk(self, message, rnum, snum, flag):
        """
        find pubkey Y from message m, signature (r,s)
        Y = (R*s-G*m)/r

        note that there are 2 pubkeys related to a signature
        """
        m = self.GFn.value(message)
        r = self.GFn.value(rnum)
        s = self.GFn.value(snum)

        R = self.ec.decompress(r, flag)

        #return (R*s - self.G * m)*(1//r)
        return R*(s//r) - self.G * (m//r)

    def findpk2(self, r1, s1, r2, s2, flag1, flag2):
        """
        find pubkey Y from 2 different signature on the same message
        sigs: (r1,s1) and (r2,s2)

        returns  (R1*s1-R2*s2)/(r1-r2)
        """
        R1 = self.ec.decompress(r1, flag1)
        R2 = self.ec.decompress(r2, flag2)

        rdiff = self.GFn.value(r1-r2)

        return (R1*s1-R2*s2)*(1//rdiff)

    def crack2(self, r, s1, s2, m1, m2):
        """
        find signsecret and privkey from duplicate 'r'

        signature (r,s1) for message m1
        and signature (r,s2) for message m2

        s1 = (m1 + x*r)/k
        s2 = (m2 + x*r)/k

        subtract -> (s1-s2) = (m1-m2)/k  ->  k = (m1-m2)/(s1-s2)

        -> privkey =  (s1*k-m1)/r  .. or  (s2*k-m2)/r
        """
        sdelta = self.GFn.value(s1-s2)
        mdelta = self.GFn.value(m1-m2)

        secret = mdelta // sdelta
        x1 = self.crack1(r, s1, m1, secret)
        x2 = self.crack1(r, s2, m2, secret)

        if x1 != x2:
            print("x1=%s" % x1)
            print("x2=%s" % x2)

        return (secret, x1)

    def crack1(self, rnum, snum, message, signsecret):
        """
        find privkey, given signsecret k, message m, signature (r,s)

        x = (s*k-m)/r
        """
        m = self.GFn.value(message)
        r = self.GFn.value(rnum)
        s = self.GFn.value(snum)
        k = self.GFn.value(signsecret)
        return (s*k-m)//r

def secp256k1():
    """
    create the secp256k1 curve
    """
    GFp = FiniteField(2**256 - 2**32 - 977)
    ec = EllipticCurve(GFp, 0, 7)
    generator = ec.point( 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798, 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8  )
    grouporder = 2**256 - 432420386565659656852420866394968145599
    return ECDSA(ec, generator, grouporder)

def verifytest(calced, expected, descr):
    """
    verifytest is used to verify test results
    """
    if type(calced)!=type(expected):
        if type(expected)==str:
            calced = "%s" % calced
    if calced!=expected:
        print("ERROR in %s: got %s, expected %s" % (descr, calced, expected))

def test_gfp():
    gfp = FiniteField(37)
    verifytest(gfp.value(2)**16,  "0x9", "<2>^16")
    verifytest(2**gfp.value(16),  "0x9", "2^<16>")
    verifytest(gfp.value(1)//4,  "0x1c", "<1>/4")
    verifytest(1//gfp.value(4),  "0x1c", "1/<4>")
    verifytest(gfp.value(1)==38,  "True", "<1>==38")
    verifytest(1==gfp.value(38),  "True", "1==<38>")
    verifytest(gfp.value(1)!=38,  "False", "<1>!=38")
    verifytest(gfp.value(16)*16,  "0x22", "<16>*16")
    verifytest(16*gfp.value(16),  "0x22", "16*<16>")

def test_ec():
    GFp = FiniteField(2**256 - 2**32 - 977)
    GFn = FiniteField(2**256 - 432420386565659656852420866394968145599)
    ec = EllipticCurve(GFp, 0, 7)
    G = ec.point( 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798, 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8  )
    verifytest(G.isoncurve(), True, "oncurve")

    a = 901283019823
    b = 100980987879
    P = G*a
    Q = G*b

    verifytest(G*a+G*b, G*(a+b), "p(a+b)=pa+pb")
    verifytest((P+Q)*a, P*a+Q*a, "(p+q)a=pa+qa")

    verifytest(P//GFn.value(a), G, "P/a")


def test_dsa():
    dsa = secp256k1()

    priv = 0x1234
    pubkey = dsa.calcpub(priv)
    signsecret1 = 0x1111
    signsecret2 = 0x2222

    msg1 = 0x1234123412341234123412341234123412341234123412341234123412341234
    (r1,s1) = dsa.sign(msg1, priv, signsecret1)
    check1 = dsa.verify(msg1, pubkey, r1, s1)

    verifytest(pubkey, "(0x37a4aef1f8423ca076e4b7d99a8cabff40ddb8231f2a9f01081f15d7fa65c1ba,0xb96ced90a1b8f9b43a18fc900ff55af2be0e94b90a434fca5b9e226b835024cd)", "pubkey")

    verifytest(priv, dsa.crack1(r1, s1, msg1, signsecret1), "crack1")

    verifytest(r1, "0x7592aab5d43618dda13fba71e3993cd7517a712d3da49664c06ee1bd3d1f70af", "r1")
    verifytest(s1, "0x8e578a508331374bcbb1618ea3a8c9c63d49e9d42e0ed605b8c74910cfa50c11", "s1")
    verifytest(check1, True, "verify1")

    msg2 = 0x1111111111111111111111111111111111111111111111111111111111111111
    (r2,s2) = dsa.sign(msg2, priv, signsecret1)
    check2 = dsa.verify(msg2, pubkey, r2, s2)

    verifytest(r2, "0x7592aab5d43618dda13fba71e3993cd7517a712d3da49664c06ee1bd3d1f70af", "r2")
    verifytest(s2, "0x351a200607d9aae72e3fb30fe41cf92dcd0b22117f57123df005974290d9429b", "s2")
    verifytest(check2, True, "verify2")

    verifytest(r1, r2, "rvalues")
    (crackedsecret, crackedprivkey) = dsa.crack2(r1, s1, s2, msg1, msg2)

    verifytest(crackedsecret, signsecret1, "crackedsecret")
    verifytest(crackedprivkey, "0x%x" % priv, "crackedpriv")

    """
    Note that we can find public keys, which can validate a signature, but did not sign it!
    """
    for flag in range(2):
        pk1 = dsa.findpk(msg1, r1, s1, flag)
        check1 = dsa.verify(msg1, pk1, r1, s1)
        print("%d -> %s pk1=%s" % (flag, check1, pk1))

        pk2 = dsa.findpk(msg2, r2, s2, flag)
        check2 = dsa.verify(msg2, pk2, r2, s2)
        print("%d -> %s pk2=%s" % (flag, check2, pk2))

    """
    only one of these 4 verifies our msg correctly
    """
    (r3,s3) = dsa.sign(msg1, priv, signsecret2)
    for flag1 in range(2):
        for flag2 in range(2):
            pk = dsa.findpk2(r1,s1,r3,s3,flag1,flag2)
            check1 = dsa.verify(msg1, pk, r1,s1)
            check3 = dsa.verify(msg1, pk, r3,s3)
            print("%d,%d : %s %s -> %s" % (flag1, flag2, check1, check3, pk))

def test_crack():
    """
    Demonstrate cracking an actual bitcoin key, see https://gist.github.com/nlitsme/f3c9953a420012bd413a684068a770ff
    """
    pk=(0xdbd0c61532279cf72981c3584fc32216e0127699635c2789f549e0730c059b81, 0xae133016a69c21e23f1859a95f06d52b7bf149a8f2fe4e8535c8a829b449c5ff)
    r=0xd47ce4c025c35ec440bc81d99834a624875161a26bf56ef7fdc0f5d52f843ad1
    s1=0x44e1ff2dfd8102cf7a47c21d5c9fd5701610d04953c6836596b4fe9dd2f53e3e
    s2=0x9a5f1c75e461d7ceb1cf3cab9013eb2dc85b6d0da8c3c6e27e3a5a5b3faa5bab
    m1=0xc0e2d0a89a348de88fda08211c70d1d7e52ccef2eb9459911bf977d587784c6e
    m2=0x17b0f41c8c337ac1e18c98759e83a8cccbc368dd9d89e5f03cb633c265fd0ddc

    # k:0x7a1a7e52797fc8caaa435d2a4dace39158504bf204fbe19f14dbb427faee50ae
    # x:0xc477f9f65c22cce20657faa5b2d1d8122336f851a508a1ed04e479c34985bf96

    E = secp256k1()
    k, x = E.crack2(r, s1, s2, m1, m2)
    print("k = %s" % k)
    print("x = %s" % x)

def main():
    test_gfp()
    test_ec()
    test_dsa()
    test_crack()

if __name__ == '__main__':
    main()

And update the script to work with python3

corrected an incorrect 'gcd' implementation, now the script works with the real world example from https://bitaps.com/9ec4bc49e828d924af1d1029cacf709431abbde46d59554b62bc270e3b29c4b1

And the latest update: using bool and int operators to get a bit more readable code.

I made a new repository: https://github.com/nlitsme/bitcoinexplainer doing the above, and more in javascript.