Revisiting $(W^2 + X^2 + Y^2 + Z^2) = (A^2 + B^2)$ This query presents a narrow question that was broached by some of the responses
to the following query: Can a sum of $n$ squares be expressed as the sum of $n/2$ squares?
I'm unsure whether the question that I am specifically asking was intended by the
OP of the referenced query above.  Anyway...
Given any two fixed positive integers $A$ and $B$,
find all positive integer solutions to 
$(W^2 + X^2 + Y^2 + Z^2) = (A^2 + B^2)$
What I am looking for is analysis + an elegant algorithm to be applied in the
general case, as opposed to specific cases [e.g. where $(A^2 + B^2)$ is below
a certain number].  Shown below is the only (ugly) algorithm that I can come up with.
It seems to me that this algorithm requires computer support.
Can this algorithm be improved upon?  Especially pertinent, can sufficient
analysis be used to solve the general problem without computer support.
Alternatively, assume that there is no elegant algorithm.  In that case can : 
(1) necessary conditions be identified for a solution to exist 
(2) sufficient conditions be identified for a solution to exist 
(3) both (1) and (2)
$\underline{\textbf{my ugly algorithm}}$
Let $C = (A^2 + B^2)$. 
Let $S = \{(i,j) : i,j \in \mathbb{Z^+}, i \leq j, (i + j) = C\}.$ 
Go through each element of $S$ individually, identifying all solutions
$(W,X,Y,Z)$ for each element.  By a solution, I intend one of the three below:
(a)
$W^2 = i$ and $(X^2 + Y^2 + Z^2) = j.$
(b)
$W^2 = j$ and $(X^2 + Y^2 + Z^2) = i.$
(c)
$(W^2 + X^2) = i$ and $(Y^2 + Z^2) = j.$
When inspecting an individual element from $S$, follow the algorithm in the
"Inspect Individual Element" section below.  This query closes with an
"Analysis Around Pythagorean Triplets" section.  Perhaps the best that can be
offered is to add more special cases similar to the closing section of this query.
Inspect Individual Element
If neither $i$ nor $j$ happens to be a perfect square then solutions of types
(a) and (b) are immediately eliminated.  Assume that one of $i$ and $j$ is a
perfect square.  I know of no analysis that considers the question of whether
(for example) the Diophantine $(X^2 + Y^2 + Z^2) = j$ (with $j$ a fixed number)
is solvable, and if so,
how to identify all possible solutions.  Assuming that there is no such analysis,
then I see no way to attack the Diophantine $(X^2 + Y^2 + Z^2) = j$ other than
by a computer search.
The remainder of this section assumes that regardless of whether $i$ or $j$ is
a perfect square, solutions of type (c) are going to be sought.  I have seen
(elementary number theory) analysis around generating all possible pythagorean
triplets, but it is unclear how this analysis could be used against (for example)
the Diophantine
$(W^2 + X^2) = i$, 
where $i$ is a fixed positive integer that may or may not be a perfect square.
As in the search for solutions of types (a) and (b), if no analysis can be
brought to bear on the Diophantines involved in type (c), then (again) I think
that a computer search is inevitable.
Analysis Around Pythagorean Triplets
Suppose (for example) that when considering a specific element of $S$, that
$i$ is a perfect square (i.e. $i = k^2$).
Then analysis around
generating pythagorean triplets could be invoked.  That is, in this case,
any solution of type (c) that involves this particular element from $S$ must
involve a pythagorean triplet : 
$W^2 + X^2 = k^2.$
This means that if you could establish that no such pythagorean triplet is
possible tnat involves $k^2$, then you have established that no type (c) solution
is possible for this particular element of S.  Further, if you could enumerate
all possible solutions to 
$W^2 + X^2 = k^2$, 
then any type (c) solution for this particular element from $S$ must involve
one of these pythagorean triplet solutions.
However, even here, assuming that $j$ is not a perfect square you would still
be left with solving $Y^2 + Z^2 = j.$
Addendum
While exploring the "Lagrange's four-square theorem" link in John Omielan's answer,
I (think that I have) discovered a way to moderately refine the algorithm, based on the
Sum of two squares theorem.  Alternatively, perhaps I'm misreading the theorem.
Suppose, for a specific element $(i,j) \in S,$ neither $i$ nor $j$ is a perfect square.
Then the only solutions possible that utilize this element will be of type (c).  If
either $i$ or $j$ fails the criteria of the "Sum of two squares" theorem, then this
specific element of S can not be used to generate a solution.
Alternatively, if $i$ and $j$ both pass the criteria of the "Sum of two squares" theorem,
then at least one type (c) solution (and perhaps more than one) will be possible with
this particular element of $S$.
Addendum-2
The "Sum of two squares theorem" linked in the addendum also (indirectly) affects
the discussion in the "Analysis Around Pythagorean Triplets" section.  Suppose
that for a specific element $(i,j) \in S,$ that $i$ is a perfect square (i.e. $i = k^2$).
Further suppose that (for the moment) you are only exploring (possible) type (c) solutions
for this element [as opposed to type (a) solutions].
Further suppose that you have used "Pythagorean triplet analysis" to determine that
there is at least one solution to $W^2 + X^2 = k^2.$
As indicated in the "Analysis Around Pythagorean Triplets" section, this still
leaves the $Y^2 + Z^2 = j$ Diophantine to be solved.
If $j$ also happens to be a perfect square (i.e. $j = l^2$), then "Pythagorean triplet analysis"
can (also) be used to attack the $Y^2 + Z^2 = l^2$ Diophantine.
Alternatively, if $j$ is not a perfect square, then the "Sum of two squares theorem"
can be used to attack the $Y^2 + Z^2 = j$ Diophantine.
Addendum-3
Although the algebra in my Addendum and Addendum-2 sections seem valid, in light of the
comment posted by Steven Stadnicki following John Omielan's answer, I question the
utility of these two sections.  Given fixed integers $A$ and $B$ that are both
large, nothing in the Addendum or the Addendum-2 sections
permits the specific problem that involves the specific values for $A$ and $B$
to be comprehensively attacked without computer support.
Consequently, to evaluate the utility of the Addendum and/or Addendum-2 sections,
one would first have to examine the computer efficiency of the algorithm(s)
referred to by Steven Stadnicki (which I do not understand).
Then, one would have to consider the cost of identifying each element $(i,j) \in S,$
and (perhaps) applying the "Sum of Two Squares" theorem to this element.  This would
require that for each element $(i,j)$ you would (usually) have to compute the prime factorizations
of both $i$ and $j$.  Among other considerations, this means that you need a table of
prime numbers for $\{1, 2, 3, 4, \cdots, [A^2 + B^2]\}.$
To some extent, the cost of
construction of the prime number table might be "amortized", if you were solving
more than one problem at a time.  Anyway, my (blind) instinct is that my Addendum and
Addenum-2 analysis might not have accomplished anything.
Note that although I am questioning the utility of the analysis that I added in my Addendum and Addendum-2 sections, I am not similarly questioning the utility of the analysis provided by John Omielan's answer.  This is because his analysis (e.g. determining the residue mod 16 of a number) seems computer non-intensive to me.
Finally, I must confess that this (Addendum-3) section is motivated by meta-cheating.
Steven Stadnicki's comment suggests that "rabbits" are unlikely.  Am I to believe
that someone with an undergraduate level comprehension of number theory (me), who
just happened upon the "Sum of Two Squares" theorem, has suddenly re-invented the wheel?
Addendum-4 Response to Case A in Arnold's answer.
Case A in this answer seems to beg some questions:
(1)
Given fixed positive integers $(p,q)$ with $p^2 + q^2 = r,$ 
can any non-zero integer values for $(a,b,c,d)$ be found 
so that $a^2 + b^2 + c^2 + d^2 = r$
where simultaneously the pattern of case A does not apply.
Situations where the pattern of case A does not apply would be that
either:
(1a)
The elements in $\{a,b,c,d\}$ can not be paired up so that
(for example) 
$(|a|+|b|) = p$ and $(|c|+|d|) = q.$
(1b)
The elements in $\{a,b,c,d\}$ can not be paired up so that
(for example) 
$(ab + cd) = 0.$
(2)
Ignoring the difficulties around (1) directly above, can the analysis
in Case A be applied to large values of $(p,q)$ manually
(i.e. where no computer support is used)?
This means that (for example) although you could (arguably) compute the residue
of $p$ or $q$ modulo some positive integer, you could not compute the
prime factorizations of $p$ or $q.~~$ Arguably, you (also) would not
have access to a list of prime numbers.
This means, that with $(p,q)$ fixed, you have to analytically identify
all satisfying sets $[a,b,c,d]$ of non-zero integers such that
$(|a|+|b|) = p, ~ (|c|+|d|) = q, ~ \text{and} ~ (ab + cd) = 0.$
(3)
Ignoring the difficulties in (1) above, and assuming that
the approach taken in Case A is used
[as described in (2) directly above], and further assuming that this approach
does require computer support,
then the (perhaps unanswered) question is:
How does the computer effiency of the case A approach compare with
the computer efficiency of the algorithm(s) referenced by
Steven Stadnicki in his comment to John Omielan's answer?
 A: I don't know of any elegant algorithm to find all of the solutions in the general case, or any particular way to improve your algorithm. Nonetheless, note Lagrange's four-square theorem states

