The aim of this tutorial is to implement the X25519 key exchange protocol in Magma, learning some useful Magma techniques along the way. X25519, defined in RFC7748, is essentially a formalization of the techniques used in Bernstein's Curve25519 Elliptic Curve Diffie--Hellman software.

Useful references:

1 Montgomery arithmetic

Our first step is to define Montgomery curves and their arithmetic in Magma.

Start by downloading the Magma file Montgomery.m. This contains some useful pre-defined code, which will let you get straight to the interesting part of the tutorial more quickly.

Montgomery.m defines a function MontgomeryModel(A,B) which, given elements \(A \neq \pm 2\) and \(B \neq 0\) of some field, constructs

  1. a projective plane curve \(C: BY^2Z = X(X^2 + AXZ + Z^2)\), using Magma's intrinsic Curve, which returns an object of type CrvEll;
  2. the morphism \(x: C \to \mathbb{P}^1\) mapping \((X:Y:Z)\) to \((X:Z)\);
  3. the Weierstrass model curve \(E: y^2 = x(x^2 + ABx + B^2)\) (as a Magma CrvEll object, created using EllipticCurve) and the isomorphism \(i: E \to C\);
  4. the x-line pseudo-addition xADD operation for \(C\); and
  5. the x-line pseudo-doubling xDBL operation for \(C\).

It then stores all of these associated items inside the curve \(C\) as attributes (to avoid re-computing them all the time), and returns \(C\)... At least, that's what it will do once you've finished this exercise!

The Weierstrass model and isomorphism will be important for testing our code against the equivalent built-in operations for Weierstrass models in Magma.

1.1 Montgomery models

Our first task in completing MontgomeryModel is to construct the underlying curve \(C = C_{A,B}\). This also allows us to practice defining algebraic curves in Magma.

Start by looking for // EXERCISE 1.1 in the Montgomery.m source file, and completing the defining equation of \(C\).

Now test your code: in the Magma interpreter we should be able to replicate the following (after running load "Montgomery.m"):

> // Constructing Curve25519:
> p := 2^255 - 19;
> IsPrime(p);
true
> k := FiniteField(p);
> Curve25519 := MontgomeryModel(k!486662,k!1);
> Curve25519;
Curve over GF(57896044618658097711785492504343953926634992332820282019728792003\
956564819949) defined by
57896044618658097711785492504343953926634992332820282019728792003956564819948*X\
    ^3 + 5789604461865809771178549250434395392663499233282028201972879200395656\
    4333287*X^2*Z + Y^2*Z + 578960446186580977117854925043439539266349923328202\
    82019728792003956564819948*X*Z^2

The defining equation looks a bit hairy like that, but that's what happens when you negate small quantities modulo medium-sized primes:

> -DefiningEquation(Curve25519);
X^3 + 486662*X^2*Z + 5789604461865809771178549250434395392663499233282028201972\
8792003956564819948*Y^2*Z + X*Z^2

To get the curve parameters \(A\) and \(B\) back, we can use AParameter and BParameter (defined near the top of Montgomery.m):

> // The Montgomery curve parameters
> A := AParameter(Curve25519);
> A;
486662
> B := BParameter(Curve25519);
> B;
1

1.2 The map to the x-line

The arithmetic of Montgomery models is mostly computed using only \(x\)-coordinates, so our next step is to fix the projection to the \(x\)-line. This also allows us to practice defining a morphism in Magma using the map< ... -> ... | ... > constructor.

Extend your MontgomeryModel function so that it

  1. constructs the projective line over the same base field as the curve \(C\), then
  2. defines the map \(x: C \to \mathbb{P}^1\) sending \((X:Y:Z) \mapsto (X:Z)\), and
  3. stores \(x\) as the attribute xMap of \(C\) (we do this bit for you).

The relevant lines of Montgomery.m are

    // EXERCISE 1.2: Construct PP^1 and x: C -> PP^1
    PP1 := ProjectiveSpace(k,1);
    /* uncomment and complete the following lines: */
    // x := map< C -> PP1 | /* FILL IN IMAGE HERE */ >;
    // C`xMap := x;  // set associated attributed here

Note that we can access the map to the \(x\)-line using the function xMap defined near the top of Montgomery.m. Continuing the example above:

> // try out x-mapping:
> C := Curve25519;  // curve defined as above
> x := xMap(C);
> pt := C![0,0,1];
> pt_PP1 := x(pt);
> pt_PP1;
(0 : 1)

Note: the \(x\) map that we have defined behaves slightly unexpectedly on the point at infinity:

