[curves] Another try at point compression

[Michael Hamburg]

Better late than never.

I’m still working on this, but I think the results are now basically acceptable.  Unfortunately, the code requires constant-time conditional select and negate, but at least that way it works with 5-mod-8 primes.  (And with Curve/Ed25519, but it only reduces the cofactor to 2 in that case).

TLDR:
* Algo basically works.
* All operations tested in SAGE, both with high-level code (no z coordinates, cost is no object) and low-level code.
* Uses no more “expensive” field operations than any other compression scheme.
* Could still be simpler.
* Needs checking to make sure there are no special cases.
* Should probably check if inverting one of the coordinates simplifies things at all.

#################
# The Jacobi Quartic
#################

Define the modified Jacobi Quartic form as

s^2 + (2-4d) + 1/s^2 = t^2

Identity: (inf,inf)
Negation map: (s,t) <-> (s,-t)
2-torsion: (s,t) -> (s,t), (-s,t), (1/s,-t), (-1/s,-t)

A point on J is encoded as s such that s is even and finite, and t is even or infinite.  (We could replace t with 1/t here to keep it finite, but I don’t think that actually makes things simpler.)  The encoding of the identity is 0.

###############
# Edwards curves
###############

J is 2-isogenous to the Edwards curve E_d by:
x = 2 / (s + 1/s)
y = (1/s - s) / t

and to the Twisted Edwards curve E_(-1,d-1) by:
x = 2 / (1/s - s)
y = (s + 1/s) / t

To encode from these curves, invert the isogeny.  Define local_d = d-1 on the twisted curve, and let a=+-1 be the curve’s twist coefficient.  Compute:

sy = sqrt(a*local_d-1) / sqrt(a*(y^2-1)) # sqrt(a*local_d-1) is constant
jac_t = -2*a*sy
if is_odd(jac_t):
    sy = -sy
return make_even(a/x + y*sy)

Surprisingly, if the input point is given in (X:Y:Z:T) extended homogeneous coordinates, this can be done without the “inverse square root trick", though of course it still requires an inverse square root to compute sy.  This takes something like S + 7M + isqrt.

Decoding to affine requires the isqrt trick to compute 1 / t and 1 / (1 + as^2) as a batch.  It costs something like 4S + 9M + isqrt.  Then again, Edwards decompression is usually on the complicated side and requires the isqrt trick, so we’re not that far above par.  Decoding to projective or extended is a little simpler.

To test point equality, check if X1:Y1 == X2:Y2, i.e. if Y1 * X2 == X2 * Y1.

There was some concern (I think from RR?) about batching signatures or ephemeral generation.  This can be done, but only if you’re encoding 2P instead of P.  You have to use the dual isogeny instead of the inverse isogeny.  I haven’t coded up the formula for this, though.

##################
# Montgomery curves
##################

The Montgomery case is more complicated, but not absolutely terrible.  It requires the odd-ladder trick, but it doesn’t require the scalar to actually be odd.  Store at every step sqrt(odd) : sqrt(oddz), where odd:oddz is the odd side of what came out of the ladder.  It’s computed as an intermediate in the ladder operation, so you just have to save it.  The encoding algorithm also needs to know whether the scalar is odd or not.

Write the Montgomery curve as u + 1/u + (2-4d) = t^2.  It’s 2-isogenous to the Edwards curve by u = s^2.

The identity we’ll use this that if the ladder state is u:v:w with W=U+V, then y_u = (1 + vw - uv - uw) / sqrt(4uvw).  It could also be -that, of course, but if you compute y_u and y_v this way (and y_w negatively) then they are consistent.

At the end of the ladder, you have the base sqrt(u) and u; V:Z_V and W:Z_W, and sqrt(odd):sqrt(oddz), where either V,Z_V or W,Z_W = odd^2, oddz^2.  Let’s say you’ve done the final cswap already, so the target point is V:Z_V.

So if you begin by computing R := 1 / sqrt(4 * u * V * Z_V * W * Z_W), and compute y_u, y_v with the above formula.  Then you can tell if you want sqrt(v) or 1/sqrt(v) by checking if y_u and y_v have the same parity.

Finally, you need sqrt(v) = sqrt(V / Z_V) or its reciprocal.  You don’t need another sqrt operation to get this.  Note that 2R * sqrt(u) = 1/sqrt(V * Z_V * W * Z_W).  So if the scalar is even, then odd,oddz = sqrt(W,Z_W) and 2R*sqrt(u)*(V or Z_V)*odd*oddz is the answer.  If it’s odd, then odd,oddz = sqrt(V,Z_V), and (2R*sqrt(u))^2 * R*Z_R * (V or Z_V)*odd*odd_z is the answer.

Overall, this is a pretty complicated algorithm, but the total cost comes out to something like 1S + 14M + isqrt + two conditional selects.  So it’s at least not horribly complicated.  But it does require some finesse to avoid breaking the batched division in the case that w = 0 or 1/0, which can happen.  (If u or v is 1 or 1/0, then the answer is 0 anyway.)

Cheers,
— Mike