It's been a while since I read Elsholtz's article, but after doing so I felt none the wiser. Below I have translated Heath-Brown's proof into the language of binary quadratic forms; Zagier's proof looks more interesting from this point of view (the connections to Gauss reduction are much closer), but when working out the details I got stuck in the middle.

One essential ingredient for the proofs by Heath-Brown and Zagier was pointed out already by Frick in 1918, who showed that if $p = a^2 + 4b^2$ is an odd prime number, then the indefinite binary quadratic form $Q = (-b,a,b)$ with discriminant $p$ is Gauss reduced and is contained in the principal cycle.

For proving that such a form exists without assuming that $p$ is a sum of two squares, we consider all forms $(A,B,C)$ with discriminant $p$ such that $A < 0$ and $C > 0$. From $p = B^2 - 4AC$ it then follows that the set
$$ S = \{(A,B,C): B^2 - 4AC = p, A < 0, C > 0\} $$
is finite. The obvious map
$$ \mu: S \to S, \quad (A,B,C) \to (-C,B,-A) $$
is an involution; if $S$ had odd cardinality, it would follow that $\mu$ has a fixed point, say $(A,B,-A)$, from which we would get $p = B^2 + 4A^2$. Unfortunately, $S$ has even cardinality since the involution
$$ \nu: S \to S, \quad (A,B,C) \to (A,-B,C) $$
has no fixed points: this is because $B = 0$ implies $p = 4AC$, which is impossible for prime numbers $p$.

We now would like to find a subset $U \subset S$ of $S$ with odd cardinality on which $\mu$ is still defined. The most natural idea would be considering the forms with $B > 0$. For showing that this set of forms has odd cardinality, we have to define an involution $(A,B,C) \to (A',B',C')$ on this subset that has exactly one fixed point. To find such an involution, we start with $(A,B,C) \to (A,-B,C)$ and then apply reduction by changing the middle coefficient modulo $2A$ and then adjusting the last coefficient so that the discriminant is $p$. This gives
$$ (A,-B,C) \to (A',B',C') = (A,-2A-B,A+B+C). $$
Now we are facing the problem that it is not clear at all that $B' = -2A-B > 0$, or that $C' = A+B+C > 0$. But if we set
$$ U = \{(A,B,C) \in S: A+B+C > 0 \}, $$
then the map
$$ \gamma: (A,B,C) \to (A,-2A-B,A+B+C) $$
actually is an involution on $U$. Moreover, $(A,B,C)$ is a fixed point if and only if $-2A-B = B$ and $A+B+C = C$, which is equivalent to $A = -B$. Since $p = B^2 - 4AC = B^2 + 4BC = B(B+4C)$ is prime, we must have $|A| = |B| = 1$. Since $A < 0$, this implies that the fixed point is $(-1,1,\frac{p-1}4)$; this form is equivalent to the principal form $(1,1,\frac{p-1}4)$.

The involution $\gamma$ on $U$ shows that $U$ has odd cardinality; the map
$$ (A,B,C) \to (-C,-B,-A) $$
is an involution on $S$ sending $U$ to $S \setminus U$, which impliesthat $|S| = 2 |U|$. The involution $\nu$ on $S$ sends elements with $B > 0$ to elements with $B < 0$, hence
$$ T = \{(A,B,C) \in S: B > 0\} $$
has the same number of elements as $U$, and in particular, it has odd cardinality. Finally, $\mu$ is an involution on $T$, and now the Two-Squares Theorem follows.

References

- H. Frick,
* Über den Zusammenhang der Perioden quadratischer Formen
positiver Determinante mit der Zerlegung einer Zahl in die
Summe zweier Quadrate*, Diss. ETH Zürich, 1918