Computing $x$ and $y$

Suppose $p$ is a prime of the form $4m+1$. There is a construction of Legendre of $x$ and $y$ that is explained on pages 120-121 of Davenport. I'm unconvinced that it is any more efficient than the following naive algorithm: compute $\sqrt{p-x^2}$ for $x=1,2,\ldots$ until it's an integer. This takes at most $\sqrt{p}$ steps. Here's a simple PARI program which implements this algorithm.

{sumoftwosquares(n) =
   local(y);
   for(x=1,floor(sqrt(n)),
      y=sqrt(n-x^2); 
      if(y-floor(y)==0, return([x,floor(y)]))
   );
   error(n," is not a sum of two squares.")
}