How do I multiply a Twisted Edwards point in Montgomery space?

Question

EdDSA (and ed25519) signatures require a scalar multiplication. Currently, I do this directly in Twisted Edwards space. (The code can be found in my crypto library.) My research and my tests suggest it would be quite a bit faster to do that multiplication in Montgomery space instead. This would mean 4 steps:

1. Convert the Twisted Edwards point to Montgomery space.
2. Perform the Montgomery ladder.
3. Recover the Y coordinate.
4. Convert back to Twisted Edwards space.

Steps 1, 2, and 4 seem pretty simple (even straightforward). Step 3 is more complex, but is documented in this paper, which I have followed. Here is my attempt, which unfortunately gives the wrong result:

typedef i32 fe[10]
typedef struct { fe X; fe Y; fe Z; fe T; } ge;

sv ge_scalarmult(ge *p, const ge *q, const u8 scalar[32])
{
    // convert q to montgomery format
    fe x1, y1, z1, x2, z2, x3, z3, t1, t2, t3, t4;
    fe_sub(z1, q->Z, q->Y);  fe_mul(z1, z1, q->X);  fe_invert(z1, z1);
    fe_add(t1, q->Z, q->Y);
    fe_mul(x1, q->X, t1  );  fe_mul(x1, x1, z1);
    fe_mul(y1, q->Z, t1  );  fe_mul(y1, y1, z1);
    fe_1(z1); // coherence

    // montgomery scalarmult
    // The same ladder is used for x25519,
    // it comes from the ref10 implementation.
    // Field elements are modified in-place
    x25519_ladder(x1, x2, z2, x3, z3, scalar);

    // recover the y1 coordinate
    fe_mul(t1, x1, z2);  // t1 = x1 * z2
    fe_add(t2, x2, t1);  // t2 = x2 + t1
    fe_sub(t3, x2, t1);  // t3 = x2 − t1
    fe_sq (t3, t3);      // t3 = t3^2
    fe_mul(t3, t3, x3);  // t3 = t3 * x3
    fe_mul973324(t1, z2);// t1 = 2a * z2
    fe_add(t2, t2, t1);  // t2 = t2 + t1
    fe_mul(t4, x1, x2);  // t4 = x1 * x2
    fe_add(t4, t4, z2);  // t4 = t4 + z2
    fe_mul(t2, t2, t4);  // t2 = t2 * t4
    fe_mul(t1, t1, z2);  // t1 = t1 * z2
    fe_sub(t2, t2, t1);  // t2 = t2 − t1
    fe_mul(t2, t2, z3);  // t2 = t2 * z3
    fe_add(t1, y1, y1);  // t1 = y1 + y1
    fe_mul(t1, t1, z2);  // t1 = t1 * z2
    fe_mul(t1, t1, z3);  // t1 = t1 * z3
    fe_mul(x1, t1, x2);  // x1 = t1 * x2
    fe_sub(y1, t2, t3);  // y1 = t2 − t3
    fe_mul(z1, t1, z2);  // z1 = t1 * z2

    // convert back to twisted edwards
    fe_sub(t1  , x1, z1);
    fe_add(t2  , x1, z1);
    fe_mul(p->X, x1, t2);
    fe_mul(p->Y, y1, t1);
    fe_mul(p->Z, y1, t2);
    fe_mul(p->T, x1, t1);
}

Can someone tell me what I did wrong?

Answer 1

I finally got the answer from Hacker News. My mistake came from the conversion code. I used a formula from Wikipedia, which is slightly different from the correct one. Simply put, I missed a multiplication by a constant (specifically, the square root of -486664, or 6853475219497561581579357271197624642482790079785650197046958215289687604742 —the other square root also works.)

The correct code for converting to Montgomery format is this (note the inversion, needed to ensure that z1 == 1, which allows a faster Montgomery ladder):

// sqrt(-486664), in 26/25 bit limbs.
fe K = { 54885894, 25242303, 55597453,  9067496, 51808079,
         33312638, 25456129, 14121551, 54921728,  3972023 };

// Actual conversion
fe x1, y1, z1, x2, z2, x3, z3, t1, t2, t3, t4;
fe_sub(z1, q->Z, q->Y);  fe_mul(z1, z1, q->X);  fe_invert(z1, z1);
fe_add(t1, q->Z, q->Y);
fe_mul(x1, q->X, t1  );  fe_mul(x1, x1, z1);
fe_mul(y1, q->Z, t1  );  fe_mul(y1, y1, z1);
fe_mul(y1, K, y1);  // missing multiplication
fe_1(z1);

(Even though we use projective coordinates, the Montgomery ladder is faster when it can assume z1 is equal to one. This compensates the cost of the inversion.)

The correct code for converting back to Twisted Edwards format is this:

fe_sub(t1  , x1, z1);
fe_add(t2  , x1, z1);
fe_mul(x1  , K , x1);  // missing multiplication
fe_mul(p->X, x1, t2);
fe_mul(p->Y, y1, t1);
fe_mul(p->Z, y1, t2);
fe_mul(p->T, x1, t1);

Answer 2

For Ed448, this issue has happened for me. I design my scheme to work in twisted Edward space, and all results are verified by given test vector in RFC 8032.

Now I wanna work over Montgomery space. The first three mentioned steps, i.e., convert the base point to Mont., Mont. ladder execution, and y-coordinate recovery can be performed simply. However, the back transformation has a problem. Based on the equation described in RFC 7748, the map between Mont. space and Ed space are as follows: \begin{multline} (x, y) = \frac{(4*v*(u^2 - 1)}{(u^4 - 2*u^2 + 4*v^2 + 1)},- \frac{(u^5 - 2*u^3 - 4*u*v^2 + u)}{(u^5 - 2*u^2*v^2 - 2*u^3 - 2*v^2 + u))} \end{multline}

This is my attempt (in sage) to achieve the point in Ed448-Goldilocks:

X2=mod(X*X,p)
Z2=mod(Z*Z,p)
Y2=mod(Y*Y,p)
X3=mod(X2*X,p)
Z3=mod(Z2*Z,p)
X4=mod(X2*X2,p)
Z4=mod(Z2*Z2,p)
X5=mod(X4*X,p)

a1=mod(4*(X2-Z2),p)
a1=mod(a1*Y,p)
a1=mod(a1*Z,p)

a2=Integer(mod(2*X2*Z2,p))
a2=Integer(mod(X4-a2,p))
a2=Integer(mod(a2+4*Y2*Z2,p))
a2=Integer(mod(a2+Z4,p))

b1=mod(4*X*Y2,p)
b1=mod(b1*Z2,p)
b1=mod(X5-b1,p)
b1=mod(b1-2*X3*Z2,p)
b1=mod(b1+X*Z4,p)
b1=mod(-b1,p)

b2=Integer(mod(2*X2*Y2,p))
b2=Integer(mod(b2*Z,p))
b2=Integer(mod(X5-b2,p))
b2=Integer(mod(b2-2*X3*Z2,p))
b2=Integer(mod(b2-2*Y2*Z3,p))
b2=Integer(mod(b2+X*Z4,p))

a2_inv=Integer(inverse_mod(a2,p))
x=Integer(mod(a1*a2_inv,p))

b2_inv=Integer(inverse_mod(b2,p))
y=Integer(mod(b1*b2_inv,p))

The results are not correct. what mistake I did?