Algorytm Shora

Teraz skupimy naszą uwagę na problemie faktoryzacji liczb całkowitych i zobaczymy, jak można go efektywnie rozwiązać na komputerze kwantowym, używając estymacji fazy. Algorytm, który otrzymamy, to algorytm Shora faktoryzacji liczb całkowitych. Shor nie opisał swojego algorytmu konkretnie w kategoriach estymacji fazy, ale jest to naturalny i intuicyjny sposób wyjaśnienia, jak on działa.

Zaczniemy od omówienia pośredniego problemu znanego jako problem znajdowania rzędu i zobaczymy, jak estymacja fazy dostarcza rozwiązania tego problemu. Następnie zobaczymy, jak efektywne rozwiązanie problemu znajdowania rzędu daje nam efektywne rozwiązanie problemu faktoryzacji liczb całkowitych. (Gdy rozwiązanie jednego problemu zapewnia rozwiązanie innego problemu w ten sposób, mówimy, że ten drugi problem redukuje się do pierwszego — więc w tym przypadku redukujemy faktoryzację liczb całkowitych do znajdowania rzędu.) Ta druga część algorytmu Shora w ogóle nie korzysta z obliczeń kwantowych; jest całkowicie klasyczna. Obliczenia kwantowe są potrzebne tylko do rozwiązania problemu znajdowania rzędu.

Problem znajdowania rzędu

Kilka podstawowych pojęć z teorii liczb

Aby wyjaśnić problem znajdowania rzędu i sposób, w jaki można go rozwiązać za pomocą estymacji fazy, pomocne będzie rozpoczęcie od kilku podstawowych pojęć z teorii liczb oraz wprowadzenie po drodze przydatnej notacji.

Na początek, dla dowolnej dodatniej liczby całkowitej $N,$ zdefiniujmy zbiór $\mathbb{Z}_N$ w następujący sposób.

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

Na przykład, $\mathbb{Z}_1 = \{0\},\;$ $\mathbb{Z}_2 = \{0,1\},\;$ $\mathbb{Z}_3 = \{0,1,2\},\;$ i tak dalej.

Są to zbiory liczb, ale możemy myśleć o nich jako o czymś więcej niż tylko zbiorach. W szczególności możemy rozważać operacje arytmetyczne na $\mathbb{Z}_N,$ takie jak dodawanie i mnożenie — a jeśli umówimy się, że zawsze bierzemy nasze wyniki modulo $N$ (czyli dzielimy przez $N$ i bierzemy resztę jako wynik), zawsze pozostaniemy wewnątrz tego zbioru podczas wykonywania tych operacji. Te dwie konkretne operacje dodawania i mnożenia, obie brane modulo $N,$ czynią $\mathbb{Z}_N$ pierścieniem, który jest fundamentalnie ważnym typem obiektu w algebrze.

Na przykład, $3$ i $5$ są elementami $\mathbb{Z}_7,$ a jeśli pomnożymy je przez siebie, otrzymamy $3\cdot 5 = 15,$ co daje resztę $1$ po podzieleniu przez $7.$ Czasami wyrażamy to w następujący sposób.

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

Ale możemy również po prostu napisać $3 \cdot 5 = 1,$ o ile jasne jest, że pracujemy w $\mathbb{Z}_7,$ po to tylko, aby utrzymać naszą notację tak prostą, jak to możliwe.

Jako przykład, oto tabliczki dodawania i mnożenia dla $\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}$$

Spośród $N$ elementów $\mathbb{Z}_N,$ szczególne są elementy $a\in\mathbb{Z}_N$ spełniające $\gcd(a,N) = 1.$ Często zbiór zawierający te elementy oznacza się gwiazdką, w następujący sposób.

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

Jeśli skupimy naszą uwagę na operacji mnożenia, zbiór $\mathbb{Z}_N^{\ast}$ tworzy grupę — konkretnie grupę abelową — która jest kolejnym ważnym typem obiektu w algebrze. Podstawowym faktem dotyczącym tych zbiorów (i skończonych grup w ogóle) jest to, że jeśli wybierzemy dowolny element $a\in\mathbb{Z}_N^{\ast}$ i będziemy wielokrotnie mnożyć $a$ przez siebie, to w końcu zawsze otrzymamy liczbę $1.$

Jako pierwszy przykład weźmy $N=6.$ Mamy, że $5\in\mathbb{Z}_6^{\ast},$ ponieważ $\gcd(5,6) = 1,$ a jeśli pomnożymy $5$ przez siebie, otrzymamy $1,$ co potwierdza powyższa tabela.

