量子アルゴリズムの基本：Shorのアルゴリズムの確認（その２）

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

はじめに

前回の（その１）で、アルゴリズムの説明をしたので、自作の量子計算シミュレータqlazyで動作確認をします。処理の流れ全体を改めて整理しておきます。

• [STEP.0] 素数判定
• [STEP.1] 因数が自明な場合、計算してアルゴリズム終了
• [STEP.2] aをランダムに選ぶ
• [STEP.3] a,Nの最大公約数gcd(a,N)を計算する
• [STEP.4] a,Nに関する位数を求める(量子位数発見)
  • [STEP.4-1] 第1レジスタを準備
  • [STEP.4-2] 第2レジスタを準備
  • [STEP.4-3] べき剰余の計算
  • [STEP.4-4] 第1レジスタを逆量子フーリエ変換
  • [STEP.4-5] 第1レジスタを測定して位数推定
• [STEP.5] 位数から因数を求める

実装

ということで、早速Pythonで実装したコード全体を示します。

【2021.9.5追記】qlazy最新版でのソースコードはここに置いてあります。

コード全体

shor.py

import sys
import random
import math
import numpy as np
from qlazypy import QState

#
# custom gates
#

def swap(self,q0,q1):

    self.cx(q0,q1).cx(q1,q0).cx(q0,q1)
    return self

def hadamard(self,id=None):

    for i in range(len(id)):
        self.h(id[i])
    return self

def iqft(self,id=None):

    dim = len(id)

    for i in range(dim):
        phase = -1.0/2**i
        for j in range(i):
            self.cp(id[j],id[i],phase=phase)
            phase *= 2.0
        self.h(id[i])
            
    i = 0
    while i < dim-1-i:
        self.swap(id[i], id[dim-1-i])
        i += 1
        
    return self

def mod_exp(self,a,N,id_up,id_dn):

    # make unitary matrix: |j>|k>->|j>|a^j k mod N>
    
    bitnum_up = len(id_up)
    bitnum_dn = len(id_dn)
    
    matdim_up = 2**bitnum_up
    matdim_dn = 2**bitnum_dn
    matdim = matdim_up * matdim_dn

    mat = np.array([[0]*matdim]*matdim)
    for j in range(matdim_up):
        for k in range(matdim_dn):
            col = matdim_dn*j + k
            if k < N:
                row = matdim_dn*j + ((a**j)*k)%N
            else:
                row = col
            mat[row][col] = 1.0

    # apply unitary matrix(mat) to the quantum state

    self.apply(mat)
    
    return self

#
# functions
#

def create_register(num):

    return [0]*num

def init_register(*args):

    idx = 0
    for i in range(len(args)):
        for j in range(len(args[i])):
            args[i][j] = idx
            idx += 1
    return idx

def approx_frac(number):

    iter = 10
    eps = 0.0001

    # if input number is integer
    if number == int(number):
        return (int(number),1)
    
    # expand continued fraction
    x = number
    confra = []
    for i in range(iter):
        a = int(x)
        confra.append(a)
        x -= a
        if x < eps:
            break
        x = 1/x

    # approximate fraction (get denominator,numerator)
    denom_0, numer_0 = confra[0], 1
    denom_1, numer_1 = confra[1] * denom_0 + 1, confra[1]
    denom, numer = denom_1, numer_1
    for i in range(2,len(confra)):
        denom, numer = confra[i]*denom_1+denom_0, confra[i]*numer_1+numer_0
        denom_0, denom_1 = denom_1, denom
        numer_0, numer_1 = numer_1, numer

    return (denom, numer)
    
def discover_order(a,N):

    bitnum_dn = len(str(bin(N-1)))-2
    # bitnum_up = 2*bitnum_dn + 1
    bitnum_up = 4
    id_up = create_register(bitnum_up)
    id_dn = create_register(bitnum_dn)
    qubit_num = init_register(id_up,id_dn)
    
    qs = QState(qubit_num)
    qs.hadamard(id_up)
    qs.x(id_dn[0])
    qs.mod_exp(a,N,id_up,id_dn)

    qs.iqft(id_up)
    md = qs.m(id=id_up,shots=100)

    order = 1
    for i in range(len(md.freq_list)):
        if md.freq_list[i] != 0:
            s_r = i/(2**bitnum_up)  # s/r
            (s,r) = approx_frac(s_r)
            if r != 0 and r < N and (a**r)%N == 1 and r%2 == 0:
                order = r
                break
    qs.free()
    
    return r  # return 1 if order discovery failure

def discover_order2(a,N):
    
    for i in range(1,N):
        x = (a**i)%N
        if x == 1:
            break
    return i

def is_prime(number):

    if number == 1:
        return False
    for i in range(2,int(math.sqrt(number)+1)):
        if number%i == 0:
            return False
    return True
        
