ショアのアルゴリズム

次に、整数因数分解問題に注目し、位相推定を用いて量子コンピュータ上でこの問題を効率的に解く方法を見ていきましょう。 導出されるアルゴリズムが 整数因数分解のためのショアのアルゴリズム です。 ショアは自身のアルゴリズムを特に位相推定の言葉で説明していませんでしたが、その仕組みを説明するうえで自然かつ直感的なアプローチです。

まず、位数発見問題 と呼ばれる中間問題について説明し、位相推定がこの問題の解を与えることを確認します。 次に、位数発見問題の効率的な解が整数因数分解問題の効率的な解をもたらすことを見ていきます。 （ある問題の解が別の問題の解を与えるとき、後者の問題は前者に 帰着する と言います。ここでは整数因数分解を位数発見に帰着しています。） ショアのアルゴリズムのこの第二部分は量子計算をまったく使いません。完全に古典的です。 量子計算が必要なのは位数発見を解くためだけです。

位数発見問題

初等整数論の基礎

位数発見問題と、それを位相推定で解く方法を説明するために、まずいくつかの基本的な整数論の概念を確認し、便利な記法を導入しておきましょう。

はじめに、任意の正整数 $N$ に対して、集合 $\mathbb{Z}_N$ を次のように定義します。

$$\mathbb{Z}_N = \{0,1,\ldots,N-1\}$$

例えば、 $\mathbb{Z}_1 = \{0\},\;$ $\mathbb{Z}_2 = \{0,1\},\;$ $\mathbb{Z}_3 = \{0,1,2\},\;$ などとなります。

これらは数の集合ですが、単なる集合以上のものとして考えることができます。 特に、$\mathbb{Z}_N$ 上の 演算（加算や乗算など）を考えることができ、答えを常に $N$ で割った余り（すなわち $N$ を法とした剰余）として取ることで、これらの演算を施しても常にこの集合の中に留まります。 $N$ を法とした加算と乗算という二つの演算によって、$\mathbb{Z}_N$ は代数において基本的に重要なオブジェクトである 環 になります。

例えば、$3$ と $5$ は $\mathbb{Z}_7$ の元であり、これらを掛け合わせると $3\cdot 5 = 15$ となり、これを $7$ で割ると余りが $1$ になります。 これを次のように表すことがあります。

$$3 \cdot 5 \equiv 1 \; (\textrm{mod } 7)$$

ただし、$\mathbb{Z}_7$ 内で計算していることが明確な場合は、記法を簡潔にするために単に $3 \cdot 5 = 1$ と書くこともできます。

例として、$\mathbb{Z}_6$ の加算表と乗算表を示します。

$$\begin{array}{c|cccccc} + & 0 & 1 & 2 & 3 & 4 & 5 \\\hline 0 & 0 & 1 & 2 & 3 & 4 & 5 \\ 1 & 1 & 2 & 3 & 4 & 5 & 0 \\ 2 & 2 & 3 & 4 & 5 & 0 & 1 \\ 3 & 3 & 4 & 5 & 0 & 1 & 2 \\ 4 & 4 & 5 & 0 & 1 & 2 & 3 \\ 5 & 5 & 0 & 1 & 2 & 3 & 4 \\ \end{array} \qquad \begin{array}{c|cccccc} \cdot & 0 & 1 & 2 & 3 & 4 & 5 \\\hline 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 2 & 3 & 4 & 5 \\ 2 & 0 & 2 & 4 & 0 & 2 & 4 \\ 3 & 0 & 3 & 0 & 3 & 0 & 3 \\ 4 & 0 & 4 & 2 & 0 & 4 & 2 \\ 5 & 0 & 5 & 4 & 3 & 2 & 1 \\ \end{array}$$

$\mathbb{Z}_N$ の $N$ 個の元のうち、$\gcd(a,N) = 1$ を満たす元 $a\in\mathbb{Z}_N$ は特別です。 これらの元からなる集合は、しばしば次のようにアスタリスクを用いて表されます。

$$\mathbb{Z}_N^{\ast} = \{a\in \mathbb{Z}_N : \gcd(a,N) = 1\}$$

乗算という演算に注目すると、集合 $\mathbb{Z}_N^{\ast}$ は 群、とりわけ アーベル群 を形成します。これも代数において重要なオブジェクトの一つです。 これらの集合（および有限群一般）についての基本的な事実として、任意の元 $a\in\mathbb{Z}_N^{\ast}$ を選んで $a$ を繰り返し自身に掛けていくと、いずれ必ず $1$ が得られます。

最初の例として、$N=6$ を取りましょう。 $\gcd(5,6) = 1$ なので $5\in\mathbb{Z}_6^{\ast}$ であり、$5$ を自身に掛けると $1$ が得られます。これは上の表からも確認できます。