> infty := C![0,1,0];
> infty;
(0 : 1 : 0)
> x(infty);
(1 : 0)

Why do you think this is?

1.3 Isomorphism to a Weierstrass model

Before going much further, we are going to need a means of effectively testing our implementation of Montgomery arithmetic. We will do this by defining an isomorphism between \(C\) and a Weierstrass model, which is the form used by Magma's built-in elliptic curve machinery. This also lets us practice defining isomorphisms between geometric objects in Magma using the iso< ... -> ... | ... > constructor.

Extend your MontgomeryModel function so that it

  1. computes an isomorphism to a Weierstrass model for \(C\), and
  2. stores the isomorphism in the WeierstrassIsomorphism attribute of \(C\) (we do this bit for you).

The relevant lines of Montgomery.m are

    // EXERCISE 1.3: Construct E and the isomorphism C <-> E
    E := EllipticCurve(Polynomial([0,B^2,A*B,1]));
    /* uncomment and fill in the following lines: */
    // i := iso< E -> C | /* FILL IN IMAGES HERE */ , /* FILL IN PREIMAGES HERE */ >;
    // C`WeierstrassIsomorphism := i;

For example, The isomorphism from the Weierstrass model can be accessed using WeierstrassModel:

> // Using the isomorphism to a Weierstrass model:
> C := Curve25519; // the curve above
> E, E_to_C := WeierstrassModel(C);
> infty := E_to_C(E!0);
> infty;
(0 : 1 : 0)
> // Let's build a more interesting point...
> // If E is defined by y^2 + h(x)y = f(x) for some f and h, then
> // we can recover f and h from E using HyperellipticPolynomials(E).
> f, h := HyperellipticPolynomials(E);
> assert h eq 0;
> xcoord := k!9;
> ycoord := Sqrt(Evaluate(f,xcoord));
> base_point_E := E![xcoord,ycoord];
> r := 2^252 + 27742317777372353535851937790883648493;
> IsPrime(r);
true
> r*base_point_E;
(0 : 1 : 0)
> base_point_C := E_to_C(base_point_E);
> base_point_C;
(9 : 14781619447589544791020593568409986887264606134616475288964881837755586237401 : 1)

1.4 Pseudo-doubling

The first operation to define for Montgomery curves is the pseudo-doubling operation \(x(P) \mapsto x([2]P)\). For this we will use Magma's func< ... | ... > constructor, which conveniently constructs a function object given arguments and an expression in terms of those arguments.

Recall that if \(x(P) = (X_P:Z_P)\) and \((X_{[2]P}:Z_{[2]P}) = x([2]P)\), then \[(X_{[2]P}:Z_{[2]P}) = (Q\cdot R : S\cdot(R+\tfrac{A+2}{4}S))\] where \(Q = (X_P+Z_P)^2\), \(R = (X_P-Z_P)^2\), and \(S = Q - R\) (so in fact \(S = 4X_PZ_P\)).

Extend your MontgomeryModel function so that it

  1. constructs a function xDBL(xP), which takes one point xP = x(P) on \(\mathbb{P}^1\) and returns x2P = x([2]P) and xDBL on \(\mathbb{P}^1\), and
  2. stores the result in the attribute xDBLOperation (we do this for you).