Lagrange's four-square theorem, also known as Bachet's conjecture, states that every natural number can be represented as the sum of four integer squares.

However, this also includes cases where one or more of those squares are $0$. Since you state all $4$ squares are to be positive, note the Uniqueness section states the only positive integers which cannot be expressed as the sum of four non-zero squares are

... the eight odd numbers $1, 3, 5, 9, 11, 17, 29, 41$ and all numbers of the form $2(4^{k}),6(4^{k})$ or $14(4^{k})$.

Thus, as long as $A^2 + B^2$ is not one of these numbers (e.g., $A = 1$ and $B = 2$ has no solution since $1 + 4 = 5$, and also $A = B = 2^{n}$ for any $n \ge 0$ have no solution either), there is at least one solution. This gives both the necessary and sufficient conditions you asked about for a solution to exist.
Update: In case you weren't already aware of these, there are several fairly simple & efficient ways you can cut down on the number of checks your algorithm, as stated in the question, needs to make. For example, note that all perfect squares are congruent to an element of $\{0, 1, 4, 9\}$ modulo $16$. You can use this to check on $i$ in case (a) and $j$ in your case $b$.
Also, the sum of $2$ squares can only be congruent to an element of $\{0, 1, 2, 4, 5, 8, 9, 10, 13\}$ modulo $16$, which you can use to check $i$ and $j$ in your case (c).
Finally, as shown in Dickson's "Modern Elementary Theory of Numbers" book's Table 5, the expression $ax^2 + by^2 + cz^2$, where $a = b = c = 1$, cannot be equal to any positive integers of the form $A = 4^{k}(8n + 7)$. You can use this to check $j$ in your case (a) and $i$ in your case (b).
A: We can always find $4$ squares that equal $2$ squares if all terms are parts of Pythagorean triples i.e. where $A_1^2+B_1^2+A_2^2+B_2^2=C_1^2+C_2^2.$
Let us start with Euclid's formula $\quad A=m^2-k^2\quad B=2mk\quad C=m^2+k^2\quad$ and find triples, if they exist, for any two $C$-values. Any value of C that returns an integer within the limits for $m$ provides the $m,k$ values needed to generate a Pythagorean triple.
$$C=m^2+k^2\implies k=\sqrt{C-m^2}\qquad\text{for}\qquad \bigg\lfloor\frac{ 1+\sqrt{2C-1}}{2}\bigg\rfloor \le m \le \lfloor\sqrt{C-1}\rfloor$$
The lower limit ensures $m>k$ and the upper limit ensures $k\in\mathbb{N}$.
For example:
$$C=65\implies \bigg\lfloor\frac{ 1+\sqrt{130-1}}{2}\bigg\rfloor=6 \le m \le \lfloor\sqrt{65-1}\rfloor=8\quad\land \quad m\in\{7,8\}\Rightarrow k\in\{4,1\}\\$$
$$F(7,4)=(33,56,65)\qquad \qquad F(8,1)=(63,16,65) $$
This example provides two pairs of $A^2,B^2$ that add up to $65^2$ but we need choose only one of the pairs unless we want the sum to be  $2C^2$. Now, using the same technique for another known $C$-value:
$$C=85\implies \bigg\lfloor\frac{ 1+\sqrt{170-1}}{2}\bigg\rfloor=7 \le m \le \lfloor\sqrt{85-1}\rfloor=9\quad\land \quad m\in\{7,9\}\Rightarrow k\in\{6,2\}\\$$
$$F(7,6)=(13,84,85)\qquad \qquad F(9,2)=(77,36,85)$$
Now we have two sets of $A^2,B^2$ (choose one) that add up to $65^2$, likewise for $85^2$, and there are four ways to combine these. Some $C$-values will have one, or many pairs for $A,B$ but if no value of $m$ yields an integer $k$, then no triple exists for that $C$-value.
A: We have the identity's:
Case A
$(a,b,c,d)^2=(p,q)^2$
Where, $(p,q)=[(a+b),(c+d)]$
And Condition is: $(ab+cd=0)$
For, $(a,b,c,d)=(9,7,-21,3)$ we get:
$(p,q)=(16,18)$
Case B:
$(a,b,c,d,e,f)^2=(p,q,r)^2$
Where, $(p,q)=[(a+b),(c+d),(e+f)]$
And Condition is: $(ab+cd+ef=0)$
For, $(a,b,c,d)=(3,2,4,1,-10,1)$ we get:
$(p,q,r)=(5,5,9)$
Case C:
$(a,b,c,d,e,f,g,h)^2=(p,q,r,s)^2$
Where, $(p,q)=[(a+b),(c+d),(e+f),(g+h)]$
And Condition is: $(ab+cd+ef+gh=0)$
For, $(a,b,c,d)=(15,13,12,9,8,4,-67,5)$ we get:
$(p,q,r,s)=(28,21,12,62)$
Similarly one can apply the method for $(2n)$ squares
on LHS & get $(n)$ squares on the RHS for any 'n' &
suitable values of $(a,b,c,d,-----)$.
A: I think what you might be missing is the gap in efficiency. The size of your set $S$ is linear in $N$, where $N$ is the number simultaneously being expressed as a sum of two and a sum of four squares. This means that any algorithm that relies on iterating through the members of $S$ will take time at best linear in $N$ and may take longer depending on how many operations have to be done for each element.
This might not seem bad on the surface, but a shift in perspective might help understand why it's 'slow'. Since a number can be written 'efficiently' using the standard base-2 or base-10 notation, for number-theoretic algorithms what's usually thought of as 'efficient' is based on the logarithm of the input numbers, or equivalently the size ('in memory') of them. You couldn't count to 31415926 without taking a very long time, but you know how to add it to 27182818 in just a handful of operations. What's more, using just a few more steps you could multiply the two numbers together. In particular, the addition of $n$-bit numbers — that is, numbers of size $\leq N=2^n$ — can be done in $O(n)$ time. This is $O(\log N)$ when represented in terms of the actual numbers being manipulated. Similarly, naive multiplication can be done in $O(n^2)$ time, and it can be shown that division can be done in comparable time. One of the most foundationally important results in Computer Science in the last few decades was the result that primality testing can be done in time polynomial in the length of the number tested — that is to say, in $O(n^k)$ for some exponent $k$.
Thought of in these terms, the algorithm for solving the four-squares problem that was mentioned in my comment takes time $O(n^2)$ (maybe up to some much smaller factors, I suspect); your algorithm, by contrast, takes time roughly $2^n$; it's exponentially worse than the probabilistic algorithm.
I hope this helps make sense of the notion!
