We ran into a similar problem with the complex numbers. If p is 1 mod 4, it is the sum of two squares. Using the gcd algorithm in the complex plane, we developed an efficient algorithm for finding a2+b2, even if p is very large. Call this the pab algorithm; we'll need it below.
If n is even, remove the factors of 2. We can always multiply the answer by 1+i as necessary. In fact you may want to do this for other small factors of n.
By a theorem in analytic number theory (which is not yet on-line), every arithmetic progression contains infinitely many primes. Let q be a prime equal to 1 mod 4 and -1 mod n. Use the pab algorithm to find a2+b2 = q. Then let w = 1+aj+bk, and note that |w| is divisible by n.
If |w| = n we are done, but if not, the gcd algorithm comes to our rescue, just as it did in the complex plane.
Let g be the right gcd of w and n. Note that g is not an integer, else g would have to divide the coefficients of w. This means g is not n, or an associate of n. Rather, g is a quaternion that divides n, with |g| < n2.
At the same time, |g| divides |w|. Most of the time |w|/n will be coprime to n. If not, select another q and build another w. With |w|/n coprime to n, |g| divides n.
Suppose |g| < n. In other words, g spans the integer k, where k is a proper divisor of n. That means w and n also span k. If x is a linear combination of w and n, what is the norm of x mod n? We examined this question on the previous page. The generator n contributes nothing, so |x| is a multiple of |w|, which is a multiple of n. The norm of everything spanned by w and n must be divisible by n. The norm of any left multiple of g is divisible by n. So g cannot generate k. Therefore |g| = n, and we have split n into the sum of four squares.
There are many ways to tweak this algorithm. We could have set w = 1+i+aj+bk, where q is a prime that is 1 mod 4 and -2 mod n. Or, w could be 2+aj+bk, where q is 1 mod 4 and -4 mod n. In fact, q could be -p mod n, where p is any prime equal to 1 mod 4. Find s2+t2 = p, and set w = s+ti+aj+bk.
If n is hundreds of digits, or even tens of digits, you'll probably find a suitable prime q by considering n-k2-l2, for various values of k and l. In this case there's no need for the quaternion gcd algorithm, because the norm is already n.