フェルマーテストを用いてa^2-b^2=Nとなる組み合わせをpythonで探す。

久しぶりに素因数分解の話。

フェルマーテストとは、

$${a^{\tfrac{n-1}{2}}mod \quad n=1}$$

をnが素数の時に満たせば、aはnに対する平方剰余である。つまり、nを底とする剰余で2乗するとaになるものが存在することになることを利用して、ある数が素数を底として、また別のある数に対する平方数であるかどうかを判定する方法です。
ちなみに、全ての平方数は、全ての素数を底とする平方剰余であることは簡単に示すことができます。

今回は、平方数のこの性質を用いて、a^2-b^2=Nとなるabの組み合わせを探すプログラムを作ってみたので紹介します。

N=717651*131651
a=13537924+3

ls=[0 for i in range(len(pr))]
ls2=[0 for i in range(len(pr))]
ls3=[[0]+[j for j in range(1,pr[i])if fermat(j,pr[i])==1]for i in range(len(pr))]
print(ls3) 
def fermat(a,n):
    if a>0 and bekiyolists(n,a%n,[(n-1)//2])[0]==1:
        return 1
    else:
        if ugojo(a,n)>1:
            #print(a)
            ug=ugojo(a,n)
            while a%ug**2==0:
                a//=ug**2
            if bekiyolists(n,a,[(n-1)//2])[0]==1:
                return 1
        return 0 
def combin(a,d=[],c=[]):
    dd=[]
    if len(a)==0:
        dd+=[c]
    else:
        for j in range(len(a[0])):
            dd+=combin(a[1:len(a)],d,c+[a[0][j]])
    return dd
for i in range(10):
    
    ls=[fermat(a,pr[j]) for j in range(len(pr))]
    ls2=[fermat(a+N,pr[j]) for j in range(len(pr))]
    print(ls,ls2,a,a+N)
    if sum(ls)==len(ls) and sum(ls2)==len(ls2):
        print('is sqrt',np.sqrt(a),np.sqrt(a+N))
        a+=0
        ls=[fermat(a,pr[j]) for j in range(len(pr))]
        ls2=[fermat(a+N,pr[j]) for j in range(len(pr))]
        print(ls,ls2,a)
        print(np.where(np.array(bbmx)==len(ls))[0])
        break
    b=1
    bg=1
    bbls=[[] for j in range(len(pr))if ls[j]==1 and ls2[j]==1]
    bbor=[pr[j] for j in range(len(pr))if ls[j]==1 and ls2[j]==1]
    jj=0
    for j in range(len(pr)):
        if ls[j]==1 and ls2[j]==1:
            bb=a%pr[j]
            bb2=(a+N)%pr[j]
            bb3=ls3[j].index(bb)
            for k in range(len(ls3[jj])+1):
                bb4=(ls3[j][(bb3+k)%len(ls3[j])]-bb)%pr[j]
                if (bb2+bb4)%pr[j] not in ls3[j]:
                    continue
                bbls[jj]+=[bb4]
            jj+=1
    print(bbls)
    bbcm=combin(bbls)
    #print(bbcm)
    bbmx=[]
    bbb=[]
    bgt=np.prod(np.array(pr))
    for j in range(len(bbcm)):
        bg=1
        btm=1
        for k in range(len(bbcm[j])):
                #print(bb4,'bb4',pr[j],k)
                
                gc=xgcd(bg,bbor[k])
                    #print(gc,b,bb4,bb,bb2)
                btm=(bbcm[j][k]*bg*gc[1]+btm*bbor[k]*gc[2])%(bg*bbor[k])
                    #print(b)
                bg*=bbor[k]
        lstm=[fermat(a+btm,pr[j]) for j in range(len(pr))]
        lstm2=[fermat(a+btm+N,pr[j]) for j in range(len(pr))]
        #print(lstm,lstm2,np.sqrt(a+btm),np.sqrt(a+btm+N))
        bbmx+=[sum(np.array(lstm)*np.array(lstm2))]
        bbb+=[btm]
        if bbmx[-1]==len(ls):
            sq=np.sqrt((a+bbb[-1]))
            sq2=np.sqrt((a+bbb[-1]+N))
            if abs(sq-int(sq))<1.0e-5 and abs(sq2-int(sq2))<1.0e-5:
                print('sqrt',bbb[-1]%N,sq,sq2,bbb[-1],ugojo(abs(int(sq)**2-int(sq2)**2),N))
        
    #b=np.prod([pr[j]**(ls[j]*ls2[j]) for j in range(len(pr)-2)])
    #print(b)
    a+=bbb[bbmx.index(max(bbmx))]
        
    #print(a)
        
print(a)

いろいろと余分な物が入っていますが、とりあえず上記のプログラムで動かすことができます。

初期値からの加算値を特定の順序に従って決定します。まずaとa+nを用意して、小さないくつかの素数を用意します。次に、aとa+nの両方でフェルマーテストを満たす素数をリストアップします。リストアップした素数のそれぞれについて、aとa+nの共通する加算値を調べて、両者のフェルマーテストを両方満たす剰余のリストを生成します。
リストアップした全ての剰余について、それぞれひとつを取った組み合わせを作ります。全ての組み合わせについて中国の剰余定理で、全ての剰余を満たす自然数を生成します。
こうして作った加算値は全て、リストアップした全ての素数について、aとa+nの両方でフェルマーテストを満たすことがわかります。この中のいくつかの数の組み合わせは、リストアップされていない他の素数についてフェルマーテストを両者満たす可能性があります。

フェルマーテストを満たす素数の数が1番多いを採用してしばらく更新を続けていれば、$${a^2-b^2=N}$$を満たす数が範囲内にあればこれを満たす数値の組み合わせになる可能性が高くなります。

$${a^2-b^2=N}$$を満たす組み合わせでは、必ずaとbがそれぞれフェルマーテストを満たすことになります。aとbの組み合わせがわかれば、高い確率で、nが合成数であれば素因数分解をすることができます。

小さな合成数では早く動きますが、大きな合成数になると、範囲に掛からなくなる可能性が高く、また、調べる素数を大きくすればそれだけ、組み合わせ数が増えるため、全ての場合を調べるのに時間がかかる可能性が出てきます。
これが、大きな合成数の素因数分解にこの手法を使う場合の課題になります。

加算値の制約をnの倍数に制限すれば、さらに効果的なプログラムになりそうな気がします。

(pythonの言語としての計算の制約に注意。あまりにも大き過ぎる数値は掛け算さえままならなくなる。)

↓いつものモジュール群

import numpy as np

def bekiyolists(bunkaiN,tei0,bekilist):
    #bunkaiN=497
    #tei=4
    #beki=12
    beki=len(bekilist)
    bekiyolist=[0]*beki
    blist=[]
    tei=tei0
    l=0
    bekiloglist=[[]]*beki
    mx=0
    #print(bekilist)
    for j in bekilist:
        #print(j)
        beki2=j
        #if l%100==1:
            #print(l)
        i=0
        listes=[]
        while beki2>0:
            listes.append(beki2%2)
            #blist.append(beki2%2)
            beki2=beki2//2
            i+=1
        bekiloglist[l]=listes
        if mx<i:
            mx=i
        l+=1
    #print(bekiloglist)
    result=[1 ]*beki
    for i in range(mx):
        for j in range(beki):
            if len(bekiloglist[j])>i:
                result[j]=(result[j]*tei**bekiloglist[j][i])%bunkaiN
        tei=(tei**2)%bunkaiN
    bekiyolist=result
    #print(result)
    #print(beki,blist)
#blist.reverse()
        #result=1

        #for bv in blist :
            #if bv==1:
                #result=(result*tei)%bunkaiN
            #tei=(tei*tei)%bunkaiN
        #bekiyolist[k]=result
        #k+=1

        #print(tei,result)
    return bekiyolist
    
    
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 ugojod(a,b):
    if a<0:
        a*=-1
    if b<0:
        b*=-1
    if a>b:
        a,b=b,a
    while a>0:
        c=int(b)//int(a)
        a,b=b-c*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):
    g,x,y=xgcd(a,m)
    if g!=1:
        return 'not exist'
    else:
        return x%m