$$5^2 = 1 \quad \text{(working within $\mathbb{Z}_6$)}$$

Jako drugi przykład weźmy $N = 21.$ Jeśli przejdziemy przez liczby od $0$ do $20,$ to te, które mają NWD równe $1$ z $21,$ są następujące.

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

Dla każdego z tych elementów możliwe jest podniesienie tej liczby do dodatniej potęgi całkowitej, aby otrzymać $1.$ Oto najmniejsze potęgi, dla których to działa:

$$\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}$$

Naturalnie pracujemy w $\mathbb{Z}_{21}$ dla wszystkich tych równań, czego nie zawracaliśmy sobie głowy zapisywaniem — przyjmujemy to za domyślne, aby uniknąć zaśmiecania. Będziemy to kontynuować przez resztę lekcji.

Sformułowanie problemu i związek z estymacją fazy

Możemy teraz sformułować problem znajdowania rzędu.

Znajdowanie rzędu

Wejście: dodatnie liczby całkowite $N$ i $a$ spełniające $\gcd(N,a) = 1$
Wyjście: najmniejsza dodatnia liczba całkowita $r$ taka, że $a^r \equiv 1$ $(\textrm{mod } N)$

Alternatywnie, używając notacji, którą właśnie wprowadziliśmy powyżej, mamy dane $a \in \mathbb{Z}_N^{\ast},$ i szukamy najmniejszej dodatniej liczby całkowitej $r$ takiej, że $a^r = 1.$ Tę liczbę $r$ nazywamy rzędem elementu $a$ modulo $N.$

Aby powiązać problem znajdowania rzędu z estymacją fazy, rozważmy operację zdefiniowaną na układzie, którego stany klasyczne odpowiadają $\mathbb{Z}_N,$ gdzie mnożymy przez ustalony element $a\in\mathbb{Z}_N^{\ast}.$

$$M_a \vert x\rangle = \vert ax \rangle \qquad \text{(dla każdego $x\in\mathbb{Z}_N$)}$$

Dla jasności: mnożenie wykonujemy w $\mathbb{Z}_N,$ więc niejawnie zakładamy, że bierzemy iloczyn modulo $N$ wewnątrz ketu po prawej stronie równania.

Na przykład, jeśli weźmiemy $N = 15$ i $a=2,$ to działanie $M_2$ na bazie standardowej $\{\vert 0\rangle,\ldots,\vert 14\rangle\}$ wygląda następująco.

$$\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}$$

Jest to operacja unitarna, pod warunkiem że $\gcd(a,N)=1;$ permutuje ona elementy bazy standardowej $\{\vert 0\rangle,\ldots,\vert N-1\rangle\},$ więc jako macierz jest macierzą permutacji. Z jej definicji jasno wynika, że operacja ta jest deterministyczna, a prostym sposobem, by zobaczyć, że jest odwracalna, jest rozważenie rzędu $r$ elementu $a$ modulo $N$ i zauważenie, że odwrotnością $M_a$ jest $M_a^{r-1}.$

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

Istnieje inny sposób myślenia o odwrotności, który nie wymaga żadnej wiedzy o $r$ (które w końcu jest tym, co próbujemy obliczyć). Dla każdego elementu $a\in\mathbb{Z}_N^{\ast}$ zawsze istnieje jednoznaczny element $b\in\mathbb{Z}_N^{\ast}$ spełniający $ab=1.$ Ten element $b$ oznaczamy przez $a^{-1}$ i można go wydajnie obliczyć; rozszerzenie algorytmu GCD Euklidesa robi to z kosztem kwadratowym względem $\operatorname{lg}(N).$ A zatem

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

Tak więc operacja $M_a$ jest zarówno deterministyczna, jak i odwracalna. Oznacza to, że jest opisana przez macierz permutacji, a więc jest unitarna.

Zastanówmy się teraz nad eigenvectorami i eigenvalue'ami operacji $M_a,$ zakładając, że $a\in\mathbb{Z}_N^{\ast}.$ Jak właśnie uzasadniono, to założenie mówi nam, że $M_a$ jest unitarne.

Istnieje $N$ eigenvalue'ów $M_a,$ być może zawierających tę samą eigenvalue powtórzoną wielokrotnie, i ogólnie istnieje pewna dowolność w wyborze odpowiadających eigenvectorów — ale nie musimy się martwić wszystkimi możliwościami. Zacznijmy prosto i zidentyfikujmy tylko jeden eigenvector $M_a.$

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

Liczba $r$ to rząd $a$ modulo $N,$ tutaj i w całej pozostałej części lekcji. Eigenvalue związana z tym eigenvectorem to $1,$ ponieważ nie zmienia się on, gdy mnożymy przez $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$$

