モジュラ逆数の冪剰余演算に関連する検証結果

モジュロ演算を主に使う素数の判定法であるフェルマーテスト。

これをちょこっと調べていたのですが、意味合いをすぐには判断できない結果を見つけたので、記録しておきます。

ある数nと、互いに素な数aについて、nを法とするaのモジュラ逆数をa^-1としたとき、

$${a^{-1*b}(mod n)=a^{b} (mod n)}$$のとき、2bは法をnとするaの冪剰余の、周期の倍数である。

一応、カーマイケル数についても、周期発見ができるっぽい。
大きな数の素因数分解も、モジュラ逆数からそれを基にしてさらに順に小さいモジュラ逆数を求めて、小さいモジュラ逆数から順に周期発見を繰り返せば、大きな合成数の周期発見ができそうな気がしなくもない…。

これから検証予定の素因数分解アルゴリズムの手順(4087の因数を求める手順を例に考える。)
①求めたい数と互いに素な自然数をひとつ決める。(今回は5)
②求めたい数を法とする、決めた自然数に対するモジュラ逆数を求める。(4087のモジュラ逆数は1635)
③求めたモジュラ逆数を新しく法にして、5に対するモジュラ逆数を、モジュラ逆数が十分に小さくなるまで順に求めて、数列として記録しておく。(モジュラ逆数は、拡張ユークリッド互除法を使って求めることができます。)ただし、決めた自然数、今回では5と割り切れる場合にはモジュラ逆数が求められなくなるので、その場合は法とする数の方を5で割って、新しく法とする。(大きい方から順に、4087、1635→807、323、194、39、8となる。)
④それぞれのモジュラ逆数の数列について、小さい方から順にひとつ隣の大きい数値を法に、冪剰余の周期を求めて、この数字を大きくさせる。ここで予想では、小さい方から順に求めると、後の数字が全て、先の数字の倍数になっているというものだけど、ほんとに上手くいくのかしらん？？

8 39 : 4 (×4)
39 194 : 96 (×24)
194 323 : 144 (×3/2)
323 807 : 402 (67/24)
1635 4087 : 330 (55/67)

うーん、微妙…関係なくもないとは思うけど…。隣合う数同士が最小公倍数を共有してる気もする…。

ある数の倍数であることと、冪剰余がその数の冪剰余の周期の倍数の周期を持つこととは等価であることと考えられるため、モジュラ逆数の関係にある数同士の冪剰余の周期には、全て元の数の周期(この場合では5の周期と同じ周期)の倍数であるという共通項が見られると考えられるのでしょう。

ちなみに、周期を見つけた後は、周期が偶数であれば、素因数分解は簡単にできます。
$${a^{b/2}-1,a^{b/2}+1}$$の2つの数に高い確率で因数が分かれて入ることがわかっているので、ユークリッドの互除法を使って、元の数との最小公倍数を求めれば、2つの因数を得ることができます。

追記、ただいまプログラムを作成中なのですが、ある大きな数についての因数を探す問題を、この方法にてアプローチする方法を試しています。
ある互いに素な数nとaについて、aとモジュラ逆数の関係にある数a^-1の、nを法とする冪剰余の結果が一致するための数は、常にある数の倍数のみに限られる。この数をbとすると、a^b(mod n)≡a^-b(mod n)≡a-1であり、a^2b(mod n)≡a^-2b(mod n)≡1である。
また、この時のbを使って、a^-b(mod n)-1またはa^-b(mod n)+1はnをふたつの因数に分けることができる。
今、nがある数mに対するaのモジュラ逆数a^-1(mod m)≡nであったと仮定すると、nの周期を発見することにより、mを2つの因数に分けることができるはず。
ある大きい数について、適当な数に対するモジュラ逆数を求める。その数が仮に素因数分解が全てできたとすると、

…さらに追記。詳しくは別に記事を書きますが、うまく周期発見を進めることが出来そうなプログラムを作ることができました。とりあえずは宣伝のため、プログラムのみはっ付けます。

import numpy as np
import cmath
import math
import random
import time
def ugojo(a,b):
    if a>b:
        a,b=b,a
    while a>0:
        a,b=b%a,a
        #print(a,b)
    return b
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*tei)%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 bekiyolist(bunkaiN,tei,bekimax):
    #bunkaiN=497
    #tei=4
    #beki=12
    beki=bekimax
    #bekiyolist=[0]*(beki+1)
    bekiyolist=[0]*1
    blist=[]
    j=0
    result=1
    bekiyolist[0]=1
    pre=0
    while j<beki:

        #if j%10000000==1:
            #print(j)
        result=(result*tei)%bunkaiN
        #bekiyolist[1+j]=result
        if result>2**700 and j%50==1:
            print(result-pre)
        if result<2**900:
            bekiyolist.append(result)
        j+=1
        if result==1:
            print(j)
        #print(tei,result)
    return bekiyolist
def modinv(a,n):#tukaimasenn
    #a>n
    inv=0
    for i in range(1,n):
        if a*i%n==1:
            inv=i
            break

    return inv
def xgcd(a,b):#net yori inyou kansha
    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):#net yori inyou kansha
    g,x,y=xgcd(a,m)
    if g!=1:
        return 'not exist'
    else:
        return x%m

def shukisuitei(n,tei):#main program!!
    obj=n
    inv=modinv2(tei,obj)
    jougen=3
    beki=bekiyolists(obj,tei,[obj])[0]
    #sab=(inv*beki)%obj
    numlis=[]
    itilis=[]
    modlis=[]
    ook=20000
    lists=bekiyolist(obj,tei,jougen)
    picks=[]
    iss=[]
    for i in range(19960,ook):
        #if lists[i]<ook:
        if ugojo(i,tei)==1 and ugojo(i,obj)==1:
            picks.append(i)
            iss.append(i)
    print(inv)
    for k in range(len(picks)):
        i=picks[k]
        iti=1
        if ugojo(obj,i)>1 or i%2==0 or ugojo(i,tei)>1:
            continue
        for j in range(2,i):
            if bekiyolists(i,tei,[j])[0]==1:
                iti=j
                break
        #sabmodi=sab%i
        bekii=bekiyolists(i,tei,[obj%iti])[0]
        bekimodi=beki%i
        if ugojo(bekimodi,i)>1:
            continue
        ch=(inv)
        chi=ch%i
        chn=ch%obj
        chin=(ch*i)%(obj*i)
        #print(ch,chi,chn,'ch',i,iss[k])
        for j in range(iti):
            samp=(i*bekiyolists(i*obj,tei,[j])[0])%(i*obj)
            if samp!=chin:
                continue
            else:
                numlis.append(i)
                itilis.append(iti)
                modlis.append(j)
                print('find!')
                break
        chin=(ch)%(i)                
        #print(ch,chi,chn,'ch',i,iss[k])
        for j in range(i):
            samp=(bekiyolists(obj,tei,[j])[0])%(i)
            if samp!=chin:
                continue
            else:
                numlis.append(i)
                itilis.append(iti)
                modlis.append(j)
                print('find! ex!')
    tmpbai=1
    print(numlis)
    print(itilis)
    print(modlis)
    return    
objnum=18757
tei=3
print(shukisuitei(objnum,tei))

実行結果は以上のスクショのようになります。説明は記事を分けます。