Recently I needed to convert a Montgomery form elliptic curve (Curve25519) to Weierstrass form. This post captures some of that work and code for arbitrary Montgomery form curves as well. I am going to do a lot of math hand waving, but will try to provide references as necessary.
Curve25519 is a Montgomery form curve discovered by Daniel Bernstein. It is fast and has some interesting security properties. It is a Montgomery curve as well, as opposed to the Weierstrass curve, which are what most readers will be familiar with. The NIST elliptic curves are Weierstrass curves.
Montgomery curves are of the form y^2 = x^3 + a*x^2 + x (mod p).
Weierstrass curves are of the form y^2 = x^3 + a*x + b (mod p).
For Curve25519, the Weierstrass parameters are:
P: 0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffedL A: 0x2aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa984914a144L B: 0x7b425ed097b425ed097b425ed097b425ed097b425ed097b4260b5e9c7710c864L Gx: 0x2aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaad245a Gy: 0x20ae19a1b8a086b4e01edd2c7748d14c923d4d7e6d7c61b229e9c5a27eced3d9
There is a formula for conversion on Wikipedia already, I just went ahead and implemented it. Anyways, here is code to print out the Weierstrass parameters for Curve25519 as well as the generator point coordinates.
# Python code def egcd(a, b): if a == 0: return (b, 0, 1) else: g, y, x = egcd(b % a, a) return (g, x - (b // a) * y, y) def modinv(a, m): g, x, y = egcd(a, m) if g != 1: raise Exception('modular inverse does not exist') else: return x % m def modular_sqrt(a, p): """ Find a quadratic residue (mod p) of 'a'. p must be an odd prime. Solve the congruence of the form: x^2 = a (mod p) And returns x. Note that p - x is also a root. 0 is returned is no square root exists for these a and p. The Tonelli-Shanks algorithm is used (except for some simple cases in which the solution is known from an identity). This algorithm runs in polynomial time (unless the generalized Riemann hypothesis is false). """ # Simple cases # if legendre_symbol(a, p) != 1: return 0 elif a == 0: return 0 elif p == 2: return p elif p % 4 == 3: return pow(a, (p + 1) / 4, p) # Partition p-1 to s * 2^e for an odd s (i.e. # reduce all the powers of 2 from p-1) # s = p - 1 e = 0 while s % 2 == 0: s /= 2 e += 1 # Find some 'n' with a legendre symbol n|p = -1. # Shouldn't take long. # n = 2 while legendre_symbol(n, p) != -1: n += 1 # Here be dragons! # Read the paper "Square roots from 1; 24, 51, # 10 to Dan Shanks" by Ezra Brown for more # information # # x is a guess of the square root that gets better # with each iteration. # b is the "fudge factor" - by how much we're off # with the guess. The invariant x^2 = ab (mod p)gx_w = (9 + a_m/3)%p # is maintained throughout the loop. # g is used for successive powers of n to update # both a and b # r is the exponent - decreases with each update # x = pow(a, (s + 1) / 2, p) b = pow(a, s, p) g = pow(n, s, p) r = e while True: t = b m = 0 for m in xrange(r): if t == 1: break t = pow(t, 2, p) if m == 0: return x gs = pow(g, 2 ** (r - m - 1), p) g = (gs * gs) % p x = (x * gs) % p b = (b * g) % p r = m def legendre_symbol(a, p): """ Compute the Legendre symbol a|p using Euler's criterion. p is a prime, a is relatively prime to p (if p divides a, then a|p = 0) Returns 1 if a has a square root modulo p, -1 otherwise. """ ls = pow(a, (p - 1) / 2, p) return -1 if ls == p - 1 else ls # Montgomery form Curve25519 parameters (b*y^2 = x^3 + a*x^2 + x (mod p)) # Change these for a a different Montgomery curve p = pow(2,255)-19 a_m = 486662 b_m = 1 # Here be dragons! This works, so don't really change it interim_a_numerator = 3-pow(a_m,2) interim_a_denominator = 3*pow(b_m,2) interim_a_denominator = modinv(interim_a_denominator,p) a_w = interim_a_numerator * interim_a_denominator # This IS multiplcation, since its modular inverse a_w = a_w % p interim_b_numerator = 2*pow(a_m,3)-9*a_m interim_b_denominator = 27 * pow(b_m,3) interim_b_denominator = modinv(interim_b_denominator, p) b_w = interim_b_numerator * interim_b_denominator # This IS multiplcation, since its modular inverse b_w = b_w % p if __name__ == "__main__": x = 9 # This is Curve25519 specific x = (x + a_m * modinv(3,p))%p # This is needed, since we are doing a change of variables, so this is the analogous change. y2 = pow(x,3) + a_w*x + b_w y2 = y2 % p y = modular_sqrt(y2,p) print "P: " + hex(p) print "A: " + hex(a_w) print "B: " + hex(b_w) print "Gx: " + hex(x) print "Gy: " + hex(y)
Some of the preceding code is from Eli Bendersky’s page. I highly recommend visiting it!
Hopefully this code sample helps you in converting Montgomery curves to Weierstrass form. Do be aware that Weierstrass form elliptic curve math is typically not as fast as Montgomery form, due to not being able to use some tricks. If you’d like to know a little bit more about performance comparisons, a nice overview is available on this Wikipedia page. At a quick glance though, Short Weierstrass curves (what we’re using) take doubling a point takes 11 operations, compared to the 4 of a Montgomery curve. If you can get away with it, use Montgomery curves.
EDIT:
I received a comment (it’s below) that mentioned that since the Montgomery curve was using a change of variables to get the Weierstrass form, the generator point would also need to be updated. I overlooked this and it is now in the code. Thank you very much Brian!
Another note that I had was that it would be much easier to do this math with Sage, rather than Python directly. It will handle all the field math, such as division, relatively transparently and has many operations built into it directly, rather than requiring us to implement it ourself.
I think you’ve made a mistake in applying the formula from Wikipedia; one step in the derivation is to replace x (essentially) with t – A/3. This means the x coordinate in the new equation (t) is shifted from the x coordinate in the original equation (x), so Gx in the Weierstrass form is actually 9 + A/3.
I’ve found a useful check is to plug x back into the original Curve25519 equation, and seeing if I get the same value for y^2 (mod p) (which should be the same if the equations really are equivalent).
Thanks Brian! I have updated the entry to reflect your comment.
How did you calculate Gy: 0x20ae19a1b8a086b4e01edd2c7748d14c923d4d7e6d7c61b229e9c5a27eced3d9?
All the other values I get with sage, too, but not this one
I tried within GF(P)
y2 = Gx^3 + A* Gx + b
and looked on both roots, but get
Gy = 1c0ee62f128259ba9ae4279367a19badb11bc7061367be4696529a289de7a18c
Gy = 63f119d0ed7da645651bd86c985e64524ee438f9ec9841b969ad65d762185e61
Got it!
please remove my comments and sorry for the noise
Glad you were able to figure it out fxa! What ended up being the trouble? What are you working on by chance?
My trouble was to mix a with A. After adding prefixes m_ for montgommery and w_ for weierstrass I did the job. If you are interested in, I can post the sage script
And please remove my comments and your reply, because it does not help any reader. We can communicate by email. I will add a new complete reply with sage-script explaining my pain — that will be informative for all readers
Thanks for the post. It is very useful since many libraries such as MIRACL do not support Montgomery curves and they need to be converted.
gmpy2 makes modular arithmetic much easier:
import gmpy2
mod = 0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed
mA = 486662
mB = 1
xm= 9
ym= 0x20ae19a1b8a086b4e01edd2c7748d14c923d4d7e6d7c61b229e9c5a27eced3d9
a = (((3- (mA**2))%mod)*gmpy2.invert(3*mB**2,mod)) % mod
b = (((2*(mA**3) - 9*mA)%mod)*gmpy2.invert(27*(mB**3),mod)) % mod
x=(xm*gmpy2.invert(mB,mod) + mA*gmpy2.invert(3*mB,mod))%mod
y = (ym*gmpy2.invert(mB,mod)) % mod
print('P: {0:x}'.format(mod))
print('A: {0:x}'.format(a))
print('B: {0:x}'.format(b))
print('Gx: {0:x}'.format(x))
print('Gy: {0:x}'.format(y))
Hi,
But , Can i get the order of the G point?
I believe the line “x = (x + a_m/3)%p” must be “x = (x + a_m*modinv(3, p))%p”
Great point! Updated the post with the change.
Preventing errors like this is one of the reasons I’m a big fan of Sage these days – if you do division with field elements, it understands you really mean “multiply by modular inverse”.