「量子コンピュータで世界中の暗号が破られる？」ショアのアルゴリズムを疑似実装したプログラムを書いた話。

目次。

目次。
はじめに。
きっかけ（海外大学院に向けて）。
参考書籍。
ショアのアルゴリズムの疑似実装。

この文章を読んで、面白い！役に立った！...と思った分だけ、投げ銭していただけると嬉しいです。

ofuse.me

【宣伝】ギターも歌も下手だけど、弾き語りをやっているので、よければ聴いてください。

www.youtube.com

はじめに。

帰省途中にパソコンが壊れて、前回の文章は「一太郎」というワープロソフトを使って書いたが、普段、Wordに慣れすぎているせいで、ものすごく使いづらい...。一太郎を使うのは、断念してこの文章はワードパットで書いている。

（途中で、数式やプログラムのコードが出てくるので、結局ワードパットも断念して、はてなブログの投稿ページに直書き）

この文章は、途中から一般向けではなくなる。悪しからず...。量子コンピュータとかショアのアルゴリズムの話をきちんとするのであれば、「量子力学」の知識は絶対に必要であるが、量子コンピュータが暗号をどうやって破るのかについて、ゆるふわな感じで全体像を理解したいのであれば、特に量子力学をはじめとした前提知識はいらないと思う。

近いうちに『「量子コンピュータで世界中の暗号が破られる？」ショアのアルゴリズムを疑似実装したプログラムを書いた話。』を前提知識がなくても分かるようにかみ砕いた文章を書きたいと思う。ただ、今は修論を書かないといけないので少し待ってください...。

⇓書いた。

blog.sun-ek2.com

量子コンピュータで破られるといわれている暗号の名前は、「RSA暗号」。

「RSA暗号...そんなもの聞いたことないから、私にはきっと関係ない話！」と思われる方がいるかもしれないが、これは大きな間違い。スマホやパソコンでTwitter, Facebook, InstagramをはじめとするSNSを見たり、ネットサーフィンしたり、メールの送受信をしているのであれば、知らず知らずのうちにRSA暗号を使っている。少なくとも、僕の文章を読んでいる時点で、既にRSA暗号を使っているのである。

RSA暗号は、インターネットで安全に情報を送るために必要不可欠な暗号。

きっかけ（海外大学院に向けて）。

僕の学部の専門は、生物工学。詳しくは、以下の文章を参照のこと。

blog.sun-ek2.com

年末に海外大学院出願作業を終えたのだが、出願したのはすべて物理学専攻の量子情報。結局、生物系の専攻には一切出願しなかった。

合格率を上げるためにやらなければならないことは、生物工学が専門の僕でも物理学専攻の量子情報でちゃんとやっていけることを向こうの大学院の先生方に納得してもらうこと。

「専門は生物工学ですが、量子情報の勉強も別でやっています」と言ってもきっと説得力がない。そのため、量子情報の分野でものすごく有名な「ショアのアルゴリズム」を疑似実装するプログラムを書いて、向こうの先生方を説得しようと思い立ったのである。

「学部の専門は生物工学やけど、量子情報の勉強をちゃんとやんじょるけん、おたくの大学院に入れていた」

...って言うよりも

「学部の専門は生物工学やけど、量子情報の勉強もちゃんとやんじょって、おたくの分野でものっそ有名なアルゴリズムのプログラム書いたけん見てんまい」

って言った方が、100倍説得力があるのである。

（高松に帰省しているので、讃岐弁で。大阪に住むようになって、関西弁にはできるだけ染まらないようにしていたのだが、そうしていたら逆に讃岐弁がきつくなったような気がする）

参考書籍。

参考にしたのは、「量子情報理論第3版」という本。量子力学の知識があれば、すらすら読める。本の中で詳しく説明されているので、情報理論の知識は持っていなくても大丈夫。第1章は、量子力学の復習みたいになっているが、こちらは量子力学の知識がない状態で読み進めるのはしんどいと思う。

posted with カエレバ

ショアのアルゴリズムの疑似実装。

使用した言語は、Python。ショアのアルゴリズムは、ある整数を2つの素数に素因数分解するアルゴリズム。この文章の目的は、ショアのアルゴリズムの説明ではないので、詳細を知りたい方は、「量子情報理論第3版」を読んでください。

こういったアルゴリズムを実際に実装してみるのは、ものすごく重要であるということがよく分かった。「本当は理解していないもの」を「理解した」と錯覚していることは多々あって、実際にアルゴリズムの実装を試みて初めて、自分に理解が足りていないということが明らかになるのである。

僕は、ショアのアルゴリズムを完全に理解したと思っていたが、いざ実装を試みると色々と壁にぶつかり、自分の理解がいかに甘かったか、痛いほど思い知った。

よく分かっていないこと（誰か教えてください...）。

ショアのアルゴリズムは、第1レジスタと第2レジスタ、2つの量子ビット列を使う。

