Converting Montgomery Elliptic Curves to Weierstrass Form

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)
		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')
		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:
			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.

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.

About samkerr

I'm an eclectic person. I like to dabble in a multitude of things. I'm sure you'll find my blog reflects that.
This entry was posted in Uncategorized. Bookmark the permalink.

13 Responses to Converting Montgomery Elliptic Curves to Weierstrass Form

  1. Brian says:

    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).

  2. samkerr says:

    Thanks Brian! I have updated the entry to reflect your comment.

  3. fxa says:

    How did you calculate Gy: 0x20ae19a1b8a086b4e01edd2c7748d14c923d4d7e6d7c61b229e9c5a27eced3d9?
    All the other values I get with sage, too, but not this one

  4. fxa says:

    I tried within GF(P)
    y2 = Gx^3 + A* Gx + b
    and looked on both roots, but get
    Gy = 1c0ee62f128259ba9ae4279367a19badb11bc7061367be4696529a289de7a18c
    Gy = 63f119d0ed7da645651bd86c985e64524ee438f9ec9841b969ad65d762185e61

  5. fxa says:

    Got it!
    please remove my comments and sorry for the noise

  6. samkerr says:

    Glad you were able to figure it out fxa! What ended up being the trouble? What are you working on by chance?

  7. fxa says:

    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

  8. fxa says:

    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

  9. absinthe says:

    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.

  10. panax says:

    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))

  11. liu says:

    But , Can i get the order of the G point?

  12. Sven Anderson says:

    I believe the line “x = (x + a_m/3)%p” must be “x = (x + a_m*modinv(3, p))%p”

  13. samkerr says:

    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”.

Leave a Reply

Your email address will not be published. Required fields are marked *