An elliptic curve is made up of the following:
- A field, F_p.
- What this means is all arithmetic on scalars is modulo p.
- Modern ECC have p as some large prime.
- For curve25519, p = 2^255 - 19, a prime.
- An equation and it's parameters:
- Curve25519 chose to use a Montgomery curve for performance reasons.
- Montogomery curves have the form B y^2 = x^3 + A x^2 + x, where x and y are members of F_p
- For curve25519, A = 486662, B = 1
- A base point / A generator.
- The generator is what is used to derive public and private keys.
- For curve25519, the generator is the point on the curve for which x = 9.
There are security considerations for the above selection of parameters, but since we're implementors we don't care about any of that.
An elliptic curve has the following properties:
- There exist points on the curve, P = (x, y). (x, y) MUST statisfy the curve equation above, and x and y MUST be members of the field F_p.
- Elliptic curves has the property where you have an operation where you can add points together: In other words, P + Q = R. This IS NOT simply R.x = P.x + Q.x.
- Scalar multiplication of a point is implemented through repeated addition. So 3 P = P + P + P.
The security of an elliptic curve relies on the conjecture that inverting scalar multiplication of points is hard.
If I needed to look up the addition formulae of a curve, I would take a look at https://hyperelliptic.org/EFD/. There are may be several different curve representations I could use.
Since our curve is a Montgomery curve, we look that up here: https://hyperelliptic.org/EFD/g1p/auto-montgom.html
We see that the naive implementation of P + Q = R is:
R.x = b*(Q.y-P.y)^2 / (Q.x-P.x)^2 - a - P.x - Q.x
R.y = (2*P.x+Q.x+a) * (Q.y-P.y) / (Q.x-P.x) - b*(Q.y-P.y)^3 / (Q.x-P.x)^3 - P.y
This won't be the formulas we use as faster formulas exist.
One such formula is the (X, Z) representation. In this representation, we can retrieve our original (x, y) coordinates by X/Z = x. The y coordinate can be recalculated from x and the curve equation, but curve25519 was designed so that the y coordinate is unnecessary and can be ignored. Equations for this representation can be found here: https://hyperelliptic.org/EFD/g1p/auto-montgom-xz.html#diffadd-dadd-1987-m-3 . Almost all curve25519 implementations will use this representation.
This is the formula most implementations use: https://hyperelliptic.org/EFD/g1p/auto-montgom-xz.html#ladder-ladd-1987-m-3
Addition and doubling:
A = X2+Z2
AA = A^22
B = X2-Z2
BB = B^2
E = AA-BB
C = X3+Z3
D = X3-Z3
DA = D*A
CB = C*B
X5 = Z1*(DA+CB)^2
Z5 = X1*(DA-CB)^2
X4 = AA*BB
Z4 = E*(BB+a24*E)
An example implementation from libsodium:
// (X2, Z2) = P
// (X3, Z3) = Q
fe_sub(tmp0,x3,z3); // D = X3 - Z3
fe_sub(tmp1,x2,z2); // B = X2 - Z2
fe_add(x2,x2,z2); // A = X2 + Z2
fe_add(z2,x3,z3); // C = X3 + Z3
fe_mul(z3,tmp0,x2); // DA = D * A
fe_mul(z2,z2,tmp1); // CB = C * B
fe_sq(tmp0,tmp1); // BB = B^2
fe_sq(tmp1,x2); // AA = A^2
fe_add(x3,z3,z2); // DA + CB
fe_sub(z2,z3,z2); // DA - CB
fe_mul(x2,tmp1,tmp0); // X4 = AA * BB (output)
fe_sub(tmp1,tmp1,tmp0); // E = AA - BB
fe_sq(z2,z2); // X5 = Z1 * (DA - CB) ^ 2 (output; Z1 = 1)
fe_mul121666(z3,tmp1); // a24 * E
fe_sq(x3,x3); // (DA + CB) ^ 2
fe_add(tmp0,tmp0,z3); // BB + a24 * E
fe_mul(z3,x1,z2); // Z5 = X1 * (DA - CB) ^ 2 (output)
fe_mul(z2,tmp1,tmp0); // Z4 = E * (BB + a24 * E) (output)
// 2P = (X4, Z4)
// P+Q = (X5, Z5)
Another example would be agl's implementation: https://github.com/agl/curve25519-donna/blob/master/curve25519-donna-c64.c#L293-L324
Scalar multiplication is implemented by repeated addition and doublings.
The idea behind this is that by breaking down the scalar into its bits all you need to do is double P enough times and add those doublings together. For example, if you want to multiply P by 13: 13 * P = (8 + 4 + 1) * P = 8 * P + 4 * P + 1 * P.
Say you want to multiply the point represented by x = x1 by the scalar s:
X2, Z2, X3, Z3 = 1, 0, X1, 1
for i = 254 ..0:
bit = (s >> i) & 1
PX = cmov(X2, X3, bit)
PZ = cmov(Z2, Z3, bit)
QX = cmov(X2, X3, !bit)
QZ = cmov(Z2, Z3, !bit)
X4, Z4 = point_doubling(PX, PZ)
X5, Z5 = point_addition(PX, PZ, QX, QZ, X1)
X2 = cmov(X4, X5, bit)
Z2 = cmov(Z4, Z5, bit)
X3 = cmov(X4, X5, !bit)
Z3 = cmov(Z4, Z5, !bit)
return X2 / Z2
where cmov is a constant-time conditional move. Most implementations use a constant-time swap (cswap) instead of a cmov here.
CONSTANT-TIME IMPLEMENTATION CONSIDERATION:
djb requires private keys to always have the 255th bit set in order to allow for the above constant-time implementation. ENSURE ALL PRIVATE KEYS HAVE THIS BIT SET.
IMPLEMENTATION NOTE:
The division X2 / Z2 is normally implemented as X2 * Z2^(p-2).
public key = private key * B
where B is the point on the curve for which B.x = 9
You need to implement a scalar type to represent scalars modulo p = 2^255 - 19.
djb's implementation is rather performance-oriented in this aspect, and might be difficult to follow unless you've read his paper.
Two reference implementation that are easier to read and don't use assembly are below: