ある数を法とする平方剰余から元の平方根を求めるプログラム

def ugojo(a,b):
    if a<0:
        a*=-1
    if b<0:
        b*=-1
    if a>b:
        a,b=b,a
    while a>0:
        a,b=b%a,a
        #print(a,b)
    return b
def xgcd(a,b):
    x0,y0,x1,y1=1,0,0,1
    while b!=0:
        q,a,b=a//b,b,a%b
        x0,x1=x1,x0-q*x1
        y0,y1=y1,y0-q*y1
    return a,x0,y0
def modinv2(a,m):
    a=a%m
    g,x,y=xgcd(a,m)
    if g!=1:
        if m==0:
            return 'not exist'
        d=ugojo(a,m)
        a,m=a,m//d
        return d*modinv2(a,m)
    else:
        return x%m
if 1:
    def pmoti(p,sisu,a):
        f=[]
        rev=[]
        for i in range(p):
            if (i**2-a)%p==0 or (i**2+a)%p==0:
                if (i**2-a)%p==0:
                    f+=[i]
                    rev+=[False]
                else:
                    f+=[i]
                    rev+=[True]
        f=list(set(f))
        #print(f)
        if sisu==1:
            return f
        else:
            b=[modinv2(2*f[j],p) for j in range(len(f))]
            #print(b,'b')
            for i in range(len(f)):
                for j in range(2,sisu+1):
                    if rev[i]:
                        f[i]=(f[i]-(f[i]**2+a)*b[i])%p**j
                    else:
                        f[i]=(f[i]-(f[i]**2-a)*b[i])%p**j
            f=list(set(f))
            return f
    def kon2(p,sisu,a,num,rawmode=False):
        d=pmoti(p,1,a%p)
        d+=pmoti(p,1,(num-a)%p)
        d=list(range(p))
        rev=[]
        for i in range(len(d)):
            if (d[i]**2*(num)-a)%(p)==0:
                rev+=[False]
            else:
                rev+=[True]
        e=[]
        nnp=1
        aap=a%p
        for k in range(sisu):
            f=[]
            app=a
            jj=0
            if num%p==0:
                nump=num
                while nump%p==0:
                    nump//=p
                    jj+=1
            nnnp=(num//p**jj)%p**(k+1)
            if nnnp<p**(k) and nnnp!=num and k>0:
                nnnp+=p**(k)
            ug=ugojo(nnnp,app)
            if ug>1:
                np=(nnnp)//ug
                ap=app//ug
                if k>1:
                    f+=kon2(p,k+1,ap,np)
                else:
                    f+=kon2(p,k+1,ap,np)
                    f+=pmoti(p,1,ap)
                    f+=pmoti(p,1,np-ap)
            else:
                ug=1
                ap=app
                np=nnnp
            for q in range(p):
                for i in range(len(d)):
                    j=0
                    ks=1
                    pp=0
                    r=d[i]
                    while j<p*2:
                        t=d[i]+j*p**k
                        if abs(t**2//nnp*nnp-t**2//np*np-q*p**k)==0:
                            f+=[t%(p**(k+1)*np)*ug]
                            
                        r+=p**pp*ks
                        
                        if ((r**2-ap%p**(k+1)-q*p**k))%np==0 or ((r**2+ap%p**(k+1))%p**k-q*p**k)%np==0:
                            f+=[r%(p**(k+1)*np)]
                            break
                        if (((r**2-ap-q*p**k))%p**(pp)==0 or ((r**2+ap-q*p**k))%p**(pp)==0) and pp<k:
                            pp+=1
                        if (((r**2-ap-q*p**k))%np==(ap-aap) or ((r**2+ap-q*p**k))%np==(ap-aap))and ap!=aap:
                            u=ap-aap
                            while u%p==0 and u>1:
                                u//=p
                            if ks%u!=0:
                                ks*=u
                        if (((r**2-ap%p**(k+1)-q*p**k))%np)==0 or (((r**2+ap%p**(k+1)-q*p**k))%np)==0:
                            l=0
                            while l<p*2:
                                s=r+j*l*p**(k)
                                if ((s**2-ap%p**(k+1)-q*p**k))%(np)==0:
                                    f+=[(s)%(p**(k+1)*(np))*ug]
                                    break
                                if ((s**2+ap%p**(k+1)-q*p**k))%(np)==0:
                                    f+=[(s)%(p**(k+1)*(np))*ug]
                                    break
                                if ((s**2-ap%p**(k+1)-q*p**k))%(np)==num:
                                    f+=[(s)%(p**(k+1)*(np))*ug]
                                    break
                                if ((s**2+ap%p**(k+1)-q*p**k))%(np)==num:
                                    f+=[(s)%(p**(k+1)*(np))*ug]
                                    break
                                l+=1
                            else:
                                j+=1
                                continue
                            break
                        j+=1
                    
            d=list(set(f))
            nnp=np
            aap=ap
            e=[d[i]%np for i in range(len(d))]+[(np-d[i])%np for i in range(len(d))]
            d=[d[i]%num for i in range(len(d))]+[(np-d[i])%num for i in range(len(d))]
            d=list(set(e+d))
        if rawmode:
            return d
        d=[d[i]%num for i in range(len(d))]
        d=list(set(d))
        for i in range(len(d)):
            d+=[(num-d[i])%num]
        d=list(set(d))
        i=0
        while i <len(d):
            if (d[i]**2-a)%num!=0 and (d[i]**2+a)%num!=0:
                del d[i]
            else:
                i+=1
        return d
        print(d)
    acou=0
    bcou=0
    lislis=list(range(1,365))
    for i in range(364):
        lis=kon2(19,3,1+i,365)
        for j in range(len(lis)):
            if lis[j] in lislis:
                del lislis[lislis.index(lis[j])]
        if len(lis)>0:
            acou+=1
            bcou+=len(lis)
        print(lis,1+i)
    print(acou,bcou,len(lislis),lislis)

大きな数を法とする平方剰余の平方根を求める問題は一般に、素因数分解と同じくらい難しい問題と見られているようです。また、素因数分解の問題とどちらかが効率的に解くことができればもう一方も同じく効率的に解くことができるという関係性にあるらしいです。

今回は、ヘンゼルの補題を用いた根の持ち上げを使って合成数の平方剰余の平方根を求めるプログラムを作ってみました。ただし、4桁くらいまでの数までしか求められないようです。

このプログラムの特徴
①平方剰余と、平方剰余を足したらmod(合同式)が0になるという意味で、平方剰余と対になる数を一緒に求めています。
②求めたい大きな数とは互いに素となる、11や17などの小さな素数の合同式を法に計算しています。求めたい数の合同式との関係性は以下の式で計算することができます。

$${a \hspace{2pt} mod \hspace{2pt} N \equiv([a/N(aのNによる商:整数)] *N )\hspace{2pt} mod \hspace{2pt} p\hspace{2pt}+a \hspace{2pt} mod \hspace{2pt} p}$$

平方剰余(とその対になる数)の分布(365を法にした例で1から20まで)

[1, 27, 74, 173, 192, 291, 338, 364] 1
[] 2
[]3
[2, 19, 54, 148, 217, 311, 346, 363] 4
[] 5
[33, 113, 131, 161, 204, 234, 252, 332]6
[] 7
[] 8
[3, 81, 143, 154, 211, 222, 284, 362] 9
[]10
[]11
[]12
[] 13
[] 14
[] 15
[4, 38, 69, 108, 257, 296, 327, 361]16
[] 17
[]18
[47, 101, 172, 174, 191, 193, 264, 318]19
[]20