Shank -tonelli

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