【できました】素因数分解界のゲームチェンジャー(になるかもしれない)アルゴリズム

こないだの記事から派生させて作っていたプログラムがついに素因数分解をしました！

ということで素案のプログラムを書き出しておきます！

n=103*983のそのままの素因数分解の流れ

def kyoritan(pa,uug,ko=0):
    #ko=0
    #print(pa,uug)
    for i in range(27,-1,-1):
        if abs(ko*(ko+uug)+uug-pa)>abs((ko+2**i)*(ko+2**i+uug)+uug-pa):
            ko+=2**i
        elif ko>2**i and abs(ko*(ko+uug)+uug-pa)>abs((ko-2**i)*(ko-2**i+uug)+uug-pa):
            ko-=2**i
        #print(uug,ko,i,abs(ko*(ko+uug)+uug-pa),abs((ko+2**i)*(ko+2**i+uug)+uug-pa))
    return ko,ko*(ko+uug)+uug-pa 
def apbasb(apb,asb,cut1,cut2,mins=-1):
    resls=[]
    kols=[]
    cutls=[]
    matchls=[]
    for i in range(-800,800):
        res=kyoritan(apb-i*cut1,asb+i*cut1*cut2,1)
        if res[0]==0:
            continue
        resls+=[res]
        kols+=[abs(res[1]-mins)]
        cutls+=[i]
        if res[1]==mins:
            print(res,i)
            matchls+=[[res,i,i*cut1,i*cut2*cut1]]
            #break
        #print(res,i)
    ind=kols.index(min(kols))
    #print(kols)
    return resls[ind][1],cutls[ind],cutls[ind]*cut1,cutls[ind]*cut2*cut1,matchls
ti=time.time()
n=103*983
#print(n,apbasb(n,40,1,-1,1))
#print(kyoritan(n-20,40-20))

#print(xs)
nn=n
bls=[]
while nn>0:
    bls+=[nn%2]
    nn//=2
#bls[0]-=1
bbls=bls.copy()
k=0
rem=0
kake1=0
kake2=0
#del bls[-1]
ll=len(bls)//2
prekake1=0
prekake2=0
prekake3=0
pre2k=1
kake3=1
hg=1
ndif=n
res=[0,0]
print(bls,n)
while ll>0:
    k+=1
    #if sum(bls[2*ll:len(bls)])==0:
#        ll-=1
#        k=0
#        continue
    ndif=n-kake1*kake2*2**(2*ll)
    apb=kake1*kake2#+res[1]
    #apb=0
    cut1=1#kake1-kake2
#    if kake1==kake2:
#        cut1=kake1
    cut2=2**k
    asb=ndif//2**(2*(ll)-k)+1
    #asb=(ndif/2**(2*(ll-k)))%2**k+rem
    #print('k',k,'ll',ll,rem,ndif)
    ret=apbasb(apb,asb+2,cut1,cut2,1)
    res=kyoritan(apb-ret[2],asb+2+ret[3],1)
    #cut3=kake2-ret[1]
#    ret2=apbasb(apb-ret[2],asb+ret[3],cut3,cut2,1)
#    res2=kyoritan(apb-ret[2]-ret2[2],asb+ret[3]+ret2[3],1)
    #print(ret,'asb',asb,'apb',apb,'cut1',cut1,'cut2',cut2,bls,asb+ret[3])
    atai1,atai2=res[0]+1,res[0]+asb+ret[3]+1
    als=[1,-1]
    #kake1,kake2=2**k*(res[0])+1,2**k*(res[0]+asb+ret[1])+1
    kake1,kake2=2**k*(atai1)+1,2**k*(atai2)+1
    for i in range(len(als)):
        for j in range(len(als)):
            if atai1*als[i]+atai2*als[j]==asb+ret[3]:
                #kake1,kake2=2**k*(atai1)+als[j]*1,2**k*(atai2)+als[i]*1
                kake1,kake2=2**k*(atai1)+als[j]*1,2**k*(atai2)+als[i]*1
                print('match',als[i],als[j])
    #kake1,kake2=2**k*(kake2-ret[1])+1,2**k*(kake1-ret2[1])-1
    print('asb',asb,'apb',apb,'kake',kake1,kake2,'atai',atai1,atai2,'zure',apb,-ret[2],',',asb,ret[3],'asb-apb,pair',apb,asb,'apb-n',apb-ret[2],'asb+2n',asb+ret[3])
    kake1,kake2=abs(kake1),abs(kake2)
    #if kake1==kake2 and kake2>1:
        #kake1,kake2=kake1-1,kake1+1
    kake3=kake1+asb+ret[3]-1
    #ndif=n-kake1*kake2*2**(2*(ll-k)-2)
    ndif=n-kake1*kake2*2**(2*(ll-k))
    if kake1>kake2:
        kake1,kake2=kake2,kake1
