https://en.wikipedia.org/wiki/Cipolla%27s_algorithm describes the steps to compute a modular square root, assuming a prime modulus.
The converse of this algorithm could be: considering a quadratic residue, a failure to compute its square root would indicate that the modulus is not prime.
This looks like a possible primality test.
Steps to check the primality of a odd number p > 3
get the the smallest integer i larger than $ \sqrt p $, which satisfies $ a = i^2 \mod p $ is not a perfect square (a is a (non-trivial) quadratic residue mod p)
get the first integer n larger than i where $ u = (n^2 - a) \mod p $ is not a quadratic residue (checked $ jacobi(u|q) = -1 $)
Now in the field $ F_{p^2} = F_p(\sqrt u) $
let $ w = (n + \sqrt u) $
Try to compute the square root of w
$$ z = w^{(q+1)/2} $$
Then verify z, and curiously there are only 2 possibilities
$ z = (x + y * \sqrt u )$ with either $ y \not = 0 $ , either $ x^2 \not = a \mod p $ and curiously, p seems to always be composite
$ z = (x + 0 * \sqrt u )$ with $ x^2 = a \mod p $ and curiously, p seems to always be prime
I tested all integers p up to $ 2^{48} $ against either a false positive, either a false negative. I checked against all spsps base 2 up to $ 2^{64} $ (http://www.cecm.sfu.ca/Pseudoprimes) I can't find a counter-example for that Cipolla primality test
I did not find many papers related to this, the closest being https://arxiv.org/abs/1307.7920
How can I find at least one pseudo-prime to this test ?
For the sake of convenience, I will redefine variables in this answer. Let $n$ be a composite number. $z^2=b\pmod n$ has a solution ($z$) if and only if for each prime $p | n$, there is a solution to $z^2=b\pmod p$. In the case that the Kronecker symbol $K(a | n)=1$, but there is a prime $p|n$ where $K(a | p)=-1$, then $n$ will be proven composite by this algorithm.
The solution $z$ if it exists, is $(a + \sqrt{a^2-b})^{(n+1)/2} \pmod n$ where Jacobi$(a^2-b|n)=-1$.
The polynomial $x^2 - 2ax + b$ has the roots $x = a ± \sqrt{a^2-b}$.
Therefore, we have, $z=x^{(n+1)/2} = c_1x+c_2 \pmod {x^2-2ax+b,n}$.
These follow from properties of Lucas Sequences. $c_1$ must be $0$ assuming $a^2-b$ is not a square $\pmod n$. The rest of the work seems trivial enough, so this amounts to counterexamples being Lucas pseudo primes for the parameters $(a,b)$.