The relevant piece of code in Montgomery.m is

    // EXERCISE 1.4: Pseudo-doubling
    /* uncomment and fill in the following lines: */
    xDBL := func< xP | /* FILL IN IMAGE HERE */ >;
    C`xDBLOperation := xDBL;  // set associated attribute

Hint: look at the type of what your xDBL function is supposed to be taking and returning: points on \(\mathbb{P}^1\). If PP1 is the projective line, then you can create explicit points \((a:b)\) on it using PP1![a,b].

Hint: you can use the where keyword to break complicated expressions down into a sequence of more manageable sub-expressions. For example, the following Magma code defines a function that computes a root of a quadratic polynomial (assuming it has one):

qroot := func< quadratic | (-b + Sqrt(Delta))/(2*a) 
                           where Delta is (b^2 - 4*a*c) 
                           where a is Coefficient(quadratic,2)
                           where b is Coefficient(quadratic,1)
                           where c is Coefficient(quadratic,0) >;

You can test your code as follows (this will run without any problems if your code is correct):

load "Montgomery.m";
k := GF(1009);
A := Random(k);
B := k!1;
C := MontgomeryModel(A,B);
E, E_to_C := WeierstrassModel(C);
x := xMap(C);
xDBL := xDBLOperation(C);

for trial in [1..100] do
    P_E := Random(E);
    P_C := E_to_C(P_E);
    xP := x(P_C);
    x2P := xDBL(xP);
    assert x2P eq x(E_to_C(2*P_E));
end for;

(If your code is correct, then the above will run quietly without any problems...)

1.5 Pseudo-addition

Now we will define pseudo-addition, which is slightly more complicated than pseudo-doubling. As before, we use Magma's func< ... | ... > constructor (though this time it will have more than one argument).

Recall that if \(x(P) = (X_P:Z_P)\), \(x(Q) = (X_Q:Z_Q)\), \(x(P\oplus Q) = (X_\oplus: Z_\oplus)\), and \(x(P\ominus Q) = (X_\ominus: Z_\ominus)\), then \[(X_\oplus: Z_\oplus) = (Z_\ominus\cdot(U+V)^2:X_\ominus\cdot(U-V)^2)\] where \(U = (X_P-Z_P)(X_Q+Z_Q)\) and \(V = (X_P+Z_P)(X_Q-Z_Q)\).

Extend your MontgomeryModel function so that it

  1. computes a function xADD(xP,xQ,xD) which takes three points xP, xQ, and xD on \(\mathbb{P}^1\), where we assume xP \(= x(P)\), xQ \(= x(Q)\), and xD \(= x(P-Q)\) for some \(P\) and \(Q\) on \(C\), and returns xS \(= x(P+Q)\); and
  2. stores that function in the attribute xADDOperation of \(C\) (we do this bit for you).

The relevant lines of Montgomery.m are:

    // EXERCISE 1.5: Pseudo-addition
    /* uncomment and fill in the following lines: */
    // xADD := func< xP, xQ, xD | /* FILL IN IMAGE HERE */ >;
    // C`xADDOperation := xADD;  // set associated attribute

Hint: you should definitely use where here, both to simplify your code and to make the evaluation more efficient by re-using common subexpressions!

You can test your code as follows:

load "Montgomery.m";
k := GF(1009);
A := Random(k);
B := k!1;
C := MontgomeryModel(A,B);
E, E_to_C := WeierstrassModel(C);
x := xMap(C);
xADD := xADDOperation(C);

for trial in [1..100] do
    P_E := Random(E);
    Q_E := Random(E);
    xP := x(E_to_C(P_E));
    xQ := x(E_to_C(Q_E));
    if xP eq xQ then // we should use xDBL, not xADD, in this case
        continue;
    end if;
    xD := x(E_to_C(P_E - Q_E));
    xS := xADD(xP,xQ,xD);
    assert xS eq x(E_to_C(P_E + Q_E));
end for;

(If your code is correct, this should run quietly without any problems.)

2 Scalar multiplication on Montgomery models

Now we can implement scalar multiplication for Montgomery curves \(C\). Our basic tool is the Montgomery ladder, which implements the pseudoscalar multiplication on \(\mathbb{P}^1\) associated with \(C\). We can extend this to full scalar multiplication on \(C\) using Okeya and Sakurai's \(y\)-coordinate recovery trick.

2.1 The Montgomery ladder

First, recall Montgomery's ladder algorithm (here in Python):

def ladder(m,xP):
    # The Montgomery laddder: compute x([m]P) given x(P) and m
    # Assumes m is positive
    reg_0 = [1,0]  # Image on PP^1 of pt at infinity
    reg_1 = xP
    for b in bits(m): # from most significant down to least significant bit
        if b == 0:
            (reg_0,reg_1) = (xDBL(reg_0),xADD(reg_0,reg_1,xP))
        else:
            (reg_0,reg_1) = (xADD(reg_0,reg_1,xP),xDBL(reg_1))
    end for
    return reg_0

Define a new function Ladder(C,m,xP) implementing this algorithm in Montgomery.m. Your function should take a Montgomery model \(C\), a positive integer \(m\), and a point xP \(= x(P)\) in \(\mathbb{P}^1(k)\) for some \(P\) in \(C(k)\), and return the point \(x([m]P)\) in \(\mathbb{P}^1(k)\).

Hint: to get a sequence containing the bits of \(m\) in the right order, from most to least significant, you can use Reverse(IntegerToSequence(m,2)).

Test your code using the following examples:

