I'm trying to understand group operations on elliptic Montgomery curves in affine representation. Let's say the curve I use is Curve25519, i.e.:
$$y^2 = x^3 + A\,x^2 + x\quad\text{where}\quad A=486662 \quad\text{over the finite field}\quad \mathbb{F}_{2^{255}-19}$$
Now, let's take two points on this curve in affine representation: $$\begin{align} P &= ({\scriptstyle 9,\ 43114425171068552920764898935933967039370386198203806730763910166200978582548})\\ Q &= ({\scriptstyle 14847277145635483483963372537557091634710985132825781088887140890597596352251}, \\ & \qquad{\scriptstyle 48981431527428949880507557032295310859754924433568441600873610210018059225738}) \end{align}$$
Their sum, according to SAGE, should be $$\begin{align} P + Q &= (\scriptstyle 12697861248284385512127539163427099897745340918349830473877503196793995869202, \\ & \qquad{\scriptstyle 39113539887452079713994524130201898724087778094240617142109147539155741236674})\end{align}$$
Now I'm trying to calculate this $P + Q$ in affine space (projective coordinates lose the $y$ coordinate, which I happen to be interested in). Montgomery gives in his paper "Speeding the Pollard and Elliptic Curve Methods of Factorization" from 1987 the formula on page 19:
$$x_3 = \left(\frac{y_1-y_2}{x_1-x_2}\right)^2 - A - x_1 - x_2$$
This works great in SAGE:
x1 = F(9)
y1 = F(43114425171068552920764898935933967039370386198203806730763910166200978582548)
x2 = F(14847277145635483483963372537557091634710985132825781088887140890597596352251)
y2 = F(48981431527428949880507557032295310859754924433568441600873610210018059225738)
print(((y1 - y2) / (x1 - x2)) ** 2 - F(486662) - x1 - x2)
Which gives me the correct $x$ coordinate of $P+Q$. But how do I get the $y$ coordinate of the point addition? The Montgomery paper does not mention it and all references that I've found so far only use projective coordinates (where $y$ is discarded).
A solution I thought of was to recover $y$ from the calcuated $x$ by doing a square root -- but this would be ambiguous (give me two possible solutions) and I don't know which one I should accept as the correct one. Also I'd like to avoid the performance penalty for the Tonelli-Shanks square rooting operation.
Setting $x=X-\frac{A}{3}$ should transform the Montgomery form to short Weierstrass form. Then you can use the known Weierstrass formulae.
However, to save work, helpful people have prepared an explicit formulas database which might interest you, particularly this entry. It should give you $$y_3 = (2x_1+x_2+A)\frac{y_2-y_1}{x_2-x_1}-\left(\frac{y_2-y_1}{x_2-x_1}\right)^3-y_1$$ Combined with the equation for the $x_3$ coordinate that you have given, this can be simplified to $$y_3 = (x_1-x_3)\frac{y_2-y_1}{x_2-x_1}-y_1$$ This makes use of the $x_3$ computed before. The formula is plausible: You can easily interpret it as calculating the $y$ coordinate corresponding to the line through $P$ and $Q$ at $x=x_3$ and then flipping its sign.