A simple implementation of the Tonelli-Shanks algorithm to compute a square root in Z/pZ where p is prime. It could probably be made quite faster by using a faster pow_mod function instead of the recursive one and also by trying to avoid some of the modulus calculations.

tonelli_shanks.c

long pow_mod(long x, long n, long p) {
  if (n == 0) return 1;
  if (n & 1)
    return (pow_mod(x, n-1, p) * x) % p;
  x = pow_mod(x, n/2, p);
  return (x * x) % p;
}

/* Takes as input an odd prime p and n < p and returns r
 * such that r * r = n [mod p]. */
long tonelli_shanks(long n, long p) {
  long s = 0;
  long q = p - 1;
  while ((q & 1) == 0) { q /= 2; ++s; }
  if (s == 1) {
    long r = pow_mod(n, (p+1)/4, p);
    if ((r * r) % p == n) return r;
    return 0;
  }
  // Find the first quadratic non-residue z by brute-force search
  long z = 1;
  while (pow_mod(++z, (p-1)/2, p) != p - 1);
  long c = pow_mod(z, q, p);
  long r = pow_mod(n, (q+1)/2, p);
  long t = pow_mod(n, q, p);
  long m = s;
  while (t != 1) {
    long tt = t;
    long i = 0;
    while (tt != 1) {
      tt = (tt * tt) % p;
      ++i;
      if (i == m) return 0;
    }
    long b = pow_mod(c, pow_mod(2, m-i-1, p-1), p);
    long b2 = (b * b) % p;
    r = (r * b) % p;
    t = (t * b2) % p;
    c = b2;
    m = i;
  }
  if ((r * r) % p == n) return r;
  return 0;
}

Thanks!

I rewrited this to js BigInt -> https://github.com/gkucmierz/utils/blob/main/src/tonelli-shanks.mjs#L8

But why inner pow_mod in this line ?
long b = pow_mod(c, pow_mod(2, m-i-1, p-1), p);

2 ** (m - i - 1) is a bit shift, where m is initially amount of trailing zero bits in p (except the lowest bit), so 2 ** (m - something) can never exceed p - 1 (m stays same or decreases on each iteration).