def base_exponent(number):

    size = len(str(bin(number)))-2
    for exp in reversed(range(1,size)):
        bas = int(math.pow(number,1/exp))
        if bas**exp == number:
            break
    return bas,exp
    
def shor(N):

    # if N <= 1
    if N <= 1:
        print("N must be bigger than 1")
        sys.exit()
        
    # if N is prime
    if is_prime(N) == True:
        print("N is prime")
        sys.exit()
    
    # if N is even, then factor->2
    if N%2 == 0:
        print("N is even")
        return 2
    
    # if N = bas^exp, then factor->bas
    bas,exp = base_exponent(N)
    if exp != 1:
        print("N is {0:}^{1:}, factor -> {0:}".format(bas,exp))
        return bas

    # otherwise
    for _ in range(10):
        a = random.randint(2,N-1)
        print(">> try a = ", a)
        fac = math.gcd(a,N)
        print("gcd({0:},{1:}) -> {2:}".format(a,N,fac))
        if fac != 1:
            print("factor -> ",fac)
            break
        else:
            print("order discovery start")
            r = discover_order(a,N)
            print("order -> ", r)

            if r%2 != 0:
                print("order discovery failure")
                continue

            fac = math.gcd(int(a**(r/2))-1,N)
            if fac != 1 and fac != N:
                print("factor -> ",fac)
                break
            fac = math.gcd(int(a**(r/2))+1,N)
            if fac != 1 and fac != N:
                print("factor -> ",fac)
                break
            print("factor estimation failure")

    return fac

if __name__ == '__main__':

    QState.swap = swap
    QState.hadamard = hadamard
    QState.iqft = iqft
    QState.mod_exp = mod_exp

    print("== input number ==")
    args = sys.argv
    N = int(args[1])
    print("N = ",N)

    print("== Shor's algorithm start ==")
    fac_0 = shor(N)
    fac_1 = N // fac_0

    print("== factorization ==")
    print("N = {0:} = {1:} * {2:}".format(N,fac_0,fac_1))

main処理部の、

fac_0 = shor(N)

が、Shorのアルゴリズムの本体を呼び出しているところです。この関数shorの中身が本質なので、その中身を見ながら、順に説明していきます。

大まかな流れ

[STEP.0] 素数判定

まず、

# if N <= 1
if N <= 1:
    print("N must be bigger than 1")
    sys.exit()

で、与えられた整数Nが1以下の場合は対象外なのではじいた上で、素数判定です。

# if N is prime
if is_prime(N) == True:
    print("N is prime")
    sys.exit()

is_prime関数がTrueを返したら、「素数だよ」と言って終了します。素数判定の関数is_primeは、

def is_prime(number):

    if number == 1:
        return False
    for i in range(2,int(math.sqrt(number)+1)):
        if number%i == 0:
            return False
    return True

です。2から順番に割り算していくのですが、$\sqrt{N}+1$まで試せば十分です。というわけで、おバカなアルゴリズムで実装しました。

[STEP.1] 因数が自明な場合、計算してアルゴリズム終了

# if N is even, then factor->2
if N%2 == 0:
    print("N is even")
    return 2

Nが偶数なら、「因数は2だよ」と言って、2をリターンします。

# if N = bas^exp, then factor->bas
bas,exp = base_exponent(N)
if exp != 1:
    print("N is {0:}^{1:}, factor -> {0:}".format(bas,exp))
    return bas

base_exponent関数で、Nをべき乗の形($N=a^b$)で書いたときの整数$a,b$を返します。$b$が1以外であれば、得られた$a$をリターンします。関数base_exponentの中身は、

def base_exponent(number):

    size = len(str(bin(number)))-2
    for exp in reversed(range(1,size)):
        bas = int(math.pow(number,1/exp))
        if bas**exp == number:
            break
    return bas,exp

です。意外とあっさり書けました。指数$b$を変化させながら、実際に$N^{1/b}$を計算して、それが整数になっているかどうかを判定していく感じです。整数になったら終了して、そのときの$a$(変数bas)と$b$(変数exp)をリターンします。最終的に指数$b$が1だとすると、$N$はべき乗の形には書けないということなので、その判断もこの関数でできます。

[STEP.2] aをランダムに選ぶ

# otherwise
for _ in range(10):
    a = random.randint(2,N-1)

forループで10回実行します(10は適当です)。何回かランダムにaを選んで、位数がわかったら、ループを抜けます。

[STEP.3] a,Nの最大公約数gcd(a,N)を計算する

fac = math.gcd(a,N)
print("gcd({0:},{1:}) -> {2:}".format(a,N,fac))
if fac != 1:
    print("factor -> ",fac)
    break
