[Rabin and Shallit](https://doi.org/10.1002%2Fcpa.3160390713) have a randomized polynomial-time algorithm to express an integer $n$ as a sum of four squares $n=a^2+b^2+c^2+d^2$ (in time $\log(n)^2$ assuming the Extended Riemann Hypothesis). 


I'm wondering why this does not give an efficient factorization algorithm? Here's what one could try: run their algorithm $m$ times, with different random steps. This should give expressions
$n=a\_l^2+b\_l^2+c\_l^2+d\_l^2, l\leq m$, presumably with many distinct representations as a sum of four squares (cf. [Jacobi's theorem](https://en.wikipedia.org/wiki/Jacobi%27s_four-square_theorem)). We can think of these as factorizations of $n$ over the Lipschitz integers,
so $n=|a+bi+cj+dk|^2=(a+bi+cj+dk)(a-bi-cj-dk)$. 


The Lipschitz integers do not admit a Euclidean algorithm, but the [Hurwitz quaternions](https://en.wikipedia.org/wiki/Hurwitz_quaternion) do. Hence one should be able to take the $\gcd$ of Hurwitz quaternions efficiently. I.e., for $N,D$ Hurwitz quaternions, there should be an efficient algorithm to find $N=QR, D=PR$, with $|R|< |N|,|D|$. 


Now, take $\gcd(a\_l+b\_li+c\_lj+d\_lk,a\_p+b\_pi+c\_pj+d\_pk)$, $1\leq l <p\leq m$, where the $\gcd$ is taken in the [Hurwitz quaternions](https://en.wikipedia.org/wiki/Hurwitz_quaternion). It should be efficient to find the $\gcd$ since the Hurwitz quaternions admit a Euclidean algorithm. In turn, this should give further factorizations of $a\_l+b\_li+c\_lj+d\_lk$ into Hurwitz quaternions, and hence the norms of these factors will give factors of $n$.


This approach will fail if it turns out that all of these Hurwitz $\gcd$ factors differ by Hurwitz units, for example if $n$ is prime. Of course, we could initially run a [polynomial-time](https://en.wikipedia.org/wiki/AKS_primality_test) primality test to make sure $n$ is not prime. 


**Question:** Are there certain composite $n$ for which one will not obtain a factorization this way with high enough probability to give a fast algorithm (i.e. $m$ has to be too large to get a pair with non-trivial $\gcd$ with high probability)? Maybe I just need to think a bit more about the proof of Jacobi's theorem...



Representing a number as a sum of four squares and factorization
I think the reason is that there are $p+1$ distinct ways of writing an odd prime $p$ as the sum of four squares up to sign changes; these correspond to the same number of elements of the Lipschitz order up to units. If you take two different Lipschitz elements of reduced norm $p$ up to units, their greatest common divisor is $1$. 


So if we take two random Lipschitz elements of reduced norm $n=pq$, then their greatest common divisor will be $1$ with probability $(1-1/(p+1))(1-1/(q+1))$, and I don't see how you win with this. (These aren't significantly different odds than trying a random element modulo $n$ and hoping for a factor in common!)


