It is easy to see that if $(z,u,w)$ is a primitive pythagorean triple, then $w$ is odd and either $z$ is odd and $u$ is even or $z$ is even and $u$ is odd. WLOG suppose that $u=2x$ is even, so $z$ is odd. Then
$$u^2=4x^2=w^2-z^2=(w+z)(w-z).$$
It is plain that $\gcd(w+z,w-z)=2$, so $y_1:=(w+z)/2$ and $y_2:=(w-z)/2$ are relatively prime positive integers. Thus
$$x^2=y_1y_2,$$
and therefore $y_1$ and $y_2$ are relatively prime divisors of $x^2$. This implies that $y_1$ and $y_2$ are perfect squares, so write $y_1=a^2$ and $y_2=b^2$. Hence,
$$u=2x=2ab,\quad z=y_1-y_2=a^2-b^2,\qquad\text{and}\qquad w=y_1+y_2=a^2+b^2.$$
This proves that every primitive triple has the form $(a^2-b^2,2ab,a^2+b^2)$, as wanted.
Of course, the condition $\gcd(a,b)=1$ is needed, otherwise $(a^2-b^2,2ab,a^2+b^2)$ won't be primitive.
The computation of the Pythagorean triples of the form $T_n=(a_n,a_n+1,c_n)$ is equivalent to the computation of some convergents of the continued fraction of $\sqrt{2}$: in particular
$$ [1;\underbrace{2,2,\ldots,2,2}_{2n\text{ times}}]=\frac{2a_n+1}{c_n} $$
where
$$ c_n = \frac{(1+\sqrt{2})^{2n+1}-(1-\sqrt{2})^{2n+1}}{2\sqrt{2}},\qquad 2a_n+1 =\frac{(1+\sqrt{2})^{2n+1}+(1-\sqrt{2})^{2n+1}}{2}=d_n $$
both fulfill the recurrence $\ell_{n+2}=6\ell_{n+1}-\ell_n$. They can be expressed in terms of $D_n$ and $D_{n+1}$, where
$$ D_n = (3+2\sqrt{2})^n+(3-2\sqrt{2})^n =\sigma^n+{\bar{\sigma}}^n$$
is the trace of the $n$-th power of a $2\times 2$ matrix. This sequence fulfills
$$ D_{2n} = D_n^2-2,\qquad D_{2n+1}=D_n D_{n+1}-6 \tag{R}$$
so the couple $(D_n,D_{n+1})$ can be computed by a repeat-squaring algorithm. A concrete example will hopefully clarify how. Let us assume to want to compute $D_{23}$ and $D_{24}$. The binary representation of $23$ is $10111_2$, so we compute the couples $(D_m,D_{m+1})$ for $m=1_2,10_2,101_2,1011_2$ and finally $10111_2$ via $(R)$.
$$ (D_1,D_2)=(6,34) $$
$$ (D_2,D_3)= (34,198)$$
$$ (D_5,D_6)=(6726,39202) $$
$$ (D_{11},D_{12})=(263672646,15367968024) $$
$$ (D_{23},D_{24})=(405211279147678086,2361744410637427202)$$
This gives us $c_{23}$ and $d_{23}$, thus $T_{23}$, with no more than $3\log_2(23)$ multiplications.
Best Answer
After discovering Find all integers satisfying $m^2=n_1^2+n_1n_2+n_2^2$, the paper Abdelalim, S., & Dyani, H. (2014). The Solution of the Diophantine Equation x2+ 3y2= z2. International Journal of Algebra, 8(15), 729–732. gives an easy answer:
Theorem 2.1 Let $E : x^2 + 3 y^2 = z^2$ diophantine equation and $(x, y, z) \in \mathbb{Z}^3$ with $x \wedge y = 1$, $y$ is even and $x z \wedge 3 = 1$. Then the following properties are equivalent:
Theorem 2.2 Let $E : x^2 + 3 y^2 = z^2$ diophantine equation and $(x, y, z) \in \mathbb{Z}^3$ with $x \wedge y = 1$, $y$ is odd and $x z \wedge 3 = 1$. Then the following properties are equivalent:
Scaling the above primitive solutions $(x, y, z)$, and in our case restricting to $2 \mid z$, gives all solutions.
The following is a method that is less explicit, as it requires dividing by a $\gcd$.
Thanks to Daniel Hast's comment, I arrived at the following:
Finding such $(a, b)$ is equivalent to finding a single rational point on the ellipse $x^2 + 3 y^2 = 4$, then finding all rational lines through it. A rational point is $(1, 1)$, and a rational line through that is $t y = s x - s + t$ where $\gcd (s, t) = 1$. The other intersection with the ellipse is at:
$$ (x, y) = \left( \frac{3 s^2 - 6 s t - t^2}{3 s^2 + t^2}, - \frac{3 s^2 + 2 s t - t^2}{3 s^2 + t^2} \right) $$
So:
$$ \begin{array}{rll} \left( \frac{3 s^2 - 6 s t - t^2}{3 s^2 + t^2} \right)^2 + 3 \left( \frac{3 s^2 + 2 s t - t^2}{3 s^2 + t^2} \right)^2 & = & 4\\ (3 s^2 - 6 s t - t^2)^2 + 3 (3 s^2 + 2 s t - t^2)^2 & = & 4 (3 s^2 + t^2)^2\\ 4 (3 s^2 + t^2)^2 - 3 (3 s^2 + 2 s t - t^2)^2 & = & (3 s^2 - 6 s t - t^2)^2 \end{array} $$
Note that the numerator and denominators of $x, y$ can be scaled while the rational point stays the same. Let $d = \gcd (3 s^2 + t^2, 3 s^2 + 2 s t - t^2)$. For $k \in \mathbb{Z}$, our solutions are then:
$$ a = \frac{k}{d} \cdot (3 s^2 + 2 s t - t^2) \quad b = \frac{k}{d} \cdot (3 s^2 + t^2) $$