ttest:=T(T(10,13), T(56,101), T(1030,10009), T(44402,100049),
   T(665820697,1000000009), T(881398088036,1000000000039),
   T("41660815127637347468140745042827704103445750172002", BN(10).pow(50) + 577),
   T(1032,10009) );
foreach n,p in (ttest){ n=BN(n);
   r:=tonelli(n,p);
   assert((r*r-n)%p == 0,"(r*r-n)%p == 0 : %s,%s,%s-->%s".fmt(r,n,p,(r*r-n)%p));
   println("n=%d p=%d".fmt(n,p));
   println("   roots: %d %d".fmt(r, p-r));
}