モジュロ演算を使う冪剰余の周期を求めるプログラムを作ってみた。
モジュロ演算の冪剰余で出てくる、冪乗余の周期を求めるプログラムを作ってみましたので公開します。
モジュロ演算とは、合同演算とも呼ばれており、全ての数を、ある数で割った余り同士と考えて、それぞれの数を掛け算したり、割り算する計算方法です。式で表すとこんな感じ。
$${74 \equiv 4(mod 70)}$$
この式は、74を70で割ると、4に等しいことを表しています。
このあまりを表す数値同士は、剰余演算の等価性と呼ばれる性質のお陰で、互いに足し算したり、掛け算したりしても、答えが変わらないことが保証されています。
$${a(mod N)+b(mod N) \equiv a+b(mod N)}$$
$${(a(mod N)×b(mod N))(mod N)\equiv a+b(mod N)}$$
代表的な式が以上の式です。これを使って、同じ数で割った余り同士では、四則演算をすることができます。
モジュラ逆数
モジュロ演算で余り同士の計算をする時に使う大事な概念があります。それが、モジュラ逆数です。
すなわち、モジュロ演算で言うところの自然数の逆数、ある数とその数のモジュラ逆数を掛けると、1(mod N)になる数です。これの何が良いかというと、四則演算で逆数を掛ければ割り算になるように、モジュロ演算で言うところの割り算を行うことができます。
$${(a^{-1}(mod N)×2a(mod N))(mod N)\equiv 2 (mod N)}$$
分数にこそなりませんが、四則演算との違いに気をつけて使えば、大きくモジュロ演算の幅を広げることができます。ちなみに、モジュラ逆数は、拡張ユークリッド互除法という方法を用いて、少ない計算量で高速に求めることができます。
モジュロ演算の冪剰余
冪剰余とは、ある数を指数で冪乗した余りを求める計算です。モジュロ演算での冪剰余は、剰余演算の等価性を使って高速に計算することができます。
$${a^{x}(mod N)}$$
演算の速さは、数の大きさに対して指数関数的に速くなるので、ほとんど無視できる早さで演算することができます。フェルマーの小定理と言われる定理により、ある自然数の冪乗を、それと互いに素な自然数で割った余りは、指数がある程度以上になると元の数値に戻り、指数が大きくなる程に循環する性質があることが知られています。
$${a^{x}\equiv a(mod N)}$$
また、どの自然数の冪剰余も、0乗は1ということで、1(mod N)を経路として通ります。1(mod N)は、1を引くとNの倍数になることが自然と導かれるため、モジュロ演算では大事な数となっています。循環する周期を求めて、もしそれが偶数であった場合には、1を引いて$${(a^{x/2}-1)(a^{x/2}+1)}$$のように因数分解をすると、両者にNの元の因数が分かれて入ることがあることがわかっているため、大きな数の因数分解を行う際に、大事な要素となるとして、注目されているかと思います。
今回の記事では、この冪剰余の周期を求めるためのプログラムを書いてみました。
プログラム
いつも通りのpythonのプログラムで、以下のようなスクリプトになりました。
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))前々から作っている冪剰余の関数、ユークリッド互除法の関数に加え、拡張ユークリッド互除法を用いて、モジュラ逆数を求める関数を加えてあります。
拡張ユークリッド互除法については、このページのものを借用していると思います。参考として使わせて頂きました。
メインプログラムは一番下のプログラムです。求めたい冪乗余をa(mod N)とした時に、最初にaとNが互いに素であるかの簡単なチェックをした後、それとは別に小さな数iをいくつか生成して、Nとの積を取り、さらにいくつかの簡単な数を指数とした、a(mod N)の冪乗を生成します。一方で、a(mod N)のモジュラ逆数を生成します。($${a^{-1}(mod N)}$$)
両者に生成した数iを掛けて、さらにi×Nで割った余りを比較します。両者の余りが0になった時のみ、その時のaの指数を記録します。
この数については、求める周期を割り切れる数の場合にのみ、両者が一致して記録されるようになっているようです。
プログラムでは引き続き、実験的にモジュラ逆数とaの冪乗余をiで割るだけの計算も進めています。値は少々崩れることがあるみたいですが、いろんな値、ちいさな値で反応があることを少しだけ期待しています。
モジュラ逆数は基本的に自然数を掛けると1になる数ですので、冪乗余の周期に関係なく、周期の最後から-1した指数のmodが出力されます。そのため、モジュラ逆数が自然数の何乗で一致するのかを調べれば、その指数が冪剰余の周期と一致します。
また、ちいさな数を生成することで、大雑把に調べることができるようになります。周期の何分の1かの数で、検出することが可能になるかもしれません。
今のところ、大きな数の素因数分解はある程度早くて正確なアルゴリズムが存在しません。今回書いたプログラムで使う計算は、冪剰余の周期を発見するためのものですが、もしうまく周期を発見することができれば、そのまま合成数の因数への分解までうまくこなすことができる可能性があるため、大きな数の素因数分解を行うスピードを大きく改善できる可能性があります。今回のプログラムに期待をしています。
実行結果のスクショはこんな感じです。周期に近いと思われる数が何度かピックアップされていることがわかると思います。
この記事が参加している募集
この記事が気に入ったらサポートをしてみませんか？