A modification of Dror's comment.

This probabilistic algorithm worked for me.

The main idea is to pick some $z$, compute $m=n - z^2$, factor $m$ with trial division and express it as a sum of two squares if possible. The probability of finding prime $m=4a+1$ or $2p$ is high enough for practical purposes.

The algorithm:

1. z:=0
2. z:=z+1
3. m:=n-z^2
4. if can't trial factor m goto 2
5. if m=x^2+y^2 (the factorization is known) then x^2+y^2+z^2=n. Done
6. goto 2