ショアのアルゴリズム(量子計算)の発想を参考に、特殊関数のフーリエ変換から位数を求め、素因数を探す方法の検討

素因数分解のアルゴリズムを検討します。

import numpy as np
import cmath
import math
import random
def bekiyolist(bunkaiN,tei,bekimax):
    #bunkaiN=497
    #tei=4
    #beki=12
    beki=bekimax
    bekiyolist=[0]*(beki+1)
    blist=[]
    j=0
    result=1
    bekiyolist[0]=1
    while j<beki:

        result=((result%bunkaiN)*(tei%bunkaiN))%bunkaiN
        bekiyolist[1+j]=result
        j+=1
        #print(tei,result)
    return bekiyolist
objnum=4087#gouseisuu
sample=3000
output=1000
chstart=0


print('.')
beky=bekiyolist(objnum,3,sample)
print('.')
print(beky[0:30])

beky2=[0]*output
beky0=0
for i in range(sample+1):
    beky0+=beky[i]*pow(np.e,-2*np.pi*i*0*1.0j/objnum)
    for j in range(chstart,chstart+len(beky2)):
        beky2[j-chstart]+=beky[i]*pow(np.e,-2*np.pi*beky[i]*j*1.0j/objnum)
beky2=[abs(beky2[i]) for i in range(len(beky2))]
beky3=[0]*2000
for i in range(len(beky)):
    if beky[i]==1:
        print (i)
print(['{:.1f}'.format(abs(k)/1000)for k in beky2])
print(beky2[0])
for i in range(1,len(beky2)):
    if abs(beky2[i])>abs(beky0)/4:
        print(math.gcd(objnum,i+chstart),abs(beky2[i]),i+chstart)

print(objnum)
mod=7**110-1
print(math.gcd((1+mod),4087),math.gcd(mod-1,4087))

プログラムは上記の通りです(python)
量子計算で使うショアのアルゴリズムとは、量子コンピュータを使う素因数分解のアルゴリズムで、多項式時間で素因数分解問題を解決できる、現状唯一？のアルゴリズムとして知られています。
アルゴリズムの肝となる量子計算は、量子離散フーリエ変換です。量子計算にはいろいろな複雑な計算方法があり、未だ詳しい解説本が世の中に普及していないところも多いのですが、とりあえずここでは、量子計算ではフーリエ変換を非常に早く計算できることを利用して、計算を効率化しているということだけ覚えて頂ければ大丈夫です。
また、量子計算に依存しないコンピュータでも、巨大な整数の最大公約数を求める計算は、ユークリッドの互除法を用いて、非常に早く計算できることも、このアルゴリズムの早さを保証する大事な要素のようです。

ショアのアルゴリズムの簡単な流れです。

1. 素因数分解したい数Nと互いに素な自然数a(1<a<N)を１つ選ぶ。

2. f(x)=$${a^x(mod N)}$$を0からN-1までの自然数xについて計算して、f(x)の位数rを逆量子フーリエ変換で求める。

3. rが偶数でなければ1.に戻る。rが偶数であれば$${a^{r/2}-1,a^{r/2}+1}$$とNの最大公約数をそれぞれ計算する。

4. 求めた最大公約数が1またはNでなければ、その最大公約数をNの因数の一つとして出力する。

今回の記事では、この方法をアレンジして、f(x)=$${a^x(mod N)}$$を元に特殊な関数を定義して、この関数フーリエ変換の結果を検証することにより、大きな数の素因数分解を早めの速度で達成する方法を考えてみようと思います。

変更点

f(x)はf(x)=$${a^x(mod N)}$$を使います。
立式上の大きな変更点はフーリエ変換の所です。

$${y(t)=\int f(x)e^{-\frac{2πxt}{N}}dx → y(t)=\int f(x)e^{-\frac{2πf(x)t}{N}}dx}$$

スクリプトでのここの変更点は、次のように表現されます。

beky2[j-chstart]+=beky[i]*pow(np.e,-2*np.pi*i*(j+1)*1.0j/objnum)#フーリエ変換

↓

beky2[j-chstart]+=beky[i]*pow(np.e,-2*np.pi*beky[i]*j*1.0j/objnum)#今回使う変換

この計算式で求めた数値について、y(t)の絶対値の大きいxを求めて、求めたtについてNとの最大公約数を計算して、1でもNでもないものをNの因数のひとつとして出力します。上のショアのアルゴリズムの流れの3以降については、求められる位数が全く違う数になるため、行いません。

検証のための出力結果

このスクリーンショットでは、4087の素因数分解について、出力しています。
スクリーンショットの下は、y(t)>y(0)/4を満たすy(t)の絶対値と共に、tとNの最大公約数を出力するようにしていますが、出力される数は高頻度で、4087の素因数である、61と67を倍数に持つ数であることがわかります。

また、このスクリーンショットでは、xの範囲を0から120までの範囲にして計算しています。この計算でも、61と67の倍数をよく抽出できているようです。
一般に、この計算では、xの計算範囲を0から求めたい素因数よりも少しだけ大きい自然数の範囲まで計算すれば(範囲を連続する自然数の集合で定義する場合)、十分な計算結果が得られるようです。

素因数分解の計算量(ランダウ)は現在のところショアのアルゴリズムを除く最速で準指数計算量で求めることができて、最速で求めたい数Nの2進法での桁数nについて、

$${O(\frac{4}{\sqrt[3]{3}}n^{2/3}e^{n/3})}$$

の準指数時間で求められるみたいです。( https://dojo.qulacs.org/ja/latest/notebooks/2.4_phase_estimation_beginner.html#コラム：素因数分解と位相推定（やや難） )
ショアのアルゴリズムでは、標準的な計算量で、$${O(n^3logn)}$$で求められるようです。記事で扱ったアルゴリズムはショアのアルゴリズムよりは遅く、準指数時間くらいかな？と思います。計算を最適化したら、どれだけ早く計算ができるか、少々気になるところとなっています。

ショアのアルゴリズムは、一応1からNまでのxを求めなければいけないみたいですが、求めたい因数の大きさまで求めれば大丈夫であれば或いは…と願いたくなるような計算方法でした。

合成数565073について、xの範囲を0から1800、tの範囲を5000〜10000、ピックアップするtをy(t)>y(0)/15で計算することとして計算開始。20.8秒程で計算が完了し、797の倍数が3個抽出できた。

以上です。