$$5^2 = 1 \quad \text{($\mathbb{Z}_6$ 内での計算)}$$

第二の例として、$N = 21$ を取りましょう。 $0$ から $20$ までの数の中で、$21$ との GCD が $1$ になるものは次の通りです。

$$\mathbb{Z}_{21}^{\ast} = \{1,2,4,5,8,10,11,13,16,17,19,20\}$$

これらの元それぞれについて、その数を正整数乗すると $1$ が得られます。 その最小の乗数を示すと次の通りです。

$$\begin{array}{ccc} 1^{1} = 1 \quad & 8^{2} = 1 \quad & 16^{3} = 1 \\[1mm] 2^{6} = 1 \quad & 10^{6} = 1 \quad & 17^{6} = 1 \\[1mm] 4^{3} = 1 \quad & 11^{6} = 1 \quad & 19^{6} = 1 \\[1mm] 5^{6} = 1 \quad & 13^{2} = 1 \quad & 20^{2} = 1 \end{array}$$

これらの等式はすべて $\mathbb{Z}_{21}$ 内での計算であることは省略しています。記述を煩雑にしないよう、以後のレッスン全体でもこれを暗黙の前提とします。

問題の定式化と位相推定との関係

これで位数発見問題を定式化できます。

位数発見

入力: $\gcd(N,a) = 1$ を満たす正整数 $N$ と $a$
出力: $a^r \equiv 1$ $(\textrm{mod } N)$ を満たす最小の正整数 $r$

あるいは、先ほど導入した記法では、$a \in \mathbb{Z}_N^{\ast}$ が与えられ、$a^r = 1$ を満たす最小の正整数 $r$ を求めることになります。 この数 $r$ を、$a$ の $N$ を法とした 位数 と呼びます。

位数発見問題を位相推定と結び付けるために、古典的な状態が $\mathbb{Z}_N$ に対応するシステムに対して、固定した元 $a\in\mathbb{Z}_N^{\ast}$ で乗算する演算を考えましょう。

$$M_a \vert x\rangle = \vert ax \rangle \qquad \text{(各 $x\in\mathbb{Z}_N$ に対して)}$$

明確にしておくと、$\mathbb{Z}_N$ 内での乗算を行っているので、右辺のケットの内側の積が $N$ を法としていることは暗黙の前提です。

例えば、$N = 15$、$a=2$ とすると、$M_2$ の標準基底 $\{\vert 0\rangle,\ldots,\vert 14\rangle\}$ への作用は次の通りです。

$$\begin{array}{ccc} M_{2} \vert 0 \rangle = \vert 0\rangle \quad & M_{2} \vert 5 \rangle = \vert 10\rangle \quad & M_{2} \vert 10 \rangle = \vert 5\rangle \\[1mm] M_{2} \vert 1 \rangle = \vert 2\rangle \quad & M_{2} \vert 6 \rangle = \vert 12\rangle \quad & M_{2} \vert 11 \rangle = \vert 7\rangle \\[1mm] M_{2} \vert 2 \rangle = \vert 4\rangle \quad & M_{2} \vert 7 \rangle = \vert 14\rangle \quad & M_{2} \vert 12 \rangle = \vert 9\rangle \\[1mm] M_{2} \vert 3 \rangle = \vert 6\rangle \quad & M_{2} \vert 8 \rangle = \vert 1\rangle \quad & M_{2} \vert 13 \rangle = \vert 11\rangle \\[1mm] M_{2} \vert 4 \rangle = \vert 8\rangle \quad & M_{2} \vert 9 \rangle = \vert 3\rangle \quad & M_{2} \vert 14 \rangle = \vert 13\rangle \end{array}$$

$\gcd(a,N)=1$ のとき、この演算はユニタリです。標準基底 $\{\vert 0\rangle,\ldots,\vert N-1\rangle\}$ の元を並べ替えるため、行列として見ると置換行列となります。 定義から明らかにこの演算は決定論的であり、逆演算が存在することを確認する簡単な方法は、$a$ の $N$ を法とした位数 $r$ を考え、$M_a$ の逆演算が $M_a^{r-1}$ であることに気づくことです。

$$M_a^{r-1} M_a = M_a^r = M_{a^r} = M_1 = \mathbb{I}$$

$r$（求めようとしている量）に関する知識を必要としない逆演算の考え方もあります。 $a\in\mathbb{Z}_N^{\ast}$ のすべての元に対して、$ab=1$ を満たす一意の元 $b\in\mathbb{Z}_N^{\ast}$ が常に存在します。 この元 $b$ を $a^{-1}$ と表記し、効率的に計算できます。 ユークリッドの GCD アルゴリズムの拡張により、$\operatorname{lg}(N)$ の 2 乗オーダーのコストで計算できます。 したがって、