$\displaystyle |ψ \gt =\frac{1}{\sqrt{2^{n}}} \sum^{2^{n}-1}_{a=0} |a \gt|x^{a} \, mod \, N \gt$

$\displaystyle |ψ \gt$ は、次のように書き換えることができる。

$\displaystyle |ψ^{'} \gt =\frac{1}{\sqrt{2^{n}・r}} \sum^{2^{n}-1}_{a=0} \sum^{r-1}_{s=0} e^\frac{2πisa}{r} |a \gt|Φ_s \gt$

ちなみに $\displaystyle |Φ_s \gt$ は、

$\displaystyle |Φ_s \gt =\frac{1}{\sqrt {r}} \sum^{r-1}_{k=0} e^\frac{-2πisk}{r} |x^{k} \, mod \, N \gt$

「量子情報理論第3版」のp133によると、 $\displaystyle |ψ \gt$ を $\displaystyle |ψ^{'} \gt$ に変形することによって、第2レジスタを観測することなしに最終的な結果が得られるという。

$\displaystyle |ψ \gt$ を $\displaystyle |ψ^{'} \gt$ に変形する過程は分かる。 $\displaystyle |ψ^{'} \gt$ では、第2レジスタを観測することなしになぜ答えが得られるのかも分かる。

よく分からないのが次の2点。

数式をいじって、第2レジスタの状態ケットを $\displaystyle |x^{a} \, mod \, N \gt$ という基底から $\displaystyle |Φ_s \gt$ という基底に張り直すことができるのは分かる。けれども現実問題として、「第2レジスタを観測することなしに、つまり第2レジスタと外部の系が相互作用することなしに基底ケットを張り直すことができるのか？」って、考え始めたら頭がこんがらがってきた。 $\displaystyle |ψ^{'} \gt$ を量子フーリエ変換すると合同式の周期性を捕まえることができるが、 $\displaystyle |ψ \gt$ ではできない。けれども $\displaystyle |ψ^{'} \gt$ = $\displaystyle |ψ \gt$ 。 $\displaystyle |ψ^{'} \gt$ と $\displaystyle |ψ \gt$ は、等号で結ばれるのに前者は、第2レジスタの観測なしで答えが出て、後者は第2レジスタを観測しないと答えがでない...。結構悩んだが、よく分からんかった。

第2レジスタを観測することなしに結果が得られるメリットは？別に第2レジスタを観測したからって不都合が生じるようには思えない。「量子情報理論第3版」のp133で「第2レジスタを観測しなくてもいいからすごいだろ」って言っていたのに、p137のケーススタディではがっつり第2レジスタを観測しているじゃないですか！そして、以下の記事でも第2レジスタを観測している。ちなみに僕が書いたプログラムも第2レジスタを観測している。

slpr.sakura.ne.jp

プログラミングの勉強方法。

最近は、オンラインプログラミングスクールがどんどん誕生している。オンライン相談会や体験レッスンなどが用意されているので、色々なところを見て回って、自分に合うプログラミングスクールを見つけて、そこで勉強してみては。

ソースコード（プログラムコード）。

Pythonのインタプリタと僕がimportしているパッケージ群が入っている環境なら以下のコードをコピペすれば、ちゃんと動くと思う。

ショアのアルゴリズムが威力を発揮するのは量子コンピュータの上。古典コンピュータ（現在使われているコンピュータ）にとってみれば、ただただひたすら効率の悪いアルゴリズム。5桁の数を素因数分解するのは、やめたほうがいいかもしれない。ちなみに計算律速は、量子フーリエ変換。高速フーリエ変換の知見を取り入れたら速くなりそうな気がするが、めんどくさかったので定義式をべた書き。

時々、素因数分解に失敗することがあるが、それは量子コンピュータの仕様。例えば、第1レジスタを観測して、量子状態が $\displaystyle |0...00 \gt$ に収縮すると正しい答えは得られない。ちなみにこの状態に収縮する可能性は結構高い。（詳しくは、量子フーリエ変換後に得られるグラフを参照のこと）

スマホ表示だったらソースコードがものすごく見にくいので、パソコン表示にするか、パソコンで見た方がいいかも。

import cmath
import math
import numpy as np
import random
import matplotlib.pyplot as plt
from numba import jit, prange
import time

print('')
target_N = int(input('Please enter a number for factorization: '))
start = time.time()
n_Qbits = 2*int(math.log2(target_N)+1)+1
coprime_x = 1
n_Qbits_wave_function = np.full(2**n_Qbits, cmath.exp(0j*cmath.pi))
second_register = np.ones(2**n_Qbits)
print('')


#Preparing coprime x
for i in range(3, target_N):
    if target_N % i != 0:
        coprime_x = i
        break

two_n_mod_list = np.ones(n_Qbits)
for i in range(n_Qbits):
    two_n_mod_list[i] = coprime_x**(2**i) % target_N


#Modulo operation
@jit
def mod_func():
    print('Modulo Operation...')
    for i in range(2**n_Qbits):
        tmp = np.array(np.array(list(format(i, 'b')))[::-1], dtype=int)
        tmp = tmp * two_n_mod_list[:tmp.size]
        tmp[tmp == 0] = 1
        second_register[i] = tmp.prod() % target_N
        print('\r{0}%'.format(i/2**n_Qbits*100), eあnd='')
    print('\nDone!\n')
mod_func()

second_register_measurement = second_register[random.randint(0, 2**n_Qbits-1)]
Nth_Qbit_wave_function = np.full(n_Qbits, cmath.exp(0j*cmath.pi))
survived_states = np.where(second_register == second_register_measurement)


#Quantum fourier transform
@jit
def QFT():
    print('Quantum Fourier Transform...')
    for i in range(2**n_Qbits):
        tmp_wave_function = 1
        for j in survived_states[0]:
            tmp_wave_function += cmath.exp(2j*cmath.pi*i*j/2**n_Qbits)
        n_Qbits_wave_function[i] = tmp_wave_function
        print('\r{0}%'.format(i/2**n_Qbits*100), end='')
    print('\nDone!\n')
QFT()
elapsed_time = time.time() - start
print ('Elapsed Time: {0}'.format(elapsed_time) + '[sec]')


x = np.linspace(0, 2**n_Qbits, 2**n_Qbits)
y_tmp = (n_Qbits_wave_function*n_Qbits_wave_function.conjugate()).real
y = y_tmp/y_tmp.sum()

print('\n\nPeak Values:')
for i in range(2**n_Qbits):
    if y[i] > 1/2**(n_Qbits/2):
        print(i,'     ', end='')

plt.plot(x, y, color='black')
plt.xlabel('|x>')
plt.ylabel(r'$|<x|ψ>|^2$')
plt.show()

first_register_measurement = random.random()
measured_Qbits = 0

#Cumulative probability distribution
for i in range(2**n_Qbits):
    if  first_register_measurement < y[0:i+1].sum():
        measured_Qbits = i
        break
print('\nMeasured Qbits: ', measured_Qbits)


#Continued-fraction expansion
target = measured_Qbits/2**n_Qbits
denominator_n_minus_2 = 1
numerator_n_minus_2 = 0

if target == 0:
    a = 1
    b = 1
else:
    a = int(1/target)
    b = 1/target - a

denominator_n_minus_1 = a
numerator_n_minus_1 = 1

denominator_n = 1
numerator_n = 1

while(abs(numerator_n/denominator_n-target)-1/2/2**n_Qbits > 0):
    if b == 0:
        denominator_n = denominator_n_minus_1
        break

    a = int(1/b)
    b = 1/b - a

    denominator_n = a*denominator_n_minus_1 + denominator_n_minus_2
    numerator_n = a*numerator_n_minus_1+numerator_n_minus_2

    denominator_n_minus_2 = denominator_n_minus_1
    numerator_n_minus_2 = numerator_n_minus_1

    denominator_n_minus_1 = denominator_n
    numerator_n_minus_1 = numerator_n

print('\nContinued-fraction Expansion:')
print('Denominator: ', denominator_n)
print('Numerator: ', numerator_n, '\n')


#Euclidean algorithm
candidate = [coprime_x**int(denominator_n/2) + 1, coprime_x**int(denominator_n/2) - 1]

for i in range(2):
    tmp_target_N = target_N
    if  candidate[i] == tmp_target_N:
        print('Answer: ', max(candidate[i], tmp_target_N))
        continue
    while(True):
        if candidate[i]*tmp_target_N == 0:
            break
        if candidate[i] > tmp_target_N:
            candidate[i] %= tmp_target_N
        elif candidate[i] < tmp_target_N:
            tmp_target_N %= candidate[i]
    print('Answer',i+1 ,': ', max(candidate[i], tmp_target_N))

出力例。

量子フーリエ変換後に得られるグラフ。

「量子情報理論第3版」のp137のグラフを再現できた！

横軸; 量子ビット。
縦軸; 波動関数にその複素共役をかけたもの。

f:id:sun_ek2:20200104020430p:plain

諸々と肝心な答え。

Please enter a number for factorization: 35

Modulo Operation...
99.951171875%
Done!

Quantum Fourier Transform...
99.951171875%
Done!

Elapsed Time: 5.861655235290527[sec]

Peak Values:
0 171 341 512 683 853 1024 1195 1365 1536 1707 1877
Measured Qbits: 1195

Continued-fraction Expansion:
Denominator: 12
Numerator: 7

Answer 1 : 5
Answer 2 : 7

この文章を読んで、面白い！役に立った！...と思った分だけ、投げ銭していただけると嬉しいです。

ofuse.me

にほんブログ村