-
-
Save nlitsme/c9031c7b9bf6bb009e5a to your computer and use it in GitHub Desktop.
| 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) |
maybe a python2 compatibility issue? I have now updated it to work both with python2 and python3.
Thanks, it works now.
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?
Thank you for your replies and for your time.
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] shows that dividing by value is multiplying by '1 // value' https://github.com/nlitsme/python-bcutils/blob/master/gfp.py#L88[2] shows that 'div' = multiply by modular inverse https://github.com/nlitsme/python-bcutils/blob/master/modinv.py#L20[3] 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 (x2,y2)= pt2 I need to calculate in python both:
- pt1 / 2 or e.g. pt1 + 5
and also - 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?
Thanks again!
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.
Hello.
When I run this code, I get TypeError: pow() 3rd argument not allowed unless all arguments are integers.
Could you explain, what I'm doing wrong.