Recover Y coordinate from xz elliptic curve multiplication

Question

I have an elliptice curve in the form

y² = x³ + ax + b (mod p)

And I have a multiplication algortihm which uses only x and z coordinate

How can I recover the Y coordinate ?

I tried to use the curve equation, but because of sqrt_mod, there's 2 possibles solutions and I don't know which one to choose

def xDBLADD(P, Q, PQ):
    (X1, Z1), (X2, Z2), (X3, Z3) = PQ, P, Q
    X4 = (X2 ** 2 - a * Z2 ** 2) ** 2 - 8 * b * X2 * Z2 ** 3
    Z4 = 4 * (X2 * Z2 * (X2 ** 2 + a * Z2 ** 2) + b * Z2 ** 4)
    X5 = Z1 * ((X2 * X3 - a * Z2 * Z3) ** 2 - 4 * b * Z2 * Z3 * (X2 * Z3 + X3 * Z2))
    Z5 = X1 * (X2 * Z3 - X3 * Z2) ** 2
    X4, Z4, X5, Z5 = (c % q for c in (X4, Z4, X5, Z5))
    return (X4, Z4), (X5, Z5)
def xMUL(P, k) -> int:  # use xz coordinate
    Q, R = (1, 0), P
    for i in reversed(range(k.bit_length() + 1)):
        if k >> i & 1:
            R, Q = Q, R
        Q, R = xDBLADD(Q, R, P)
        if k >> i & 1:
            R, Q = Q, R
    return Q[0]
def mul(P, k):
    Pz = (P[0], 1)
    Qz = xMUL(Pz, k)
    return Qz[0] * pow(Qz[1], -1, q) % q

Note: a similar question has been already posted here, but for Montgomery curve, so the equation is not the same

Answer 1

Let we have short Weierstraß form (see note) $$y^2 = x^3+ a_4x + a_6$$

If one wants to find the $y$ coordinate of $[n]P$, where $P=(x_1,y_1)$ is in affine coordinates the formula is

$$y_n = \frac{2a_6+ (x_1x_n+a_4)(x_1+x_n) - (x_1-x_n)^2x_{n+1}}{2y_1}$$

This is from section 13.2.3.b of Handbook of Elliptic and Hyperelliptic Curve Cryptography

---

Note

Let field $K$ we are working and this formula must work for any Elliptic curve with short Weierstraß form $$y^2 = x^3 + a_4x + a_6.$$. Short Weierstraß form is possible if the characteristic of the field $charK\neq 2$ and $charK \neq 3$ (details of conversion). Most of the time, the curves are selected to have $a=0$ which reduces the cost of doubling.

---

Validation with SageMath

The SageMath Code ( modified from the Roberts question history ) where $a = a_4$ and $b = a_6$

from sage.all import *
from random import randint
p = 115792089210356248762697446949407573530086143415290314195533631308867097853951
a = -3
b = 41058363725152142129326129780047268409114441015993725554835256314039467401291
E = EllipticCurve(Zmod(p), [a, b])
n = randint(0, p)
P1 = E.random_point()
Pn = P1 * n
Pn1 = P1 * (n + 1)
x1, y1 = P1.xy()
xn, yn = Pn.xy()
xn1, yn1 = Pn1.xy()
test = 2 * b + (x1 * xn + a) * (x1 + xn) - ((x1 - xn) ** 2) * xn1
test = pow(2y1, -1, p)
test %= p
print(test)
print(yn)
assert test == yn

No assertion errors and the outputs

26306650830132737739245713876108232680224938874883909749780413733622205404261
26306650830132737739245713876108232680224938874883909749780413733622205404261