The equation of an Edwards curve is $x^2 + y^2 = 1 + d x^2 y^2 \,$
The addition formula is $(x_1,y_1) + (x_2,y_2) = \left( \frac{x_1 y_2 + x_2 y_1}{1 + dx_1 x_2 y_1 y_2}, \frac{y_1 y_2 - x_1 x_2}{1 - dx_1 x_2 y_1 y_2} \right) \,$
What is the intuitive explanation that $dx_1 x_2 y_1 y_2$ can't be equal to ± 1 when d is a non-square? Because I don't see why you can't have $x_1 x_2 y_1 y_2 = ± \frac{1}{d}$
On
I just wanted to leave here a comment that an Edward's curve $x^2+y^2=1+dx^2y^2$ is isomorphic to a curve in Weierstrass form of the form: $$Y^2 = X^3 + (-4d - 4)X^2 - 64dX + 256d^2 + 256d.$$ I don't think this is particularly useful here, but I was curious to see a Weierstrass model... so here it is. The torsion subgroup is generically $\mathbb{Z}/4\mathbb{Z}$, generated by $(8 ,-16d + 16)$. If $d$ is a square, then we gain another point of order $2$, as Lubin points out in another response, and we get a torsion subgroup $\mathbb{Z}/2\mathbb{Z}\oplus\mathbb{Z}/4\mathbb{Z}$.
I’m not at all familiar with this representation of an elliptic curve, but I think that when you homogenize, to the equation $x^2z^2+y^2z^2=z^4+dx^2y^2$, things may get a bit clearer. The points $(1,0,0)$ and $(0,1,0)$ on the line at infinity are singular, seems to me, and even before writing the addition-rule projectively, I’ll wager that the sum of any two nonsingular points will again be nonsingular. But then I haven’t checked this.
EDIT: After looking more closely, I realize that the situation is complicated enough that it requires careful description. For typographical and esthetic reasons, I’m changing your description somewhat. The condition that $d$ be a nonsquare means that the points at infinity are not $k$-rational, where $k$ is the constant field (characteristic $\ne2$ necessarily). This makes sense arithmetically, but geometrically it does not. So, in my description of the curve, I’ll use the equation $$x^2+y^2=1+D^2x^2y^2$$ instead. It makes all the formulas simpler. I should say that I have the real picture in mind, with $D>1$, and this is how it looks:
Start with a U-shaped curve, open upwards, centered on the positive $y$-axis, with its vertex (lowest point) at $(0,1)$, and with asymptotes $x=\pm1/D$. Then take this, and rotate it $\pi/2$ repeatedly, till you have four of these, opening out in all four directions. This is your quartic curve.
For this curve, I want to take a different but isomorphic law of combination: \begin{align} (x,y)+(x',y')&=\left(\frac{xx'-yy'}{1-D^2xx'yy'},\frac{xy'+x'y}{1+D^2xx'yy'}\right)\\ \Bbb O&=(1,0)\\ -(x,y)&=(x,-y)\,, \end{align} so that the neutral point $\Bbb O$ is the vertex of the U-shaped curve that sits straddling the positive $x$-axis. I’ve made this change only because your original law had the sine addition formula in the numerator of the $x$-coordinate, and the negative of the cosine formula in the $y$-coordinate, and that annoyed me. We can call the other three vertices $\Bbb P_{\pi/2}$, $\Bbb P_\pi$, and $\Bbb P_{3\pi/2}$.
Now, when you look at the equation for the curve, you see that its highest-degree form is just $x^2y^2$, and this tells you how the curve intersects the line at infinity: at $(1,0,0)$ and $(0,1,0)$. But the exponents tell you that these points are two-fold intersections with that line. Indeed, they’re the intersections with the four asymptotes I mentioned above.
Very crudely, we can desingularize our curve by designating four points at infinity $\Bbb I^\rightarrow_+$, $\Bbb I^\rightarrow_-$, $\Bbb I^\uparrow_+$, and $\Bbb I^\uparrow_-$. The first of these, $\Bbb I^\rightarrow_+$, can be thought of as $(\infty,1/D)$ or maybe the limit of $(t,1/D)$ as $t\to\infty$, even though these points are not on the curve. In a similar way, $\Bbb I^\rightarrow_-$ is $(\infty,-1/D)$, and you can guess the others. Let’s see how this works out by seeing what the sum $\Bbb I^\rightarrow_++\Bbb I^\rightarrow_+$ will turn out to be. The formula gives $$ \left(\frac{t^2-1/D^2}{1-D^2t^2/D^2},\frac{2t/D}{1+D^2t^2/D^2}\right) \longrightarrow (-1,0)=\Bbb P_\pi\,. $$
You easily check, using the formulas, that $[2]\bigl(\Bbb P_\pi\bigr)=\Bbb O$, which shows that $\Bbb I^\rightarrow_+$ is a four-division point on our curve. Furthermore, you also see that $[2]\bigl(\Bbb P_{\pi/2}\bigr)=\Bbb P_\pi$.
Without working (!), we’ve found eight rational points on our curve: all division points, and since there are too many $4$-division points for a cyclic group, you see that the torsion is at least $C_2\oplus C_4$, where $C_n$ means cyclic of order $n$.
Now, finally, to answer your question about infinities, it’s clear. Any time your addition formula leads you to a point on the line at infinity, it’ll be of form $(\infty,\pm1/D)$ or $(\pm1/D,\infty)$, where here, “$\infty$” refers to an expression that blows up with zero in the denominator. But if $D=\sqrt d$ is not $k$-rational, then your formulas with coefficients in $k$ can not lead you to $1/\sqrt d$.