動機
RSAといえば、某夏の戦争ですよね。
対象にしている読者
- 量子コンピュータでRSAが高速に解ける事は知っているけれど、具体的にどんな手順を踏んでいるか知りたい人。
- 量子コンピュータについて少ししか知らないけど、とりあえずShorのアルゴリズムを最短距離で理解したい人。具体的にはQuantum Native Dojoのドキュメント第0章 そもそも量子コンピュータとは？と第1章 量子情報の基礎あたりさえ読んでいるば理解できると思います。
- 理論だけ知りたい方は実装は飛ばしてもらっても構いません。
概要
タイトルの通りです。細かい理論については私なんかよりも深く理解した方々が詳細に説明しているものがネットに転がっていたのでそちらを引用する事にします。とりあえずこの記事ではShorのアルゴリズムの表層をさらうことを目的にしています。
以下のような順番で進めていきます。
- RSAの復習
- Shorのアルゴリズムについて
- 量子冪剰余の実装
- 量子位数発見の実装
- 量子計算部分以外の実装
- 『よろしくお願いしまぁす！（実行）』
RSAの復習
RSAが因数分解に関係あることしか理解していなかったので、改めて調べ直してみました。
まず暗号化を行うにあたって必要なのは、二つの素数$p$と$q$です。
$$
p, q \in \mathbb{N}. \quad p, q: Prime
$$
$p, q$は作成者しか知らない数字で秘密鍵と呼ばれるものです。ここから２つの公開鍵と$p,q$以外のもう一つの秘密鍵となる数値を計算します。
まず、公開鍵である$e, N \in \mathbb{N}$です。
$$
N=pq, \hspace{3pt}
gcd((p-1)(q-1), e)=1, \hspace{3pt}
e < N
$$
$N$はそのままですが、$e$は$(p-1)(q-1)$と互いに素な自然数の中から自由に選んで構いません。ここで『$(p-1)(q-1)$はどこから来たんだ？』と思う人もいると思いますが、いったん飲み込んでください。すぐに理由はわかります。
そして、秘密鍵$d \in \mathbb{N}$は以下の性質を満たすように定めます。
$$
ed \equiv 1 \hspace{3pt} (mod \hspace{3pt} (p-1)(q-1))
$$
これで秘密鍵$(p, q, d)$と公開鍵$(e, N)$が出揃いました。
暗号化
暗号化には$(e, N)$だけを使用します。暗号化したい数字を$A_{in},$暗号化した後の数字を$X$とすると
$$
A_{in}^e \equiv X \hspace{3pt} (mod \hspace{3pt} N)
$$
このように$X$を$e$で冪乗してから$N$でモジュラ演算を行います。
コードで言うなら以下のような感じです。
// Rust
let X = A.pow(e) % N
# Python
X = (A ** e) % N
復号
復号は$d$だけでできます。復号したい数字を$Y$とすると
$$
X^d \equiv A_{out} \hspace{3pt} (mod \hspace{3pt} N)
$$
暗号化と殆ど同じ手順です。コードだと以下のようになります。
// Rust
let A = X.pow(d) % N
# Python
A = (X ** d)
この時$A_{in}=A_{out}$となり復号ができます。
ではなぜ$A_{in}=A_{out}$となるのか。その理由には以下の定理が関係します
$a$と$m$が互いに素の時、以下の式を満たすというものです。
$$
a^{\phi(m)} \equiv 1 \hspace{3pt}(mod \hspace{3pt} m) \hspace{8pt}(オイラーの定理)
$$
ここで$\phi(m)$は$m$以下の数で$m$と互いに素なものの個数を表す関数です。
そして暗号化に用いる$N=pq$について考えると、$\phi(N)=(p-1)(q-1)$となることがわかると思います。
つまり
$$
a^{(p-1)(q-1)} \equiv 1 \hspace{3pt}(mod \hspace{3pt} N)
$$
はい、ここでやっと$(p-1)(q-1)$が出てきました。
前に述べた通り秘密鍵$d$は$ed \equiv 1 \hspace{3pt} (mod \hspace{3pt} (p-1)(q-1))$を満たすように選んでいました。これを変形すると$ed - 1 \equiv 0 \hspace{3pt} (mod \hspace{3pt} (p-1)(q-1))$という事になるので
$$
a^{ed-1} \equiv 1 \hspace{3pt}(mod \hspace{3pt} N)
$$
両辺に$a$をかければ
$$
a^{ed} \equiv a \hspace{3pt}(mod \hspace{3pt} N)
$$
暗号化は$X=A_{in}^e$、復号は$A_{out}=X^d$だったので$A_{out}$を$A_{in}$で表すと
$$
A_{out} = (A_{in}^e)^d = A_{in}^{ed} \equiv A_{in} \hspace{3pt}(mod \hspace{3pt} N)
$$
となり、$A_{out}=A_{in}$が示されました。つまり、一周して$A_{in}$に戻るように$d$を選んでいた訳です。
長々と述べた訳ですが重要なのは、公開鍵である$N=pq$を素因数分解して$p, q$がわかれば、$d$も分かり、RSA暗号は解けるという事です。
実装
理解のために実際にRSAによる暗号化と復号を行う構造体を作ってみました。
実際は暗号化と復号を行うのは別のユーザーであるため、このようにまとめて実装する事はないと思います。
struct RSA {
p: usize,
q: usize,
n: usize,
e: usize,
d: usize,
}
impl RSA {
fn new(p: usize, q: usize) -> Self {
let mut e = 2;
// eにはgcd(e, (p-1)(q-1))を満たす最小の数を使う事にする
for i in 2..((p - 1) * (q - 1)) {
if is_coprime(i, (p - 1) * (q - 1)) {
e = i;
break;
}
}
let d = mod_inv(e, (p - 1) * (q - 1));
return RSA {
p: p,
q: q,
n: p * q,
e: e,
d: d,
};
}
fn encrypt(&self, x: &Vec<usize>) -> Vec<usize> {
let mut y: Vec<usize> = Vec::new();
for code in x.iter() {
y.push(mod_pow(*code, self.e, self.n));
}
return y;
}
fn decrypt(&self, y: &Vec<usize>) -> Vec<usize> {
let mut x: Vec<usize> = Vec::new();
for code in y.iter() {
x.push(mod_pow(*code, self.d, self.n));
}
return x;
}
}
#[test]
fn test_rsa() {
let in_vec: Vec<usize> = vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
let rsa = RSA::new(5, 11);
let crypt = rsa.encrypt(&in_vec);
let decrypt = rsa.decrypt(&crypt);
assert_eq!(in_vec.len(), decrypt.len());
for i in 0..inp.len() {
assert_eq!(in_vec[i], decrypt[i]);
}
}
そしてテストを実行すれば復号したものが入力と同じことが確かめられると思います。
$ cargo test
Shorのアルゴリズム
$N$からその素因数である$p, q$が計算できれば良いと分かったところで素因数分解を行う量子アルゴリズムであるShorのアルゴリズムに移ります。
ただこのアルゴリズムでは量子計算によって直接素因数を求めるわけでなく、ランダムにとってきた$a$に対して$a^r \equiv 1(mod \hspace{3pt} N)$となる最小の$r$(位数と呼ばれる)を見つけることで素因数を求めます。
この$r$が偶数の時、
$$
(a^\frac{r}{2}-1)(a^\frac{r}{2}+1) \equiv 0(mod \hspace{3pt} N)
$$
より$(a^\frac{r}{2}-1)$、または$(a^\frac{r}{2}+1)$が$N$と共通な因数を持っているのでそれぞれと$N$でユークリッドの互除法を使えば因数を求めることができます。
ただ$a=N-1$などの時には$r=2$が出力され$(a^\frac{r}{2}+1)=N$となってしまい素因数を求められないこともあるのですが、そんな時は$a$を変えて再び試す事になります。
流れとしては
- $a(<N)$をランダムに選ぶ
- $gcd(a, N)=1$で無ければ、$a$は$N$の因数なのでaを出力
- $N$に対する$a$の位数を求める
- 位数が奇数だったら初めに戻る
- $s = gcd(a^\frac{r}{2}-1, N)$と$t = gcd(a^\frac{r}{2}+1, N)$を求める
- $1 < s < N$だったら$s$を出力、$1 < t < N$だったら$t$を出力、どちらでも無ければ初めに戻る
量子コンピュータを使うのは3.だけです。
ちなみにShorのアルゴリズムはnビットの整数に対して$O(n^3logn)$の速度で実行できるようです。$1$から$\sqrt{N}$まで割り算を繰り返すナイーブな方法が$O(N^{1/2})$、ポラード・ロー法(確率的に因数分解ができてとても早いアルゴリズム)が$O(N^{1/4})$なのを考えるとかなり高速ですね。
量子位数発見
量子位数発見は大まかに言うとモジュラ演算に対する量子位相推定です。
量子位相推定は上図のような量子回路で行われ、対象のユニタリ演算$U$とその固有状態$\ket{\psi}$を入力にします。この時固有状態$\ket{\psi}$を$U$に通すと$U\ket{\psi}=\lambda\ket{\psi}$となるのですが、この時の$\lambda$を近似するのが量子位相推定です。
量子位相推定については以下がかなり詳しく書かれています。
2-2. アダマールテスト2-2. アダマールテスト
2-3. 量子フーリエ変換
2-4. 位相推定アルゴリズム（入門編）
量子位数発見では以下のようなユニタリな演算$U_a$を考えます。
$$
U_{a}\ket{x}=\ket{ax \hspace{3pt} (mod \hspace{3pt}N)}
$$
$a$と$N$が互いに素であることから$0 \le x < N$の範囲において逆元を持つため、この演算はユニタリであることがわかります。
$$
\ket{u_s} = \frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}exp(-\frac{2\pi i s k}{r})\ket{a^k (mod \hspace{3pt}N)}
$$
$\ket{u_s}$が固有状態になる事は$U_a\ket{u_s}=exp(-\frac{2\pi i s}{r})\ket{u_s}$となることから示せます。
上記の式の証明
\begin{align}
U_a\ket{u_s} &= \frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}exp(-\frac{2\pi i s k}{r})\ket{a \cdot a^k (mod \hspace{3pt}N)} \\
&= \frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}exp(-\frac{2\pi i s k}{r})\ket{a^{k+1} (mod \hspace{3pt}N)} \\
&= \frac{1}{\sqrt{r}}\sum_{k'=1}^{r}exp(-\frac{2\pi i s (k' - 1)}{r})\ket{a^{k'} (mod \hspace{3pt}N)} \\
&= \frac{1}{\sqrt{r}}\sum_{k'=1}^{r}exp(-\frac{2\pi i s k'}{r} + \frac{2\pi i s}{r})\ket{a^{k'} (mod \hspace{3pt}N)} \\
&= \frac{exp( \frac{2\pi i s}{r})}{\sqrt{r}} \Bigl(exp(-\frac{2\pi s r}{r})\ket{a^r (mod \hspace{3pt}N)} + \sum_{k'=1}^{r-1}exp(-\frac{2\pi i s k'}{r})\ket{a^{k'} (mod \hspace{3pt}N)} \Bigr) \\
&= \frac{exp(\frac{2\pi i s}{r})}{\sqrt{r}} \Bigl(exp(-\frac{2\pi s \cdot 0}{r})\ket{a^0 (mod \hspace{3pt}N)} + \sum_{k'=1}^{r-1}exp(-\frac{2\pi i s k'}{r})\ket{a^{k'} (mod \hspace{3pt}N)} \Bigr) \\
&= exp(\frac{2\pi i s}{r})\frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}exp(-\frac{2\pi i s k}{r})\ket{a^k (mod \hspace{3pt}N)} \\
&= exp(\frac{2\pi i s}{r})\ket{u_s} \\
\end{align}
そのため$\ket{u_s}$を使って$U_a$の位相推定を行えば位相$\frac{s}{r}$を得ることができます。しかし今度は$\ket{u_s}$を得る必要があります。
そこで次のような性質を使う事にします。
\ket{1}=\frac{1}{\sqrt{r}}\sum_{s=0}^{r-1}\ket{u_s}
上記の式の証明
これを示すために以下のような性質を使用します
\sum_{s=0}^{r-1}exp(\frac{2\pi i s}{r})
単位円上の正多角形の頂点の重心が原点であることを想像してもらえばわかると思います。
そして、これを利用すれば
\begin{align}
\frac{1}{\sqrt{r}}\sum_{s=0}^{r-1}\ket{u_s}
&= \frac{1}{r}\sum_{s=0}^{r-1}\sum_{k=0}^{r-1}exp(-\frac{2\pi i s k}{r})\ket{a^k mod \hspace{3pt} N} \\
&= \frac{1}{r}\Bigl( \sum_{s=0}^{r-1}\ket{1} + \sum_{s=0}^{r-1}\sum_{k=1}^{r-1} exp(-\frac{2\pi i s k}{r})\ket{a^k mod \hspace{3pt} N} \Bigr) \\
&= \frac{1}{r}\Bigl( \sum_{s=0}^{r-1}\ket{1} + \sum_{k=1}^{r-1}\sum_{s=0}^{r-1} exp(-\frac{2\pi i s k}{r})\ket{a^k mod \hspace{3pt} N} \Bigr) \\
&= \frac{1}{r}\Bigl( \sum_{s=0}^{r-1}\ket{1} \Bigr) \\
&= \ket{1} \\
\end{align}
(証明終了)
つまり$\ket{1}$は$\ket{u_0},\ket{u_1}, ... , \ket{u_{r-1}}$を重ね合わせた状態であり、位相推定の入力として使用すれば、出力としてそれぞれの状態に対する$U_a$の位相$(\frac{0}{r}, \frac{1}{r}, \frac{2}{r}, ... ,\frac{r-1}{r})$が確率的に得られます。
もちろん、実際の出力は小数の二進数表現$j=(j_0 2^{-1} + j_1 2^{-2} + j_2 2^{-3} + ...) = 0.j_0j_1j_3...$の形式なので連分数展開などのアルゴリズムを使用して$r$を推測する必要があります。
推定した位相を表現する量子ビットが十分あれば、得られた$\frac{s}{r}$から$r$を予測できます。もちろんこの時$\frac{s}{r}$が既約分数でないこともあり得るので再度の試行が必要となることもあります。
まとめると$剰余乗算U_{a}$を$\ket{1}$で位相推定すれば位数$r$が得られるということです。
U_aの実装
ここからは対象となる演算$U_{a}$を構成していくことを目的としていきます。実装には自作のQitクレイトを使用します。
実装は以下の順番で行います。
- 加算、減算(add, sub)
- 剰余加算(mod_add)
- $U_a = $ 制御剰余乗算(controlled_mod_mul)
- 冪剰余
制御剰余乗算には剰余加算が必要で、剰余加算には通常の加算が必要なのでこの順番での実装が一番簡単だと思います。
加算の実装
今から加算の実装を行いますが、
ここで問題が発生します。
$\ket{a}$と$\ket{b}$を入力に$\ket{a+b}$を得る以下のような演算を作ろうとすると、$\ket{a}$に$m$個の量子ビットを使うと合計で$3m$量子ビットが必要になります。
$\ket{a}\ket{b}\ket{0...0}\rightarrow^{add}\ket{a}\ket{b}\ket{a+b}$
回路図は以下のようになります
a0----*--*--------*-----
| | |
b0----*--X--*--*--X-----
| | |
cin---|-----*--X--------
| |
cout--X-----X-----------
$N$が6量子ビットで表現できる程度(0~63)を考えていたので、加算だけで18量子ビットを使う事になります。これに加えて量子位相推定に使用するビットや乗算の際に必要となる量子ビットを考えると最低30量子ビット必要になります。
私の貧弱なパソコンでは25量子ビットあたりが限界なのでとても動かせるものではありません。
そこで$\ket{a}\ket{b} \rightarrow^{add} \ket{a}\ket{b + a}$と$b$のレジスタを上書きするようにすれば、加算の際に使用する量子ビットの数を$2m$まで減らせます。
…が、さらにここで問題が発生します。
これでもまだ量子ビットが足りませんでした。
全体を見直して考えればここで足し合わせる$a$というのは剰余乗算で使用する$a$と同じです。さらにいえば$N$についても量子位数発見の中では入力として与えられる定数なので、この二つに量子ビットは使わない方が良い事に気づきました。
ということでここからは$m$量子ビットでの加算を構成することを目標にします。
加算の実装 定数加算の実装
まず$n$量子ビットに対して、$2^m(m < n)$の加算を考えます。
例えば以下のような$n=5$量子ビットの数値に、$2^2=4$を加算しようとすると、一番桁の高い$b_4$が変化するかを考えると、$b_3, b_2$が共に1である時にビットが変化するので
$\ket{b_0}\ket{b_1}\ket{b_2}\ket{b_3}\ket{b_4}\rightarrow^{CCX(2, 3 \rightarrow 4)}\ket{b_0}\ket{b_1}\ket{b_2}\ket{b_3}\ket{b_2b_3 \oplus b_4}$
次に$b_3$は$b_2$の時に変化するので
$\ket{b_0}\ket{b_1}\ket{b_2}\ket{b_3}\ket{b_2b_3 \oplus b_4}\rightarrow^{CX(2 \rightarrow 3)}\ket{b_0}\ket{b_1}\ket{b_2}\ket{b_2 \oplus b_3}\ket{b_2b_3 \oplus b_4}$
最後に$b_2$については無条件で切り替わるので
$\ket{b_0}\ket{b_1}\ket{b_2}\ket{b_2 \oplus b_3}\ket{b_2b_3 \oplus b_4}\rightarrow^{X(2)}\ket{b_0}\ket{b_1}\ket{\bar{b_2}}\ket{b_2 \oplus b_3}\ket{b_2b_3 \oplus b_4}$
これを回路図にすると
b0---------------------------------
b1---------------------------------
b2--*--*--X------------------------
| |
b3--*--X---------------------------
|
b4--X------------------------------
加算するのが$2^{0}=1$の時は
b0--*--*--*--*--X------------------
| | | |
b1--*--*--*--X---------------------
| | |
b2--*--*--X------------------------
| |
b3--*--X---------------------------
|
b4--X------------------------------
実装の手間を減らすために任意の数の制御を持つXゲート(CNX)structを追加しました。
実際にこれをCXやCCXで再現しようとすると凄い数のゲートが必要になるようですが、少し邪道ですが今回の主題はShorのアルゴリズムなので考えない事にしました。
$add_{2^m}$の実装
オーバーフローを無視するwrapping_qadd_const_2_powerとオーバーフローした時にビットが立つようにしたoverflow_qadd_const_2_power関数を実装します。二つを分けているのは後々剰余加算の時にわかりやすくするためです。
use Qit::{core::Operator, gates::{OperatorVec, PushOps, X, CX, CCX, CNX, U}}
pub fn wrapping_qadd_const_2_power(b: &[usize], m: usize) -> U {
assert!(b.len() > 0);
assert!(b.len() > m);
let mut u_gates = OperatorVec::new();
// 複数のゲートを適用したい順番にu_gates.push_opsで追加する
for i in (m + 1)..b.len() {
let i = b.len() - i + m;
let controlls: Vec<usize> = (m..i).map(|x| b[x]).collect();
if controlls.len() == 1 {
u_gates.push_ops(CX::new(controlls[0], b[i]));
} else if controlls.len() == 2 {
u_gates.push_ops(CCX::new(controlls[0], controlls[1], b[i]));
} else {
u_gates.push_ops(CNX::new(controlls, b[i]));
}
}
u_gates.push_ops(X::new(b[m]));
// 名前をつけて一つのゲートにまとめる
return U::new(u_gates, String::from("w_qadd_const_2^n"));
}
pub fn overflow_qadd_const_2_power(b: &[usize], overflow: usize, m: usize) -> U {
assert!(b.len() > 0);
assert!(b.len() > m);
let b = &[b, &vec![overflow]].concat();
let mut u_gates: OperatorVec = Vec::new();
for i in (m + 1)..b.len() {
let i = b.len() - i + m;
let controlls: Vec<usize> = (m..i).map(|x| b[x]).collect();
if controlls.len() == 1 {
u_gates.push_ops(CX::new(controlls[0], b[i]));
} else if controlls.len() == 2 {
u_gates.push_ops(CCX::new(controlls[0], controlls[1], b[i]));
} else {
u_gates.push_ops(CNX::new(controlls, b[i]));
}
}
u_gates.push_ops(X::new(b[m]));
return U::new(u_gates, String::from("o_qadd_const_2^n"));
}
試しにこれを使ってみる事にします。
use Qit::core::{Applicable, Qubits};
// 量子ビットのインデックスを用意
let b = vec![0, 1, 2, 3];
// bに対して2^2=4を加算する回路を用意
let mut add_4 = wrapping_qadd_const_2_power(&b, 2);
// q_in = |0011⟩
let q_in = Qubits::from_num(4, 0b0011);
// 演算を適用
let q_out = add_4.apply(q_in);
q_out.print_probs();
// |0000⟩ : 0%
// |0001⟩ : 0%
// |0010⟩ : 0%
// |0011⟩ : 0%
// |0100⟩ : 0%
// |0101⟩ : 0%
// |0110⟩ : 0%
// |0111⟩ : 100%
// |1000⟩ : 0%
// |1001⟩ : 0%
// |1010⟩ : 0%
// |1011⟩ : 0%
// |1100⟩ : 0%
// |1101⟩ : 0%
// |1110⟩ : 0%
// |1111⟩ : 0%
3に4が加算されて7になっているのがわかると思います。
$add_a$の実装
$2^m$の加算ができるようになったのであとは、$a=2^0a_0+2^1a_1+2^2a_2+...+2^{n-1}a_{n-1}$として定数$a$の加算を実装します。
pub fn wrapping_qadd_const(b: &[usize], a_const: usize) -> U {
assert!(b.len() > 0);
assert!((a_const >> (b.len())) == 0);
let mut u_gates = OperatorVec::new();
for i in 0..(b.len()) {
if (a_const >> i) & 1 == 1 {
// add(2^m)を構成するゲートのベクトルをu_gatesへと結合する。
u_gates.extend(wrapping_qadd_const_2_power(b, i).gates);
}
}
return U::new(u_gates, String::from("w_qadd_const"));
}
pub fn overflow_qadd_const(b: &[usize], overflow: usize, a_const: usize) -> U {
assert!(b.len() > 0);
assert!((a_const >> (b.len())) == 0);
let mut u_gates = OperatorVec::new();
for i in 0..(b.len()) {
if (a_const >> i) & 1 == 1 {
u_gates.extend(overflow_qadd_const_2_power(b, overflow, i).gates);
}
}
return U::new(u_gates, String::from("o_qadd_const"));
}
今度はoverflow_qadd_constを使って試してみます
use Qit::core::{Applicable, Qubits};
// 量子ビットのインデックスを用意
let b = vec![0, 1, 2];
let overflow = 3;
// bに対して2^2=4を加算する回路を用意
let mut add_7 = overflow_qadd_const(&b, overflow, 7);
// q_in = |0011⟩
let q_in = Qubits::from_num(4, 0b0101);
// 演算を適用
let q_out = add_7.apply(q_in);
q_out.print_probs();
// |0000⟩ : 0%
// .
// .
// .
// |1011⟩ : 0%
// |1100⟩ : 100% <- 5 + 7 = 12 (10)になっている
// |1101⟩ : 0%
// |1110⟩ : 0%
// |1111⟩ : 0%
$sub_a$の実装
今度は逆に定数の減算を実装しますが、実はこれの実装はライブラリの機能で簡単にできます。
量子ゲートの性質として$UU^\dagger=U^\dagger U=I$というのがあるので$U_{add}$を構成するゲートの共役転置をとって逆の順番でかければ$U_{add}$の逆の演算である$U_{sub}$は作ることができます。
特に$X$ゲートの類は$X^\dagger=X$となるので、逆の順番で適用するだけで済みます。
pub fn wrapping_qsub_const(b: &[usize], a_const: usize) -> U {
let mut sub = wrapping_qadd_const(b, a_const);
// それぞれのゲートの共役転置をとって順番を反転
sub.inverse();
// 名前を変える
sub.rename(String::from("w_qsub_const"));
return sub;
}
pub fn overflow_qsub_const(b: &[usize], overflow: usize, a_const: usize) -> U {
let mut sub = overflow_qadd_const(b, overflow, a_const);
// それぞれのゲートの共役転置をとって順番を反転
sub.inverse();
sub.rename(String::from("o_qsub_const"));
return sub;
}
剰余加算の実装 定数剰余加算の実装
仕組み
次が少し難しいです。
$\ket{b}\ket{0}\rightarrow \ket{(b + a) \hspace{3pt} mod \hspace{3pt} N}\ket{0}$を行う演算を定義します。$\ket{0}$は変化しませんが演算の中で使用する補助量子ビットです。
以下のような回路です。補助ビットをoとしています
$b$は$n$ビットをまとめて表示しています
(1) (2) (3) (4) (5) (6)
b------| |--| |--[+N]--| |-----[+a]---
|+a | |-N | | |-a |
o------| |--| |---*----| |--X---------
bとoにわたってかかる演算はoverflowです。
一つずつみていきます。
(1: overflow_qadd_const(b, o, a))
$\ket{b}\ket{0}\rightarrow_{add_a} \ket{a+b}\ket{\mathsf{overflow}}$
(2: overflow_qadd_const(b, o, N))
$\ket{a+b}\ket{\mathsf{overflow}} \rightarrow_{sub_N} \ket{a+b-N}\ket{\mathsf{overflow}}$
この時
$$
\mathsf{overflow} = \begin{cases}
0 & (a+b \ge N) \
1 & (a+b < N)
\end{cases}
$$
となっていることが重要です。
(3: overflow -> wrapping_add_const(b, N))
ここで(3)のようにoverflowで制御した$N$の加算を行うと、overflowが1の時だけ加算する事になるので
$$
b = \begin{cases}
a+b-N & (a+b \ge N) \
a+b & (a+b < N)
\end{cases}
$$
入力が$a, b < N$を満たすという前提を考えると、この時点で$\ket{b}\rightarrow \ket{(a+b)mod \hspace{3pt} N}$が計算できています。
ここからは、$o$の中身を0に戻すための作業です。
(4: overflow_sub_const(b, o, a))
$(a+b \ge N)$の時と$(a+b < N)$の時とで分けて考える事にします。
(4-1: $(a+b \ge N)$の時)
この時$b$のレジスタには$\ket{a+b-N}$が入っています。
そして$o$のレジスタには$\ket{0}$が入っています。
$(a+b-N)-b=b-N$ですが、$b<N$より$b-N<0$なのでこの演算の結果はオーバーフローします。そのため
$\ket{a+b-N}\ket{0} \rightarrow \ket{b-N}\ket{1}$
(4-2: $(a+b < N)$の時)
この時$b$のレジスタには$\ket{a+b}$が入っています。
そして$o$のレジスタには$\ket{1}$が入っています。
$(a+b)-a=a$の演算はオーバーフローを起こさないので$o$のレジスタが切り替わる事はありません。そのため
$\ket{a+b}\ket{1} \rightarrow \ket{b}\ket{1}$
これで$o$の中身は必ず1になります。
(5: X(o))
X(Not)ゲートで$o$レジスタの中身0に戻します
$\ket{b-N}\ket{1} \rightarrow_{X(o)} \ket{b-N}\ket{0}(a+b \ge N)$
$\ket{b}\ket{1} \rightarrow_{X(o)} \ket{b}\ket{0}(a+b < N)$
(6: wrapping_qadd_const(b, a))
あとは先ほど引いた$a$を足しなおします。
$\ket{b-N}\ket{0} \rightarrow_{add_a} \ket{(a+b)mod \hspace{3pt} N}\ket{0}(a+b \ge N)$
$\ket{b}\ket{0} \rightarrow_{add_a} \ket{(a+b)mod \hspace{3pt} N}\ket{0}(a+b < N)$
これで$\ket{b}\ket{0}\rightarrow \ket{(b + a) \hspace{3pt} mod \hspace{3pt} N}\ket{0}$ができました。
実装
定数$a, N$とレジスタ$b$と$o$:overflowのインデックスを入力に剰余加算を行う量子回路を返す関数を実装します。
pub fn mod_add_const(b: &[usize], overflow: usize, a_const: usize, n_const: usize) -> U {
assert!(b.len() > 0);
let mut u_gates = OperatorVec::new();
// (1)[add_a] |b⟩ -> |a+b⟩
u_gates.extend(overflow_qadd_const(b, overflow, a_const).gates);
// (2)[sub_n] |a+b⟩ -> |a+b-N⟩
u_gates.extend(overflow_qsub_const(b, overflow, n_const).gates);
// [ovrflow] |0⟩ -> |0⟩ (a+b >= N), (a+b - N >= 0)
// -> |1⟩ (a+b < N), (a+b - N < 0, overflowed)
// (3)[cont_add_N] |a+b-N⟩ -> |a+b-N⟩ (a+b >= N)
// -> |a+b⟩ (a+b < N)
let add_n = wrapping_qadd_const(b, n_const);
let const_add_n = CU::new(overflow, add_n.gates, String::from("cu-add_N"));
u_gates.push_ops(const_add_n);
// (4)[sub_a] |a+b-N or a+b⟩ -> |b-N or b⟩
u_gates.extend(overflow_qsub_const(b, overflow, a_const).gates);
// (5)[unflag] |0 or 1⟩ -> |0⟩
u_gates.push_ops(X::new(overflow));
// (6)[add_a] |b-N or b⟩ -> |a+b-N or a+b⟩
u_gates.extend(wrapping_qadd_const(b, a_const).gates);
return U::new(u_gates, String::from("mod_add_const"));
}
剰余乗算の実装 定数剰余乗算の実装
仕組み
ここでは$\ket{x}\ket{zero:0}\ket{o}\ket{cont} \rightarrow_{cmm} \ket{x}\ket{ax \hspace{3pt} mod \hspace{3pt} N, or \hspace{3pt} x}\ket{o}\ket{cont}$を行う演算$cmm(controlled_mod_mul)$を実装します。
ちなみに$x$レジスタがn-bit、$zero$レジスタがn-bit、$o$レジスタは内部の剰余加算のoverflowに繋がっているので1-bit、contは1-bitです。
実はやる事は$add_{2^m}$から$add_a$を作った過程と殆ど変わりません。
具体的には入力するxを$x=2^0x_0+2^1x_1+2^2x_2+...+2^{n-1}x_{n-1}$と考えて
$x_0=1$の時$2^0\cdot a=a$を加算, $mod_add_a(x0 \rightarrow b)$
$x_1=1$の時$2^1\cdot a=2a$を加算, $mod_add_{2a}(x1 \rightarrow b)$
$x_1=1$の時$2^2\cdot a=4a$を加算, $mod_add_{4a}(x2 \rightarrow b)$
と繰り返します。この時$2^ma$を$N$でモジュラ演算しておかないと計算できないので注意してください。
これをすれば$zero$レジスタに$ax(mod \hspace{3pt} N)$が格納されるのはわかると思います。
剰余乗算はこれで出来上がります。
さらに乗算全体を$cont$ビットが立っている時に適用するようにします。
最後に$cont$ビットが立っていない時には$zero$レジスタに$x$を入力するようにすれば終わりです。
実装
pub fn cmm_const(
x: &[usize],
tar_reg: &[usize],
overflow: usize,
cont: usize,
a_const: usize,
n_const: usize,
) -> U {
assert!(tar_reg.len() == x.len());
assert!(a_const < (1 << x.len()));
assert!(n_const < (1 << x.len()));
let mut u_gates = OperatorVec::new();
// 剰余乗算
let mut mul = OperatorVec::new();
for i in 0..x.len() {
// add<2^i * a mod N>(zero)
let adder = mod_add_const(tar_reg, overflow, (a_const << i) % n_const, n_const);
mul.push_ops(CU::from_u(x[i], adder));
}
// 制御ユニタリ化
u_gates.push_ops(CU::new(cont, mul, String::from("cu-mmul")));
u_gates.push_ops(X::new(cont));
// cont=0の時、zeroレジスタにxをコピー
for i in 0..x.len() {
u_gates.push_ops(CCX::new(cont, x[i], tar_reg[i]));
}
u_gates.push_ops(X::new(cont));
return U::new(u_gates, String::from("cmm_const"));
}
剰余乗算の改良
実は先ほどの乗算回路は完全ではありません。
なぜかというと、実際は乗算ではないからです。上の剰余乗算はこのような演算をすると説明しました。
$\ket{x}\ket{zero:0}\ket{o}\ket{cont} \rightarrow_{cmm} \ket{x}\ket{ax \hspace{3pt} mod \hspace{3pt} N, or \hspace{3pt} x}\ket{o}\ket{cont}$
しかし、上の演算を$zero$が0でない時について考えてみると
$\ket{x}\ket{zero:m}\ket{o}\ket{cont} \rightarrow_{cmm} \ket{x}\ket{ax+m \hspace{3pt} mod \hspace{3pt} N, or \hspace{3pt} x}\ket{o}\ket{cont}$
そうです、実際は『$m$に$ax$を加算する』回路だったんです。
$U_a$はその性質から$U_aU_a\ket{x}=\ket{a^2x \hspace{3pt} mod \hspace{3pt} N}$となっている必要がありますが$U_{cmm(a)}$の場合$U_{cmm(a)}U_{cmm(a)}\ket{x}=\ket{2ax \hspace{3pt} mod \hspace{3pt} N}$となってしまうわけです。
そこで以下のような回路を考えます。
x-----|cmm|--| |--|icmm|--------
| | |swap| | |
zero--| a |--| |--|a^-1|--------
cmm(a)は上で用意した剰余乗算（モドキ）です。
swapは量子ビットを入れ替えるだけの演算です。CXゲートで実装できます。
icmm($a^-1$)はcmmの共役転置(inverse)です。そして$a^{-1}$は$a\cdot a^{-1} \equiv 1 (mod \hspace{3pt})$となる数値のことです。
ではこの回路で何が起きるか見ていきたいと思います。
なお、$cont$ビットは常に立っているものとします。
cmm(a)
$\ket{x}\ket{zero:0} \rightarrow_{cmm(a)} \ket{x}\ket{0 + a\cdot x}$
swap
$\ket{x}\ket{ax} \rightarrow_{swap} \ket{ax}\ket{x}$
icmm($a^-1$)
$\ket{ax}\ket{x} \rightarrow_{icmm(a^{-1})} \ket{ax}{x - a^{-1}\cdot ax} = \ket{ax}{x - x} = \ket{ax}\ket{0}$
これで剰余乗算$U_a$が構成できることが確認できました。
冪剰余の実装
このまま量子位数発見を実装していきたいところですが、その前に量子位相の図を思い出してください。
ここで使用するのは$U_a$ではなく$U_a^{2^m}$です。
そして$U_a\ket{x}=\ket{ax (mod\hspace{3pt} N)}$から
$U_a^{2^m}\ket{x}=\ket{a^{2^m}x (mod\hspace{3pt} N)}$ということがいえます。つまり$U_a^{2^m}=U_{a^{2^m}}$ということです。初めから$a^{2^m}$をかける剰余乗算を構成しておけば何度も$U_a$を使い回さずとも目的の回路を得ることができます。
ここで位相推定部分のレジスタを$y$として$y=2^0y_0+2^1y_1+2^2y_2+...+2^{n-1}y_{n-1}$と展開した時、
- $U_a^{2^0}=U_a$について考えると、$y_0$が立っている時に適用されるので、$\ket{y}\ket{x}\ket{zero} \rightarrow_{U_a} \ket{y}\ket{x\cdot a^{y_0} (mod\hspace{3pt} N)}$
- $U_a^{2^1}=U_{a^2}$について考えると、$y_1$が立っている時に適用されるので、$\ket{y}\ket{x\cdot a^{y_0} (mod\hspace{3pt} N)}\ket{zero} \rightarrow_{U_{a^2}} \ket{y}\ket{x\cdot a^{y_0 + 2y_1} (mod\hspace{3pt} N)}$
- $U_a^{2^2}=U_{a^4}$について考えると、$y_2$が立っている時に適用されるので、$\ket{y}\ket{x\cdot a^{y_0 + 2y_1} (mod\hspace{3pt} N)}\ket{zero} \rightarrow_{U_{a^4}} \ket{y}\ket{x\cdot a^{y_0 + 2y_1 + 4y_2} (mod\hspace{3pt} N)}\ket{zero}$
- ...
- ...
- ...
- $U_a^{2^{n-1}}=U_{a^{2^{n-1}}}$について考えると、$y_{n-1}$が立っている時に適用されるので、$\ket{y}\ket{x\cdot a^{y_0 + 2y_1 + 4y_2 + ... + 2^{n-2}y_{n-2}} (mod\hspace{3pt} N)}\ket{zero} \rightarrow_{U_{a^{2^{n-1}}}} \ket{y}\ket{x\cdot a^{y_0 + 2y_1 + 4y_2 + ... + 2^{n-1}y_{n-1}}(mod\hspace{3pt} N)}\ket{zero} = \ket{y}\ket{x\cdot a^{y} (mod\hspace{3pt} N)}\ket{zero}$
よって、全体で$a^y$と$x$の積をとる演算を行っていることがわかります。
さらには位相推定には$U_a$の固有状態$\ket{\psi}$の代わりに$\ket{1}$を用いるので$x=1$となり$\ket{y}\ket{1}\rightarrow \ket{y}\ket{a^y (mod \hspace{3pt})}$という冪乗の演算を行なっています。
そこでこれを一つの演算として実装します。
pub fn me_const(
y: &[usize],
x: &[usize],
zero: &[usize],
overflow: usize,
a_const: usize,
n_const: usize,
) -> U {
assert!(zero.len() == x.len());
assert!(x.len() >= 1);
assert!(is_coprime(a_const, n_const));
let mut u_gates = OperatorVec::new();
// a^x |0⟩ -> |1⟩
u_gates.push_ops(X::new(x[0]));
for i in 0..x.len() {
let y_i = y[i];
let const_a_yi = mod_power(a_const, 1 << i, n_const);
let _const_a_yi = mod_inv(const_a_yi, n_const);
//[cmm] |x⟩|0⟩ -> |x⟩|0 + x * a^2^x_n mod N⟩
u_gates.extend(cmm_const(x, zero, overflow, y_i, const_a_yi, n_const).gates);
u_gates.extend(swap(x, zero).gates);
//[icmm] |x⟩|x * a^2^x_n mod N⟩ -> |x - x * a^2^x_n * a^(-2^x_n)⟩|x * a^2^x_n mod N⟩
// -> |0⟩|x * a^2^x_n mod N⟩
let mut icmm = cmm_const(x, zero, overflow, y_i, _const_a_yi, n_const);
icmm.inverse();
u_gates.extend(icmm.gates);
}
return U::new(u_gates, String::from("me_const"));
}
量子位数発見の実装
量子フーリエ変換
量子位相推定には$U_a$の他に、量子フーリエ変換が必要です。
2-3. 量子フーリエ変換
理屈はこちらで解説してあるので省略しますが以下のような変換を行う回路です。
\ket{j}=\frac{1}{\sqrt{2^n}}\sum_{k=0}^{2^n-1}exp(\frac{2\pi k j}{2^n})\ket{k}
pub fn qft(x: &[usize]) -> U {
let n = x.len();
let mut u_gates = OperatorVec::new();
let (a, b): (Vec<usize>, Vec<usize>) = (
(0..(n / 2)).map(|i| x[i]).collect::<Vec<usize>>(),
(0..(n / 2)).map(|i| x[n - i - 1]).collect::<Vec<usize>>(),
);
let sw = swap(&a, &b);
u_gates.extend(sw.gates);
for i in 0..n {
// hadamard
u_gates.push_ops(H::new(x[i]));
for j in (i + 1)..n {
let angle = (-((j + 1 - i) as f64)).exp2();
let r = R::new(x[i], 2.0 * PI * angle);
u_gates.push_ops(CU::new(
x[j],
vec![Box::new(r)],
format!("r_+2^-{}", j + 1 - i),
));
}
}
return U::new(u_gates, String::from("qft"));
}
実際に使ってみます。
use Qit::circuits::qft;
use Qit::core::{Applicable, Qubits};
// 量子ビットのインデックスを用意
let b = vec![0, 1, 2];
// 量子フーリエ変換を用意する
let qft = qft(&b);
let q_in = Qubits::from_num(3, 0b011);
// 演算を適用
let q_out = qft.apply(q_in);
q_out.print_cmps();
// |000⟩ : +0.354 +0.000i
// |001⟩ : -0.250 +0.250i
// |010⟩ : -0.000 -0.354i
// |011⟩ : +0.250 +0.250i
// |100⟩ : -0.354 +0.000i
// |101⟩ : +0.250 -0.250i
// |110⟩ : +0.000 +0.354i
// |111⟩ : -0.250 -0.250i
量子位相推定に使っているのは逆量子フーリエ変換ですがフーリエ変換を構成した後にライブラリのinverseメソッドを使用して作ります。
pub fn inv_qft(x: &[usize]) -> U {
let mut iqft = qft(x);
iqft.inverse();
return iqft;
}
quantum_order_discover(量子位数発見)
いよいよ本題の関数ですが、入力は$(a: usize, N: usize)(gcd(a, N) = 1)$、出力は$a^r\equiv 1(mod \hspace{3pt}N)$を満たす$r$とします。
使用する量子ビット数は固定にします。
fn quantum_order_discover(a_const: usize, n_const: usize) -> usize {
use Qit::circuits::{inv_qft, me_const};
use Qit::core::pop_from_probs;
use Qit::core::{Applicable, Operator, Qubits};
use Qit::gates::{H, U, OperatorVec, PushOps};
assert!(n_const < 1 << 6);
assert!(a_const < 1 << 6);
// ここに(s/r)が入る
let x = vec![0, 1, 2, 3, 4, 5, 6, 7];
// U_aの計算用のレジスタ
let a_x = vec![8, 9, 10, 11, 12, 13];
let zero = vec![14, 15, 16, 17, 18, 19];
let overflow = 20;
let mut hadamards = OperatorVec::new();
for i in 0..x.len() {
hadamards.push_ops(H::new(x[i]));
}
let hadamards = U::new(hadamards, String::from("hadamrds"));
let mut pe_gates = OperatorVec::new();
pe_gates.push_ops(me_const(
&x, &a_x, &zero, overflow, a_const, n_const,
));
let me = U::new(pe_gates, String::from("mod_exponential"));
let iqft = inv_qft(&x);
let q_in = Qubits::from_num(21, 0);
// 冪剰余演算を適用
let q_out = me.apply(hadamards.apply(q_in));
// 位相推定
let q_out = iqft.apply(q_out);
// q_outの中でのxの確率の分布をVec<f64>の形で得る。
let probs = q_out._measure(&x);
loop {
// 確率分布から値を取り出す
// 実際の量子コンピュータだと1shotにあたる
let s_r = pop_from_probs(&probs, x.len());
println!("pop s_r:{:>08b}({}/{})", s_r, s_r, 1 << x.len());
if s_r == 0 {
continue;
}
// s/rの値からrの候補になる複数の値を小さい順にベクトルで出力する
let r_cands = bit_float2denoms(s_r, x.len());
println!("{r_cands:#?}");
for r in r_cands {
// rが大きすぎる時にはもう一度確率分布から取り出す
if r >= (1 << x.len()) {
break;
}
// a^r == 1 なら出力
if mod_pow(a_const, r, n_const) == 1 {
return r;
}
}
}
}
// 実数値を連分数展開する
fn cont_frac(t: f64, m: usize, th: f64) -> Vec<usize> {
let mut nums: Vec<usize> = Vec::new();
let mut x = t;
let mut b = x.floor();
for _ in 0..m {
nums.push(b as usize);
x -= b;
if x == 0.0 {
return nums;
}
if x < th {
return nums;
}
x = 1.0 / x;
b = x.floor();
}
return nums;
}
// 実数値を近似して分母の候補を出力する
fn approximate_frac_denoms(t: f64, m: usize, th: f64) -> Vec<usize> {
let mut denoms: Vec<usize> = Vec::new();
let cont = cont_frac(t, m, th);
let (mut p0, mut q0) = (cont[0], 1);
let (mut p1, mut q1) = (p0 * cont[1] + 1, cont[1]);
denoms.push(q1);
for i in 2..cont.len() {
let p2 = cont[i] * p1 + p0;
let q2 = cont[i] * q1 + q0;
denoms.push(q2);
(q0, p0) = (q1, p1);
(q1, p1) = (q2, p2);
}
return denoms;
}
// 与えられた二進数を少数にして分母の候補を出力する
fn bit_float2denoms(bit: usize, size: usize) -> Vec<usize> {
// 0.01101010 -> n / m -> [m]
let t = (bit as f64) / ((1 << size) as f64);
return approximate_frac_denoms(t, 5, 1.0 / 2.0_f64.powi(size as i32));
}
これで量子計算を使用する部分は終わりです。
量子計算部分以外の実装
後はShorのアルゴリズムの説明の通りに実装します。
Shorのアルゴリズムの入力は自然数$N$で出力はその因数です。素因数とは限らないのでご注意ください。
fn shors_algorithm(n: usize) -> usize {
//! n = pq, return p or q
let mut rng = rand::thread_rng();
loop {
// 2 ~ N の間で乱数をとってくる
let a = rng.gen_range(2..n);
println!("use a={}", a);
// とってきた数値aとnが共通の因数を持つならこの時点で終わり
if !is_coprime(a, n) {
println!("a and n is not coprime");
// 今回は量子位数発見をしたいのでそのまま続ける
continue;
}
// 量子位数発見
println!("running quantum_order_discover");
let r = quantum_order_discover(a, n);
// 位数が偶数では無いときは初めから
println!("discovered {}!", r);
if r % 2 != 0 {
println!("{} is not even", r);
continue;
}
// 位数が偶数のときは(a^(r/2) + 1)か(a^(r/2) - 1)がnと共通な因数を持つ
let (a_r0, a_r1) = (mod_pow(a, r / 2, n) - 1, mod_pow(a, r / 2, n) + 1);
println!("(a^(r/2) - 1) = {}(mod N)", a_r0);
println!("(a^(r/2) + 1) = {}(mod N)", a_r1);
let gcd_ar0_n = gcd(a_r0, n);
println!("gcd(a^(r/2) - 1, N) = {}", gcd_ar0_n);
if gcd_ar0_n != 1 && gcd_ar0_n != n {
println!("return gcd(a^(r/2) - 1, N) = {}", gcd_ar0_n);
return gcd_ar0_n;
}
let gcd_ar1_n = gcd(a_r1, n);
println!("gcd(a^(r/2) + 1, N) = {}", gcd_ar0_n);
if gcd_ar1_n != 1 && gcd_ar1_n != n {
println!("return gcd(a^(r/2) + 1, N) = {}", gcd_ar1_n);
return gcd_ar1_n;
}
println!("continue");
}
}
実行
$N$を適当に決めて実行してみます。
今回は$N=57$で実行します。素数にしか見えませんが本当に因数分解できるのでしょうか？
fn main() {
shors_algorithm(57);
}
では……よろしくおねがいしまぁす！！
use a=53
running quantum_order_discover
pop s_r:00101011(43/256)
[
5,
6,
125,
256,
]
pop s_r:11110010(242/256)
[
1,
18,
55,
73,
]
discovered 18!
(a^(r/2) - 1) = 55(mod N)
(a^(r/2) + 1) = 57(mod N)
gcd(a^(r/2) - 1, N) = 1
gcd(a^(r/2) + 1, N) = 1
continue
use a=7
running quantum_order_discover
pop s_r:00000000(0/256)
pop s_r:10101011(171/256)
[
1,
3,
253,
256,
]
discovered 3!
3 is not even
use a=40
running quantum_order_discover
pop s_r:11010101(213/256)
[
1,
5,
6,
125,
]
pop s_r:01000111(71/256)
[
3,
4,
7,
11,
]
pop s_r:00001110(14/256)
[
18,
55,
73,
128,
]
discovered 18!
(a^(r/2) - 1) = 36(mod N)
(a^(r/2) + 1) = 38(mod N)
gcd(a^(r/2) - 1, N) = 3
return gcd(a^(r/2) - 1, N) = 3
なかなか苦戦していたようですが、$a=40$の時に$r=18$を見つけることができ、最終的には3という因数を見つけることができました。
(s/r)を正確に推定するためには$N$が$n$ビットで表現できる時、$2n+1$ビットくらいを位相の表現部分に使用するのが望ましいらしいです。今回でいえば追加で7ビットくらいが必要ということですね……厳しい。
おわりに
私自身は量子コンピュータに興味はありながらも実際に何かを実装してみた事は無かったので、勉強しながらこのアルゴリズムを実装してみる経験はかなり新鮮なものでした。
あと長年の夢（暗号を解く）を擬似的に叶えられたのも結構嬉しかったです。興味のある方は是非実装してみてください。