example of bitcoin curve calculations in python

curve_example.py

from __future__ import print_function, division
"""
Example of how calculations on the secp256k1 curve work.
secp256k1 is the name of the elliptic curve used by bitcoin
see http://bitcoin.stackexchange.com/questions/25382
"""

## the prime modules used in the elliptic curve coordinate calculations
p = 2**256 - 2**32 - 977

## the group order, used with the elliptic curve scalar multiplication calculations
n = 2**256 - 0x14551231950B75FC4402DA1732FC9BEBF


def inverse(x, p):
    """
    Calculate the modular inverse of x ( mod p )
    
    the modular inverse is a number such that:
    
    (inverse(x, p) * x) % p  ==  1 
    
    you could think of this as: 1/x
    """
    inv1 = 1
    inv2 = 0
    while p != 1 and p!=0:
        inv1, inv2 = inv2, inv1 - inv2 * (x // p)
        x, p = p, x % p

    return inv2


def dblpt(pt, p):
    """
    Calculate pt+pt = 2*pt
    """
    if pt is None:
        return None
    (x,y)= pt
    if y==0:
        return None 
        
    # Calculate 3*x^2/(2*y)  modulus p
    slope= 3*pow(x,2,p)*inverse(2*y,p)
    
    xsum= pow(slope,2,p)-2*x
    ysum= slope*(x-xsum)-y
    return (xsum%p, ysum%p)

def addpt(p1,p2, p):
    """
    Calculate p1+p2
    """
    if p1 is None or p2 is None:
        return None
    (x1,y1)= p1
    (x2,y2)= p2
    if x1==x2:
        if y1==y2:
            return dblpt(p1, p)
        elif (y1+y2) % p == 0:
            return None
        else:
            raise Exception("invalid points added")

    # calculate (y1-y2)/(x1-x2)  modulus p
    slope=(y1-y2)*inverse(x1-x2, p)
    xsum= pow(slope,2,p)-(x1+x2)
    ysum= slope*(x1-xsum)-y1
    return (xsum%p, ysum%p)

def ptmul(pt,a, p):
    """
    Scalar multiplication: calculate pt*a
    
    basically adding pt to itself a times
    """
    scale= pt
    acc=None
    while a:
        if a&1:
            if acc is None:
                acc= scale
            else:
                acc= addpt(acc,scale, p)
        scale= dblpt(scale, p)
        a >>= 1
    return acc


def isoncurve(pt,p):
    """
    returns True when pt is on the secp256k1 curve
    """
    (x,y)= pt
    return (y**2 - x**3 - 7)%p == 0

# (Gx,Gy)  is the secp256k1 generator point
Gx=0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798
Gy=0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8
g= (Gx,Gy)
g2=dblpt(g, p)
print("    2*G= (%x,%x)" % g2)
print("  G+2*G= (%x,%x)" % addpt(g, g2, p))
print("2*G+2*G= (%x,%x)" % addpt(g2, g2, p))

privkey= 0xf8ef380d6c05116dbed78bfdd6e6625e57426af9a082b81c2fa27b06984c11f3
print(" -> pubkey= (%x,%x)" % ptmul(g, privkey, p))

# note that here I have to use the grouporder 'n'
halvekey = (privkey * inverse(2, n)) % n
print(" halvekey = %x" % halvekey)

halvepub = ptmul(g, halvekey, p)
print("  pub / 2 = (%x,%x)" % halvepub)
print("2 * halvepub: (%x,%x)" % addpt(halvepub, halvepub, p))

"""
for reference, the numbers printed should be:
    2*G= (c6047f9441ed7d6d3045406e95c07cd85c778e4b8cef3ca7abac09b95c709ee5,1ae168fea63dc339a3c58419466ceaeef7f632653266d0e1236431a950cfe52a)
  G+2*G= (f9308a019258c31049344f85f89d5229b531c845836f99b08601f113bce036f9,388f7b0f632de8140fe337e62a37f3566500a99934c2231b6cb9fd7584b8e672)
2*G+2*G= (e493dbf1c10d80f3581e4904930b1404cc6c13900ee0758474fa94abe8c4cd13,51ed993ea0d455b75642e2098ea51448d967ae33bfbdfe40cfe97bdc47739922)
 -> pubkey= (71ee918bc19bb566e3a5f12c0cd0de620bec1025da6e98951355ebbde8727be3,37b3650efad4190b7328b1156304f2e9e23dbb7f2da50999dde50ea73b4c2688)
 halvekey = fc779c06b60288b6df6bc5feeb73312e88f8a3f027e5ac2bf7ba6cc9b441299a
  pub / 2 = (d6c8d2ecbce68e0bf127c889be491c2a4f26a84c15b1fe8688f42728baea6c18,49a300586b09945c63e89df3c1739fb30567cec3da037d0402222c3b7f05c6e2)
2 * halvepub: (71ee918bc19bb566e3a5f12c0cd0de620bec1025da6e98951355ebbde8727be3,37b3650efad4190b7328b1156304f2e9e23dbb7f2da50999dde50ea73b4c2688)
"""

# this shows how to multiply by -1
print("      G = (%x,%x)" % g)
gminus = ptmul(g,n-1, p)
print("     -G = (%x,%x)" % gminus)
gminusminus = ptmul(gminus,n-1, p)
print("  -(-G) = (%x,%x)" % gminusminus)
gminusplus = addpt(g,gminus, p)

# this will output the point at infinity, currently represented as 'None'
print(" G+(-G) = %s" % gminusplus)

Hello nlitsme, thanks for your script.
Isn't the negative -G of a curve point G computed by ptmul(G,p-1, p)?
I tried it but it did not work (i tested it with the generator point): neither is the negative -(-G) of -G the origin G, nor is G + (-G) = NULL (infinity point).
Script appendu at the end of curve_example.py:
...

(Gx,Gy) is the secp256k1 generator point

Gx=0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798
Gy=0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8
g= (Gx,Gy)
print(" G = (%x,%x)" % g)
gminus = ptmul(g,p-1, p)
print(" -G = (%x,%x)" % gminus)
gminusminus = ptmul(gminus,p-1, p)
print(" -(-G) = (%x,%x)" % gminusminus)
gminusplus = addpt(g,gminus, p)
print(" G+(-G) = (%x,%x)" % gminusplus)

and here its result of the script:
G = (79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798,483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8)
-G = (2541d1403fc71a5d927923b20a673e769284b16e7d1f597f9413dc64e82fc48,621239e38c2af6bc145db2424cf44bb8ed12351858f67c7c6130a2f69fe4a28c)
-(-G) = (c8ec1c603c5603a266ae5c63378969b305712dd4717e2bf1787b44210fa94772,15b8bab353c0c54ef7227cdc3536c3b590c5fa5d30cd27efaf703cb1dd30eea3)
G+(-G) = (b42b34a9748238715c0b8b853b6939aabc3b5224bcdd4b4b65e902417e7af914,c1064ce5dfb9e0247fa170896d939c3b87404dca2f648560936138086ac129e9)

Is my assumption wrong that -G is computed by scalar multiplication with -1?
Cheers, Hubert

you should multiply by the grouporder minus one, this is a different large prime:

g= (Gx,Gy)
n=0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141

print(" G = (%x,%x)" % g)
gminus = ptmul(g,n-1, p)
print(" -G = (%x,%x)" % gminus)
gminusminus = ptmul(gminus,n-1, p)
print(" -(-G) = (%x,%x)" % gminusminus)
gminusplus = addpt(g,gminus, p)
print(" G+(-G) = (%x,%x)" % gminusplus)

That last calculation has the wrong result because addpt has a bug: it should have checked if y1 == y2 too before calling dblpt. and return 0 when y1==-y2

How do you print out the modular inverse?

def inverse(x, p):

Ok, I think I got it:

print(" -> inverse mod= %x" % inverse(privkey, p))

-> inverse mod= 170a6c6134c7234543b7dadad1e1a2fd79a6caf9ddfb776c41b2d63cdae6875c

Nevertheless, I am getting a different value from this:

[0xf8ef380d6c05116dbed78bfdd6e6625e57426af9a082b81c2fa27b06984c11f3 x
57896044618658097711785492504343953926418782139537452191302581570759080747169
= 0xfc779c06b60288b6df6bc5feeb73312e88f8a3f027e5ac2bf7ba6cc9b441299a]

=

[0xf8ef380d6c05116dbed78bfdd6e6625e57426af9a082b81c2fa27b06984c11f3 / 2
= 0xfc779c06b60288b6df6bc5feeb73312e88f8a3f027e5ac2bf7ba6cc9b441299a]

I added some notes about the group order in the script

Damm... This is harder to understand as I tought.

So, the halfkey of 0x800000000000000000000000000000 is 0x400000000000000000000000000000
And I am getting this value right.

How the hell is the halfkey from 0xf8ef380d6c05116dbed78bfdd6e6625e57426af9a082b81c2fa27b06984c11f3
equal to 0xfc779c06b60288b6df6bc5feeb73312e88f8a3f027e5ac2bf7ba6cc9b441299a ?

UPDATE:

Sorry, i got it. The key is odd, the reason why of this strange value.

Yes, that is a way to look at at:

if you have an odd number 2*k+1, in this case: k = 0x7C779C06B60288B6DF6BC5FEEB73312F2BA1357CD0415C0E17D13D834C2608F9

then (2*k+1)/2 = k + 1/2
and because 1/2 is the modular-inverse of 2 = 0x7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF5D576E7357A4501DDFE92F46681B20A1

the answer would be: 0x7C779C06B60288B6DF6BC5FEEB73312F2BA1357CD0415C0E17D13D834C2608F9 + 0x7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF5D576E7357A4501DDFE92F46681B20A1
= 0xFC779C06B60288B6DF6BC5FEEB73312E88F8A3F027E5AC2BF7BA6CC9B441299A

I have been able to reverse secp256k1 elliptic curve point doubling formula and point addition formula that use in the process of deriving the public key of a Bitcoin private key which I know as follow
Point Addition Formula In Python3:

s = (Qy - Gy) * pow(Qx - Gx, -1, p) % p
Rx = (s**2 - Qx - Gx) % p
Ry = (s * (Qx - Rx) - Qy) % p

Point Doubling/Multiplication Formula In Python3:

s = (3 * Qx**2) * pow(Qy*2, -1, p) % p
Rx = (s**2 - Qx*2) % p
Ry = (s * (Qx - Rx) - Qy) % p

if the original points Qx and Qy are used to derive the new points Rx and Ry, am able to use the new derive points to solve for the original points, with this how can I drive the new private key from the public key?

Hi, can this solution be added to the gist? https://math.stackexchange.com/q/4852418

So you are asking: how do I found 'Q if P = n * Q ?
You can use this

division in secp256k1 is equivalent to multiplying by the modular inverse.

So you are asking: how do I found 'Q if P = n * Q ? You can use this

division in secp256k1 is equivalent to multiplying by the modular inverse.

In that case yes. Given 257 bits long natural integer n and point P retrieving point Q.
I still think you should add the code to your gist in order to work with Koblitz curves.