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)

On Wednesday, 29 September 2021 10:23:27 CEST George Stergiopoulos wrote: > On Tuesday, 28 September 2021 19:32:47 CEST George Stergiopoulos wrote: > On Tuesday, 28 September 2021 18:34:54 CEST George Stergiopoulos wrote: Quick question: Why do you rightshift one bit at the end of the multiplication (a >>= 1)? Didn't see that in any ecdsa bitcoin guide.. > It is a very basic binary multiplication algorithm. Basically I go over the bits in 'a' one at a time, starting with the lowest bits. instead of shifting right, I could have written this as 'a //= 2' - integer divide by 2. Also, if I would not shift, 'a' would stay the same, and the loop would never stop. willem Yep found it myself, that's why I deleted the comment but you are fast! :-) On another note, do you have a short script for Point Division that I can add to extend the code above? > I assume you mean scalar division? dividing two points is very difficult, this difficulty is the basis of ellipiccurve cryptography. If you want to look into how this could be done, look for: 'discrete logarithm' this repository has a more detailed implementation of the basic elliptic curve operations: [https://github.com/nlitsme/python-bcutils/blob/master/ec.py#L99[1]](https://github.com/nlitsme/python-bcutils/blob/master/ec.py#L99%5B1%5D) shows that dividing by value is multiplying by '1 // value' [https://github.com/nlitsme/python-bcutils/blob/master/gfp.py#L88[2]](https://github.com/nlitsme/python-bcutils/blob/master/gfp.py#L88%5B2%5D) shows that 'div' = multiply by modular inverse [https://github.com/nlitsme/python-bcutils/blob/master/modinv.py#L20[3]](https://github.com/nlitsme/python-bcutils/blob/master/modinv.py#L20%5B3%5D) shows how the modular inverse is calculated willem > […](#) Thank you for your reply. Yes, my original question was aiming at python code. I am aware of the disc. log problem. My experiments with ell. curve points require both operations that you mention, for example for (x1,y1)= pt1 and (x1,y1)= pt1 I need to calculate in python both: 1) pt1 / 2 or e.g. pt1 + 5 and also 2) pt1 / pt2 If I understand correctly, division (or any similar operation) of point by element value (eg. z=5) in python is as simple as result = (pt1 / z) %p and for dividing points (d. log problem), based on your code above we need something like: result = (pt1 * inverse(pt2, p)) %p which is NP-Hard, correct?

yes, it is NP-hard. but I don't think it involves calculating the inverse of a point. since also multiplication of points is not a defined operation. willem

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.