I'm not an expert, but from memory something which it turns out is called the Brahmagupta–Fibonacci identity or occasionally Diophantus identity (used constructively in the proof of the result of Fermat) might be useful: see Wikipedia.
$$\left(a^2 + b^2\right)\left(c^2 + d^2\right) = \left(ac-bd\right)^2 + \left(ad+bc\right)^2
= \left(ac+bd\right)^2 + \left(ad-bc\right)^2$$
This allows you to strip off factors one at a time.
Where you have squares in the factorization, as in the $5\cdot 11^2$ case, note that the square can be factored out as $(p^2+q^2)m^2 = (pm)^2 + (qm)^2$
Similar to the Brahmagupta-Fibonacci two-square identity. Euler has a four square identity which involves the sum of 4 squares:
$$(a_1^2+a_2^2+a_3^2+a_4^2)(b_1^2+b_2^2+b_3^2+b_4^2) =\\
\quad(a_1b_1 - a_2b_2 - a_3b_3 - a_4b_4)^2 + (a_1b_2+a_2b_1+a_3b_4-a_4b_3)^2
+(a_1b_3 - a_2b_4 + a_3b_1 + a_4b_2)^2 + (a_1b_4 + a_2b_3 - a_3b_2 + a_4b_1)^2$$
Factor $1638$ as products of any small factors you know how to represent as sum of 4 squares. Repeat apply the formula will allow you to represent $1638$ itself as sum of 4 squares.
For example, let's say we have factored $1638$ as $2\cdot 3^2 \cdot 7 \cdot 13$, we have:
$$\begin{align}
& 2\cdot 3^2 \cdot 7 \cdot 13\\
= & (1^2+1^2+0^2+0^2)(1^2+1^2+1^2+0^2)^2(2^2+1^2+1^2+1^2)(3^2+2^2+0^2+0^2)\\
= & (0^2 + 2^2 + 1^2 + 1^2)(1^2+1^2+1^2+0^2)(2^2+1^2+1^2+1^2)(3^2+2^2+0^2+0^2)\\
= & ((-3)^2 + 1^2 + 2^2 + 2^2)(2^2+1^2+1^2+1^2)(3^2+2^2+0^2+0^2)\\
= & ((-11)^2+(-1)^2+2^2 + 0^2)(3^2+2^2+0^2+0^2)\\
= & ((-31)^2 + (-25)^2 + 6^2 + 4^2)\\
\end{align}$$
This give you a non-trivial representation of $1638$ as $31^2 + 25^2 + 6^2 + 4^2$.
In general, there are many representations of a number as a sum of 4 squares.
There is a theorem:
The total number of representations of a positive integer $n$ as the
sum of four squares, representations that differ only in order and sign
being counted as distinct, is eight times the sum of the divisors of
$n$ that are not multiple of $4$.
The above representation is only $1$ out of $8 \sum_{d\mid 1638, 4 \nmid d} d = 34944$
ways of representing $1638$ as sum of 4 squares.
Best Answer
you can get a solution with complexity $\mathcal O (1)$
there are $5$ types of solutions:
let $A,B,C,D,E$ be the number of solutions of each type.
Clearly $A=1$ if $4|N$ and $0$ otherwise.
We have $B=\lfloor (N-1)/3 \rfloor-A$
we have $C=0$ if $C$ is odd and $C=(N/2-1-A)/2$ otherwise.
We have $D=[(\sum\limits_{i=1}^{2i<N}N-2i-1)-C]/2-B$.
Finally, using stars and bars we get $E=\frac{\binom{N-1}{3}-A-4B-6C-12D}{4!}$.
And clearly you want $A+B+C+D+E$