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


Here is a pari/gp program and example:

    /*
    ? n=nextprime(10^220+30000)*nextprime(2^1000+40000000);n%8
    %136 = 3
    ? t=threesquares(n);
    
    ? ##
      ***   last result computed in 2,210 ms.
    ? round(log(n))
    %138 = 1200
    */

    pl=10^6;\\ bound for trial division, may need change
    default(primelimit,pl);
    {
    twosquares(n)=
    local(K,i,v,p,c1,c2);
    K=bnfinit(x^2+1);
    v=bnfisintnorm(K,n);
    for(i=1,#v,p=v[i];c1=polcoeff(p,0);c2=polcoeff(p,1);if(denominator(c1)==1&&denominator(c2)==1,return([c1,c2])) );
    return([]);
    }
    
    {
    threesquares(n)=
    local(m,z,i,x1,y1,j,fa,g);
    if(n/4^valuation(n,4)==7,return([]););
    for(z=1,n,
    \\forstep(z=sqrtint(n),1,-1,
    m=n-z^2;
    if(m%4==3,next);
    print1(z," ",);
    fa=factor(m,pl);
    g=1;
    for(i=1,#fa~,if(!ispseudoprime(fa[i,1])/*||!isprime(fa[i,1])*/||(fa[i,2]%2==1&&fa[i,1]%4==3),g=0;break; ));
    if(!g,next);
    print("\nfound ",z," "," m=",m,factor(m));
    j=twosquares(m);
    print("j=",j);
    x1=abs(j[1]);
    y1=abs(j[2]);
    return([x1,y1,z]);
    );
    }