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:
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
Curve
, which returns an object of type CrvEll
;EllipticCurve
) and the isomorphism \(i: E \to C\);xADD
operation for \(C\); andxDBL
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
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
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
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\), andxDBLOperation
(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
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)\); andxADDOperation
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,
Ladder
function to return the value of both of its "register" variables; andScalarMultiply
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,
if
statements) on bits of secrets, andIn 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).
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.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:
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:
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
K_A = X25519(a,9)
, which he then publishes.K_B = X25519(b,9)
.K = X25519(a,K_B)
and Bob computes K = X25519(b,K_A)
; the result is their shared secret.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