Algorytm Shora

Na potrzeby tego modułu Qiskit in Classrooms uczniowie muszą mieć działające środowisko Python z zainstalowanymi następującymi pakietami:

• qiskit w wersji 2.1.0 lub nowszej
• qiskit-ibm-runtime w wersji 0.40.1 lub nowszej
• qiskit-aer w wersji 0.17.0 lub nowszej
• qiskit.visualization
• numpy
• pylatexenc

Aby skonfigurować i zainstalować powyższe pakiety, zapoznaj się z przewodnikiem Instalacja Qiskit. Aby uruchamiać zadania na prawdziwych komputerach kwantowych, uczniowie muszą założyć konto w IBM Quantum®, postępując zgodnie z instrukcjami w przewodniku Konfiguracja konta IBM Cloud.

Ten moduł był testowany i wykorzystał trzy sekundy czasu QPU. To jedynie szacunek. Twoje rzeczywiste zużycie może być inne.

# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-ibm-runtime

# Uncomment and modify this line as needed to install dependencies
#!pip install 'qiskit>=2.1.0' 'qiskit-ibm-runtime>=0.40.1' 'qiskit-aer>=0.17.0' 'numpy' 'pylatexenc'

Wprowadzenie

Na początku lat 90. XX wieku narastał entuzjazm związany z potencjałem komputerów kwantowych do rozwiązywania problemów trudnych dla klasycznych komputerów. Kilku utalentowanych informatyków opracowało algorytmy demonstrujące moc obliczeń kwantowych dla wybranych, niszowych problemów, lecz nikomu nie udało się znaleźć jednej, przełomowej „zabójczej aplikacji" dla obliczeń kwantowych, która zrewolucjonizowałaby tę dziedzinę. Tak było aż do 1994 roku, kiedy Peter Shor wymyślił to, co dziś nazywamy algorytmem Shora – służący do faktoryzacji dużych liczb.

W tamtych czasach powszechnie wiadomo było, że znalezienie czynników pierwszych dużej liczby jest niezwykle trudne dla klasycznego komputera. Właśnie na tej trudności opierały się protokoły bezpieczeństwa w internecie. Shor odkrył sposób na wykładniczo bardziej efektywne znajdowanie tych czynników, przenosząc najtrudniejsze etapy obliczeń na teoretyczny, przyszły komputer kwantowy.

W tym module przyjrzymy się algorytmowi Shora. Najpierw zarysujemy szerszy kontekst, formalizując problem, który rozwiązuje, i wyjaśniając jego znaczenie dla cyberbezpieczeństwa. Następnie przedstawimy wprowadzenie do matematyki modularnej i pokażemy, jak zastosować ją do problemu faktoryzacji — wyjaśniając, jak faktoryzacja sprowadza się do innego zagadnienia, zwanego „znajdowaniem rzędu". Zobaczymy, jak Kwantowa Transformata Fouriera i Kwantowe Szacowanie Fazy, poznane w poprzednim module, wchodzą do gry, oraz jak ich użyć do rozwiązania problemu znajdowania rzędu.

Na koniec uruchomimy algorytm Shora na prawdziwym komputerze kwantowym! Pamiętaj jednak, że algorytm ten będzie naprawdę użyteczny dopiero wtedy, gdy będziemy dysponować dużym, odpornym na błędy komputerem kwantowym — co wciąż jest kwestią kilku najbliższych lat. Dlatego po prostu rozłożymy małą liczbę na czynniki, żeby pokazać, jak algorytm działa.

Problem faktoryzacji

Celem problemu faktoryzacji jest znalezienie czynników pierwszych liczby $N$. Dla niektórych wartości $N$ jest to dość proste. Na przykład, jeśli $N$ jest parzyste, jednym z jego czynników pierwszych będzie 2. Jeśli $N$ jest potęgą liczby pierwszej, tzn. $N=p^k$ dla pewnej liczby pierwszej $p$, znalezienie $p$ jest również stosunkowo łatwe: wystarczy przybliżyć $k$-ty pierwiastek z $N$ i poszukać pobliskich liczb pierwszych, które mogłyby być wartością $p$.

Klasyczne komputery mają jednak problemy w sytuacji, gdy $N$ jest nieparzyste i nie jest potęgą liczby pierwszej. To właśnie ten przypadek rozwiązuje algorytm Shora. Algorytm znajduje dwa czynniki $p$ i $q$ takie, że $N=pq$. Można go stosować rekurencyjnie, aż wszystkie czynniki będą pierwsze. W kolejnych sekcjach zobaczymy, jak podchodzi się do tego problemu.

Znaczenie dla cyberbezpieczeństwa

Wiele schematów kryptograficznych zostało zbudowanych w oparciu o fakt, że faktoryzacja dużych liczb jest trudna — w tym jeden powszechnie stosowany dziś schemat o nazwie RSA. W RSA klucz publiczny powstaje przez pomnożenie dwóch dużych liczb pierwszych, dając $N = p\cdot q$. Każdy może użyć tego klucza publicznego do szyfrowania danych. Jednak tylko ktoś posiadający klucz prywatny, $p$ i $q$, może je odszyfrować.

Gdyby $N$ dało się łatwo rozłożyć na czynniki, każdy mógłby ustalić wartości $p$ i $q$ i złamać szyfrowanie. Ale tak nie jest. To słynnie trudny problem. W rzeczywistości czynniki pierwsze liczby RSA1024 — mającej 1024 cyfry w zapisie binarnym i 309 cyfr w zapisie dziesiętnym — nie zostały znalezione do dziś, mimo że za jej faktoryzację jeszcze w 1991 roku zaoferowano nagrodę w wysokości $100 000.

