Diophantine equation $a^2+b^2=c^2+d^2$

Essentially you ask for a parametrization of the quadric $V(a^2+b^2=c^2+d^2)$. There is a general geometric method how to do that, which you can find here. On page 13, this example is discussed: Solutions of $a^2+b^2=c^2+d^2$ are parametrized by

$(a,b,c,d) = (p r + q s , q r - p s , p r - q s , p s + q r),$

where $p,q,r,s$ are arbitrary. But one can also derive this via complex numbers (similar to the complex numbers solution of Pythagorean triples): Let $u = p + qi$, $v = r + si \in \mathbb{C}$. Then we have

$$|u \overline{v}| = |u| |\overline{v}| = |u| |v| = |u v|.$$

Square both sides, and compute the norms of the products explicitly. This gives you

$$(p r + q s)^2 + (q r - p s)^2 = (p r - q s)^2 + (p s + q r)^2.$$

Let $N$ be a positive integer with prime power factorization $$N=2^a \prod_{i=1}^m p_i^{b_i}\prod_{j=1}^n q_j^{c_j},$$ where the $p_i$ are distinct primes congruent to $1$ modulo $4$, and the $q_j$ are distinct primes congruent to $-1$ modulo $4$. (We allow $a$ to be $0$, and $m$ or $n$ or both could be $0$.)

If at least one $c_j$ is odd, then $N$ has no representations as a sum of two squares. If all the $c_j$ are even, then the number of representations of $N$ as a sum of two squares is equal to $f(N)$, where $$f(N)=4\prod_{i=1}^m (b_i+1).$$ In this formula, in counting representations, order matters, and we allow the possibility of using negative integers.

If we want, for example, the representation $5=1^2+2^2$ to count as essentially the same as the representation $5=2^2+1^2$, and we do not want to allow negative integers, things get somewhat more complicated. Let $g(N)$ be the number of representations of $N$ as a sum of two non-negative squares, where order does not matter.

If $f(N)$ as defined above is divisible by $8$, then $g(N)=f(N)/8$.

If $f(N)$ as defined above is not divisible by $8$, then $g(N)=(f(N)+4)/8$.

The above formula gives a partial answer to your question, in that it tells us the number of solutions of the equation $a^2+b^2=c^2+d^2=N$, where order does not matter,and we allow only non-negative integers in the representation.

Your question also asks how to actually produce the solutions. The difficult part is to find representations of the $p_i$ as a sum of two squares. As you mentioned, there is, for any prime $p_i$ congruent to $1$ modulo $4$, essentially only one such representation. There are reasonably good algorithms available for producing the representation of such primes.

Once we have representations for the appropriate primes, we can obtain representations for products by repeatedly using the Brahmagupta Identity $$(s^2+t^2)(u^2+v^2)=(su+tv)^2+ (sv-tu)^2.$$ Here is a simple example. We have $13=2^2+3^2$ and $17=1^2+4^2$. So taking $s=2$, $t=3$, $u=1$, and $v=4$ we get $221=13\cdot 17= 14^2 + 5^2$.

Another essentially equivalent approach is to factor $N$ as a product of Gaussian primes. Once we have that, we can easily produce all representations of $N$ as a sum of two squares.