> C := Curve25519; // curve defined above
> x := xMap(C);
> P := C![9,Sqrt(k!9*(9^2 + 486662*9 + 1)),1]; // base point
> x(P);
(9 : 1)
> Ladder(C,1,x(P));
(9 : 1)
> Ladder(C,2,x(P));
(14847277145635483483963372537557091634710985132825781088887140890597596352251 : 1)
> Ladder(C,101,x(P));
(15872060397774487147062612687452767107535139408010832422811284850336798704295 : 1)
> r := 2^252 + 27742317777372353535851937790883648493;
> Ladder(C,r-1,x(P));
(9 : 1)
> Ladder(C,r,x(P));  
(1 : 0)

You can also generate random examples as follows:

C := Curve25519; // as above
E, E_to_C := WeierstrassModel(C);
for trial in [1..1000] do
    pt_E := Random(E);
    m := Random(4*r);
    assert Ladder(C,m,x(E_to_C(pt_E))) eq x(E_to_C(m*pt_E));
end for;

(If your code is correct, this should run quietly without any problems.)

2.2 Optional: y-coordinate recovery

Note: This part is not required for the key exchange implementation later, so you can skip over it and come back later if you get stuck.

Recall that the Montgomery ladder maintains two "register" variables: the first finally contains \(x([m]P)\), while the second contains \(x([m+1]P)\). Okeya and Sakurai, following an idea of Lopéz and Dahab, showed that if \(P = (x_P:y_P:1)\) and \(Q = (x_Q:y_Q:1)\) are points on \(C\) with \(y_P \neq 0\neq y_Q\), and \(P + Q = (x_{P+Q}:y_{P+Q}:1)\), then \[2B\cdot y_P\cdot y_Q = (x_Px_Q + 1)(x_P+x_Q + 2A) - 2A - (x_P - x_Q)^2\cdot x_{P+Q}.\] Note that \(y_{P+Q}\) never appears in this formula! The idea is to use this formula to recover \(y_Q\) given \(x_P\), \(y_P\), \(x_Q\), and \(x_{P+Q}\). (We are not really interested in the special case where \(y_P = 0\) or \(y_Q = 0\).)

If we take \(Q = [m]P\) in the above, then we can use it to recover \(y_Q = y_{[m]P}\), and hence the full point \([m]P\), given only \(P\), \(x([m]P)\), and \(x([m+1]P)\).

Now,

  1. extend your Ladder function to return the value of both of its "register" variables; and
  2. define a new function ScalarMultiply for \(C\) that computes scalar multiplication on the curve \(C\) using the new version of Ladder as a subroutine.

We suggest something like this:

    function ScalarMultiply(C,m,P)
        // Computes [m]P on the Montgomery model C.
        // Assumes m is positive.
        
        // 1. Compute xP = x(P)
        /* complete here */
    
        // 2. Apply the ladder
        xmP, xm_plus_1P := Ladder(C,m,xP);
    
        // 3. Recover the correct preimage of [m]P:
        X_mP := xmP[1];
        Z_mP := xmP[2];
        Y_mP := /* complete here */ ;
    
        return C![X_mP,Y_mP,Z_mP];
    end function;

Test your code:

C := Curve25519; // as above
E, E_to_C := WeierstrassModel(C);
for trial in [1..1000] do
    pt_E := Random(E);
    m := Random(4*r);
    assert ScalarMultiply(C,m,E_to_C(pt_E)) eq E_to_C(m*pt_E);
end for;

(If your code is correct, this should run quietly without any problems.)

2.3 Branch-free code

For cryptographic applications, it is important that the runtime execution of algorithms is always independent of secret values (otherwise we are even more vulnerable to simple side-channel attacks). In particular,

In the context of Diffie-Hellman, the secrets are the scalars that we use in calls to Ladder, so we must remove the if statement from the main loop (which branches depending on bits of the scalar).

  1. Write a conditional-swap function CSWAP(bit,val_0,val_1) which returns val_0,val_1 if bit is 0 and val_1,val_0 if bit is 1.
  2. Now replace the if statement in your Ladder function with calls to CSWAP to add a layer of side-channel protection to your code.

Don't forget to test your code using the same examples as in Exercise 2.1 above.

3 X25519 key exchange

Now that we have got a working version of Montgomery arithmetic, we can use it to implement a serious cryptographic key-exchange function: X25519, which is defined in RFC7748.

Before going any further, download the file X25519.m, which contains some more skeleton code.

3.1 Compression and decompression

So far we have been applying the ladder to elements of \(\mathbb{P}^1(k)\). For key exchange we want to avoid the ambiguity (and extra space requirements) of projective points by mapping each to a single field element, which we then represent as a 255-bit integer (ie, a 32-byte value with the most significant bit masked to 0).

