Given your constraints, I think that you can do a smart brute force.
First, note that if p = a^2 + b^2 + c^2 + d^2, each of a, b, c, d have to be less than 10^6. So just loop over a from 0 to sqrt(p). Consider q = p - a^2. It is easy to check whether q can be written as the sum of three squares using Legendre's three-square theorem. Once you find a value of q that works, a is fixed and you can just worry about q.
Deal with q the same way. Loop over b from 0 to sqrt(q), and consider r = q - b^2. Fermat's two-square theorem tells you how to check whether r can be written as the sum of two squares. Though this check requires O(sqrt(r)) time again, in practice you should be able to quickly find a value of b that works.
After this, it should be straightforward to find a (c,d) pair that works for r.
Since the loops for finding a and b and (c,d) are not nested but come one after the other, the complexity should be low enough to work in your problem.
n
as the sum of four squares (I think it's not free though): onlinelibrary.wiley.com/doi/10.1002/cpa.3160390713/… And here an explanation of the algorithm for three squares: math.stackexchange.com/questions/483101/… – Keiwan