(* Legendre symbol *)
LegendreSymbol[a_, p_] := PowerMod[a, (p - 1)/2, p]

(* Tonelli-Shanks Algorithm *)
Tonelli[n_, p_] := Module[
  {q, s, z, c, r, t, m, t2, i, b},
  (* Check that n is a quadratic residue mod p *)
  If[LegendreSymbol[n, p] =!= 1,
   Message[Tonelli::notsquare, n, p]; Abort[]];

  (* Step 1: Factor out powers of 2 from p - 1 *)
  q = p - 1;
  s = 0;
  While[EvenQ[q], q = q/2; s++];

  (* Special case when s == 1 *)
  If[s == 1,
   Return[PowerMod[n, (p + 1)/4, p]]
  ];

  (* Find a non-residue z *)
  For[z = 2, z < p, z++,
   If[LegendreSymbol[z, p] == p - 1, Break[]]
  ];

  (* Initialize variables *)
  c = PowerMod[z, q, p];
  r = PowerMod[n, (q + 1)/2, p];
  t = PowerMod[n, q, p];
  m = s;

  (* Main loop *)
  While[Mod[t - 1, p] =!= 0,
   t2 = Mod[t^2, p];
   For[i = 1, i < m, i++,
    If[Mod[t2 - 1, p] === 0, Break[]];
    t2 = Mod[t2^2, p];
   ];
   b = PowerMod[c, 2^(m - i - 1), p];
   r = Mod[r*b, p];
   c = Mod[b^2, p];
   t = Mod[t*c, p];
   m = i;
  ];
  r
]

(* Test cases *)
ttest = {
   {10, 13},
   {56, 101},
   {1030, 10009},
   {44402, 100049},
   {665820697, 1000000009},
   {881398088036, 1000000000039},
   {41660815127637347468140745042827704103445750172002, 10^50 + 577}
};

Do[
 n = ttest[[i, 1]];
 p = ttest[[i, 2]];
 r = Tonelli[n, p];
 Assert[Mod[r^2 - n, p] == 0];
 Print["n = ", n, ", p = ", p];
 Print["\troots : ", r, " ", p - r],
 {i, Length[ttest]}
]