secp256k1_gej_add_ge does the wrong thing with points P, Q with P_y = -Q_y but P≠-Q

[apoelstra]

It is possible, by using the value β from the endomorphism optimization, to obtain two points P = x : y : z and Q = βx : -y : z which are both on the curve. If you attempt to add these points using secp256k1_gej_add_ge, it incorrectly returns infinity, since it sees that the z-coordinates match and the y-coordinates are negative.

Specifically, it computes the value M = Y_1 * Z_2^3 + Y_2 * Z_1^3, and sets the returned z-coordinate to 2 * M * z_1z_2, and of course M will be zero for these two points. So it is not a simple matter of fixing how the infinity flag is set, the algebra actually does not compute the complete group law. This is derived in the paper "Eric Brier and Marc Joye, Weierstrass Elliptic Curves and Side-Channel Attacks" in Proposition 1, which explicitly assumes that the points P and Q don't happen.

The current secp256k1 API does not use this function except in one place, secp256k1_ecmult_gen, where it is impossible to get such an input to the function, so right now this bug is harmless. However #252 exposes it more directly, since when the endomorphism optimization is turned on roughly half of all scalars produce a bad pair P, Q in the wNAF ladder.

[sipa]

No other choice than to replace the function with one that computes (infinity, 2*P, P+Q) independently, and uses a conditional copy to return one of them.

[apoelstra]

Do you have an algebraic reason to believe there is no unified formula, or just that you are not aware of one?

[sipa]

I am not aware of one.

[apoelstra]

I have a borderline coherent algebraic reason that we can't find a unified formula: we want to calculate λ (the point-addition λ, not the one from the endomorphism optimization) in such a way that λ = f(y1, y2)/g(x1, x2) when x1 ≠ x2, but λ = h(x1, x2)/i(y1, y2) otherwise. That is, we are trying to unify these two equations, one of which is an x-function over a y-function and the other of which is the opposite.

Any unification scheme is therefore going to need to swap x and y somehow. The Brier/Joye paper uses the curve equation y^2 = x^3 + 7 to do this, and AFAICT this is the only way they could, since this is the only relationship we have between x and y.

However, with the offending pairs (x, y) and (βx, -y) we see that y^2 = (±y)^2 and x^3 = (βx)^3; that is, the curve equation is totally unable to distinguish these points. So since the original formula λ = (y1 - y2)/(x1 - x2) broke down for the pairs (x, y) and (x, y), no manipulation by means of the curve equation can get us a formula that never breaks down.

This isn't a very solid argument, but I think it indicates that @sipa is right and we'll have to find a way to use two equations and then cmov the final result in place. (Maybe there is an earlier point in the algebra that we can do this to avoid too much wasted work?) I might offer my calc students some cash for unifying the formula, but I can't justify spending any more of my own time on it.