$$M_{a^{-1}} M_a = M_{a^{-1}a} = M_1 = \mathbb{I}.$$

以上より、演算 $M_a$ は決定論的かつ可逆です。 これは $M_a$ が置換行列で表されることを意味し、したがってユニタリです。

では、$a\in\mathbb{Z}_N^{\ast}$ を仮定した上で、演算 $M_a$ の固有ベクトルと固有値について考えましょう。 今述べたように、この仮定から $M_a$ がユニタリであることがわかります。

$M_a$ の固有値は $N$ 個あり（同じ固有値が重複することもあります）、対応する固有ベクトルの選び方には一般に自由度があります。ただし、すべての可能性を考慮する必要はありません。 まず簡単に始めて、$M_a$ の固有ベクトルを一つだけ特定しましょう。

$$\vert \psi_0 \rangle = \frac{\vert 1 \rangle + \vert a \rangle + \cdots + \vert a^{r-1} \rangle}{\sqrt{r}}$$

$r$ は $a$ の $N$ を法とした位数であり、以後のレッスン全体でもこれを前提とします。 この固有ベクトルに対応する固有値は $1$ です。$a$ を掛けてもこのベクトルは変化しないためです。

$$M_a \vert \psi_0 \rangle = \frac{\vert a \rangle + \cdots + \vert a^{r-1} \rangle + \vert a^r \rangle}{\sqrt{r}} = \frac{\vert a \rangle + \cdots + \vert a^{r-1} \rangle + \vert 1 \rangle}{\sqrt{r}} = \vert \psi_0 \rangle$$

これは $a^r = 1$ のため、各標準基底状態 $\vert a^k \rangle$ が $k\leq r-1$ のとき $\vert a^{k+1} \rangle$ にシフトされ、$\vert a^{r-1} \rangle$ は $\vert 1\rangle$ に戻るからです。 直感的に言えば、$\vert \psi_0 \rangle$ をゆっくりかき混ぜているようなものですが、すでに完全にかき混ざっているため何も変わりません。

次に $M_a$ の別の固有ベクトルの例を示します。 こちらは位数発見と位相推定の文脈でより興味深いものです。

$$\vert \psi_1 \rangle = \frac{\vert 1 \rangle + \omega_r^{-1} \vert a \rangle + \cdots + \omega_r^{-(r-1)}\vert a^{r-1} \rangle}{\sqrt{r}}$$

あるいは、このベクトルを和の記法で次のように書けます。

$$\vert \psi_1 \rangle = \frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \omega_r^{-k} \vert a^k \rangle$$

ここで複素数 $\omega_r = e^{2\pi i/r}$ が自然に現れています。これは $N$ を法とした $a$ による乗算の仕組みによるものです。 今度の対応する固有値は $\omega_r$ です。 これを確認するために、まず次のように計算します。

$$M_a \vert \psi_1 \rangle = \frac{1}{\sqrt{r}}\sum_{k = 0}^{r-1} \omega_r^{-k} M_a\vert a^k \rangle = \frac{1}{\sqrt{r}}\sum_{k = 0}^{r-1} \omega_r^{-k} \vert a^{k+1} \rangle = \frac{1}{\sqrt{r}}\sum_{k = 1}^{r} \omega_r^{-(k - 1)} \vert a^{k} \rangle = \frac{1}{\sqrt{r}}\omega_r \sum_{k = 1}^{r} \omega_r^{-k} \vert a^{k} \rangle$$

次に、$\omega_r^{-r} = 1 = \omega_r^0$ かつ $\vert a^r \rangle = \vert 1\rangle = \vert a^0\rangle$ であることから、

$$\frac{1}{\sqrt{r}}\sum_{k = 1}^{r} \omega_r^{-k} \vert a^{k} \rangle = \frac{1}{\sqrt{r}}\sum_{k = 0}^{r-1} \omega_r^{-k} \vert a^k \rangle = \vert\psi_1\rangle,$$

よって $M_a \vert\psi_1\rangle = \omega_r \vert\psi_1\rangle$ が成り立ちます。

同様の議論により、$M_a$ の追加の固有ベクトル・固有値ペアを特定できます。 任意の $j\in\{0,\ldots,r-1\}$ に対して、

$$\vert \psi_j \rangle = \frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \omega_r^{-jk} \vert a^k \rangle$$

は $M_a$ の固有ベクトルであり、対応する固有値は $\omega_r^j$ です。

$$M_a \vert \psi_j \rangle = \omega_r^j \vert \psi_j \rangle$$