Rozwiązanie Shora

W 1994 roku Peter Shor zdał sobie sprawę, że komputer kwantowy mógłby rozkładać dużą liczbę na czynniki wykładniczo efektywniej niż klasyczny komputer. Jego spostrzeżenie opierało się na związku między problemem faktoryzacji a arytmetyką modularną. Przedstawimy krótkie wprowadzenie do arytmetyki modularnej, a następnie zobaczymy, jak można jej użyć do faktoryzacji $N$.

Arytmetyka modularna

Arytmetyka modularna to cykliczny system liczenia, w którym liczenie zaczyna się w zwykły sposób — od liczb całkowitych 0, 1, 2 itd. — lecz po pewnym czasie, po upłynięciu okresu $N$, zaczyna się od nowa. Zobaczmy, jak to działa na przykładzie. Niech nasz okres wynosi 5. Wtedy podczas liczenia, w miejscu, gdzie normalnie pojawiłoby się 5, zaczynamy od nowa od 0:

$0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, ...$

Wynika to stąd, że w „świecie modulo 5" liczba 5 jest równoważna 0. Mówimy, że $5\bmod 5 \ = 0$. W istocie wszystkie wielokrotności 5 będą równoważne $0\bmod 5$.

Sprawdź swoją wiedzę

Przeczytaj poniższe pytania, zastanów się nad odpowiedzią, a następnie kliknij trójkąt, aby odkryć rozwiązanie.

Użyj arytmetyki modularnej, aby rozwiązać następujący problem:

Wyruszasz w długą, transkontynentalną podróż pociągiem o godzinie 8:00. Podróż trwa 60 godzin. O której godzinie przyjeżdżasz?

Odpowiedź:

Okres wynosi 24, ponieważ doba ma 24 godziny. Zadanie można więc zapisać w arytmetyce modularnej jako:

$(8+60)\text{mod}(24) = 20$

Dotrzesz więc do celu o godzinie 20:00, czyli o 20:00 wieczorem.

$\mathbb{Z}_N$ i $\mathbb{Z}_N^*$

Przydatne jest wprowadzenie dwóch zbiorów: $\mathbb{Z}_N$ i $\mathbb{Z}_N^*$. $\mathbb{Z}_N$ to po prostu zbiór liczb istniejących w „świecie modulo $N$". Na przykład, gdy liczyliśmy modulo 5, zbiór ten to $\mathbb{Z}_5=\{0,1,2,3,4\}$. Inny przykład: $\mathbb{Z}_{15} = \{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14\}$. Na elementach $\mathbb{Z}_N$ można wykonywać dodawanie i mnożenie (modulo $N$), a wynik każdej z tych operacji również należy do $\mathbb{Z}_N$, co czyni $\mathbb{Z}_N$ obiektem matematycznym zwanym pierścieniem.

Dla algorytmu Shora szczególnie interesujący jest pewien podzbiór $\mathbb{Z}_N$. Jest to podzbiór liczb z $\mathbb{Z}_N$, dla których największy wspólny dzielnik każdego elementu i $N$ wynosi 1, tzn. każdy element jest „wzajemnie pierwszym" z $N$. Jeśli weźmiemy zbiór tych liczb wraz z operacją mnożenia modularnego, otrzymamy inny obiekt matematyczny zwany grupą. Grupę tę oznaczamy $\mathbb{Z}_N^*$. Okazuje się, że w $\mathbb{Z}_N^*$ (i skończonych grupach w ogólności), jeśli wybierzemy dowolny element $a \in \mathbb{Z}_N^*$ i będziemy wielokrotnie mnożyć $a$ przez siebie, zawsze w końcu otrzymamy liczbę $1$. Minimalna liczba mnożeń $a$ przez siebie potrzebna do uzyskania $1$ nosi nazwę rzędu elementu $a$. Fakt ten będzie bardzo ważny w naszej dyskusji o tym, jak faktoryzować liczby.

Sprawdź swoją wiedzę

Przeczytaj poniższe pytania, zastanów się nad odpowiedzią, a następnie kliknij trójkąt, aby odkryć rozwiązanie.

Czym jest $\mathbb{Z}_{15}^*$?

Odpowiedź:

$\mathbb{Z}_{15}^* = \{1,2,4,7,8,11,13,14\}$

Wykluczyliśmy następujące liczby:

$$\begin{aligned} 3: GCD(3,15)=3 \\ 5: GCD(5,15)=5 \\ 6: GCD(6,15)=3 \\ 9: GCD(9,15)=3 \\ 10: GCD(10,15)=5 \\ 12: GCD(12,15)=3 \\ \end{aligned}$$

Jaki jest rząd każdego z elementów $\mathbb{Z}_{15}^*$?

Odpowiedź:

Rząd $r$ to najmniejsza liczba taka, że $a^r\text{mod}(15)=1$ dla każdego elementu $a$.

$$\begin{aligned} 1^1\text{mod}(15) = 1, r=1 \\ 2^4\text{mod}(15) = 1, r=4 \\ 4^2\text{mod}(15) = 1, r=2 \\ 7^4\text{mod}(15) = 1, r=4 \\ 8^4\text{mod}(15) = 1, r=4 \\ 11^2\text{mod}(15) = 1, r=2 \\ 13^4\text{mod}(15) = 1, r=4 \\ 14^2\text{mod}(15) = 1, r=2 \\ \end{aligned}$$

Zauważ, że choć byliśmy w stanie wyznaczyć rząd liczb w $\mathbb{Z}_{15}^*$, w ogólności — dla większych $N$ — NIE jest to łatwe zadanie. To sedno problemu faktoryzacji i powód, dla którego potrzebujemy komputera kwantowego. Zobaczymy dlaczego, gdy będziemy przechodzić przez dalszą część tego notesu.

