def residuo_non_quadratico(p):
m=2
while(pow(m,(p-1)/2)%p==1):
m=m+1
return m
def Shank_Tonelli(p,a):
if(pow(a,int((p-1)/2),p)==1):
m=p-1
s=0
while(m%2==0):
m=m/2
s=s+1
z=residuo_non_quadratico(p)
c=pow(z,int(m),p)
t=pow(a,int(m),p)
r=pow(a,int((m+1)/2),p)
while(True):
if(t==0):
return 0
elif(t==1):
return r
else:
i=1
while(pow(t,int(pow(2,i,p)),p)!=1):
i=i+1
b=pow(c,int(pow(2,s-i-1)),p)
s=i
c=(b*b)%p
t=(t*c)%p
r=(r*b)%p
(Credo di aver risolto, non so come eliminare il topic)