else:
    xxx

最大公約数は、mathモジュールを使って計算できます。facが1以外だった場合は、「それが因数だよ」と言って、ループを抜けます。

[STEP.4] a,Nに関する位数を求める(量子位数発見)

それ以外の場合(facが1の場合)、else節を実行します。というわけで、Shorのアルゴリズムの肝である「量子位数発見」です。

print("order discovery start")
r = discover_order(a,N)
print("order -> ", r)

が、大事なところなので、一旦後回しにします。rが得られたとして、先に進みます。

if r%2 != 0:
    print("order discovery failure")
    continue

rが奇数になることがあるので、その場合、continueして、別のaを選んで再度同じことをやります。偶数のrが得られたら、次に進みます。

[STEP.5] 位数から因数を求める

以下、見ての通りです。

fac = math.gcd(int(a**(r/2))-1,N)
if fac != 1 and fac != N:
    print("factor -> ",fac)
    break
fac = math.gcd(int(a**(r/2))+1,N)
if fac != 1 and fac != N:
    print("factor -> ",fac)
    break
print("factor estimation failure")

gcd(int(a**(r/2))-1,N)かgcd(int(a**(r/2))+1,N)のどちらか一方が1でもなくNでもなければ、それを因数だとして、ループを抜けます。そして、得られた因数をリターンします。

量子位数発見

それでは、先程説明を省略した「量子位数発見」について見ていきます。関数discover_orderの中身を見ながら、順に説明します。

[STEP.4-1] 第1レジスタを準備 / [STEP.4-2] 第2レジスタを準備

まず、必要なビット数を設定します。bitnum_upが第1レジスタのビット数、bitnum_dnが第2レジスタのビット数を表しています。

bitnum_dn = len(str(bin(N-1)))-2
# bitnum_up = 2*bitnum_dn + 1
bitnum_up = 4

第2レジスタは、Nのビット長に等しくします。第1レジスタの方は、Nのビット長の2倍よりも大きくすべし、ということらしいのですが、今回量子ビット節約のため4ビットにしました(後で動作確認するくらいの整数であれば、4ビットで大丈夫でした)。

id_up = create_register(bitnum_up)
id_dn = create_register(bitnum_dn)
qubit_num = init_register(id_up,id_dn)

各レジスタのビット数に応じて、レジスタ(実態は量子番号のリスト)を設定します。関数create_registerは、以下のようなもので、とりあえず中身が0のリストをリターンします。

def create_register(num):

    return [0]*num

中身を入れるのは、以下のinit_register関数です。

def init_register(*args):

    idx = 0
    for i in range(len(args)):
        for j in range(len(args[i])):
            args[i][j] = idx
            idx += 1
    return idx

ご覧の通り、引数として与えるリストの数はいくつでも良いです。指定されたリストの順番に番号を順序よく入れていきます。そして全体の量子ビット数をリターンします。

qs = QState(qubit_num)
qs.hadamard(id_up)
qs.x(id_dn[0])

第1レジスタ全体にアダマールをかけ、第2レジスタの最初のビットを反転させます。これで、はじめに用意しておく状態の準備ができました。

[STEP.4-3] べき剰余の計算

qs.mod_exp(a,N,id_up,id_dn)

ここで、問題の「べき剰余」が登場します。本来であれば、以前の記事で「べき剰余」を実現した量子回路で実行できれば良かったのですが、とても小さな整数の計算にも多くの量子ビットを使ってしまうので、今回は諦めました。代わりにちょっと卑怯なのですが、べき剰余を実行した結果の量子状態を人為的に用意しました。

前回の記事で説明したべき剰余後の状態、

\frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} \ket{k} \ket{a^k \space mod \space N} \tag{1}

をつくるために、まず以下のようなユニタリ変換を考えます。

U_{a,N}: \ket{k} \ket{t} \rightarrow \ket{k} \ket{a^{k} t \space mod \space N} \tag{2}

このユニタリ変換を状態$\ket{k} \ket{1}$に適用すれば、式(1)の状態が得られます。なので、今回、このユニタリ変換を直接、行列の形にして、それを量子状態に適用するようにしました。

それを、以下のように実装しました。

def mod_exp(self,a,N,id_up,id_dn):

    # make unitary matrix: |k>|t>->|k>|a^k t mod N>
    
    bitnum_up = len(id_up)
    bitnum_dn = len(id_dn)
    
    matdim_up = 2**bitnum_up
    matdim_dn = 2**bitnum_dn
    matdim = matdim_up * matdim_dn

    mat = np.array([[0]*matdim]*matdim)
    for j in range(matdim_up):
        for k in range(matdim_dn):
            col = matdim_dn*j + k
            if k < N:
                row = matdim_dn*j + ((a**j)*k)%N
            else:
                row = col
            mat[row][col] = 1.0

    # apply unitary matrix(mat) to the quantum state

    self.apply(mat)
    
    return self