Zastosowanie arytmetyki modularnej do problemu faktoryzacji

Kluczem do znalezienia czynników $p$ i $q$ takich, że $N=pq$, jest znalezienie pewnej innej liczby całkowitej $x$ spełniającej

$x^2 \equiv 1 \bmod N$ i $x \not\equiv \pm 1 \bmod N.$

Jak znalezienie $x$ pomaga nam wyznaczyć czynniki $p$ i $q$? Przeprowadźmy teraz to rozumowanie. Ponieważ $x^2 \equiv 1 \bmod N$, oznacza to, że $x^2 - 1 \equiv 0 \bmod N$. Innymi słowy, $x^2 - 1$ jest wielokrotnością $N$. A zatem dla pewnej liczby całkowitej $l$:

$x^2 - 1 = l N$

Możemy rozłożyć $x^2 - 1$ na czynniki:

$(x+1)(x-1) = l N$

Z naszych początkowych założeń wiemy, że $x \not\equiv \pm 1 \bmod N$, więc $N$ nie dzieli równo ani $x+1$, ani $x-1$. Zatem dwa czynniki $N$, czyli $p$ i $q$, muszą każdy dzielić odpowiednio $x-1$ i $x+1$. Albo $p$ jest czynnikiem $x-1$, a $q$ czynnikiem $x+1$, albo odwrotnie. Dlatego obliczając największe wspólne dzielniki (NWD) między $N$ a $x-1$ i $x+1$, otrzymamy czynniki $p$ i $q$. Wyznaczenie NWD dwóch liczb jest klasycznie łatwym zadaniem, które można wykonać np. za pomocą algorytmu Euklidesa.

Sprawdź swoją wiedzę

Przeczytaj poniższe pytania, zastanów się nad odpowiedzią, a następnie kliknij trójkąt, aby odkryć rozwiązanie.

Zrozumienie każdego kroku powyższego rozumowania może być niełatwe, dlatego spróbuj prześledzić je na przykładzie. Użyj $N=15$ i $x=11$. Najpierw zweryfikuj, że $x^2 \equiv 1 \text{mod}(N)$ i $x \not\equiv \pm 1 \bmod N$. Następnie sprawdź każdy kolejny krok. Na koniec oblicz $\text{GCD}(11\pm1,15)$ i zweryfikuj, że są to czynniki liczby $15$.

Odpowiedź:

$11^2 = 121$, co wynosi $15*8 + 1$, zatem $11^2\bmod 15 = 1$. $\checkmark$

$11 - 1 = 10$, co nie jest równoważne $0\bmod 15$. $\checkmark$

$11 + 1 = 12$, co nie jest równoważne $0\bmod 15$. $\checkmark$

Wiemy teraz, że $(x+1)(x-1) = l N$ dla pewnej liczby całkowitej $l$. Weryfikujemy to, podstawiając $x$ i $N$: $(12)(10) = l 15$ dla $l = 8$. $\checkmark$

Teraz musimy obliczyć $\text{GCD}(12,15)$ i $\text{GCD}(10,15)$.

$$\begin{aligned} \text{GCD}(12,15) = 3 \\ \text{GCD}(10,15) = 5 \end{aligned}$$

Znaleźliśmy czynniki liczby $15$!

The algorithm

Teraz, gdy wiemy już, jak znalezienie pewnej liczby całkowitej $x$ spełniającej warunek $x^2 \equiv 1\bmod N$ pomaga nam rozłożyć $N$ na czynniki, możemy przejść przez algorytm Shora. Sprowadza się on w zasadzie do znalezienia $x$:

1. Wybierz losową liczbę całkowitą Wybierz losową liczbę całkowitą $a$ taką, że $1 < a < N$.

• Oblicz $\text{GCD}(a, N)$ klasycznie.
  • Jeśli $\text{GCD}(a, N) > 1$, znalazłeś już czynnik. Zakończ.
  • W przeciwnym razie kontynuuj.

2. Znajdź rząd $r$ liczby $a$ modulo $N$ Znajdź najmniejszą dodatnią liczbę całkowitą $r$ spełniającą $a^r \equiv 1 \pmod N$.

3. Sprawdź, czy rząd jest parzysty

• Jeśli $r$ jest nieparzyste, wróć do kroku 1 i wybierz nowe $a$.
• Jeśli $r$ jest parzyste, przejdź do kroku 4.

4. Oblicz $x = a^{r/2} \bmod N$

• Sprawdź, czy $x \not\equiv 1 \pmod N$ oraz $x \not\equiv -1 \pmod N$.
  • Jeśli $x \equiv \pm 1 \pmod N$, wróć do kroku 1 i wybierz nowe $a$.
• W przeciwnym razie oblicz NWD, aby wyłuskać czynniki:

$$p = \text{GCD}(x-1, N), \quad q = \text{GCD}(x+1, N)$$

Będą to nietrywialnie czynniki $N$.

5. W razie potrzeby rozłóż rekurencyjnie

• Jeśli $p$ i/lub $q$ nie są liczbami pierwszymi, zastosuj algorytm rekurencyjnie, aby je w pełni rozłożyć.
• Gdy wszystkie czynniki są pierwsze, rozkład jest kompletny.

Na podstawie tej procedury może nie być oczywiste, dlaczego do wykonania tego zadania potrzebny jest komputer kwantowy. Jest to konieczne, ponieważ krok 2 — znalezienie rzędu $a$ modulo $N$ — jest klasycznie bardzo trudnym problemem. Złożoność rośnie wykładniczo względem $N$. Jednak przy użyciu komputera kwantowego wystarczy skorzystać z Kwantowej Estymacji Fazy (Quantum Phase Estimation), aby go rozwiązać. Krok 4, czyli znalezienie NWD dwóch liczb całkowitych, jest natomiast klasycznie dość prostym zadaniem. Mówimy zatem, że problem rozkładu na czynniki „redukuje się" do problemu wyznaczania rzędu.

