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:**

- C. Costello and B. Smith,
*Montgomery curves and their arithmetic: the case of large characteristic fields.*; - A. Langley, M. Hamburg, and S. Turner, RFC7748:
*Elliptic Curves for Security*. - The slides from this morning's talk

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

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

; - the morphism \(x: C \to \mathbb{P}^1\) mapping \((X:Y:Z)\) to \((X:Z)\);
- 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\); - the x-line pseudo-addition
`xADD`

operation for \(C\); and - 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.*

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

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

- constructs the projective line over the same base field as the curve \(C\), then
- defines the map \(x: C \to \mathbb{P}^1\) sending \((X:Y:Z) \mapsto (X:Z)\), and
- 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?

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

- computes an isomorphism to a Weierstrass model for \(C\), and
- 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)
```

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

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

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

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

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.

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

*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,

**extend**your`Ladder`

function to return the value of*both*of its "register" variables; and**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.)

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,

- we cannot have any branching (
`if`

statements) on bits of secrets, and - we cannot have any array indexing (or memory access) based on bits of secrets.

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

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

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.

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:

- "affine" points \((X:Z)\) with \(Z \neq 0\) map to \(X/Z\);
- the point "at infinity", \((X:0)\), maps to \(0\).

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

`Compressed`

, mapping points in \(\mathbb{P}^1(\mathbb{F}_p)\) to 256-bit integers as above,`Decompressed`

, mapping 256-bit integers to elements of \(\mathbb{F}_p\) and then to \(\mathbb{P}^1(\mathbb{F}_p)\) as above.

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*:

- the value is truncated (or filled with zeroes if necessary) to a 256-bit value \(x_0x_1\ldots x_{255}\);
- bit 255 is masked to 0;
- bit 254 is set to 1;
- 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.

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`

.

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

- Alice chooses her secret \(a\), a 32-byte (256-bit) value, and computes his public key
`K_A = X25519(a,9)`

, which he then publishes. - Bob chooses his secret \(b\), a 32-byte (256-bit) value, and computes his public key
`K_B = X25519(b,9)`

. - Alice computes
`K = X25519(a,K_B)`

and Bob computes`K = X25519(b,K_A)`

; the result is their shared secret. - To convert K into a cryptographic secret key, Alice and Bob should both apply a
*key derivation function*to \(K\), \(K_A\), and \(K_B\).

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`