変数matにその行列を表すnumpy配列が格納されます。それをqlazyのQStateクラスのapplyメソッドによって演算するようにしています(applyメソッドはv0.0.18で追加しました)。量子状態の次元は常に2のべき乗なので、行列の次元も2のべき乗にしたいのですが、Nの値によってはそうならないので、若干工夫をしています。一旦、2のべき乗の次元を用意して、上のユニタリ変換で定義されない数については恒等置換するようにしました。これで、べき剰余を実行した結果の状態が得られることになります。

discover_order関数の中身の説明に戻ります。

[STEP.4-4] 第1レジスタを逆量子フーリエ変換

次は、

qs.iqft(id_up)

ということで、逆量子フーリエ変換です。関数の中身は見ての通りということで説明省略します。

[STEP.4-5] 第1レジスタを測定して位数推定

次に、

md = qs.m(id=id_up,shots=100)

order = 1
for i in range(len(md.freq_list)):
    if md.freq_list[i] != 0:
        s_r = i/(2**bitnum_up)  # s/r
        (s,r) = approx_frac(s_r)
        if r != 0 and r < N and (a**r)%N == 1 and r%2 == 0:
            order = r
            break

ということで、第1レジスタを測定した結果に基づき、位数を推定します。連分数展開を使うのですが、それをapprox_frac関数で実行しています(アルゴリズムはいろいろな記事で紹介されていると思うので、中身の説明は省略します)。得られた分母が0でなくて、本当に位数になっているのかをチェックして、さらにここで偶数かどうかもチェックするようにしました。OKとなったら、その位数をリターンします。これで位数発見の処理は完了です。

動作確認

プログラムは、因数分解したい整数をコマンドラインから指定するようになっています。例えば、

python shor.py 123

のようにやれば計算してくれます。では、いくつかの整数でやってみます。上で説明したようにaはランダムに選ばれるので、場合によっては(aがたまたま因数だった場合とか)、量子位数発見せずに結果を返してしまいます。それでは面白くないので、何度か実行して量子位数発見を通った場合の結果のみを以下では示すようにします。

15の因数分解

== input number ==
N =  15
== Shor's algorithm start ==
>> try a =  4
gcd(4,15) -> 1
order discovery start
order ->  2
factor ->  3
== factorization ==
N = 15 = 3 * 5

となりました。aとして4が選ばれて、15と共通の約数がない(互いに素)ので、位数発見します。結果は2でした。これに基づき因数を推定した結果、3でした。ということで、15=3*5のように因数分解できました。

57の因数分解

==input number ==
N =  57
== Shor's algorithm start ==
>> try a =  35
gcd(35,57) -> 1
order discovery start
order ->  16
factor ->  3
== factorization ==
N = 57 = 3 * 19

となり、こちらも正しく因数分解できました。

119の因数分解

最後に119ではどうでしょうか。

== input number ==
N =  119
== Shor's algorithm start ==
>> try a =  108
gcd(108,119) -> 1
order discovery start
order ->  16
factor ->  17
== factorization ==
N = 119 = 17 * 7

ちゃんとできました。というわけで、Nのビット長が7ビットまで確認できました。2から127までのすべての数で一応やってみましたが、正しく動作することが確認できました。

おまけ

今回、量子位数発見を量子計算シミュレータで確認しましたが、この程度の整数であれば、古典コンピュータのおバカなアルゴリズムでも十分できます(もちろん)。例えば、

def discover_order2(a,N):
    
    for i in range(1,N):
        x = (a**i)%N
        if x == 1:
            break
    return i

こんな関数で行けます。上のプログラムのdiscover_orderをdiscover_order2に置き換えれば、量子計算をやらないバージョンのShorのアルゴリズムがお試しいただけます。意味ないですが、Shorシミュレータのデバッグには役立ちました。

おわりに

Shorのアルゴリズムは全体通すととてもややこしいですが、一歩ずつ勉強しながら、何とか理解できたと思います。やはりプログラムを実装しながらというのは、時間はかかるのですが、確実に理解するための１つの良い方法だと改めて感じました。理解が曖昧なところでは実装の手は止まりますし、間違って理解しているところでは確実にバグが出ますので。

今回の本質は、discover_order関数に集約されています。というわけで、このdiscover_orderのところをみなさまお好きな環境(シミュレータ)で置き換えていただけると、勉強になりますし、またShorのアルゴリズム全体を通してお楽しみ(謎)いただけると思いますので、チャレンジしてみてはいかがでしょうか。でっかい整数もこんな感じで因数分解できたよ、ということがあれば、後で教えて下さい！

以上