$M_a$ には他にも固有ベクトルが存在しますが、それらは気にしなくて構いません。ここでは特定した固有ベクトル $\vert\psi_0\rangle,\ldots,\vert\psi_{r-1}\rangle$ にのみ注目します。

位相推定による位数発見

与えられた $a\in\mathbb{Z}_N^{\ast}$ に対する位数発見問題を解くには、演算 $M_a$ に位相推定手順を適用します。

そのためには、$M_a$ だけでなく $M_a^2$、$M_a^4$、$M_a^8$ など、位相推定手順から十分精度の高い推定値を得るために必要なだけ量子回路で効率的に実装する必要があります。 ここではその方法を説明し、必要な精度の度合いについては後ほど詳しく確認します。

まず $M_a$ 単体の演算から始めましょう。 量子回路モデルで作業しているため、当然ながら $0$ から $N-1$ までの数を2進法で符号化します。 符号化すべき最大の数は $N-1$ なので、必要なビット数は

$$n = \operatorname{lg}(N-1) = \lfloor \log(N-1) \rfloor + 1$$

となります。

例えば、$N = 21$ の場合は $n = \operatorname{lg}(N-1) = 5$ です。 $\mathbb{Z}_{21}$ の元を長さ $5$ の2進数文字列として符号化すると次のようになります。

$$\begin{gathered} 0 \mapsto 00000\\[1mm] 1 \mapsto 00001\\[1mm] \vdots\\[1mm] 20 \mapsto 10100 \end{gathered}$$

そして、$M_a$ を $n$ 量子ビット演算として定義した正確な定義を示します。

$$M_a \vert x\rangle = \begin{cases} \vert ax \; (\textrm{mod}\;N)\rangle & 0\leq x < N\\[1mm] \vert x\rangle & N\leq x < 2^n \end{cases}$$

ここで重要なのは、$M_a$ がどのように機能するかは $\vert 0\rangle,\ldots,\vert N-1\rangle$ に対してのみ関心がありますが、残りの $2^n - N$ 個の標準基底状態に対しても動作を指定する必要があり、それがユニタリ演算となるようにしなければならないことです。 残りの標準基底状態に対して $M_a$ が何もしないと定義することでこれを実現できます。

前のレッスンで説明した整数乗算・除算のアルゴリズムと、それらのリバーシブルかつガベージフリーな実装の方法論を用いることで、任意の $a\in\mathbb{Z}_N^{\ast}$ に対して $O(n^2)$ のコストで $M_a$ を実行する量子回路を構築できます。 その実現方法の一つを以下に示します。

1. 演算

   $$\vert x \rangle \vert y \rangle \mapsto \vert x \rangle \vert y \oplus f_a(x)\rangle$$

   を実行する回路を構築します。ここで

   $$f_a(x) = \begin{cases} ax \; (\textrm{mod}\;N) & 0\leq x < N\\[1mm] x & N\leq x < 2^n \end{cases}$$

   前のレッスンで説明した方法を使います。これにより $O(n^2)$ のサイズの回路が得られます。

2. $n$ 個のスワップゲートを使って2つの $n$ 量子ビットシステムを入れ替えます（量子ビットを個別にスワップします）。

3. 最初のステップと同様に、演算

   $$\vert x \rangle \vert y \rangle \mapsto \vert x \rangle \bigl\vert y \oplus f_{a^{-1}}(x)\bigr\rangle$$

   のための回路を構築します。ここで $a^{-1}$ は $\mathbb{Z}_N^{\ast}$ における $a$ の逆元です。

下側の $n$ 量子ビットを初期化して3つのステップを合成すると、次の変換が得られます。

$$\vert x \rangle \vert 0^n \rangle \stackrel{\text{ステップ1}}{\mapsto} \vert x \rangle \vert f_a(x)\rangle \stackrel{\text{ステップ2}}{\mapsto} \vert f_a(x)\rangle \vert x \rangle \stackrel{\text{ステップ3}}{\mapsto} \vert f_a(x)\rangle \bigl\vert x \oplus f_{a^{-1}}(f_a(x)) \bigr\rangle = \vert f_a(x)\rangle\vert 0^n \rangle$$

この方法はワークスペース量子ビットを必要としますが、最終的に初期化された状態に戻されるため、位相推定にこれらの回路を使用できます。 得られる回路の総コストは $O(n^2)$ です。

$M_a^2$、$M_a^4$、$M_a^8$ などを実行するには、全く同じ方法を使えます。ただし、$a$ を $\mathbb{Z}_N^{\ast}$ の元として $a^2$、$a^4$、$a^8$ などに置き換えます。 すなわち、任意の冪 $k$ に対して、$M_a$ の回路を $k$ 回繰り返すのではなく、$b = a^k \in \mathbb{Z}_N^{\ast}$ を計算して $M_b$ の回路を使うことで $M_a^k$ の回路を作成できます。

