We know that an integer $n$ is the sum of two squares if and only if all its prime divisors $p$ of the form $p \equiv 3 \pmod4$ have an even exponent in the prime factor decomposition of $n$.
My question is, if I know beforehand that $n$ is representable as a sum of two squares, how can I find its representation as a sum of two squares without having to find its prime divisors. If that is not possible (or impractical), then how can I, given a prime $p$ such that $p \equiv 1 \pmod4$ decompose it as a sum of two squares?
Best Answer
Rabin has given a randomized algorithm for expressing prime $p = 4k+1$ as a sum of two squares in expected time $O(\log p)$. This previously published method is described in the first section of Rabin and Shallit's 1986 paper, "Randomized Algorithms in Number Theory" (Comm. in Pure and Applied Math v.39(supplement):S239-S256). See also the previous Question Efficiently finding two squares which sum to a prime.
For a general non-prime $n$ which you "know beforehand ... is representable as a sum of two squares", I'm not aware of an equally fast method. Perhaps this knowledge comes with some sort of "certificate" that would be useful in finding the decomposition. A somewhat naive approach would be perform a saddleback search for $a \ge b \ge 0$ such that $a^2 + b^2 = n$ starting at a guess $a = \lfloor \sqrt n \rfloor$.
The well-known identity:
$$ (a^2+b^2)(c^2+d^2)=(ac+bd)^2+(ad-bc)^2 $$
allows us (among other things) to focus on the case $n$ odd by eliminating factors of $2=1^2 + 1^2$.