#        print('change')
    kake3=abs(kake3)
    rem=ndif//(2**(2*(ll-k)))+1
    nnn=n-kake1*kake2*2**(2*(ll-k))
    hg=np.sign(nnn)
    if hg==1:
        nnn=abs(nnn)
    if hg==-1:
        nnn=abs(nnn)+1
    bls=[]
    while nnn>0:
        bls+=[nnn%2]
        nnn//=2
    bls+=[0]*(len(bbls)-len(bls))
    ll-=k
    k=0
    print(bls,ret,res,'rem',rem,(n//(2**(2*ll-2*k)))-kake1*kake2)
    print('k',k,'ll',ll,'hg',hg,'kake1',kake1,'kake3',kake3,'kake2',kake2,'ret',ret[1],'nnn',kumi(bls,2),'ndif',ndif,'rem',rem,kake1*kake2)
    print('2*',atai1,'-+1,','2*(',atai2,'+',asb,'+',ret[3],')-+1','amari',ndif//2**(2*(ll-k)),'*',2**(2*(ll-k)))
    #print(res[0],'+1,(',res[0],'+',asb,'+',ret[1],')-1')
    #print('2*(',kake2,'-',ret[1],')+1,','2*(',kake1,'-',ret2[1],')-1')
    prekake1,prekake2,prekake3,pre2k=kake1,kake2,kake3,2**k
    print(ret,res,atai1,atai2)
    print(kake1+kake2,kake1*2-kake2,kake2-kake1,ugojo(n,kake1+kake2),ugojo(n,kake1*2-kake2),ugojo(n,kake1-kake2),ugojo(n,atai1),ugojo(n,atai2))
print(time.time()-ti)
print(xs)

いまできたばかりで余計なコメントが多いですが、形が整ったら修正します。このコードは2進法をベースに作っています。

基本的な流れ

①与えられた数nを越えない最大の4の累乗数Lを調べる。(103*983=101249の場合は$${4^8}$$。)

②nの$${2^{15}}$$で割った商aを求める。

③aを$${a=2bc+b+c+1}$$を満たすbとcを求める。これは後で説明する特殊関数で大雑把に求めています。

④上で求めたbとcは2倍したaに対して$${2a=4bc+2b+2c+2=(2b+1)(2c+1)+1}$$を満たすので2b+1と2c+1を因数を生成するための候補p、qとする。

⑤nから$${(2b+1)×(2c+1)4^{L-1}}$$を引いた数を、次の候補を調べるためのnとする。

⑥Lの次数を一つづつ下げて調べていき、nと因数を生成するための値pとqを更新する。

順々にLを下げてLが0になった時に余りが0となる数の組み合わせが見つかれば、それが元の101249の因数である。

③のbとcの組み合わせを見つける方法だけ難しく、2次式の因数分解とか他の技術を使わなければならないはずのものですが、総当たりに近いですが、少し効率化する方法を考えました。

bcとb+cは(x-b)(x-c)=0の解として求めることができますが、この式について、$${x^2+αx+β=0}$$について、αに変数を足し込み、βに同じ変数の-2倍を足し込み、この式の解を求めることによって解決してます。

2次方程式の計算に持ち込み、これが成り立つbとcの組み合わせがひとつでも見つかれば次の計算へ進むことができます。

Lの次数が下がるほどに計算は時間がかかるようになり、現実問題どれだけ時間をかける必要があるかわかりませんが、もしかしたら、計算しやすくなるかもしれないと思いますので、記録として残しておきます。

この流れに先日話した他の数の因数分解から導く方法を考えてみようかと思います。③の工程は先日記事にしたような隣数から導く因数分解と同じ考え方であり、そちらの方で置換できそうな気がしますので。

この記事の続きです。