How do I multiply a Twisted Edwards point in Montgomery space?
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-&gt;Z, q-&gt;Y);  fe_mul(z1, z1, q-&gt;X);  fe_invert(z1, z1);
fe_add(t1, q-&gt;Z, q-&gt;Y);
fe_mul(x1, q-&gt;X, t1  );  fe_mul(x1, x1, z1);
fe_mul(y1, q-&gt;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-&gt;X, x1, t2);
fe_mul(p-&gt;Y, y1, t1);
fe_mul(p-&gt;Z, y1, t2);
fe_mul(p-&gt;T, x1, t1);