The hard part: order finding

Teraz omówimy, jak można użyć komputera kwantowego do wyznaczania rzędu. Najpierw wyjaśnijmy dokładniej, co rozumiemy przez „rząd". Oczywiście podałem już jego matematyczne znaczenie: to pierwsza niezerowa liczba całkowita $r$ taka, że $a^r = 1 \pmod N$. Spróbujmy jednak nieco lepiej wyczuć intuicję kryjącą się za tym pojęciem.

Dla wystarczająco małego $N$ możemy po prostu wyznaczyć rząd, obliczając kolejne potęgi $a$, biorąc modulus $N$ każdego wyniku, a następnie zatrzymując się, gdy znajdziemy potęgę $r$ spełniającą $a^r = 1 \text{mod}(N)$. Tak właśnie zrobiliśmy w powyższym przykładzie z $N=15$. Przyjrzyjmy się wykresom tych potęg modularnych dla kilku przykładowych wartości $a$ i $N$:

Zauważasz coś? To są funkcje periodyczne! A rząd $r$ jest równy okresowi! Wyznaczanie rzędu jest zatem równoważne z wyznaczaniem okresu.

Komputery kwantowe doskonale nadają się do znajdowania okresu funkcji. W tym celu możemy skorzystać z podprocedury algorytmicznej zwanej Kwantową Estymacją Fazy (Quantum Phase Estimation). Omawialiśmy QPE i jej związek z Kwantową Transformatą Fouriera w poprzednim module. Szczegółowe przypomnienie znajdziesz w module QFT lub w lekcji Johna Watrouca dotyczącej Kwantowej Estymacji Fazy w jego kursie algorytmów kwantowych. Teraz omówimy ogólny zarys procedury:

W Kwantowej Estymacji Fazy (QPE) zaczynamy od unitarnego operatora $U$ oraz własnego stanu tego unitarnego $|\psi\rangle$. Następnie używamy QPE do przybliżenia odpowiedniej wartości własnej, która — ponieważ operator jest unitarny — ma postać $e^{2\pi i \theta}$. Znalezienie wartości własnej jest zatem równoważne ze znalezieniem wartości $\theta$ w funkcji periodycznej. Układ (Circuit) wygląda następująco:

gdzie liczba Qubitów sterujących (górne $m$ Qubitów na powyższym rysunku) wyznacza dokładność przybliżenia.

W algorytmie Shora stosujemy QPE na unitarnym operatorze $M_a$:

$M_a|y\rangle \equiv |ay \mod N \rangle .$

Tutaj $|y\rangle$ oznacza stan z bazy obliczeniowej wieloQubitowego rejestru, gdzie binarna wartość Qubitów odpowiada liczbie całkowitej $y$. Na przykład, jeśli $N=15$ i $y = 2$, to $|y\rangle$ jest reprezentowane przez czteroQubitowy stan bazowy $|0010\rangle$, ponieważ do zakodowania liczb do 15 potrzebne są cztery Qubity. (Jeśli to pojęcie jest Ci nieznane, zajrzyj do wstępnego modułu Qiskit w Classrooms, gdzie znajdziesz przypomnienie o binarnym kodowaniu stanów kwantowych.)

Musimy teraz znaleźć stan własny tego operatora unitarnego. Jeśli zaczęlibyśmy w stanie $|1\rangle$, każde kolejne zastosowanie $U$ przemnożyłoby stan naszego rejestru przez $a \pmod N$, a po $r$ zastosowaniach wrócimy do stanu $|1\rangle$. Na przykład dla $a = 3$ i $N = 35$:

$$\begin{aligned} M_3|1\rangle &= |3\rangle & \\ M_3^2|1\rangle &= |9\rangle \\ M_3^3|1\rangle &= |27\rangle \\ & \vdots \\ M_3^{(r-1)}|1\rangle &= |12\rangle \\ M_3^r|1\rangle &= |1\rangle \end{aligned}$$

Zatem superpozycje stanów w tym cyklu ($|\psi_j\rangle$) postaci:

$|\psi_j\rangle = \tfrac{1}{\sqrt{r}}\sum_{k=0}^{r-1}{e^{\frac{2 \pi i j k}{r}} |a^k \rangle}$

są stanami własnymi $M_a$. (Istnieje więcej stanów własnych niż tylko te. Ale interesują nas wyłącznie stany powyższej postaci.)

Check your understanding

Przeczytaj poniższe pytanie(-a), zastanów się nad odpowiedzią, a następnie kliknij trójkąt, aby ujawnić rozwiązanie.

Znajdź stan własny operatora unitarnego odpowiadającego $a=2$ i $N = 15$.

Odpowiedź:

$$\begin{aligned} M_2|1\rangle &= |2\rangle & \\ M_2^2|1\rangle &= |4\rangle \\ M_2^3|1\rangle &= |8\rangle \\ M_2^4|1\rangle &= |1\rangle \\ \end{aligned}$$

Zatem rząd $r=4$. Interesujące nas stany własne będą równą superpozycją wszystkich stanów, przez które przechodził powyższy cykl, z różnymi fazami:

$$\begin{aligned} |\psi_0\rangle &= \frac{1}{2}(|1\rangle+|2\rangle+|4\rangle+|8\rangle) \\ |\psi_1\rangle &= \frac{1}{2}(e^{2 \pi i \frac{0}{4}}|1\rangle+e^{2 \pi i \frac{1}{4}}|2\rangle+e^{2 \pi i \frac{2}{4}}|4\rangle+e^{2 \pi i \frac{3}{4}}|8\rangle) \\ &= \frac{1}{2}(|1\rangle+i|2\rangle-|4\rangle-i|8\rangle) \\ |\psi_2\rangle &= \frac{1}{2}(e^{2 \pi i \frac{0}{4}}|1\rangle+e^{2 \pi i \frac{2}{4}}|2\rangle+e^{2 \pi i \frac{4}{4}}|4\rangle+e^{2 \pi i \frac{6}{4}}|8\rangle) \\ &= \frac{1}{2}(|1\rangle-|2\rangle+|4\rangle-|8\rangle) \\ |\psi_3\rangle &= \frac{1}{2}(e^{2 \pi i \frac{0}{4}}|1\rangle+e^{2 \pi i \frac{3}{4}}|2\rangle+e^{2 \pi i \frac{6}{4}}|4\rangle+e^{2 \pi i \frac{9}{4}}|8\rangle) \\ &= \frac{1}{2}(|1\rangle-i|2\rangle-|4\rangle+i|8\rangle) \\ \end{aligned}$$

Powiedzmy, że bylibyśmy w stanie zainicjować stan naszego Qubita w jednym z tych stanów własnych (uwaga — nie jesteśmy. Przynajmniej nie w prosty sposób. Za chwilę wyjaśnimy, dlaczego i co można zamiast tego zrobić). Wtedy moglibyśmy użyć QPE do oszacowania odpowiedniej wartości własnej, $\omega_j = e^{2 \pi i \theta_j}$, gdzie $\theta_j = \frac{j}{r}$. Następnie będziemy mogli wyznaczyć rząd $r$ z prostego równania:

$r = \frac{j}{\theta_j}.$

Pamiętaj jednak, że QPE jedynie przybliża $\theta_j$ — nie podaje dokładnej wartości. Potrzebujemy, aby przybliżenie było wystarczająco dobre, by odróżnić $r$ od $r+1$. Im więcej Qubitów sterujących $m$, tym lepsze będzie przybliżenie. W zadaniach na końcu lekcji zostaniesz poproszony(-a) o wyznaczenie minimalnej wartości $m$ potrzebnej do rozłożenia liczby $N$ na czynniki.

Musimy teraz rozwiązać pewien problem. Całe powyższe wyjaśnienie, jak znaleźć $r$, zaczyna się od przygotowania stanu własnego $|\psi_j\rangle = \tfrac{1}{\sqrt{r}}\sum_{k=0}^{r-1}{e^{\frac{2 \pi i j k}{r}} |a^k \rangle}$. Nie wiemy jednak, jak to zrobić bez uprzedniej znajomości $r$. Rozumowanie jest kołowe. Potrzebujemy sposobu na oszacowanie wartości własnej bez inicjowania stanu własnego.

Zamiast zaczynać od stanu własnego $M_a$, możemy przygotować stan początkowy jako $n$-Qubitowy stan odpowiadający $|1\rangle$ w zapisie binarnym (tzn. $|000...01\rangle$). Choć ten stan sam w sobie oczywiście nie jest stanem własnym $M_a$, jest superpozycją wszystkich stanów własnych $|\psi_k\rangle$:

$|1\rangle = \frac{1}{\sqrt{r}} \sum\limits_{k=0}^{r-1}{|\psi_k\rangle}$

Check your understanding

Przeczytaj poniższe pytanie(-a), zastanów się nad odpowiedzią, a następnie kliknij trójkąt, aby ujawnić rozwiązanie.

Sprawdź, że $|1\rangle$ jest równoważne superpozycji stanów własnych, które znalazłeś(-aś) dla $N=15$ i $a=2$ w poprzednim pytaniu sprawdzającym.

Odpowiedź:

Cztery stany własne to:

$$\begin{aligned} |\psi_0\rangle &= \frac{1}{2}(|1\rangle+|2\rangle+|4\rangle+|8\rangle) \\ |\psi_1\rangle &= \frac{1}{2}(|1\rangle+i|2\rangle-|4\rangle-i|8\rangle) \\ |\psi_2\rangle &= \frac{1}{2}(|1\rangle-|2\rangle+|4\rangle-|8\rangle) \\ |\psi_3\rangle &= \frac{1}{2}(|1\rangle-i|2\rangle-|4\rangle+i|8\rangle) \\ \end{aligned}$$

Zatem:

$$\begin{aligned} \frac{1}{\sqrt{r}} \sum\limits_{k=0}^{r-1}{|\psi_k\rangle} &= \frac{1}{2}(|\psi_0\rangle + |\psi_1\rangle + |\psi_2\rangle + |\psi_3\rangle ) \\ &= \frac{1}{4}(|1\rangle+|2\rangle+|4\rangle+|8\rangle+|1\rangle+i|2\rangle-|4\rangle-i|8\rangle+|1\rangle-|2\rangle+|4\rangle-|8\rangle + |1\rangle-i|2\rangle-|4\rangle+i|8\rangle) \\ &= \frac{1}{4}(4|1\rangle) = |1\rangle \end{aligned}$$

Jak to pozwala nam znaleźć rząd $r$? Ponieważ stan początkowy jest superpozycją wszystkich stanów własnych powyższej postaci, algorytm QPE jednocześnie szacuje każde z $\theta_k$ odpowiadających tym stanom własnym. Pomiar $m$ Qubitów sterujących na końcu da zatem przybliżenie wartości $k/r$, gdzie $k \in \{0,1,2,...,r-1\}$ jest losowo wybraną wartością własną. Jeśli kilkukrotnie powtórzymy ten obwód i uzyskamy kilka próbek z różnymi wartościami $k$, szybko będziemy mogli wyznaczyć $r$.

