Find Ed25519 Y Coordinate from X Coordinate

Question

Are there any available formulas to determine a $Y$ coordinate given only an $X$ for the Ed25519 curve? The closest thing I've come to find is the recover_x(y, sign) function in the EdDSA RFC at https://www.rfc-editor.org/rfc/rfc8032.html#section-6

---

Edit: Using @kelalaka 's steps, I implemented this python code, but it seems to be missing something as its incorrect:

p = 57896044618658097711785492504343953926634992332820282019728792003956564819949 # 2**255 -19
d=(121665 * pow(121666, -1, p)) # had to mod inverse parts of d

def get_y(x):
   a = 1 + x**2 # top line
   b = 1 + d*(x**2) # bottom line
   c = a * pow(b, -1, p) # calculation
   return math.isqrt(c) % p # finding y

f = G * 2

print(f.y) # this is correct
# Output is 15549675580280190176352668710449542251549572066445060580507079593062643049417
print(get_y(f.x)) # this is incorrect
# Output is 25882070986792651040465013186999776942933444728837758367190487200702647909051

Answer 1

1. Get the curve equation $$-x^2 + y^2 = 1 - \frac{121665}{121666}x^2y^2$$

2. Call $d = \frac{121665}{121666}$ $$-x^2 + y^2 = 1 - dx^2y^2$$

3. Take $y^2$ to left $$y^2 + dx^2y^2 = 1 + x^2$$

4. Take $y^2$ out of parenthesis $$y^2( 1+ dx^2) = 1 + x^2$$

5. Divide $$y^2 = \frac{1 + x^2}{( 1+ dx^2)}$$ this is possible since $1+ dx^2$ is not zero.

   take $-1 = dx^2$. Since $p \equiv 1 \pmod 4$ then $-1$ is a $QR$. $d$ is chosen to be non-square $(QNR)$ therefore equality is impossible.

6. Plug $x$ and $d$ to find $y^2$

7. Find $y$ by taking square roots modulo $2^{255}-19$. this will leave out one or two possible solutions.

Square root when modulus congruent to 5 modulo 8

Since $2^{255}-19 \equiv 5 \pmod{8} $ instead of Tonelli-Shanks we can use special formula for $r^2 = a$;

1. compute the square root of $-1$

   • $u \equiv (2a)^{(p-5)/8} \pmod{2^{255}-19}$
   • $i \equiv 2au^2 \pmod{2^{255}-19}$

2. compute the square root

   • $ r = \pm au(i-1)\pmod{2^{255}-19}$

---

Sagemath for calculation $Y$ given $X$ coordinate

p = 2^255-19

#d = 121665 * 121666
dd =  inverse_mod(121666, p)
d = mod(121665 * dd,p)

x = mod(15112221349535400772501151409588531511454012693041857206046113283949847762202,p)
y = mod(46316835694926478169428394003475163141307993866256225615783033603165251855960,p)

a = mod( (1 + x^2) / (1 +d*x^2), p)


v = mod((p-5)/8,p-1)
u = mod( (2*a)^(v) ,p) 
i = mod( 2*a*u^2, p)

yp = mod( a*u*(i-1),p)
yn = mod( -yp ,p)

print("yp= ", yp)
print("yn= ", yn)
print("y = ", y)

This code calculates for the base point, and the output is

yp=  46316835694926478169428394003475163141307993866256225615783033603165251855960
yn=  11579208923731619542357098500868790785326998466564056403945758400791312963989
y =  46316835694926478169428394003475163141307993866256225615783033603165251855960

The corrected Python code ( version >= 3.8 )

import math

def get_y(x,p,d):
        
   a = (1 + x**2 ) * pow((1 + d*(x**2) ), -1, p)
   
   %Modulo square root part
   v = (p-5)//8 % (p-1)

   u = pow( 2*a,v, p) 
   i = 2*a*u**2 % p

   yp = a*u*(i-1) %p
   yn = -yp %p

   return (yp,yn)# finding y


x = 15112221349535400772501151409588531511454012693041857206046113283949847762202

p = 2**255 -19

d =(121665 * pow(121666, -1, p)) % p# had to mod inverse parts of d

print("X = " , get_y(x,p,d))