フーリエ分析の知識をベースに、素因数分解のための量子ユニット(？)シミュレーションを自作してみる。
前回の記事で、素因数分解のためのショアのアルゴリズムを説明して、それとは少し違った形の因数を調べる方法を考えていました。
今回は、いくらか勉強して試してみて、量子ユニット(？)っぽい動きをする何かを古典的な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
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
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
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)
return bekiyolist
def teibai(lists,tei):
#jouken: ishhu ga len(lists)ni onaji
list1=lists
lens1=len(list1)
lens2=len(list1[0])
teib=tei
tei1=teib%lens1
tei2=teib%lens2
list2=np.array([[list1[(j*tei1)%lens1][(i*tei2)%lens2] for i in range(lens2)]for j in range(lens1)],dtype=complex)
return list2
def teibai2(lists,tei1,tei2):
#jouken ishhu ga len(lists) ni onaji
list1=lists
#print(list1,'lists')
lens1=len(list1)
lens2=len(list1[0])
list2=np.array([[list1[(j*tei1)%lens1][(i*tei2)%lens2] for i in range(lens2)]for j in range(lens1)],dtype=complex)
return list2
def toridasi(lists,div1,div2,objn,count):
#jouken ishhu ga len(lists) ni onaji
list1=np.array(lists)
#print(list1,'lists')
lens1=div1
lens2=div2
#print(lens1,lens2)
#teib=tei
objnum=objn
tbunkatu=count
kaiten=pow(np.e,2j*np.pi*(1/div1/div2))
list2=np.array([[0 for i in range(lens2)]for j in range(lens1)],dtype=complex)
list2a=np.array([[pow(np.e,2j*np.pi*(j*i/div1)) for i in range(lens1)]for j in range(lens1)],dtype=complex)
list2b=np.array([[pow(np.e,2j*np.pi*(i*j/div2)) for i in range(lens2)]for j in range(lens2)],dtype=complex)
lit=np.array([0 for i in range(tbunkatu)])
#print(list2b)
list2=list2a@list1@list2b/lens1/lens2
for i in range(lens1):
for j in range(lens2):
if abs(list2[i][j])>0.0001:
a=cmath.phase(list2[i][j])*objnum/np.pi/2
if a<-0.5:
a+=objnum
if np.round(a)<tbunkatu:
lit[int(np.round(a))]+=1
#
return lit
def modinv(a,n):
#a>n
inv=0
for i in range(1,n):
if a*i%n==1:
inv=i
break
return inv
tei=3
objnum=4087
i2sisu=4
kenshou=106
num1=59
num2=71
check1=59
check2=71
for k in solist:
dft1=np.array([[pow(np.e,-2j*np.pi*((i+0)/num2+(j+0)/num1))for i in range(num2)]for j in range(num1)],dtype=complex)
dftk1=np.array([[pow(np.e,-2j*np.pi*((j+0)*i/num1))for i in range(num1)]for j in range(num1)],dtype=complex)
dftk2=np.array([[pow(np.e,-2j*np.pi*((i+0)*j/num2))for i in range(num2)]for j in range(num2)],dtype=complex)
dft1mat=np.array(dftk1**(-1))
dft2mat=np.array(dftk2**(-1))
ang=pow(np.e,2j*np.pi/objnum)
dft1=dft1
dft=dft1
kenshoudft=teibai(kenshoudft,tei**kenshou)
kenshoudft=1.0*kenshoudft*ang**(tei**kenshou%objnum)
#print(ugojo(3*71-1,objnum%k),1,'modinv')
dft=kenshoudft
if i2sisu>0:
dftmat=np.array(dft)
dftmat=dft2mat@(dftmat.T)@dft1mat/num1/num2
dft0=np.array(dftk1)@(dftmat**(tei)).T@np.array(dftk2)
dft=teibai(dft0,tei)+dft
print('len',k,'hanni',kenshou,'~',kenshou+2**(i2sisu)-1)
for i in range(1,i2sisu):
print(i)
dftmat=np.array(dft)
dftmat=dft2mat@(dftmat.T)@dft1mat/num1/num2
dftmat=dftmat**((tei**2**i)%objnum)
dft0=np.array(dftk1)@(dftmat).T@np.array(dftk2)
dft=teibai(dft0,(tei**(2**i))%objnum)+dft
print(toridasi(dft,num1,num2,objnum,20))
dfti=np.array([[0 for i in range(check2)]for j in range(check1)],dtype=complex)
dfti=dft1mat@np.array(dft)@dft2mat/num1/num2
if abs(dfti[1][1])>0.9:
count1+=abs(dfti[1][1])
for i in range(check1):
for j in range(check2):
if abs(dfti[i][j])>0.1:#pow(np.e,-2)*0.1:
a=-cmath.phase(dfti[i,j])
if a<0:
a+=2*np.pi
if i+j-1 or[num1,num2]:
print ('check',i,j)
print(i,j,abs(dfti[i][j]),math.log(abs(dfti[i][j]),1.2),objnum-a/np.pi*objnum/2,ugojo(objnum,num1*j+i))
if i==1 and j==1 and False:
print(abs(dfti[1][1]))
print('insuu score(1.0de positive)',count1/len(solist)/2**(i2sisu)*2)冪剰余の関数はオマケ。実行結果はこんな感じです。
まずこれは何かということから説明しますと、メインで動かしているユニットが59×71個の２次元のユニットの集合で、それぞれに複素数の数値を入れることができます。それぞれのユニットに、整数の周波数の波を表す数値を入れることにより、２次元の波形データを作り出すことになります。
波形としてプロットするのは
$${f(x,y)=e^{-2πi(x/59+y/71)}}$$(x、yはユニットの中での番地)を基本単位とする周波数が自然数の波の波形です。この波形を最初に入力して、適当に倍数倍したり周波数を増やしたりして、最後に逆フーリエ変換に相当する操作、
$${F(x,y)=f(x,y)e^{2πi(x/59+y/71)}}$$を算出すると、それぞれどの周波数の波がどれだけの大きさで、位相がどのくらいであるのかを全て知ることができます。
基本的な特徴をリストアップしてみます。
ユニットの数以上の周波数を表現することはできないため、基本はモジュロ演算(掛けた時の余りが出力される演算)です。
基本の周波数(周波数1の波とする)と同じ波同士を掛け合わせると、周波数2の波を作ることができます。
周波数nの波と周波数mの波を掛け合わせると、周波数n＋mの波ができます。
周波数nの波と周波数mの波を足し合わせると、周波数nの波と周波数mの波の合成波になります。このふたつの波は、逆フーリエ変換に相当する操作にて分離することができます。
1つ飛びでユニットが順に並ぶように並び替えると全ての波の周波数が2倍されます。同じように、n個飛ばしでユニットを取って順に並べると、全ての波の周波数がn+1倍になります。(恐らくこれはスワップと呼ばれる操作です。)
以上の操作の中で全て、波の振幅と位相(波の角度が実数軸からどれだけ離れているか)は保存されます。
for i in range(1,i2sisu):
print(i)
dftmat=np.array(dft)
dftmat=dft2mat@(dftmat.T)@dft1mat/num1/num2
dftmat=dftmat**((tei**2**i)%objnum)
dft0=np.array(dftk1)@(dftmat).T@np.array(dftk2)
dft=teibai(dft0,(tei**(2**i))%objnum)+dft
以上の部分のように、逆フーリエ変換した上で、全ての波の位相をn倍すると言った操作も、可能ではあります。
周波数をn倍する操作に、モジュラ逆数$${m^{-1}}$$を用いることによって、どうやら割り算のモジュロ演算を達成することができるようです。
$${e^{-2πi(1/n)}}$$を掛けると、全ての波の位相を角度1/n×360度分までずらすことができます。
これらの性質は恐らく量子計算とある程度似た動作をするみたいですので、素因数分解をするアルゴリズムを理解する一助になるかもしれません。
ということで、このプログラムでは、位相と周波数の両方で3を底とする累乗のモジュロ計算をしています。指数関数的に数値を増やしており、この条件では3の106乗から122乗までのmod59、71の周波数と、mod4087の位相を全て、見ることができます。toridasiで逆フーリエ変換に相当する操作をして、位相の小さな数値だけを取り出して計上しています。
いろいろと工夫する余地がありそうですね。