Implementacja w Qiskit

Jak wspomnieliśmy wcześniej, nasze urządzenia sprzętowe nie są jeszcze w stanie faktoryzować ogromnych liczb, takich jak RSA1024. Sfaktoryzujemy więc małą liczbę, aby zademonstrować, jak działa algorytm. W tym demo posłużymy się uproszczoną wersją kodu przedstawionego w samouczku dotyczącym algorytmu Shora. Jeśli chcesz poznać więcej szczegółów, odwiedź tamten samouczek.

Uruchomimy algorytm przy użyciu naszego standardowego frameworku do rozwiązywania problemów kwantowych, zwanego frameworkiem wzorców Qiskit (ang. Qiskit patterns framework). Składa się on z czterech kroków:

1. Odwzorowanie problemu na obwód kwantowy
2. Optymalizacja Circuit do uruchomienia na sprzęcie kwantowym
3. Wykonanie Circuit na komputerze kwantowym
4. Przetwarzanie końcowe pomiarów

1. Odwzorowanie

Sfaktoryzujemy $N=15$, wybierając $a=2$ jako naszą liczbę pierwszą względem $N$.

Najpierw musimy zbudować Circuit implementujący unitarny operator mnożenia modularnego $M_a$. To faktycznie najtrudniejsza część całej implementacji i może być bardzo kosztowna obliczeniowo, w zależności od sposobu realizacji. W tym celu trochę „oszukamy": wiemy, że startujemy ze stanu $|1\rangle$, a z wcześniejszego pytania sprawdzającego wiemy, że:

$$\begin{aligned} M_2|1\rangle &= |2\rangle & \\ M_2|2\rangle &= |4\rangle \\ M_2|4\rangle &= |8\rangle \\ M_2|8\rangle &= |1\rangle \\ \end{aligned}$$

Zbudujemy więc operator unitarny, który wykonuje prawidłowe operacje na tych czterech stanach, a wszystkie pozostałe stany pozostawia bez zmian. To jest właśnie „oszustwo", bo korzystamy z naszej wiedzy o rzędzie $2\bmod 15$, aby uprościć operator unitarny. Gdybyśmy faktycznie próbowali faktoryzować liczbę, której czynniki są nam nieznane, nie moglibyśmy tego zrobić.

Sprawdź swoje rozumienie

Przeczytaj poniższe pytanie/pytania, zastanów się nad odpowiedzią, a następnie kliknij trójkąt, aby ujawnić rozwiązanie.

Znając działanie operatora $M_2$ na powyższe stany, skonstruuj ten operator z szeregu bramek SWAP, które zamieniają stany dwóch Qubitów. (Wskazówka: pomaga zapisanie każdego stanu $|i\rangle$ w reprezentacji binarnej.)

Odpowiedź:

Przepiszmy działanie $M_2$ na stanach w reprezentacji binarnej:

$$\begin{aligned} M_2|0001\rangle &= |0010\rangle \\ M_2|0010\rangle &= |0100\rangle \\ M_2|0100\rangle &= |1000\rangle \\ M_2|1000\rangle &= |0001\rangle \\ \end{aligned}$$

Każde z tych działań można zrealizować za pomocą prostej bramki SWAP. $M_2|0001\rangle$ uzyskujemy, zamieniając stany Qubitu $0$ i $1$. $M_2|0010\rangle$ uzyskujemy, zamieniając stany Qubitu $1$ i $2$. I tak dalej. Możemy więc rozłożyć macierz $M_2$ na następującą serię bramek SWAP:

$$M_2 = SWAP(0,1)SWAP(1,2)SWAP(2,3)$$

Pamiętając, że operatory działają od prawej do lewej, sprawdźmy, że daje to oczekiwany efekt na każdym ze stanów:

$$\begin{aligned} M_2|0001\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|0001\rangle \\ &= SWAP(0,1)SWAP(1,2)|0001\rangle \\ &= SWAP(0,1)|0001\rangle \\ &=|0010\rangle \checkmark \\ M_2|0010\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|0010\rangle \\ &= SWAP(0,1)SWAP(1,2)|0010\rangle \\ &= SWAP(0,1)|0100\rangle \\ &=|0100\rangle \checkmark \\ M_2|0100\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|0100\rangle \\ &= SWAP(0,1)SWAP(1,2)|1000\rangle \\ &= SWAP(0,1)|1000\rangle \\ &=|1000\rangle \checkmark \\ M_2|1000\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|1000\rangle \\ &= SWAP(0,1)SWAP(1,2)|0100\rangle \\ &= SWAP(0,1)|0010\rangle \\ &=|0001\rangle \checkmark \\ \end{aligned}$$

Możemy teraz zakodować Circuit równoważny temu operatorowi w Qiskit.

Najpierw importujemy niezbędne pakiety:

# Import necessary packages

import numpy as np
from fractions import Fraction
from math import floor, gcd, log

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFTGate
from qiskit.transpiler import generate_preset_pass_manager
from qiskit.visualization import plot_histogram

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

Następnie tworzymy operator $M_2$:

def M2mod15():
    """
    M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)

    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)

    U = U.to_gate()
    U.name = f"M_{b}"

    return U

# Get the M2 operator
M2 = M2mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

Algorytm QPE używa bramki kontrolowanej-$U$. Skoro mamy już Circuit z $M_2$, musimy zbudować Circuit kontrolowany-$M_2$:

def controlled_M2mod15():
    """
    Controlled M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)

    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)

    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()

    return c_U

