F legendre(a, p)
   R pow(a, (p - 1) I/ 2, p)

F tonelli(n, p)
   assert(legendre(n, p) == 1, ‘not a square (mod p)’)
   V q = p - 1
   V s = 0
   L q % 2 == 0
      q I/= 2
      s++
   I s == 1
      R pow(n, (p + 1) I/ 4, p)
   V z = 2
   L
      I p - 1 == legendre(z, p)
         L.break
      z++
   V c = pow(z, q, p)
   V r = pow(n, (q + 1) I/ 2, p)
   V t = pow(n, q, p)
   V m = s
   V t2 = BigInt(0)
   L (t - 1) % p != 0
      t2 = (t * t) % p
      V i = 1
      L(ii) 1 .< m
         I (t2 - 1) % p == 0
            i = ii
            L.break
         t2 = (t2 * t2) % p
      V b = pow(c, Int64(1 << (m - i - 1)), p)
      r = (r * b) % p
      c = (b * b) % p
      t = (t * c) % p
      m = i
   R r

V ttest = [(BigInt(10), BigInt(13)), (BigInt(56), BigInt(101)), (BigInt(1030), BigInt(10009)), (BigInt(44402), BigInt(100049)),
           (BigInt(665820697), BigInt(1000000009)), (BigInt(881398088036), BigInt(1000000000039)),
           (BigInt(‘41660815127637347468140745042827704103445750172002’), BigInt(10) ^ 50 + 577)]
L(n, p) ttest
   V r = tonelli(n, p)
   assert((r * r - n) % p == 0)
   print(‘n = #. p = #.’.format(n, p))
   print("\t  roots : #. #.".format(r, p - r))