冪 $a^k \in \mathbb{Z}_N$ の計算は、前のレッスンで言及した モジュラー冪乗 問題です。 この計算は前のレッスンで言及したモジュラー冪乗アルゴリズム（計算的数論では 冪アルゴリズム とも呼ばれる）を用いて 古典的に 実行できます。 実際に必要なのは $a$ の 2の冪乗 のみ、具体的には $a^2, a^4, \ldots a^{2^{m-1}} \in \mathbb{Z}_N^{\ast}$ であり、これらは $m-1$ 回の反復的な平方演算によって得られます。 各平方演算は $O(n^2)$ のサイズのブール回路で実行できます。

本質的に、ここで行っているのは、$M_a$ を最大 $2^{m-1}$ 回繰り返す問題を効率的な古典計算にオフロードすることです。 そしてこれが可能であることは幸運なことです！ 位相推定問題における量子回路の任意の選択では、これが可能でない場合が多く、その場合の位相推定のコストは制御量子ビット数 $m$ に対して 指数的 に増大します。

都合のよい固有ベクトルが与えられた場合の解法

位相推定を使って位数発見問題をどのように解けるかを理解するために、まず操作 $M_a$ に固有ベクトル $\vert\psi_1\rangle$ を用いて位相推定の手順を実行する場合を考えてみましょう。 実はこの固有ベクトルを手に入れることは容易ではないため、これが話の終わりではありませんが、ここから始めるのが理解の助けになります。

固有ベクトル $\vert \psi_1\rangle$ に対応する $M_a$ の固有値は

$$\omega_r = e^{2\pi i \frac{1}{r}}$$

です。つまり、$\theta = 1/r$ として $\omega_r = e^{2\pi i \theta}$ となります。 したがって、固有ベクトル $\vert\psi_1\rangle$ を使って $M_a$ に対して位相推定を実行すると、$1/r$ の近似値が得られます。 その逆数を計算することで、近似が十分に精度よければ $r$ を求めることができます。

さらに詳しく説明すると、$m$ 個の制御量子ビットを使って位相推定の手順を実行した場合、得られるのは $y\in\{0,\ldots,2^m-1\}$ という数です。 次に $y/2^m$ を $\theta$ の推定値とします。今の場合、$\theta = 1/r$ です。 この近似から $r$ を求めるには、近似の逆数を取って最近傍の整数に丸めるのが自然な方法です。

$$\left\lfloor \frac{2^m}{y} + \frac{1}{2} \right\rfloor$$

例えば、$r = 6$ の場合に固有ベクトル $\vert\psi_1\rangle$ を用いて $m = 5$ 個の制御ビットで $M_a$ に対する位相推定を行うとしましょう。 $1/r = 1/6$ に対する最良の $5$ ビット近似は $5/32$ であり、位相推定から結果 $y=5$ が得られる確率はかなり高く（この場合約 $68\%$）なります。 このとき

$$\frac{2^m}{y} = \frac{32}{5} = 6.4$$

となり、最近傍の整数に丸めると $6$ が得られ、これは正しい答えです。

一方、精度が不十分な場合は正しい答えが得られないこともあります。 例えば、位相推定で $m = 4$ 個の制御量子ビットを使った場合、$1/r = 1/6$ に対する最良の $4$ ビット近似である $3/16$ が得られる可能性があります。 その逆数を取ると

$$\frac{2^m}{y} = \frac{16}{3} = 5.333 \cdots$$

となり、最近傍の整数に丸めると誤った答えの $5$ が得られてしまいます。

では、正しい答えを得るためにはどれほどの精度が必要でしょうか？ 位数 $r$ は整数であることがわかっており、直感的には $1/r$ を近傍の候補（$1/(r+1)$ や $1/(r-1)$ を含む）と区別できるだけの精度が必要です。 $1/r$ に最も近い、注意すべき数は $1/(r+1)$ であり、これら2つの数の距離は

$$\frac{1}{r} - \frac{1}{r+1} = \frac{1}{r(r+1)}$$

です。

したがって、$1/r$ を $1/(r+1)$ と間違えないようにするには、最良近似 $y/2^m$ が $1/r$ に対して $1/(r+1)$ よりも近くなるような精度を保証すれば十分です。 精度が十分で

$$\left\vert \frac{y}{2^m} - \frac{1}{r} \right\vert < \frac{1}{2 r (r+1)}$$

すなわち誤差が $1/r$ と $1/(r+1)$ の距離の半分未満であれば、$y/2^m$ は $1/(r+1)$ や $1/(r-1)$ を含む他のどの候補よりも $1/r$ に近くなります。

これを次のように確認できます。

