Solving for $z^2 = x^2 -xy + y^2$ Recently, I came across the following solution to finding integer solutions for $z^2 = x^2 - xy + y^2$: 


*

*$x = k(-n^2 -2mn)$

*$y = k(m^2 - n^2)$

*$z = k(mn + m^2 + n^2)$


I've been scratching my head trying to figure out how to derive this solution.
I initially thought that this was derived in the same way as the solution to the Pythagorean Triples based on:
$$(p^2 + q^2)^2 = (2pq)^2 + (p^2 - q^2)^2$$
So that, I get:
$(p^2 + q^2)^2 = z^2 + (2pq)(p^2 - q^2) = (2pq)^2 + (p^2 -q^2)^2$
But I don't see how to solve for $z$ based on the above.
Can anyone provide me the outline of how one figures out the general solution for $z^2 = x^2 -xy + y^2$?

Edit:  Added clarification about solutions in integers

Here is a way to find a parametrization.
First, notice that an integer solution means a rational solution $1=\left(\frac{x}{z}\right)^2-\frac{x}{z}\cdot\frac{y}{z}+\left(\frac{y}{z}\right)^2$ (when $z=0$ the only possible solution is $x=y=z=0$).  So, we can reduce this problem to looking for rational points on the ellipse $1=a^2-ab+b^2$.
The point $(1,0)$ is on the ellipse, and every rational point $(a,b)$ determines a line through $(a,b)$ and $(1,0)$ with rational slope.  If we solve the system $b-1=ma$ and $1=a^2-ab+b^2$ for $a$ in terms of $m$, we get $a=\frac{1-2m}{m^2-m+1}$ (with $m=\frac{1}{2}$ giving $a=0$), and from the equation of the line, $b=\frac{1-m^2}{m^2-m+1}$.  So a rational slope also determines a rational $(a,b)$.
When $m=\frac{p}{q}$ with integers $p,q$ we have $a=\frac{q^2-2pq}{p^2-pq+q^2}$ and $b=\frac{q^2-p^2}{p^2-pq+q^2}$.  With this, we can write $x=k(q^2-2pq)$, $y=k(q^2-p^2)$, and $z=k(p^2-pq+q^2)$ for arbitrary integers $k,p,q$ for solutions to $z^2=x^2-xy+y^2$.  The $k$ comes from the fact we may scale solutions.
This is very similar to your parametrization.  Perhaps if we selected a different point for the rational-sloped lines we would find yours.