Dzieje się tak, ponieważ $a^r = 1,$ więc każdy stan bazy standardowej $\vert a^k \rangle$ zostaje przesunięty do $\vert a^{k+1} \rangle$ dla $k\leq r-1,$ a $\vert a^{r-1} \rangle$ zostaje przesunięty z powrotem do $\vert 1\rangle.$ Mówiąc nieformalnie, to tak, jakbyśmy powoli mieszali $\vert \psi_0 \rangle,$ ale jest już całkowicie wymieszany, więc nic się nie zmienia.

Oto kolejny przykład eigenvectora $M_a.$ Ten jest bardziej interesujący w kontekście znajdowania rzędu i estymacji fazy.

$$\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}}$$

Alternatywnie, możemy zapisać ten wektor za pomocą sumowania w następujący sposób.

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

Tutaj widzimy, że liczba zespolona $\omega_r = e^{2\pi i/r}$ pojawia się w naturalny sposób, ze względu na to, jak działa mnożenie przez $a$ modulo $N.$ Tym razem odpowiadająca eigenvalue to $\omega_r.$ Aby to zobaczyć, możemy najpierw obliczyć w następujący sposób.

$$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$$

Następnie, ponieważ $\omega_r^{-r} = 1 = \omega_r^0$ oraz $\vert a^r \rangle = \vert 1\rangle = \vert a^0\rangle,$ widzimy, że

$$\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,$$

a zatem $M_a \vert\psi_1\rangle = \omega_r \vert\psi_1\rangle.$

Stosując to samo rozumowanie, możemy zidentyfikować dodatkowe pary eigenvector/eigenvalue dla $M_a.$ Dla dowolnego wyboru $j\in\{0,\ldots,r-1\}$ mamy, że

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

jest eigenvectorem $M_a,$ którego odpowiadająca eigenvalue to $\omega_r^j.$

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

Istnieją inne eigenvectory $M_a,$ ale nie musimy się nimi zajmować — skupimy się wyłącznie na eigenvectorach $\vert\psi_0\rangle,\ldots,\vert\psi_{r-1}\rangle,$ które właśnie zidentyfikowaliśmy.

Znajdowanie rzędu za pomocą estymacji fazy

Aby rozwiązać problem znajdowania rzędu dla danego wyboru $a\in\mathbb{Z}_N^{\ast},$ możemy zastosować procedurę estymacji fazy do operacji $M_a.$

W tym celu musimy wydajnie zaimplementować nie tylko $M_a$ za pomocą obwodu kwantowego, lecz także $M_a^2,$ $M_a^4,$ $M_a^8$ i tak dalej, posuwając się na tyle daleko, na ile to konieczne, aby uzyskać wystarczająco precyzyjne oszacowanie z procedury estymacji fazy. Wyjaśnimy tutaj, jak można to zrobić, a później ustalimy dokładnie, jaka precyzja jest potrzebna.

Zacznijmy od samej operacji $M_a.$ Naturalnie, ponieważ pracujemy z modelem obwodów kwantowych, użyjemy notacji binarnej do zakodowania liczb pomiędzy $0$ a $N-1.$ Największą liczbą, którą musimy zakodować, jest $N-1,$ zatem liczba potrzebnych bitów wynosi

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

Na przykład, jeśli $N = 21,$ mamy $n = \operatorname{lg}(N-1) = 5.$ Oto jak wygląda kodowanie elementów $\mathbb{Z}_{21}$ jako ciągów binarnych długości $5.$

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

A teraz oto precyzyjna definicja tego, w jaki sposób $M_a$ jest zdefiniowana jako operacja $n$-kubitowa.

$$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}$$

Chodzi o to, że chociaż interesuje nas tylko, jak $M_a$ działa dla $\vert 0\rangle,\ldots,\vert N-1\rangle,$ musimy określić, jak działa dla pozostałych $2^n - N$ stanów bazy standardowej — i musimy to zrobić w taki sposób, aby nadal otrzymać operację unitarną. Zdefiniowanie $M_a$ tak, aby nic nie robiła z pozostałymi stanami bazy standardowej, spełnia ten warunek.

Korzystając z algorytmów mnożenia i dzielenia liczb całkowitych omówionych w poprzedniej lekcji, wraz z metodologią ich odwracalnych, bezśmieciowych implementacji, możemy zbudować obwód kwantowy realizujący $M_a,$ dla dowolnego wyboru $a\in\mathbb{Z}_N^{\ast},$ o koszcie $O(n^2).$ Oto jeden ze sposobów, w jaki można to zrobić.