$$\frac{y}{2^m} = \frac{1}{r} + \varepsilon$$

ここで $\varepsilon$ は

$$\vert\varepsilon\vert < \frac{1}{2 r (r+1)}$$

を満たすとします。逆数を取ると

$$\frac{2^m}{y} = \frac{1}{\frac{1}{r} + \varepsilon} = \frac{r}{1+\varepsilon r} = r - \frac{\varepsilon r^2}{1+\varepsilon r}$$

となります。分子を最大化し分母を最小化することで、$r$ からどれだけ離れているかを以下のように評価できます。

$$\left\vert \frac{\varepsilon r^2}{1+\varepsilon r} \right\vert \leq \frac{ \frac{r^2}{2 r(r+1)}}{1 - \frac{r}{2r(r+1)}} = \frac{r}{2 r + 1} < \frac{1}{2}$$

$r$ から $1/2$ 未満の距離にありますので、予想どおり丸めると $r$ が得られます。

残念ながら、$r$ がまだわからないため、必要な精度を $r$ から直接求めることはできません。 代わりに、$r$ は $N$ より小さいという事実を利用して、十分な精度を確保することができます。 具体的には、$1/r$ に対する最良近似 $y/2^m$ が

$$\left\vert \frac{y}{2^m} - \frac{1}{r} \right\vert \leq \frac{1}{2N^2}$$

を満たすだけの精度を確保すれば、逆数を取ったときに $r$ を正しく決定できます。 $m = 2\operatorname{lg}(N)+1$ とすることで、前述の方法を使ってこの精度の推定が高い確率で得られるようになります。 （成功確率の下限が 40% でも問題なければ、$m = 2\operatorname{lg}(N)$ で十分です。）

一般的な解法

先ほど見たように、$M_a$ の固有ベクトル $\vert \psi_1 \rangle$ があれば、十分な数の制御量子ビットを使って精度を確保することで、位相推定を通じて $r$ を求めることができます。 しかし、固有ベクトル $\vert\psi_1\rangle$ を手に入れることは容易ではないため、どのように進めるかを考える必要があります。

一時的に、$\vert\psi_1\rangle$ の代わりに任意の $k\in\{0,\ldots,r-1\}$ に対する固有ベクトル $\vert\psi_k\rangle$ を用いて同様に進めると仮定しましょう。 位相推定の手順から得られる結果は、近似値

$$\frac{y}{2^m} \approx \frac{k}{r}$$

になります。

$k$ も $r$ もわからないという前提で考えると、これが $r$ の特定につながる場合もあればそうでない場合もあります。 例えば、$k = 0$ の場合は $0$ の近似値 $y/2^m$ が得られますが、残念ながらこれからは何もわかりません。 ただし、これは例外的なケースであり、$k$ が非ゼロの値であれば、少なくとも $r$ に関する何らかの情報を得ることができます。

連分数アルゴリズムと呼ばれるアルゴリズムを使うと、近似値 $y/2^m$ を近傍の分数に変換できます。近似が十分に精度よければ $k/r$ も求めることができます。 ここでは連分数アルゴリズムの説明は省略します。 代わりに、このアルゴリズムに関する既知の事実を述べます。

事実

整数 $N\geq 2$ と実数 $\alpha\in(0,1)$ が与えられたとき、$v\neq 0$ かつ $\gcd(u,v)=1$ を満たし、$\vert \alpha - u/v\vert < \frac{1}{2N^2}$ を満たす整数 $u,v\in\{0,\ldots,N-1\}$ の選択は高々1つしか存在しません。 $\alpha$ と $N$ が与えられると、連分数アルゴリズムは $u$ と $v$ を求めるか、それらが存在しないことを報告します。 このアルゴリズムは $O((\operatorname{lg}(N))^3)$ のサイズを持つブール回路として実装できます。

$k/r$ に対する非常に精度の高い近似 $y/2^m$ があり、$N$ と $\alpha = y/2^m$ に対して連分数アルゴリズムを実行すると、事実に記載された $u$ と $v$ が得られます。 この事実の分析から、次の結論が導かれます。

$$\frac{u}{v} = \frac{k}{r}.$$

特に注意すべき点として、必ずしも $k$ と $r$ を個別に学習するわけではなく、既約分数の形で $k/r$ のみを学習することになります。

例えば、すでに述べたように、$k=0$ からは何も学べません。 ただし、それが起きるのは $k=0$ の場合だけです。 $k$ が非ゼロの場合、$k$ と $r$ の間に公約数があるかもしれませんが、連分数アルゴリズムで得られる $v$ は少なくとも $r$ の約数になります。