# Get the controlled-M2 operator
controlled_M2 = controlled_M2mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M2, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

Mamy teraz naszą bramkę kontrolowaną-$U$. Aby jednak uruchomić algorytm Quantum Phase Estimation, będziemy potrzebować kontrolowanego-$U^2$, kontrolowanego-$U^4$, aż do kontrolowanego-$U^{2^{m-1}}$, gdzie $m$ to liczba Qubitów użytych do estymacji fazy. Im więcej Qubitów, tym dokładniejsza będzie estymacja fazy. Do naszej procedury estymacji fazy użyjemy $m=8$ Qubitów kontrolnych. Potrzebujemy więc:

$$M_{a^{2^k}}|y\rangle \equiv |a^{2^k} y \bmod N \rangle$$

gdzie indeks $k$, dla $0 \le k \le m-1 = 7$, odpowiada Qubitowi kontrolnemu. Obliczmy teraz $a^{2^k}\bmod N$ dla każdej wartości $k$:

def a2kmodN(a, k, N):
    """Compute a^{2^k} (mod N) by repeated squaring"""
    for _ in range(k):
        a = int(np.mod(a**2, N))
    return a

k_list = range(8)
b_list = [a2kmodN(2, k, 15) for k in k_list]

print(b_list)

[2, 4, 1, 1, 1, 1, 1, 1]

Ponieważ $a^{2^k} \bmod N = 1$ dla $k \ge 2$, wszystkie odpowiadające im operatory ($M_8$ i wyższe) są równoważne operatorowi tożsamości. Musimy więc skonstruować tylko jedną dodatkową macierz, $M_4.$

Uwaga: To uproszczenie działa tutaj tylko dlatego, że rząd $2 \bmod 15$ wynosi $4$. Gdy $k=2$ (czyli $2^k = 4$), każda kolejna potęga operatora jest operatorem tożsamości. Ogólnie rzecz biorąc, dla większych $N$ lub innych wyborów $a$ nie można pominąć konstruowania wyższych potęg. To właśnie jeden z powodów, dla których jest to uważane za przykład zabawkowy: małe liczby pozwalają na skróty, które nie działałyby w przypadkach o większej skali.

def M4mod15():
    """
    M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)

    U.swap(1, 3)
    U.swap(0, 2)

    U = U.to_gate()
    U.name = f"M_{b}"

    return U

# Get the M4 operator
M4 = M4mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M4, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

I podobnie jak poprzednio, tworzymy z niego operator kontrolowany-$M_4$:

def controlled_M4mod15():
    """
    Controlled M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)

    U.swap(1, 3)
    U.swap(0, 2)

    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()

    return c_U

# Get the controlled-M4 operator
controlled_M4 = controlled_M4mod15()

# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M4, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

Teraz możemy złożyć to wszystko razem, aby znaleźć rząd $2\bmod 15$ za pomocą Circuit kwantowego, korzystając z estymacji fazy:

# Order finding problem for N = 15 with a = 2
N = 15
a = 2

# Number of qubits
num_target = floor(log(N - 1, 2)) + 1  # for modular exponentiation operators
num_control = 2 * num_target  # for enough precision of estimation

# List of M_b operators in order
k_list = range(num_control)
b_list = [a2kmodN(2, k, 15) for k in k_list]

# Initialize the circuit
control = QuantumRegister(num_control, name="C")
target = QuantumRegister(num_target, name="T")
output = ClassicalRegister(num_control, name="out")
circuit = QuantumCircuit(control, target, output)

# Initialize the target register to the state |1>
circuit.x(num_control)

