If $xy=z^2+1$ for natural $x$, $y$, $z$, prove that $x=a^2+b^2$, $y=c^2+d^2$, $z=ac+bd$ for some integers $a$, $b$, $c$, $d$ A contest-math number theory problem

If $x,y,z \in \mathbb{N}$ and $xy=z^2+1$, prove that there exist integers $a$, $b$, $c$, $d$ such that $x=a^2+b^2$, $y=c^2+d^2$, $z=ac+bd$.


After working for 20 hours, I learned that there is a solution using complex numbers.  But I have no idea about the solution.
I don't know what will come out of this.
$$xy=(z-i)(z+i)$$
I guess if I could prove that $x = (a-bi)(a+bi)$ and that $y=(c-di)(c+di)$ then maybe I would get a result.
But I don't know how to relate this to the fact that $a, b, c, d$ are integers.

Question 1: (Main)
Can we find a solution using pure number theory without using a complex numbers?
Question 2:
How can a solution be made using complex numbers?

I would love to get a detailed answer.  Thank you!
 A: The method below does not require knowing all the primes that can be expressed as the sum of two squares.
We have determinant $1$ in
$$
A =
\left(
\begin{array}{cc}
x & z \\
z & y
\end{array}
\right)
$$
As the entries and trace are positive, it is positive definite.
Gauss reduction is the finite process of multiplying successively in the form $A \mapsto P^T A P,$ with the choice of
$$
P =
\left(
\begin{array}{cc}
1 & n \\
0 & 1
\end{array}
\right)
$$
or
$$
P =
\left(
\begin{array}{cc}
0 & -1 \\
1 & 0
\end{array}
\right)
$$
until the final output is the identity matrix. This comes from the inequalities of Gauss reduction. So $I = P_j^T A P_j $ for an integer matrix $P_j$ of determinant $1.$   By taking $Q = P_j^{-1}$ we get
$$ A = Q^T Q  $$
$$
Q =
\left(
\begin{array}{cc}
a & c \\
b & d
\end{array}
\right)
$$
NOTE: Gauss reduced means $\langle r,2s,t \rangle$
with $1 \leq r \leq t,$ then $-r < 2s \leq r,$ in the matrix below. A little fiddling with inequalities is needed.
$$
\left(
\begin{array}{cc}
r & s \\
s & t
\end{array}
\right)
$$
A: Hint. Notice that every odd prime factor $p$ dividing $z^2+1$ satisfies $p\equiv 1\bmod 4$, since it is the only way that $-1$ is a quadratic residue (we will treat the case $2\mid z^2+1$ at the end). Besides, Fermat's theorem establishes that for every $p\equiv 1\bmod 4$ there exist integers $r,s$ such that $p=r^2+s^2$ (this can be proven, e.g., with Thue's Lemma). Thus, both $x,y$ are a product of odd numbers (and possibly powers of these odd numbers).
Yet, we also have from Brahmagupta's Identity that $$(r^2 + s^2)(t^2 + u^2) = (rt + su)^2 + (ru - st)^2 $$ You can, hence, take two odd primes - say $p_1$ and $p_2$ - , and write them as a sum of squares: $p_1\cdot p_2=(r^2+s^2)(t^2+u^2)=(rt+su)^2+(ru-st)^2$. Use this identity repeatedly. Can you finish?
For the case $2\mid z^2+1$, notice that $2=1^2+1^2$.
Now, we just need to prove that these $a,b,c,d$ satisfy $ad - bc=\pm1$. This looks harder using elementary techniques; I might try it later...