In X25519 we follow Bernstein's suggestion:

If we are working over \(\mathbb{F}_p\), then we can compute all of this in "constant time" via \((X:Z) \mapsto XZ^{(p-2)}\) in \(\mathbb{F}_p\); we can then get a 256-bit integer representation using Integers( )!value, where value is the result of the computation in \(\mathbb{F}_p\).

Going in the opposite direction, we map any element \(x\) of \(\mathbb{F}_p\) to the point \((x:1)\) in \(\mathbb{P}^1(\mathbb{F}_p)\).

Open X25519.m and complete the functions

3.2 Scalar clamping

The scalars used in X25519 have some special properties to ensure that the scalar multiplication algorithm runs in a regular way, and to avoid some possible attacks involving the use of small-order points.

To convert a random integer to an X25519 scalar, we use a procedure that has become known as clamping:

  1. the value is truncated (or filled with zeroes if necessary) to a 256-bit value \(x_0x_1\ldots x_{255}\);
  2. bit 255 is masked to 0;
  3. bit 254 is set to 1;
  4. bits 0, 1, and 2 are masked to 0.

This means that the set of all X25519 scalars is \(\{ 2^{254} + 8x : x \in [0,2^{251}) \}\).

Complete the function Clamped(x) in X25519.m, which takes any integer \(x\) and returns its "clamped" value as defined above.

Hint: in Magma, this will be much easier to do using integer operations than by working on individual bits.

3.3 The X25519 function

You can now define the X25519 function, ready to carry out key exchange:

function X25519(m,u)
    // The X25519 function, as defined in RFC7748.
    // The arguments m and u are non-negative integers
    // to be interpreted as 256-bit values (truncating if necessary).
    
    // Step 1: Clamp m
    m := Clamped(m);

    // Step 2: Decompress the point

    xP := Decompressed(u);

    // Step 3: Apply the ladder
    xmP := Ladder(Curve25519,m,xP);

    // Step 4: Compress the result
    v := Compressed(xmP);

    return v;
end function;

We are going to test your code using the test vectors from RFC7748.

First,

> test_scalar := 0xc49a44ba44226a50185afcc10a4c1462dd5e46824b15163b9d7c52f06be346a5
> test_base := 0x4c1cabd0a603a9103b35b326ec2466727c5fb124a4c19435db3030586768dbe6
> expected := 0x5285a2775507b454f7711c4903cfec324f088df24dea948e90c6e99d3755dac3
> X25519(test_scalar,test_base) eq expected;
true

The second kind of test involves iterating X25519 calls.

    function iterated_test(n)
        m := 9;
        u := 9;
        for i in [1..n] do
            new_m := X25519(m,u);
            u := m;
            m := new_m;
        end for;
        return m;
    end function;

What should happen is this:

> expected_1 := 0x7930ae1103e8603c784b85b67bb897789f27b72b3e0b35a1bcd727627a8e2c42;
> iterated_test(1) eq expected_1;
true
> expected_1000 := 0x684cf59ba83309552800ef566f2f4d3c1c3887c49360e3875f2eb94d99532c51;
> iterated_test(1000) eq expected_1000;
true

There is a further standard test with 1000000 iterations, but given the overhead of the Magma function calls, this will take a long time to run! But if you have an hour or two to spare and some cycles to burn, then you can check that iterated_test(1000000) eq 0x2454664fd2d24d5fdf303c88c001c63b6f5e577e29974486fd8625abe011397c.

3.4 Key exchange

To complete a Diffie-Hellman key exchange using X25519, we proceed as follows. First, we work in the subgroup of Curve25519 generated by a point with \(x\)-coordinate 9: that is, \(9\) is the compressed form of the public "generator".

Now

Check your code by simulating this example:

Alice's private key, a:

0x2a2cb91da5fb77b12a99c0eb872f4cdf4566b25172c1163c7da518730a6d0777

Alice's public key, X25519(a, 9):

0x6a4e9baa8ea9a4ebf41a38260d3abf0d5af73eb4dc7d8b7454a7308909f02085

Bob's private key, b:

0xebe088ff278b2f1cfdb6182629b13b6fe60e80838b7fe1794b8a4a627e08ab5d

Bob's public key, X25519(b, 9):

0x4f2b886f147efcad4d67785bc843833f3735e4ecc2615bd3b4c17d7b7ddb9ede

Their shared secret, K:

0x4217161e3c9bf076339ed147c9217ee0250f3580f43b8e72e12dcea45b9d5d4a