# Add the Hadamard gates and controlled versions of the
# multiplication gates
for k, qubit in enumerate(control):
    circuit.h(k)
    b = b_list[k]
    if b == 2:
        circuit.compose(
            M2mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    elif b == 4:
        circuit.compose(
            M4mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    else:
        continue  # M1 is the identity operator

# Apply the inverse QFT to the control register
circuit.compose(QFTGate(num_control).inverse(), qubits=control, inplace=True)

# Measure the control register
circuit.measure(control, output)

circuit.draw("mpl", fold=-1)

2. Optymalizacja

Teraz, gdy zmapowaliśmy nasz obwód, kolejnym krokiem jest jego optymalizacja pod kątem uruchomienia na konkretnym komputerze kwantowym. Najpierw musimy wczytać Backend.

service = QiskitRuntimeService()

backend = service.backend("ibm_marrakesh")

Jeśli nie masz dostępnego czasu na swoim koncie lub z jakiegokolwiek powodu chcesz użyć symulatora, możesz uruchomić poniższą komórkę, aby skonfigurować symulator naśladujący wybrane przez nas urządzenie kwantowe:

pm = generate_preset_pass_manager(optimization_level=2, backend=backend)

transpiled_circuit = pm.run(circuit)

print(f"2q-depth: {transpiled_circuit.depth(lambda x: x.operation.num_qubits==2)}")
print(f"2q-size: {transpiled_circuit.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circuit.count_ops()}")
transpiled_circuit.draw(output="mpl", fold=-1, style="clifford", idle_wires=False)

2q-depth: 188
2q-size: 281
Operator counts: OrderedDict({'sx': 548, 'rz': 380, 'cz': 281, 'measure': 8, 'x': 6})

3. Wykonanie

# Sampler primitive to obtain the probability distribution
sampler = Sampler(backend)

# Turn on dynamical decoupling with sequence XpXm
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XpXm"
# Enable gate twirling
sampler.options.twirling.enable_gates = True

pub = transpiled_circuit
job = sampler.run([pub], shots=1024)

result = job.result()[0]
counts = result.data["out"].get_counts()

plot_histogram(counts, figsize=(35, 5))

Widzimy cztery wyraźne piki przy 00000000, 01000000, 10000000 i 11000000, z pewną liczbą wyników w innych bitstrings spowodowanych szumem komputera kwantowego. Zignorujemy je i zachowamy jedynie dominującą czwórkę, narzucając próg: tylko zliczenia powyżej tego progu są uznawane za rzeczywisty sygnał ponad szumem.

# Dictionary of bitstrings and their counts to keep
counts_keep = {}
# Threshold to filter
threshold = np.max(list(counts.values())) / 2

for key, value in counts.items():
    if value > threshold:
        counts_keep[key] = value

print(counts_keep)

4. Post-processing

W algorytmie Shora znaczna część algorytmu jest wykonywana klasycznie. Dlatego resztę umieścimy w kroku „post-processingu", po uzyskaniu pomiarów z komputera kwantowego. Każdy z powyższych pomiarów można przekształcić w liczby całkowite, które po podzieleniu przez $2^m$ stanowią nasze przybliżenia $\frac{k}{r}$, gdzie $k$ jest za każdym razem losowe.

a = 2
N = 15

FACTOR_FOUND = False
num_attempt = 0

while not FACTOR_FOUND:
    print(f"\nATTEMPT {num_attempt}:")
    # Here, we get the bitstring by iterating over outcomes
    # of a previous hardware run with multiple shots.
    # Instead, we can also perform a single-shot measurement
    # here in the loop.
    bitstring = list(counts_keep.keys())[num_attempt]
    num_attempt += 1
    # Find the phase from measurement
    decimal = int(bitstring, 2)
    phase = decimal / (2**num_control)  # phase = k / r
    print(f"Phase: theta = {phase}")

    # Guess the order from phase
    frac = Fraction(phase).limit_denominator(N)
    r = frac.denominator  # order = r
    print(f"Order of {a} modulo {N} estimated as: r = {r}")

    if phase != 0:
        # Guesses for factors are gcd(a^{r / 2} ± 1, 15)
        if r % 2 == 0:
            x = pow(a, r // 2, N) - 1
            d = gcd(x, N)
            if d > 1:
                FACTOR_FOUND = True
                print(f"*** Non-trivial factor found: {x} ***")

ATTEMPT 0:
Phase: theta = 0.0
Order of 2 modulo 15 estimated as: r = 1

ATTEMPT 1:
Phase: theta = 0.75
Order of 2 modulo 15 estimated as: r = 4
*** Non-trivial factor found: 3 ***

Wnioski

Po przejściu przez ten moduł być może ogarnął cię nowy podziw dla geniuszu Petera Shora, który wpadł na tak sprytny algorytm. Mamy jednak nadzieję, że osiągnąłeś też nowy poziom zrozumienia jego zwodniczej prostoty. Choć algorytm może wydawać się imponująco (lub przytłaczająco) złożony, jeśli rozłożysz go na poszczególne kroki i przejdziesz przez nie powoli, ty również będziesz w stanie uruchomić algorytm Shora.

Choć jesteśmy daleko od wykorzystania tego algorytmu do faktoryzacji liczb takich jak RSA1024, nasze komputery kwantowe poprawiają się z dnia na dzień, a gdy zostanie osiągnięty próg zwany odpornością na błędy, takie algorytmy będą wkrótce gotowe do użycia. To ekscytujący czas, by uczyć się o obliczeniach kwantowych!

Zadania

Kluczowe pojęcia:

• Nowoczesne systemy kryptograficzne opierają się na klasycznej trudności faktoryzacji dużych liczb całkowitych.
• Arytmetyka modularna — w tym struktury $\mathbb{Z}_N$ i $\mathbb{Z}_N^*$ — stanowi matematyczny fundament algorytmu Shora.
• Problem faktoryzacji liczby całkowitej $N$ można sprowadzić do problemu znalezienia rzędu liczby modulo $N$.
• Kwantowe wyznaczanie rzędu wykorzystuje techniki kwantowej estymacji fazy do wyznaczenia okresu funkcji $a^x \mod N$.
• Algorytm Shora składa się z klasyczno-kwantowego przepływu pracy, który wybiera podstawę, wykonuje kwantowe wyznaczanie rzędu, a następnie klasycznie oblicza czynniki z uzyskanego wyniku.

Prawda/Fałsz:

1. P/F Wydajność algorytmu Shora zagraża bezpieczeństwu szyfrowania RSA.
2. P/F Algorytm Shora można efektywnie wykonać na dowolnym współczesnym komputerze kwantowym.
3. P/F Algorytm Shora wykorzystuje kwantową estymację fazy (QPE) jako kluczową procedurę pomocniczą.
4. P/F Klasyczna część algorytmu Shora polega na obliczeniu największego wspólnego dzielnika (NWD).
5. P/F Algorytm Shora działa wyłącznie dla faktoryzacji liczb parzystych.
6. P/F Pomyślne uruchomienie algorytmu Shora zawsze gwarantuje znalezienie właściwych czynników.

Krótkie odpowiedzi:

1. Dlaczego algorytm Shora jest uważany za potencjalne przyszłe zagrożenie dla szyfrowania RSA?
2. Dlaczego wyznaczenie okresu, czyli rzędu, modularnej funkcji wykładniczej pomaga w faktoryzacji liczby w algorytmie Shora?

Zadania dla zaawansowanych:

1. Ile kubitów sterujących $m$ potrzebujemy dla danej liczby $N$, którą chcemy faktoryzować, aby uzyskać w QPE precyzję wystarczającą do znalezienia właściwej wartości rzędu $r$?

2. Stosując opisaną tutaj procedurę faktoryzacji liczby 15, spróbuj teraz sfaktoryzować liczbę 21.