Algorithm for finding the representation of an integer as a sum of two squares

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?

• Sure, I will edit it – proofromthebook Jun 16 '16 at 16:34
• Wikipedia has an article collecting existence proofs for the one-prime case. I haven't read it closely enough to see if any of them are constructive. It also mentions the identity $$(a^2+b^2)(p^2+q^2) = (ap+bq)^2 + (aq-bp)^2$$ which shows how to combine solutions for the one-prime case into a solution for the general case. – Henning Makholm Jun 16 '16 at 16:44
• There are definitely constructive proofs for primes; see the book of Niven/Zuckerman/Montgomery for one. – Greg Martin Jun 16 '16 at 17:10
• @proofromthebook, I have a solution I will be posting it later (maybe not here ). It's really simple. – user25406 Oct 17 '16 at 14:11
• Here it is, my solution. math.stackexchange.com/questions/1972771/… – user25406 Oct 17 '16 at 16:14

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$.
$$(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$.