In Zagier's one-sentence proof, why is S defined to be {(x,y,z)∈ℕ^3:x^2+4yz=p,p prime}? I've looked at a very clear explanation of Zagier's proof (specifically, it can be found here:http://danielmath.wordpress.com/2012/12/26/one-sentence-proof/) but the first step still eludes me: why is S defined to be {(x,y,z)∈ℕ^3:x^2+4yz=p,p prime}? How does one arrive at that definition of S from the equation p=4k+1? I do not see where the three variables x,y,z came from, why they have to be a part of ℕ that is cubed,such that x^2+4yz=p - in essence, I really have no idea how this definition of S came to be.
Many, many thanks! 
 A: He wants to solve $p=x^2+y^2$.
If there is a solution, $p$ is odd so one of the squares has to be even.  So instead, look at $p=x^2+(2y)^2=x^2+4yy$.
The clever bit is to look at $p=x^2+4yz$ instead of $p=x^2+4yy$.  Most solutions come in pairs, because $x^2+4yz=x^2+4zy$.  Only if $y=z$ does a solution not have a partner.
He counts the number of solutions, and because of the other formula finds there is an odd number.
So they don't all come in pairs, so there is a solution with $y=z$ - that is, $p=x^2+(2y)^2$
A: Michael's response is as clear as you will get without getting too involved.  But still, in "the clever bit", why not write $p=xy+4z^2$ or even $p=xy + 4 wz$?  And in Zagier's original proof, although you don't mention it, the complicated piecewise-defined involution should be at least as mysterious as consideration of the set S.
More understanding can be gained through several different perspectives.  Many are found in Elsholtz's survey here:
http://www.math.tugraz.at/~elsholtz/WWW/papers/papers30nathanson-new-address3.pdf
(He also gives a beautiful alternative one-sentence proof that uses geometry in the torus $\mathbb{Z}_p \times \mathbb{Z}_p$.) He shows, in particular, that if one searches for a linear involution to accomplish what Zagier's proof does, then Zagier's is the unique one with suitably "nice" properties.
Some of Elsholtz's survey is drawn from work of Dijkstra:
http://www.cs.utexas.edu/users/EWD/ewd11xx/EWD1154.PDF
On replacing the relation $p=x^2 + 4y^2$ by $p=x^2 + 4yz$, Dijkstra says: "we establish a one-to-one correspondence between the solutions of $p=x^2 + 4y^2$ and the fixed points of an involution. For the construction of that involutions, we do something with which every computer scientist is very familiar: replacing in a target relation -- here $p=x^2 + 4y^2$ -- something by a fresh variable”. 
A perspective in terms of binary quadratic forms can be found here:
http://www.rzuser.uni-heidelberg.de/~hb3/publ/bf.pdf
where we see $p=x^2+4yz$ arise naturally when fixing the discriminant of a general binary quadratic form.
