Is computing a square root of a number and having more than 2 roots a reliable way to decide primality? I was reading CLRS and it mentions that if $p$ is a prime of the form $4k+3$ and $a$ is a quadratic residue, then $a^{k+1}$ is a square root of $a$.  One can also easily show that $a^{-k}$ is a square root of $a$ as well.

Can we use this to get a primality test, for numbers $N$ of the form $N=4k+3$?

After reviewing Miller-Rabin (and after asking the related question Is there any efficient algorithm for primality testing for numbers that are of the form $4k+3$ using the square root function? ) and thinking about it a little more I thought of an algorithm that seems promising.  The intuition is that if $N=pq$ (or some other composite) there will be more than 2 square roots. Therefore, if we can find more than 2 square roots, we have evidence that $N$ is a composite.  With that idea, here is the algorithm:

PrimalityTest$(N)$: (where $N = 4k + 3$)


Choose $x \in \mathbb{Z}_N$ and compute $a = x^2 \pmod N$.
Compute $y = a^{(N+1)/4} \pmod N$.  (Notice that if $N$ is prime, this is just computing a square root of $a$.)
If $y^2 \equiv a \pmod N$ but $y \not \equiv \pm x \pmod N$ then there are more than 2 distinct square roots of $a$. Therefore, return "composite".
If $y^2 \not \equiv a \pmod N$ then return "composite".  (If $N$ were prime, then step 2 would have computed a valid square root of $a$, so by the contrapositive, if $y$ is not a valid square root, then $N$ is not prime.)
If $y^2 \equiv a \pmod N$ and $y \equiv \pm x \pmod N$ we have no idea what $N$ is. Therefore, go back to step 1 and choose a different $x$.  If this step happens too often, simply return "prime".


I believe this will succeed with probability 1/2, since it has half chance of choosing a "bad" square root that gives us no information. The only step I am not sure how to analyze probabilistically is step 4, but it seems that step 2 might do nonsense when $N$ is composite, it seems the probability that 4 catches is is quite high. I'd guess around $\frac{N-2}{N}$.
There are composite numbers where this primality test fails with probability 1: it always outputs "prime", even though the number is actually composite.  These numbers are basically the Carmichael numbers, with a slight twist.  Therefore, this is unfortunately not a correct primality testing algorithm.

Recall that $N$ is a Carmichael number if $x^{N-1}\equiv 1 \pmod N$ holds for all $x \in \mathbb{Z}_N^*$.  All Carmichael numbers are odd.  I'll define a super-Carmichael number to be a number $N$ such that $x^{(N-1)/2}\equiv 1 \pmod N$ holds for all $x \in \mathbb{Z}_N^*$.

Now consider what happens in your algorithm if $N$ is a super-Carmichael number.  We pick $x$ randomly, then compute

$$y=a^{(N+1)/4}=(x^2)^{(N+1)/4}=x^{(N+1)/2} = x^{(N-1)/2} \cdot x = x \pmod N,$$

since $x^{(N-1)/2} = 1 \pmod N$.  It follows that we also have $y^2 = x^2 \pmod N$.  Thus, no matter what value of $x$ you choose, you'll always find that $y^2 \equiv a \pmod N$ is true and that $y \equiv \pm x \pmod N$ (when $N$ is a super-Carmichael number).  So, your algorithm will always give an incorrect answer, when you run it on a super-Carmichael number.

For instance, $1729 = 7 \times 13 \times 19$ is a super-Carmichael number: we have $x^{36} \equiv 1 \pmod{1729}$ for all $x \in \mathbb{Z}_{1729}^*$, and $36$ divides $(1729-1)/2 = 864$.  Therefore, your algorithm gives an incorrect answer on 1729.  I suspect there are infinitely many super-Carmichael numbers.  I also suspect that your algorithm might actually give an incorrect answer on all Carmichael numbers (not just the super-Carmichael ones), but I haven't tried to check the details to see if that's actually true.