1. Zbuduj obwód realizujący operację

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

   gdzie

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

   przy użyciu metody opisanej w poprzedniej lekcji. Daje to obwód o rozmiarze $O(n^2).$

2. Zamień dwa $n$-kubitowe układy, używając $n$ bramek swap do indywidualnej zamiany kubitów.

3. Podobnie jak w pierwszym kroku, zbuduj obwód dla operacji

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

   gdzie $a^{-1}$ jest odwrotnością $a$ w $\mathbb{Z}_N^{\ast}.$

Inicjalizując dolne $n$ kubitów i składając te trzy kroki, uzyskujemy następującą transformację:

$$\vert x \rangle \vert 0^n \rangle \stackrel{\text{krok 1}}{\mapsto} \vert x \rangle \vert f_a(x)\rangle \stackrel{\text{krok 2}}{\mapsto} \vert f_a(x)\rangle \vert x \rangle \stackrel{\text{krok 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$$

Metoda wymaga kubitów roboczych, ale są one na końcu przywracane do stanu początkowego, co pozwala nam wykorzystać te obwody do estymacji fazy. Całkowity koszt uzyskanego obwodu wynosi $O(n^2).$

Aby wykonać $M_a^2,$ $M_a^4,$ $M_a^8$ i tak dalej, możemy zastosować dokładnie tę samą metodę, z tą różnicą, że zastępujemy $a$ przez $a^2,$ $a^4,$ $a^8$ i tak dalej, jako elementy $\mathbb{Z}_N^{\ast}.$ To znaczy, dla dowolnej potęgi $k,$ którą wybierzemy, możemy utworzyć obwód dla $M_a^k$ nie przez $k$-krotne iterowanie obwodu dla $M_a,$ lecz obliczając $b = a^k \in \mathbb{Z}_N^{\ast}$ i następnie używając obwodu dla $M_b.$

Obliczanie potęg $a^k \in \mathbb{Z}_N$ to problem potęgowania modularnego wspomniany w poprzedniej lekcji. Obliczenie to można wykonać klasycznie, używając algorytmu potęgowania modularnego wspomnianego w poprzedniej lekcji (często nazywanego algorytmem potęgowania w obliczeniowej teorii liczb). W rzeczywistości potrzebujemy jedynie potęg $a$ o wykładnikach będących potęgami dwójki, a mianowicie $a^2, a^4, \ldots a^{2^{m-1}} \in \mathbb{Z}_N^{\ast},$ a potęgi te możemy uzyskać poprzez iteracyjne podnoszenie do kwadratu $m-1$ razy. Każde podniesienie do kwadratu można wykonać przez obwód logiczny o rozmiarze $O(n^2).$

W istocie to, co tu faktycznie robimy, to przerzucenie problemu iterowania $M_a$ aż $2^{m-1}$ razy na wydajne obliczenia klasyczne. I mamy szczęście, że jest to możliwe! Dla dowolnego wyboru obwodu kwantowego w problemie estymacji fazy jest to mało prawdopodobne — a w takim przypadku wynikowy koszt estymacji fazy rośnie wykładniczo wraz z liczbą kubitów sterujących $m.$

Rozwiązanie przy zadanym dogodnym eigenvectorze

Aby zrozumieć, jak możemy rozwiązać problem znajdowania rzędu za pomocą estymacji fazy, zacznijmy od założenia, że uruchamiamy procedurę estymacji fazy na operacji $M_a$ używając eigenvectora $\vert\psi_1\rangle.$ Zdobycie tego eigenvectora nie jest łatwe, jak się okaże, więc nie będzie to koniec historii — ale warto od tego zacząć.

Eigenvalue operacji $M_a$ odpowiadająca eigenvectorowi $\vert \psi_1\rangle$ wynosi

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

To znaczy $\omega_r = e^{2\pi i \theta}$ dla $\theta = 1/r.$ Zatem, jeśli uruchomimy procedurę estymacji fazy na $M_a$ używając eigenvectora $\vert\psi_1\rangle,$ otrzymamy przybliżenie $1/r.$ Obliczając odwrotność będziemy w stanie poznać $r$ — pod warunkiem, że nasze przybliżenie jest wystarczająco dobre.

Bardziej szczegółowo: gdy uruchamiamy procedurę estymacji fazy z $m$ kubitami sterującymi, otrzymujemy liczbę $y\in\{0,\ldots,2^m-1\}.$ Następnie przyjmujemy $y/2^m$ jako przybliżenie $\theta,$ którym w rozważanym przypadku jest $1/r.$ Aby z tego przybliżenia odgadnąć, czym jest $r,$ naturalną rzeczą jest obliczenie odwrotności naszego przybliżenia i zaokrąglenie do najbliższej liczby całkowitej.

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

Przykładowo, załóżmy, że $r = 6$ i wykonujemy estymację fazy na $M_a$ z eigenvectorem $\vert\psi_1\rangle$ używając $m = 5$ bitów sterujących. Najlepszym $5$-bitowym przybliżeniem $1/r = 1/6$ jest $5/32,$ i mamy dość dużą szansę (około $68\%$ w tym przypadku) uzyskać wynik $y=5$ z estymacji fazy. Mamy

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

a zaokrąglenie do najbliższej liczby całkowitej daje $6,$ co jest prawidłową odpowiedzią.

Z drugiej strony, jeśli nie użyjemy wystarczającej precyzji, możemy nie otrzymać prawidłowej odpowiedzi. Na przykład, jeśli weźmiemy $m = 4$ kubity sterujące w estymacji fazy, możemy uzyskać najlepsze $4$-bitowe przybliżenie $1/r = 1/6,$ które wynosi $3/16.$ Wzięcie odwrotności daje

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

a zaokrąglenie do najbliższej liczby całkowitej daje błędną odpowiedź $5.$

Jak dużej precyzji potrzebujemy zatem, aby uzyskać prawidłową odpowiedź? Wiemy, że rząd $r$ jest liczbą całkowitą, i intuicyjnie mówiąc potrzebujemy precyzji wystarczającej do odróżnienia $1/r$ od pobliskich możliwości, w tym $1/(r+1)$ i $1/(r-1).$ Najbliższą liczbą do $1/r,$ którą musimy się przejmować, jest $1/(r+1),$ a odległość między tymi dwiema liczbami wynosi

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

Zatem, jeśli chcemy mieć pewność, że nie pomylimy $1/r$ z $1/(r+1),$ wystarczy użyć precyzji wystarczającej do zagwarantowania, że najlepsze przybliżenie $y/2^m$ liczby $1/r$ jest bliższe $1/r$ niż $1/(r+1).$ Jeśli użyjemy na tyle dużej precyzji, że

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

tak aby błąd był mniejszy niż połowa odległości między $1/r$ a $1/(r+1),$ wtedy $y/2^m$ będzie bliższe $1/r$ niż jakiejkolwiek innej możliwości, w tym $1/(r+1)$ i $1/(r-1).$

Możemy to dodatkowo zweryfikować w następujący sposób. Załóżmy, że

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

dla $\varepsilon$ spełniającego

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

Gdy weźmiemy odwrotność, otrzymamy

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

Maksymalizując licznik i minimalizując mianownik, możemy ograniczyć, jak daleko jesteśmy od $r,$ w następujący sposób.

$$\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}$$

Jesteśmy w odległości mniejszej niż $1/2$ od $r,$ więc zgodnie z oczekiwaniami otrzymamy $r$ po zaokrągleniu.

Niestety, ponieważ nie znamy jeszcze $r,$ nie możemy go użyć do określenia, jakiej dokładności potrzebujemy. Zamiast tego możemy wykorzystać fakt, że $r$ musi być mniejsze niż $N,$ aby zapewnić użycie wystarczającej precyzji. W szczególności, jeśli użyjemy dokładności wystarczającej do zagwarantowania, że najlepsze przybliżenie $y/2^m$ liczby $1/r$ spełnia

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

to będziemy mieli wystarczającą precyzję, aby poprawnie wyznaczyć $r$ po wzięciu odwrotności. Przyjęcie $m = 2\operatorname{lg}(N)+1$ zapewnia, że mamy dużą szansę uzyskać estymację z taką precyzją przy użyciu metody opisanej wcześniej. (Przyjęcie $m = 2\operatorname{lg}(N)$ jest wystarczające, jeśli akceptujemy dolne ograniczenie 40% na prawdopodobieństwo sukcesu.)

Rozwiązanie ogólne

Jak właśnie zobaczyliśmy, jeśli mamy eigenvector $\vert \psi_1 \rangle$ operatora $M_a,$ możemy poznać $r$ poprzez estymację fazy, o ile użyjemy wystarczającej liczby kubitów sterujących, aby zrobić to z wystarczającą precyzją. Niestety, nie jest łatwo zdobyć eigenvector $\vert\psi_1\rangle,$ więc musimy wymyślić, jak postępować dalej.

Załóżmy na moment, że postępujemy tak samo jak powyżej, ale z eigenvectorem $\vert\psi_k\rangle$ zamiast $\vert\psi_1\rangle,$ dla dowolnego wyboru $k\in\{0,\ldots,r-1\},$ który zechcemy rozważyć. Wynik, jaki otrzymamy z procedury estymacji fazy, będzie przybliżeniem

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

Przy założeniu, że nie znamy ani $k,$ ani $r,$ może to pozwolić nam zidentyfikować $r$ albo nie. Na przykład, jeśli $k = 0,$ otrzymamy przybliżenie $y/2^m$ do $0,$ co niestety nic nam nie mówi. To jest jednak nietypowy przypadek; dla innych wartości $k$ będziemy w stanie przynajmniej dowiedzieć się czegoś o $r.$

Możemy użyć algorytmu znanego jako algorytm ułamków łańcuchowych, aby zamienić nasze przybliżenie $y/2^m$ na pobliskie ułamki — w tym $k/r,$ jeśli przybliżenie jest wystarczająco dobre. Nie będziemy tutaj wyjaśniać algorytmu ułamków łańcuchowych. Zamiast tego, oto sformułowanie znanego faktu dotyczącego tego algorytmu.

Fakt

Dla danej liczby całkowitej $N\geq 2$ i liczby rzeczywistej $\alpha\in(0,1),$ istnieje co najwyżej jeden wybór liczb całkowitych $u,v\in\{0,\ldots,N-1\}$ spełniających $v\neq 0$ i $\gcd(u,v)=1$ oraz $\vert \alpha - u/v\vert < \frac{1}{2N^2}.$ Dla danych $\alpha$ i $N,$ algorytm ułamków łańcuchowych znajduje $u$ i $v,$ albo zgłasza, że nie istnieją. Algorytm ten można zaimplementować jako obwód boolowski o rozmiarze $O((\operatorname{lg}(N))^3).$

Jeśli mamy bardzo bliskie przybliżenie $y/2^m$ wartości $k/r$ i uruchamiamy algorytm ułamków łańcuchowych dla $N$ i $\alpha = y/2^m,$ otrzymamy $u$ i $v,$ jak opisano w powyższym fakcie. Analiza tego faktu pozwala nam wywnioskować, że

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

Zauważmy w szczególności, że niekoniecznie poznajemy $k$ i $r,$ poznajemy jedynie $k/r$ w postaci nieskracalnej.

Na przykład, jak już zauważyliśmy, niczego się nie dowiemy z $k=0.$ Ale to jedyna wartość $k,$ dla której tak się dzieje. Gdy $k$ jest niezerowe, może mieć wspólne czynniki z $r,$ ale liczba $v,$ którą otrzymujemy z algorytmu ułamków łańcuchowych, musi przynajmniej dzielić $r.$

Nie jest to oczywiste, ale jest prawdą, że jeśli mamy zdolność poznania $u$ i $v$ dla $u/v = k/r$ dla $k\in\{0,\ldots,r-1\}$ wybranego jednostajnie losowo, to z dużym prawdopodobieństwem będziemy w stanie odzyskać $r$ po zaledwie kilku próbkach. W szczególności, jeśli naszym zgadnięciem dla $r$ jest najmniejsza wspólna wielokrotność wszystkich wartości mianownika $v,$ które obserwujemy, to z dużym prawdopodobieństwem będziemy mieli rację. Intuicyjnie mówiąc, niektóre wartości $k$ nie są dobre, ponieważ mają wspólne czynniki z $r,$ a te wspólne czynniki są przed nami ukryte, gdy poznajemy $u$ i $v.$ Ale losowe wybory $k$ raczej nie będą ukrywać czynników $r$ przez długi czas, a prawdopodobieństwo, że nie zgadniemy $r$ poprawnie, biorąc najmniejszą wspólną wielokrotność obserwowanych mianowników, spada wykładniczo wraz z liczbą próbek.

Pozostaje rozwiązać problem, jak zdobyć eigenvector $\vert\psi_k\rangle$ operatora $M_a,$ na którym uruchomimy procedurę estymacji fazy. Jak się okazuje, wcale nie musimy ich tworzyć!

Zamiast tego uruchomimy procedurę estymacji fazy na stanie $\vert 1\rangle,$ przez co rozumiemy $n$-bitowe binarne kodowanie liczby $1,$ zamiast eigenvectora $\vert\psi\rangle$ operatora $M_a.$ Do tej pory mówiliśmy tylko o uruchamianiu procedury estymacji fazy na konkretnym eigenvectorze, ale nic nie stoi na przeszkodzie, aby uruchomić tę procedurę na stanie wejściowym, który nie jest eigenvectorem $M_a,$ i to właśnie robimy tutaj ze stanem $\vert 1\rangle.$ (Nie jest on eigenvectorem $M_a,$ chyba że $a=1,$ czym nie będziemy zainteresowani.)

Uzasadnieniem wyboru stanu $\vert 1\rangle$ zamiast eigenvectora $M_a$ jest to, że następujące równanie jest prawdziwe.

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

Jednym ze sposobów weryfikacji tego równania jest porównanie iloczynów skalarnych obu stron z każdym standardowym stanem bazowym, używając wzorów wspomnianych wcześniej w lekcji, aby pomóc w ocenie wyników dla prawej strony. W konsekwencji otrzymamy dokładnie te same wyniki pomiarów, jak gdybyśmy wybrali $k\in\{0,\ldots,r-1\}$ jednostajnie losowo i użyli $\vert\psi_k\rangle$ jako eigenvectora.

Bardziej szczegółowo, wyobraźmy sobie, że uruchamiamy procedurę estymacji fazy ze stanem $\vert 1\rangle$ zamiast jednego z eigenvectorów $\vert\psi_k\rangle.$ Po wykonaniu odwrotnej kwantowej transformaty Fouriera pozostaje nam stan

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

gdzie

$$\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.$$

Wektor $\vert\gamma_k\rangle$ reprezentuje stan górnych $m$ kubitów po wykonaniu na nich odwrotności kwantowej transformaty Fouriera.

Zatem, dzięki temu, że $\{\vert\psi_0\rangle,\ldots,\vert\psi_{r-1}\rangle\}$ jest zbiorem ortonormalnym, stwierdzamy, że pomiar górnych $m$ kubitów daje przybliżenie $y/2^m$ wartości $k/r,$ gdzie $k\in\{0,\ldots,r-1\}$ jest wybierane jednostajnie losowo. Jak już omówiliśmy, pozwala to nam poznać $r$ z wysokim stopniem pewności po kilku niezależnych uruchomieniach, co było naszym celem.

Całkowity koszt

Koszt implementacji każdej kontrolowanej operacji unitarnej $M_a^k$ wynosi $O(n^2).$ Mamy $m$ kontrolowanych operacji unitarnych, a ponieważ $m = O(n),$ całkowity koszt kontrolowanych operacji unitarnych wynosi $O(n^3).$ Dodatkowo mamy $m$ bramek Hadamarda (które dają wkład $O(n)$ do kosztu), a odwrotna kwantowa transformata Fouriera daje wkład $O(n^2)$ do kosztu. Zatem koszt kontrolowanych operacji unitarnych dominuje nad kosztem całej procedury — który w związku z tym wynosi $O(n^3).$

Oprócz samego obwodu kwantowego po drodze trzeba wykonać kilka klasycznych obliczeń. Obejmuje to obliczenie potęg $a^k$ w $\mathbb{Z}_N$ dla $k = 2, 4, 8, \ldots, 2^{m-1},$ które są potrzebne do skonstruowania kontrolowanych bramek unitarnych, a także algorytm ułamków łańcuchowych przekształcający przybliżenia $\theta$ na ułamki. Obliczenia te mogą być wykonane przez obwody Boole'a o łącznym koszcie $O(n^3).$

Jak to zwykle bywa, wszystkie te ograniczenia można poprawić, używając asymptotycznie szybkich algorytmów; podane ograniczenia zakładają stosowanie standardowych algorytmów dla podstawowych operacji arytmetycznych.

Faktoryzacja przez znajdowanie rzędu

Ostatnia rzecz, którą musimy omówić, to jak rozwiązanie problemu znajdowania rzędu pomaga nam w faktoryzacji. Ta część jest całkowicie klasyczna — nie ma w niej nic specyficznie związanego z obliczeniami kwantowymi.

Oto podstawowa idea. Chcemy rozłożyć liczbę $N$ na czynniki, i możemy to zrobić rekurencyjnie. Dokładniej, możemy skupić się na zadaniu rozszczepienia $N,$ co oznacza znalezienie dowolnych dwóch liczb całkowitych $b,c\geq 2,$ dla których $N = bc.$ Nie jest to możliwe, jeśli $N$ jest liczbą pierwszą, ale możemy najpierw efektywnie sprawdzić, czy $N$ jest pierwsze, używając algorytmu testowania pierwszości, a jeśli $N$ nie jest liczbą pierwszą, spróbujemy ją rozszczepić. Po rozszczepieniu $N$ możemy po prostu zastosować rekurencję na $b$ i $c$, aż wszystkie nasze czynniki będą pierwsze i otrzymamy rozkład $N$ na czynniki pierwsze.

Rozszczepianie liczb parzystych jest łatwe: po prostu zwracamy $2$ i $N/2.$

Łatwo jest również rozszczepiać potęgi doskonałe, czyli liczby postaci $N = s^j$ dla liczb całkowitych $s,j\geq 2,$ po prostu poprzez przybliżanie pierwiastków $N^{1/2},$ $N^{1/3},$ $N^{1/4},$ i tak dalej, oraz sprawdzanie pobliskich liczb całkowitych jako kandydatów na $s.$ Nie musimy posuwać się dalej niż $\log(N)$ kroków w tym ciągu, ponieważ w tym momencie pierwiastek spada poniżej $2$ i nie ujawni dodatkowych kandydatów.

Dobrze, że potrafimy zrobić obie te rzeczy, ponieważ znajdowanie rzędu nie pomoże nam w faktoryzacji liczb parzystych ani potęg liczb pierwszych, gdzie liczba $s$ akurat jest pierwsza. Jeśli jednak $N$ jest nieparzysta i nie jest potęgą liczby pierwszej, znajdowanie rzędu pozwala nam rozszczepić $N.$

Algorytm probabilistyczny rozszczepiający nieparzystą, złożoną liczbę całkowitą N, która nie jest potęgą liczby pierwszej

1. Wybierz losowo $a\in\{2,\ldots,N-1\}.$

2. Oblicz $d=\gcd(a,N).$

3. Jeśli $d > 1,$ to zwróć $b = d$ i $c = N/d$ i zatrzymaj się. W przeciwnym razie przejdź do następnego kroku, wiedząc, że $a\in\mathbb{Z}_N^{\ast}.$

4. Niech $r$ będzie rzędem $a$ modulo $N.$ (Tutaj potrzebujemy znajdowania rzędu.)

5. Jeśli $r$ jest parzyste:

   5.1 Oblicz $x = a^{r/2} - 1$ modulo $N$
   5.2 Oblicz $d = \gcd(x,N).$
   5.3 Jeśli $d>1,$ to zwróć $b=d$ i $c = N/d$ i zatrzymaj się.

6. Jeśli ten punkt został osiągnięty, algorytm nie znalazł czynnika $N.$

Pojedyncze uruchomienie tego algorytmu może nie znaleźć czynnika $N.$ Konkretnie, dzieje się tak w dwóch sytuacjach:

• Rząd $a$ modulo $N$ jest nieparzysty.
• Rząd $a$ modulo $N$ jest parzysty oraz $\gcd\bigl(a^{r/2} - 1, N\bigr) = 1.$

Korzystając z podstawowej teorii liczb, można udowodnić, że dla losowego wyboru $a,$ z prawdopodobieństwem co najmniej $1/2$ żadne z tych zdarzeń nie zachodzi. W rzeczywistości prawdopodobieństwo, że którekolwiek ze zdarzeń zachodzi, wynosi co najwyżej $2^{-(m-1)},$ gdzie $m$ jest liczbą różnych czynników pierwszych $N,$ co jest powodem, dla którego potrzebne jest założenie, że $N$ nie jest potęgą liczby pierwszej. (Założenie, że $N$ jest nieparzyste, jest również wymagane, aby ten fakt był prawdziwy.)

Oznacza to, że każde uruchomienie ma co najmniej 50% szans na rozszczepienie $N.$ Dlatego, jeśli uruchomimy algorytm $t$ razy, losowo wybierając $a$ za każdym razem, uda nam się rozszczepić $N$ z prawdopodobieństwem co najmniej $1 - 2^{-t}.$

Podstawowa idea algorytmu jest następująca. Jeśli wybierzemy $a,$ dla którego rząd $r$ liczby $a$ modulo $N$ jest parzysty, to $r/2$ jest liczbą całkowitą i możemy rozważyć liczby

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

Korzystając ze wzoru $Z^2 - 1 = (Z+1)(Z-1),$ wnioskujemy, że

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

Wiemy, że $a^r \; (\textrm{mod}\; N) = 1$ z definicji rzędu — co jest innym sposobem powiedzenia, że $N$ dzieli bez reszty $a^r - 1.$ Oznacza to, że $N$ dzieli bez reszty iloczyn

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

Aby tak było, wszystkie czynniki pierwsze $N$ muszą być również czynnikami pierwszymi $a^{r/2} - 1$ lub $a^{r/2} + 1$ (lub obu) — a dla losowego wyboru $a$ okazuje się, że jest mało prawdopodobne, aby wszystkie czynniki pierwsze $N$ dzieliły jeden z wyrazów, a żaden nie dzielił drugiego.