PARI code analyzing whether given integer is presentable as sum of two integer squares

[Maximilian Hasler]

I found the PARI function I wrote some time ago.
It gives the list of all {a,b} such that a²+b²=n.


(In fact it's a collection of 2 functions, the second one being a
helper function treating the case of prime n where at most 1 solution
is possible.)
It should be clear how the function could be simplified and
accelerated if you only need to know if there is at least one
solution.
hope that helps
Maximilian

sum2sqr(n)={ local( t, f, p = 1, s=[] ); if( n<=1, return( if( n<0,
[], [[0, n]] )));

	if( isprime(n), return( if( n=sum2sqr_prime(n),[n],[])));

	for( i=1,matsize(f = factor(n))[1],

	  if( f[i,1] % 4 == 1, s = concat( s, [[sum2sqr_prime( f[i,1]
)*[1,I]~, f[i,2]]] ),

	    if( f[i,1] == 2, p = (1+I)^f[i,2] /* always comes first */,

	      if( f[i,2] % 2 == 1, return([]), p *= f[i,1]^(f[i,2]>>1))

	) ) );

	s = vector( #s, i, f=s[i][1]; vector( t=s[i][2]+1, j, f^(t-j)*conj(f)^(j-1) ));

	/* run over all elements of the cartesian product = n-uples with i-th
component in s[i] */

	t=[];

	forvec( T = vector(#s,i,[1,#s[i]]), f = prod( j=1,#T, s[j][T[j]], p);

		t = setunion( t, [ vecsort( abs([ real(f),imag(f) ]))] )

	);

	vecsort( eval( t ),1 )

}

sum2sqr_prime(p)={ local(c,d); if( p==2, [1,1], if( p%4 != 3,
 if( p%8==5, c=2, if( p%3==2, c=3, c=5; while( Mod(p,c)^(c\2)+1,
c=nextprime( c+1 ))));
 c=lift(Mod(c,p)^(p\4)); d=sqrtint(p); while(c > d, c=(p+0)%(p=c)); [c,p%c] ))}