自明ではありませんが、$k\in\{0,\ldots,r-1\}$ を一様ランダムに選んで $u/v = k/r$ の $u$ と $v$ を学習できる能力があれば、数回のサンプリングで高確率に $r$ を復元できることは事実です。 具体的には、観測した分母 $v$ のすべての値の最小公倍数を $r$ の推定値とすれば、高い確率で正解になります。 直感的に言えば、$k$ の値によっては $r$ との公約数があるために都合が悪く、$u$ と $v$ を学習しても $r$ の因数が隠されてしまうことがあります。 しかし、ランダムな $k$ の選択は長期間 $r$ の因数を隠し続けることはほぼなく、観測した分母の最小公倍数による推定が失敗する確率は、サンプル数の増加に対して指数関数的に減少します。

残る問題は、位相推定を実行するための $M_a$ の固有ベクトル $\vert\psi_k\rangle$ をどのように手に入れるかです。 実は、固有ベクトルを実際に準備する必要はないのです！

代わりに行うのは、固有ベクトル $\vert\psi\rangle$ の代わりに、数 $1$ の $n$ ビット二進数表現を意味する状態 $\vert 1\rangle$ を用いて位相推定を実行することです。 これまでは特定の固有ベクトルに対して位相推定を実行することについてのみ説明しましたが、$M_a$ の固有ベクトルでない入力状態に対して手順を実行することも問題ありません。ここではまさにその状態 $\vert 1\rangle$ を用いています。 （$a=1$ でない限り、これは $M_a$ の固有ベクトルではありませんが、$a=1$ は興味の対象となる選択ではありません。）

$M_a$ の固有ベクトルの代わりに状態 $\vert 1\rangle$ を選ぶ根拠は、次の等式が成り立つことです。

$$\vert 1\rangle = \frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \vert \psi_k\rangle$$

この等式を確認する一つの方法は、両辺と各標準基底状態の内積を比較することです。右辺の評価には、本レッスンで前述した公式を活用します。 この結果として、$k\in\{0,\ldots,r-1\}$ を一様ランダムに選んで固有ベクトル $\vert\psi_k\rangle$ を使った場合とまったく同じ測定結果が得られます。

さらに詳しく説明すると、固有ベクトル $\vert\psi_k\rangle$ の1つの代わりに状態 $\vert 1\rangle$ を用いて位相推定を実行するとしましょう。 逆量子フーリエ変換が実行された後、残る状態は

$$\frac{1}{\sqrt{r}} \sum_{k = 0}^{r-1} \vert \psi_k\rangle \vert \gamma_k\rangle$$

となります。ここで

$$\vert\gamma_k\rangle = \frac{1}{2^m} \sum_{y=0}^{2^m - 1} \sum_{x=0}^{2^m-1} e^{2\pi i x (k/r - y/2^m)} \vert y\rangle$$

です。ベクトル $\vert\gamma_k\rangle$ は、逆量子フーリエ変換が適用された後の上位 $m$ 量子ビットの状態を表しています。

したがって、$\{\vert\psi_0\rangle,\ldots,\vert\psi_{r-1}\rangle\}$ が正規直交集合であるという事実から、上位 $m$ 量子ビットの測定により、$k\in\{0,\ldots,r-1\}$ が一様ランダムに選ばれた値 $k/r$ に対する近似値 $y/2^m$ が得られることがわかります。 すでに述べたように、これにより複数回の独立した実行を経て高い確信度で $r$ を求めることができ、これが目標でした。

総コスト

各制御ユニタリ $M_a^k$ を実装するコストは $O(n^2)$ です。 制御ユニタリ演算は $m$ 個あり、$m = O(n)$ ですので、制御ユニタリ演算の総コストは $O(n^3)$ となります。 さらに、$m$ 個のアダマールゲート（コストへの寄与は $O(n)$）と、$O(n^2)$ を寄与する逆量子フーリエ変換があります。 したがって、制御ユニタリ演算のコストが手順全体のコストを支配し、全体コストは $O(n^3)$ となります。

量子回路自体に加えて、途中でいくつかの古典的計算も必要です。 これには、制御ユニタリゲートの作成に必要な $\mathbb{Z}_N$ における冪 $a^k$（$k = 2, 4, 8, \ldots, 2^{m-1}$）の計算、および $\theta$ の近似値を分数に変換する連分数アルゴリズムが含まれます。 これらの計算は、総コスト $O(n^3)$ のブール回路で実行できます。

通常のケースと同様に、漸近的に高速なアルゴリズムを使用することでこれらの上限をすべて改善できます。上記の上限は基本的な算術演算に標準的なアルゴリズムを使用していると仮定しています。

位数発見による素因数分解

最後に議論すべきことは、位数発見問題を解くことが素因数分解にどのように役立つかです。 この部分は完全に古典的なものであり、量子計算とは直接関係しません。

