Expected runtime analysis for sums of four squares (Rabin and Shallit)

I've been reading a paper of Rabin and Shallit ("Randomized Algorithms in Number Theory"), which gives a brief sketch of an ERH-conditional algorithm to compute a representation of a positive integer $$n$$ as a sum of four squares, originally presented in some unpublished notes. The claimed expected time complexity of the algorithm is $$O(\log^2 n)$$. I understand why the algorithm is correct (conditional on ERH), but I'm confused about why the time complexity is as claimed. In particular, the key step is the computation of a prime $$p \leq n^5$$ satisfying a certain congruence condition mod $$4n$$. ERH guarantees that if you keep choosing numbers satisfying this congruence condition in $$[1, n^5]$$ uniformly and independently at random, then you'll come across a prime in expected $$O(\log n)$$ time. But still, at every trial you must test whether the candidate $$p$$ is prime (in fact I was considering the possibility of seeing what happens later in the algorithm if $$p$$ is composite, hoping that it will terminate quickly with a wrong answer or something, but this doesn't seem to be the case; also Rabin--Shalit specifically writes on page S243 to test it for primality).

So in order to satisfy the desired time complexity, this means that the amortized time complexity of the primality test needs to be $$O(\log n)$$. As far as I know, even if you allow randomization there isn't any known algorithm for primality testing with this complexity (especially in 1986 when this paper was published) unless you allow some constant error probability. Can anyone point out what I am missing?

• I don't have access to that article. In this article, an improved version of the algorithm is presented, and here it isn't required that the number found is actually a prime; the density of primes guaranteed by ERH is merely used to prove that the required kind of number, which will typically be a prime, can be found with a certain probability (see p. $6$). Perhaps the same is true for the original algorithm? Dec 10, 2019 at 12:27
• I thought this might work as well, but here's the reason I think it's unlikely: in the improved algorithm of Pollack--Trevino, once you have your prime candidate p, nothing you do with it really relies on it being prime. However, in Rabin's original algorithm (according to the description given in Pollack--Trevino as well), you need to find integers a, b such that p = a^2 + b^2. As far as I know the solution to that one really does rely on p being prime, since it's equivalent to finding a square root of -1 in F_p. Dec 10, 2019 at 12:38
• In fact the standard solution to that problem (finding a square root of -1 in F_p for p = 1 mod 4) uses the Euclidean algorithm in F_p[X] and facts about the Legendre symbol mod p, which makes me think there might be some problems with generalizing to composite p Dec 10, 2019 at 12:42
• But the algorithm of Pollack–Treviño also relies on finding a square root of $-1$. They do this by forming $u^{(p-1)/4}$ for random $u\in[1,p-1]$, using the fact that if $p$ is prime, this is a square root of $-1$ for the quadratic nonresidues, i.e. half the values. Dec 10, 2019 at 13:18
• Okay, I think you are right. In fact, there's no reason you can't just substitute this this step of Pollack--Trevino for the analogous one in Rabin--Shallit. In fact even though Z/n isn't even an integral domain I think the Euclidean algorithm procedure will still terminate (though possibly with a wrong answer or just because it has found a zero-divisor coefficient of a polynomial) in the desired amount of time. Thanks! Dec 11, 2019 at 2:08

As discussed in the comments, apparently the number used need to be a prime. We only need to find a square root of $$-1$$ modulo that number, and the prime density from ERH is used only to show that we have a certain probability of doing so.