Representing a number as a sum of four squares and factorization

Question

Rabin and Shallit 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). 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 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. 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 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...

Answer 1

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!)