Q: Rabin and Shallit Algorithm I want to implement Rabin and Shallit algorithm for representing a positive integer as a sum of three squares. Can anyone give me a rough sketch of the method? I searched through the internet but couldn't find any clues.
Thanks in advance.

A: Rabin and Shallit's 1986 paper, "Randomized Algorithms in Number Theory" (Comm. in Pure and Applied Math v.39(supplement):S239-S256), uses "a number of Diophantine problems involving sums of squares" to illustrate efficiencies that are possible with random choices, reducing the expected number of operations while still always producing correct answers (assuming some reasonable conjectures in number theory).
Section 1 revisits a previously published randomized algorithm of Rabin for representing a prime $p=4k+1$ as the sum of two squares, requiring an expected number $O(\log p)$ of operations.  This conclusion doesn't rest on any unproven conjectures.
Sections 2 and 3 deal with three randomized algorithms to represent a positive integer $n$ as a sum of four squares.  The first one has lowest expected operation/time complexity $O(log^2 n)$, but proving this complexity depends on an ERH assumption.  The other two also exhibit random polynomial time, but without dependence on unproven conjectures.
Section 4 treats expressing $n \gt 0$ as a sum of three integer squares, which is possible if and only if $n$ is not of the form $4^a (8k+7)$ for $a,k \ge 0$.  By a usual argument, we may remove any factors of $4$ from $n$.
They briefly argue that any prime $p \not\equiv 7 \pmod{8}$ may be represented by a sum of three (or fewer) squares in random polynomial time, without assuming any unproven conjecture.  For $p=2$ and $p \equiv 1 \pmod{4}$, representation as the sum of two squares is possible "in a random polynomial time using the algorithms of Section 1."
That leaves, among the prime cases, where $p \equiv 3 \pmod{8}$.  Appealing to $\mathbb{Z}[\sqrt{-2}]$ being a (Euclidean) domain where GCD's can be computed "in polynomial time" and to $-2$ being a quadratic residue for primes of the form $8k+3$, they say "the techniques of Section 2" suffice to get $p = x^2 + 2y^2$, and thus a sum of three squares, in random polynomial time.  "The details are left to the reader."  However I suspect they actually mean here to apply the techniques of Section 1.
In cases where $n$ is not prime (or itself a square), they propose to try writing $n$ as the sum of a prime $p$ and a square if $n \equiv 1 \text{ or } 2 \pmod{4}$, or as a square plus twice a prime $2p$ if $n \equiv 3 \pmod{8}$.  In the former cases expressing $p$ as a sum of two squares (resp. in the latter cases expressing $2p$ as a sum of two squares) then gives the desired representation of $n$ as a sum of three squares.
A detail whose discussion they defer to the paper's last page (excepting the bibliography) is helpful at this point.  Suppose $n \equiv 3 \pmod{8}$ and $n = x^2 + 2p$, as in the latter cases of the previous paragraph.  Since squares modulo $8$ are congruent to $0$ or $1$, $2p \equiv 2 \pmod{8}$ is compelled, and $p \equiv 1 \text{ or } 5 \pmod{8}$.  Either way $p \equiv 1 \pmod{4}$ and is thus a sum of two squares, say $p = y^2 + z^2$.  Then also $2p = (y+z)^2 + (y-z)^2$ is a sum of two squares.
Similarly, in cases where $n \equiv 1 \text{ or } 2 \pmod{4}$ and $n = x^2 + p$, taking both sides modulo $4$ proves that $p \equiv 0,1,\text{ or } 2 \pmod{4}$, so unless $p=2$, we have $p \equiv 1 \pmod{4}$.  "In either case, we can express $p$ as the sum of two squares."
The adequacy of the proposed algorithm then rests on the chance of success in getting $n = x^2 + p$ if $n \equiv 1 \text{ or } 2 \pmod{4}$ (resp. $n = x^2 + 2p$ if $n \equiv 3 \pmod{8}$) with prime $p$ for a randomly chosen integer $0 \le x \le \sqrt{n}$.
For the first of these the authors restate Conjecture H of Hardy and Littlewood, specifically:
Every sufficiently large number n is either a square or the sum of a prime and a square.
together with a supposition about the number $N(n)$ of such representations.  They propose for the second set of cases a similar, apparently novel conjecture about expressions "as the sum of a square and twice a prime."
The details of these conjectures do not affect the implementation of the algorithm, at least as given in their pseudo-code presentation, but I'm puzzled by their intended handling of twenty-two known small exceptions to Conjecture H when $n \equiv 1 \text{ or } 2 \pmod{4}$:
$$ 5,10,13,34,37,58,61,85,130,214,226,370,526,706,730,829,1414,1549,1906,2986,7549,9634 $$
which they say are "probably all the exceptions."  While most of these are already sums of two squares, I found three which are not so: $214,526,1414$.