基本的な考え方は以下のとおりです。 数 $N$ を素因数分解したいと考えており、これを再帰的に行うことができます。 具体的には、$N$ の分割（任意の2つの整数 $b,c\geq 2$ で $N = bc$ を満たすものを見つけること）に集中できます。 $N$ が素数の場合はこれが不可能ですが、まず素数判定アルゴリズムで $N$ が素数かどうかを効率的に確認でき、素数でなければ分割を試みます。 $N$ を分割したら、$b$ と $c$ に対して再帰を繰り返し、すべての因数が素数になって $N$ の素因数分解が得られるまで続けます。

偶数の分割は簡単で、$2$ と $N/2$ を出力するだけです。

完全冪数、すなわち整数 $s,j\geq 2$ に対して $N = s^j$ の形の数の分割も容易です。 $N^{1/2},$ $N^{1/3},$ $N^{1/4}$ などの根を近似し、$s$ の候補として近傍の整数を確認することで行えます。 根が $2$ を下回って追加の候補が現れなくなるため、この数列の $\log(N)$ ステップ先まで調べれば十分です。

これらの両方ができることは重要です。位数発見は偶数の因数分解や、$s$ が素数である素数冪の因数分解には役立たないからです。 しかし、$N$ が奇数かつ素数冪でない場合は、位数発見を使って $N$ を分割することができます。

素数冪でない奇数の合成数 N を分割する確率的アルゴリズム

1. $a\in\{2,\ldots,N-1\}$ をランダムに選ぶ。

2. $d=\gcd(a,N)$ を計算する。

3. $d > 1$ ならば $b = d$ と $c = N/d$ を出力して終了する。そうでなければ $a\in\mathbb{Z}_N^{\ast}$ であることがわかっているので次のステップへ進む。

4. $r$ を $a$ の $N$ に対する位数とする（ここで位数発見が必要になります）。

5. $r$ が偶数の場合：

   5.1 $x = a^{r/2} - 1$（$N$ を法として）を計算する
   5.2 $d = \gcd(x,N)$ を計算する。
   5.3 $d>1$ ならば $b=d$ と $c = N/d$ を出力して終了する。

6. この地点に達した場合、アルゴリズムは $N$ の因数を見つけることに失敗しています。

このアルゴリズムの実行は $N$ の因数を見つけることに失敗する場合があります。 具体的には、次の2つの状況で失敗します。

• $a$ の $N$ に対する位数が奇数である。
• $a$ の $N$ に対する位数が偶数であり、$\gcd\bigl(a^{r/2} - 1, N\bigr) = 1$ である。

基本的な整数論を用いると、ランダムな $a$ の選択に対して、少なくとも $1/2$ の確率でこれらのいずれの事象も起きないことが証明できます。 実際、どちらかの事象が起きる確率は $2^{-(m-1)}$ 以下であり、ここで $m$ は $N$ の異なる素因数の数です。 これが $N$ が素数冪でないという仮定が必要な理由です。 （$N$ が奇数であるという仮定も、この事実が成立するために必要です。）

つまり、各実行で $N$ を分割できる確率は少なくとも 50% です。 したがって、$a$ をそれぞれランダムに選びながらアルゴリズムを $t$ 回実行すれば、少なくとも $1 - 2^{-t}$ の確率で $N$ の分割に成功します。

アルゴリズムの基本的な考え方は次のとおりです。 位数 $r$ が偶数となる $a$ の選択があれば、$r/2$ は整数となり、次の数を考えることができます。

$$a^{r/2} - 1\; (\textrm{mod}\; N) \quad \text{and} \quad a^{r/2} + 1\; (\textrm{mod}\; N).$$

公式 $Z^2 - 1 = (Z+1)(Z-1)$ を使うと、

$$\bigl(a^{r/2} - 1\bigr) \bigl(a^{r/2} + 1\bigr) = a^r - 1$$

が導かれます。

位数の定義から $a^r \; (\textrm{mod}\; N) = 1$ であることがわかります。つまり $N$ は $a^r - 1$ を整除するということです。 したがって $N$ は次の積を整除します。

$$\bigl(a^{r/2} - 1\bigr) \bigl(a^{r/2} + 1\bigr).$$

これが成り立つためには、$N$ のすべての素因数が $a^{r/2} - 1$ または $a^{r/2} + 1$（あるいは両方）の素因数でなければなりません。ランダムな $a$ の選択では、$N$ のすべての素因数が一方の項のみを割り切り、もう一方は割り切らないことはほぼ起こりません。 そうでなければ、$N$ の素因数の一部が第一項を割り切り、残りが第二項を割り切る限り、第一項との GCD を計算することで $N$ の非自明な因数を見つけることができます。