אלגוריתם שור

הערכת שימוש: שלוש שניות על מעבד Eagle r3 (הערה: זוהי הערכה בלבד. זמן הריצה שלך עשוי להשתנות.)

אלגוריתם שור, שפותח על ידי פיטר שור בשנת 1994, הוא אלגוריתם קוונטי פורץ דרך לפירוק גורמים של מספרים שלמים בזמן פולינומי. משמעותו טמונה ביכולתו לפרק גורמים של מספרים שלמים גדולים במהירות אקספוננציאלית לעומת כל אלגוריתם קלאסי ידוע, תוך איום על אבטחת מערכות קריפטוגרפיות בשימוש נרחב כמו RSA, המסתמכות על הקושי של פירוק גורמים של מספרים גדולים. באמצעות פתרון יעיל של בעיה זו על מחשב קוונטי חזק מספיק, אלגוריתם שור עשוי לחולל מהפכה בתחומים כגון קריפטוגרפיה, אבטחת סייבר ומתמטיקה חישובית, תוך הדגשת הכוח הטרנספורמטיבי של חישוב קוונטי.

מדריך זה מתמקד בהדגמת אלגוריתם שור באמצעות פירוק הגורמים של 15 במחשב קוונטי.

ראשית, אנחנו מגדירים את בעיית מציאת הסדר ובונים מעגלים תואמים מפרוטוקול הערכת הפאזה הקוונטי. לאחר מכן, אנחנו מריצים את מעגלי מציאת הסדר על חומרה אמיתית תוך שימוש במעגלים בעומק הקצר ביותר שאנו יכולים לתרגם. החלק האחרון משלים את אלגוריתם שור על ידי חיבור בעיית מציאת הסדר לפירוק גורמים של מספרים שלמים.

אנחנו מסיימים את המדריך בדיון על הדגמות אחרות של אלגוריתם שור על חומרה אמיתית, תוך התמקדות הן ביישומים גנריים והן באלו המותאמים לפירוק גורמים של מספרים שלמים ספציפיים כגון 15 ו-21. הערה: מדריך זה מתמקד יותר ביישום ובהדגמה של המעגלים הנוגעים לאלגוריתם שור. למשאב חינוכי מעמיק על החומר, ראה בקורס יסודות האלגוריתמים הקוונטיים מאת ד"ר ג'ון ווטרוס, ומאמרים בסעיף הפניות.

דרישות

לפני תחילת מדריך זה, ודא שיש לך את הדברים הבאים מותקנים:

• Qiskit SDK v2.0 ומעלה, עם תמיכה ב-ויזואליזציה
• Qiskit Runtime v0.40 ומעלה (pip install qiskit-ibm-runtime)

הכנה

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

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

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFT, UnitaryGate
from qiskit.transpiler import CouplingMap, generate_preset_pass_manager
from qiskit.visualization import plot_histogram

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

שלב 1: מיפוי קלט קלאסי לבעיה קוונטית

רקע

אלגוריתם שור לפירוק גורמים של מספרים שלמים משתמש בבעיה מתווכת המכונה בעיית מציאת הסדר. בחלק זה, אנחנו מדגימים כיצד לפתור את בעיית מציאת הסדר באמצעות הערכת פאזה קוונטית.

בעיית הערכת פאזה

בבעיית הערכת הפאזה, ניתן לנו מצב קוונטי $\ket{\psi}$ של $n$ קיוביטים, יחד עם מעגל קוונטי יוניטרי הפועל על $n$ קיוביטים. ניתנת לנו ההבטחה ש-$\ket{\psi}$ הוא וקטור עצמי של המטריצה היוניטרית $U$ המתארת את פעולת המעגל, והמטרה שלנו היא לחשב או לקרב את ערך העצמי $\lambda = e^{2 \pi i \theta}$ שאליו $\ket{\psi}$ מתאים. במילים אחרות, המעגל צריך להוציא קירוב למספר $\theta \in [0, 1)$ המקיים $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ מטרת מעגל הערכת הפאזה היא לקרב את $\theta$ ב-$m$ ביטים. מבחינה מתמטית, אנחנו רוצים למצוא $y$ כך ש-$\theta \approx y / 2^m$, כאשר $y \in {0, 1, 2, \dots, 2^{m-1}}$. התמונה הבאה מציגה את המעגל הקוונטי שמעריך את $y$ ב-$m$ ביטים על ידי ביצוע מדידה על $m$ קיוביטים. במעגל שלמעלה, $m$ הקיוביטים העליונים מאותחלים במצב $\ket{0^m}$, ו-$n$ הקיוביטים התחתונים מאותחלים ב-$\ket{\psi}$, שמובטח שהוא וקטור עצמי של $U$. המרכיב הראשון במעגל הערכת הפאזה הם הפעולות היוניטריות המבוקרות האחראיות לביצוע בעיטת פאזה לקיוביט הבקרה המתאים שלהן. יוניטריים מבוקרים אלו מועלים בחזקה בהתאם למיקום של קיוביט הבקרה, החל מהביט המשמעותי ביותר ועד לביט הפחות משמעותי. מכיוון ש-$\ket{\psi}$ הוא וקטור עצמי של $U$, המצב של $n$ הקיוביטים התחתונים אינו מושפע מפעולה זו, אך מידע הפאזה של ערך העצמי מתפשט ל-$m$ הקיוביטים העליונים. מסתבר שלאחר פעולת בעיטת הפאזה באמצעות יוניטריים מבוקרים, כל המצבים האפשריים של $m$ הקיוביטים העליונים אורתונורמליים זה לזה עבור כל וקטור עצמי $\ket{\psi}$ של היוניטרי $U$. לכן, מצבים אלו ניתנים להבחנה באופן מושלם, ואנו יכולים לסובב את הבסיס שהם יוצרים בחזרה לבסיס החישובי כדי לבצע מדידה. ניתוח מתמטי מראה שמטריצת סיבוב זו מתאימה לטרנספורמציה של פורייה הקוונטית (QFT) ההופכית במרחב הילברט בעל מימד $2^m$. האינטואיציה מאחורי זה היא שהמבנה המחזורי של אופרטורי ההעלאה בחזקה המודולרית מקודד במצב הקוונטי, וה-QFT ממיר מחזוריות זו לפסגות ניתנות למדידה בתחום התדר.

להבנה מעמיקה יותר של מדוע נעשה שימוש במעגל QFT באלגוריתם שור, אנחנו מפנים את הקורא לקורס יסודות האלגוריתמים הקוונטיים. כעת אנחנו מוכנים להשתמש במעגל הערכת הפאזה למציאת סדר.

בעיית מציאת הסדר

כדי להגדיר את בעיית מציאת הסדר, אנחנו מתחילים עם כמה מושגי תורת המספרים. ראשית, עבור כל מספר שלם חיובי נתון $N$, נגדיר את הקבוצה $\mathbb{Z}_N$ כ-$\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ כל הפעולות האריתמטיות ב-$\mathbb{Z}_N$ מבוצעות מודולו $N$. בפרט, כל האלמנטים $a \in \mathbb{Z}_n$ שהם זרים זה לזה עם $N$ הם מיוחדים ומהווים את $\mathbb{Z}^*_N$ כ-$\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ עבור אלמנט $a \in \mathbb{Z}^*_N$, המספר השלם החיובי הקטן ביותר $r$ כך ש-$a^r \equiv 1 \; (\mathrm{mod} \; N)$ מוגדר כ-סדר של $a$ מודולו $N$. כפי שנראה מאוחר יותר, מציאת הסדר של $a \in \mathbb{Z}^*_N$ תאפשר לנו לפרק גורמים של $N$. כדי לבנות את מעגל מציאת הסדר ממעגל הערכת הפאזה, אנחנו זקוקים לשני שיקולים. ראשית, אנחנו צריכים להגדיר את היוניטרי $U$ שיאפשר לנו למצוא את הסדר $r$, ושנית, אנחנו צריכים להגדיר וקטור עצמי $\ket{\psi}$ של $U$ כדי להכין את המצב ההתחלתי של מעגל הערכת הפאזה.

כדי לחבר את בעיית מציאת הסדר להערכת פאזה, אנחנו שוקלים את הפעולה המוגדרת על מערכת שמצביה הקלאסיים מתאימים ל-$\mathbb{Z}_N$, שבה אנחנו מכפילים באלמנט קבוע $a \in \mathbb{Z}^*_N$. בפרט, אנחנו מגדירים את אופרטור הכפל $M_a$ כך ש-$M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ עבור כל $x \in \mathbb{Z}_N$. שימו לב שזה משתמע שאנו לוקחים את המכפלה מודולו $N$ בתוך ה-ket בצד ימין של המשוואה. ניתוח מתמטי מראה ש-$M_a$ הוא אופרטור יוניטרי. יתרה מכך, מסתבר ש-$M_a$ יש זוגות וקטור עצמי וערך עצמי המאפשרים לנו לחבר את הסדר $r$ של $a$ לבעיית הערכת הפאזה. באופן ספציפי, עבור כל בחירה של $j \in \{0, \dots, r-1\}$, יש לנו $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ הוא וקטור עצמי של $M_a$ שערך העצמי המתאים לו הוא $\omega^{j}_{r}$, כאשר $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$ על ידי תצפית, אנחנו רואים שזוג וקטור עצמי/ערך עצמי נוח הוא המצב $\ket{\psi_1}$ עם $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. לכן, אם היינו יכולים למצוא את הוקטור העצמי $\ket{\psi_1}$, היינו יכולים להעריך את הפאזה $\theta=1/r$ עם המעגל הקוונטי שלנו ולכן לקבל הערכה של הסדר $r$. עם זאת, זה לא קל לעשות זאת, ואנחנו צריכים לשקול אלטרנטיבה.

בואו נבחן מה יהיה תוצאת המעגל אם נכין את המצב החישובי $\ket{1}$ כמצב ההתחלתי. זה לא מצב עצמי של $M_a$, אבל הוא הסופרפוזיציה האחידה של מצבים עצמיים שתיארנו לעיל. במילים אחרות, היחס הבא מתקיים. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ המשמעות של המשוואה לעיל היא שאם נקבע את המצב ההתחלתי ל-$\ket{1}$, נקבל בדיוק את אותה תוצאת מדידה כאילו בחרנו $k \in \{ 0, \dots, r-1\}$ באופן אחיד אקראי והשתמשנו ב-$\ket{\psi_k}$ כוקטור עצמי במעגל הערכת הפאזה. במילים אחרות, מדידה של $m$ הקיוביטים העליונים מניבה קירוב $y / 2^m$ לערך $k / r$ כאשר $k \in \{ 0, \dots, r-1\}$ נבחר באופן אחיד אקראי. זה מאפשר לנו ללמוד את $r$ עם רמת ביטחון גבוהה לאחר מספר הרצות עצמאיות, שזו הייתה המטרה שלנו.

אופרטורי העלאה בחזקה מודולריים

עד כה, קישרנו את בעיית הערכת הפאזה לבעיית מציאת הסדר על ידי הגדרת $U = M_a$ ו-$\ket{\psi} = \ket{1}$ במעגל הקוונטי שלנו. לכן, המרכיב האחרון שנותר הוא למצוא דרך יעילה להגדיר חזקות מודולריות של $M_a$ כ-$M_a^k$ עבור $k = 1, 2, 4, \dots, 2^{m-1}$. כדי לבצע חישוב זה, אנחנו מוצאים שעבור כל חזקה $k$ שאנו בוחרים, אנחנו יכולים ליצור מעגל עבור $M_a^k$ לא על ידי איטרציה $k$ פעמים של המעגל עבור $M_a$, אלא על ידי חישוב $b = a^k \; \mathrm{mod} \; N$ ולאחר מכן שימוש במעגל עבור $M_b$. מכיוון שאנו זקוקים רק לחזקות שהן בעצמן חזקות של 2, אנחנו יכולים לעשות זאת באופן יעיל קלאסית באמצעות העלאה בריבוע איטרטיבית.

שלב 2: אופטימיזציה של הבעיה לביצוע על חומרה קוונטית

דוגמה ספציפית עם $N = 15$ ו-$a=2$

אנחנו יכולים לעצור כאן כדי לדון בדוגמה ספציפית ולבנות את מעגל מציאת הסדר עבור $N=15$. שימו לב ש-$a \in \mathbb{Z}_N^*$ לא-טריוויאליים אפשריים עבור $N=15$ הם $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. בדוגמה זו, אנחנו בוחרים $a=2$. אנחנו נבנה את האופרטור $M_2$ ואת אופרטורי ההעלאה בחזקה המודולריים $M_2^k$. פעולת $M_2$ על מצבי הבסיס החישובי היא כדלקמן. $M_2 \ket{0} = \ket{0} \quad M_2 \ket{5} = \ket{10} \quad M_2 \ket{10} = \ket{5}$ $M_2 \ket{1} = \ket{2} \quad M_2 \ket{6} = \ket{12} \quad M_2 \ket{11} = \ket{7}$ $M_2 \ket{2} = \ket{4} \quad M_2 \ket{7} = \ket{14} \quad M_2 \ket{12} = \ket{9}$ $M_2 \ket{3} = \ket{6} \quad M_2 \ket{8} = \ket{1} \quad M_2 \ket{13} = \ket{11}$ $M_2 \ket{4} = \ket{8} \quad M_2 \ket{9} = \ket{3} \quad M_2 \ket{14} = \ket{13}$ על ידי תצפית, אנחנו יכולים לראות שמצבי הבסיס מעורבבים, כך שיש לנו מטריצת תמורה. אנחנו יכולים לבנות פעולה זו על ארבעה קיוביטים עם שערי swap. להלן, אנחנו בונים את $M_2$ ואת הפעולות controlled-$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)

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)

שערים הפועלים על יותר משני קיוביטים יפורקו עוד יותר לשערים של שני קיוביטים.

circ.decompose(reps=2).draw(output="mpl", fold=-1)

כעת אנחנו צריכים לבנות את אופרטורי ההעלאה בחזקה המודולריים. כדי לקבל דיוק מספיק בהערכת הפאזה, נשתמש בשמונה קיוביטים למדידת ההערכה. לכן, אנחנו צריכים לבנות $M_b$ עם $b = a^{2^k} \; (\mathrm{mod} \; N)$ עבור כל $k = 0, 1, \dots, 7$.

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]

כפי שאנו יכולים לראות מרשימת ערכי $b$, בנוסף ל-$M_2$ שבנינו בעבר, אנחנו צריכים גם לבנות $M_4$ ו-$M_1$. שימו לב ש-$M_1$ פועל באופן טריוויאלי על מצבי הבסיס החישובי, כך שזה פשוט אופרטור הזהות.

$M_4$ פועל על מצבי הבסיס החישובי כדלקמן. $M_4 \ket{0} = \ket{0} \quad M_4 \ket{5} = \ket{5} \quad M_4 \ket{10} = \ket{10}$ $M_4 \ket{1} = \ket{4} \quad M_4 \ket{6} = \ket{9} \quad M_4 \ket{11} = \ket{14}$ $M_4 \ket{2} = \ket{8} \quad M_4 \ket{7} = \ket{13} \quad M_4 \ket{12} = \ket{3}$ $M_4 \ket{3} = \ket{12} \quad M_4 \ket{8} = \ket{2} \quad M_4 \ket{13} = \ket{7}$ $M_4 \ket{4} = \ket{1} \quad M_4 \ket{9} = \ket{6} \quad M_4 \ket{14} = \ket{11}$

לכן, תמורה זו ניתנת לבנייה עם פעולת ה-swap הבאה.

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)

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)

שערים הפועלים על יותר משני קיוביטים יפורקו עוד יותר לשערים של שני קיוביטים.

circ.decompose(reps=2).draw(output="mpl", fold=-1)

ראינו שאופרטורי $M_b$ עבור $b \in \mathbb{Z}^*_N$ נתון הם פעולות תמורה. בשל הגודל הקטן יחסית של בעיית התמורה שיש לנו כאן, מכיוון ש-$N=15$ דורש רק ארבעה קיוביטים, הצלחנו לסנתז את הפעולות הללו ישירות עם שערי SWAP על ידי בדיקה. באופן כללי, זו עשויה שלא להיות גישה ניתנת להרחבה. במקום זאת, ייתכן שנצטרך לבנות את מטריצת התמורה במפורש ולהשתמש במחלקה UnitaryGate של Qiskit ובשיטות תרגום כדי לסנתז את מטריצת התמורה הזו. עם זאת, זה יכול להביא למעגלים עמוקים משמעותית יותר. דוגמה כלהלן.

def mod_mult_gate(b, N):
    """
    Modular multiplication gate from permutation matrix.
    """
    if gcd(b, N) > 1:
        print(f"Error: gcd({b},{N}) > 1")
    else:
        n = floor(log(N - 1, 2)) + 1
        U = np.full((2**n, 2**n), 0)
        for x in range(N):
            U[b * x % N][x] = 1
        for x in range(N, 2**n):
            U[x][x] = 1
        G = UnitaryGate(U)
        G.name = f"M_{b}"
        return G

# Let's build M2 using the permutation matrix definition
M2_other = mod_mult_gate(2, 15)

# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2_other, inplace=True)
circ = circ.decompose()

# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)

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

qubits: 4
2q-depth: 94
2q-size: 96
Operator counts: OrderedDict({'cx': 45, 'swap': 32, 'u': 24, 'u1': 7, 'u3': 4, 'unitary': 3, 'circuit-335': 1, 'circuit-338': 1, 'circuit-341': 1, 'circuit-344': 1, 'circuit-347': 1, 'circuit-350': 1, 'circuit-353': 1, 'circuit-356': 1, 'circuit-359': 1, 'circuit-362': 1, 'circuit-365': 1, 'circuit-368': 1, 'circuit-371': 1, 'circuit-374': 1, 'circuit-377': 1, 'circuit-380': 1})

בואו נשווה את הספירות האלו עם עומק המעגל המקומפל של היישום הידני שלנו של שער $M_2$.

# Get the M2 operator from our manual construction
M2 = M2mod15()

# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ = circ.decompose(reps=3)

# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)

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

qubits: 4
2q-depth: 9
2q-size: 9
Operator counts: OrderedDict({'cx': 9})

כפי שאנו יכולים לראות, גישת מטריצת התמורה הביאה למעגל עמוק משמעותית אפילו עבור שער $M_2$ בודד בהשוואה ליישום הידני שלנו. לכן, נמשיך עם היישום הקודם שלנו של פעולות $M_b$. כעת, אנחנו מוכנים לבנות את מעגל מציאת הסדר המלא תוך שימוש באופרטורי ההעלאה בחזקה המודולריים המבוקרים שהגדרנו קודם לכן. בקוד הבא, אנחנו גם מייבאים את מעגל QFT מספריית המעגלים של Qiskit, שמשתמש בשערי Hadamard על כל קיוביט, בסדרה של שערי controlled-U1 (או Z, בהתאם לפאזה), ובשכבה של שערי swap.

# 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(QFT(num_control, inverse=True), qubits=control, inplace=True)

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

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

שימו לב שהשמטנו את פעולות ההעלאה בחזקה המודולריות המבוקרות משאר קיוביטי הבקרה מכיוון ש-$M_1$ הוא אופרטור הזהות. שימו לב שמאוחר יותר במדריך זה, נריץ את המעגל הזה על ה-backend ibm_marrakish. כדי לעשות זאת, אנחנו מתרגמים את המעגל בהתאם ל-backend ספציפי זה ומדווחים על עומק המעגל ומספר השערים.

service = QiskitRuntimeService()
backend = service.backend("ibm_marrakesh")
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: 187
2q-size: 260
Operator counts: OrderedDict({'sx': 521, 'rz': 354, 'cz': 260, 'measure': 8, 'x': 4})

שלב 3: ביצוע באמצעות הפרימיטיבים של Qiskit

ראשית, אנחנו דנים במה שהיינו מקבלים תיאורטית אם היינו מריצים את המעגל הזה על סימולטור אידיאלי. להלן, יש לנו סט של תוצאות סימולציה של המעגל לעיל תוך שימוש ב-1024 shots. כפי שאנו יכולים לראות, אנחנו מקבלים התפלגות אחידה בקירוב על פני ארבעה מחרוזות ביט על פני קיוביטי הבקרה.

# Obtained from the simulator
counts = {"00000000": 264, "01000000": 268, "10000000": 249, "11000000": 243}

plot_histogram(counts)

על ידי מדידת קיוביטי הבקרה, אנחנו מקבלים הערכת פאזה של שמונה ביט של האופרטור $M_a$. אנחנו יכולים להמיר את הייצוג הבינארי הזה לעשרוני כדי למצוא את הפאזה שנמדדה. כפי שאנו יכולים לראות מההיסטוגרמה לעיל, ארבעה מחרוזות ביט שונות נמדדו, וכל אחת מהן מתאימה לערך פאזה כדלקמן.

# Rows to be displayed in table
rows = []
# Corresponding phase of each bitstring
measured_phases = []

for output in counts:
    decimal = int(output, 2)  # Convert bitstring to decimal
    phase = decimal / (2**num_control)  # Find corresponding eigenvalue
    measured_phases.append(phase)
    # Add these values to the rows in our table:
    rows.append(
        [
            f"{output}(bin) = {decimal:>3}(dec)",
            f"{decimal}/{2 ** num_control} = {phase:.2f}",
        ]
    )

# Print the rows in a table
headers = ["Register Output", "Phase"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Register Output           Phase
0  00000000(bin) =   0(dec)    0/256 = 0.00
1  01000000(bin) =  64(dec)   64/256 = 0.25
2  10000000(bin) = 128(dec)  128/256 = 0.50
3  11000000(bin) = 192(dec)  192/256 = 0.75

זכרו שכל פאזה שנמדדה מתאימה ל-$\theta = k / r$ כאשר $k$ נדגמת באופן אחיד אקראי מ-$\{0, 1, \dots, r-1 \}$. לכן, אנחנו יכולים להשתמש באלגוריתם השברים המתמשכים כדי לנסות למצוא את $k$ ואת הסדר $r$. לפייתון יש פונקציונליות זו מובנית. אנחנו יכולים להשתמש במודול fractions כדי להפוך מספר עשרוני לאובייקט Fraction, למשל:

Fraction(0.666)

Fraction(5998794703657501, 9007199254740992)

מכיוון שזה נותן שברים שמחזירים את התוצאה במדויק (במקרה זה, 0.6660000...), זה יכול לתת תוצאות מסובכות כמו זו למעלה. אנחנו יכולים להשתמש במתודה .limit_denominator() כדי לקבל את השבר הדומה ביותר למספר העשרוני שלנו, עם מכנה מתחת לערך מסוים:

# Get fraction that most closely resembles 0.666
# with denominator < 15
Fraction(0.666).limit_denominator(15)

Fraction(2, 3)

זה הרבה יותר נחמד. הסדר (r) חייב להיות קטן מ-N, אז נקבע את המכנה המקסימלי להיות 15:

# Rows to be displayed in a table
rows = []

for phase in measured_phases:
    frac = Fraction(phase).limit_denominator(15)
    rows.append(
        [phase, f"{frac.numerator}/{frac.denominator}", frac.denominator]
    )

# Print the rows in a table
headers = ["Phase", "Fraction", "Guess for r"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Phase Fraction  Guess for r
0   0.00      0/1            1
1   0.25      1/4            4
2   0.50      1/2            2
3   0.75      3/4            4

אנחנו יכולים לראות ששניים מערכי העצמי שנמדדו סיפקו לנו את התוצאה הנכונה: $r=4$, ואנו יכולים לראות שלאלגוריתם שור למציאת סדר יש סיכוי לכישלון. תוצאות רעות אלו הן משום ש-$k = 0$, או משום ש-$k$ ו-$r$ אינם זרים זה לזה - ובמקום $r$, ניתן לנו גורם של $r$. הפתרון הקל ביותר לזה הוא פשוט לחזור על הניסוי עד שנקבל תוצאה מספקת עבור $r$. עד כה, יישמנו את בעיית מציאת הסדר עבור $N=15$ עם $a=2$ תוך שימוש במעגל הערכת הפאזה על סימולטור. השלב האחרון של אלגוריתם שור יהיה לקשר את בעיית מציאת הסדר לבעיית פירוק גורמים של מספרים שלמים. חלק אחרון זה של האלגוריתם הוא קלאסי לחלוטין וניתן לפתור על מחשב קלאסי לאחר שהמדידות של הפאזה התקבלו ממחשב קוונטי. לכן, אנחנו נדחה את החלק האחרון של האלגוריתם עד לאחר שנדגים כיצד אנחנו יכולים להריץ את מעגל מציאת הסדר על חומרה אמיתית.

הרצות חומרה

כעת אנחנו יכולים להריץ את מעגל מציאת הסדר שתרגמנו קודם לכן עבור ibm_marrakesh. כאן אנחנו פונים ל-dynamical decoupling (DD) לדיכוי שגיאות, ו-gate twirling למטרות צמצום שגיאות. DD כולל הפעלת רצפים של פולסי בקרה מתוזמנים במדויק להתקן קוונטי, תוך ממוצע יעיל של אינטראקציות סביבתיות לא רצויות ודקוהרנציה. Gate twirling, מצד שני, מבצע רנדומיזציה של שערים קוונטיים ספציפיים כדי להפוך שגיאות קוהרנטיות לשגיאות Pauli, שמצטברות באופן לינארי ולא ריבועי. שתי הטכניקות משולבות לעתים קרובות כדי לשפר את הקוהרנטיות והאמינות של חישובים קוונטיים.

# 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

# Assign tags before executing
sampler.options.environment.job_tags = ["TUT_SA"]

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

כפי שאנו יכולים לראות, קיבלנו את אותן מחרוזות ביט עם הספירות הגבוהות ביותר. מכיוון שלחומרה קוונטית יש רעש, יש דליפה מסוימת למחרוזות ביט אחרות, שאותן נוכל לסנן באופן סטטיסטי.

# 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)

{'00000000': 58, '01000000': 41, '11000000': 42, '10000000': 40}

שלב 4: עיבוד לאחר והחזרת תוצאה בפורמט קלאסי רצוי

פירוק גורמים של מספרים שלמים

עד כה, דנו כיצד אנחנו יכולים ליישם את בעיית מציאת הסדר באמצעות מעגל הערכת פאזה. כעת, אנחנו מקשרים את בעיית מציאת הסדר לפירוק גורמים של מספרים שלמים, מה שמשלים את אלגוריתם שור. שימו לב שחלק זה של האלגוריתם הוא קלאסי. אנחנו מדגימים זאת כעת תוך שימוש בדוגמה שלנו של $N = 15$ ו-$a = 2$. זכרו שהפאזה שמדדנו היא $k / r$, כאשר $a^r \; (\textrm{mod} \; N) = 1$ ו-$k$ הוא מספר שלם אקראי בין $0$ ל-$r - 1$. מהמשוואה הזו, יש לנו $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ מה שאומר ש-$N$ חייב לחלק את $a^r-1$. אם $r$ הוא גם זוגי, אז נוכל לכתוב $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ אם $r$ אינו זוגי, לא נוכל להמשיך והעלינו לנסות שוב עם ערך אחר עבור $a$; אחרת, יש סבירות גבוהה שהמחלק המשותף המקסימלי של $N$ ושל $a^{r/2}-1$, או $a^{r/2}+1$ הוא גורם אמיתי של $N$.

מכיוון שחלק מהריצות של האלגוריתם ייכשלו סטטיסטית, נחזור על האלגוריתם הזה עד שלפחות גורם אחד של $N$ יימצא. התא שלהלן חוזר על האלגוריתם עד שלפחות גורם אחד של $N=15$ נמצא. נשתמש בתוצאות של הריצה על החומרה למעלה כדי לנחש את הפאזה ואת הגורם המתאים בכל איטרציה.

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.25
Order of 2 modulo 15 estimated as: r = 4
*** Non-trivial factor found: 3 ***

דיון

עבודה קשורה

בחלק זה, אנחנו דנים בעבודות אבן דרך אחרות שהדגימו את אלגוריתם שור על חומרה אמיתית.

העבודה החלוצית [3] מ-IBM® הדגימה את אלגוריתם שור בפעם הראשונה, תוך פירוק גורמים של המספר 15 לגורמים הראשוניים שלו 3 ו-5 תוך שימוש במחשב קוונטי תהודה מגנטית גרעינית (NMR) בעל שבעה קיוביטים. ניסוי אחר [4] פירק גורמים של 15 תוך שימוש בקיוביטים פוטוניים. על ידי שימוש בקיוביט בודד שממוחזר מספר פעמים וקידוד רגיסטר העבודה במצבים רב-ממדיים, החוקרים צמצמו את מספר הקיוביטים הנדרש לשליש מזה שבפרוטוקול הסטנדרטי, תוך שימוש באלגוריתם מקומפל בן שני פוטונים. מאמר משמעותי בהדגמה של אלגוריתם שור הוא [5], המשתמש בטכניקת הערכת הפאזה האיטרטיבית של קיטאייב [8] כדי לצמצם את דרישת הקיוביטים של האלגוריתם. המחברים השתמשו בשבעה קיוביטי בקרה וארבעה קיוביטי מטמון, יחד עם יישום של כופלים מודולריים. עם זאת, יישום זה דורש מדידות אמצע-מעגל עם פעולות feed-forward ומיחזור קיוביטים עם פעולות איפוס. הדגמה זו נעשתה על מחשב קוונטי מלכודת יונים.

עבודה עדכנית יותר [6] התמקדה בפירוק גורמים של 15, 21 ו-35 על חומרת IBM Quantum®. בדומה לעבודה קודמת, החוקרים השתמשו בגרסה מקומפלת של האלגוריתם שהשתמשה בטרנספורמציה של פורייה קוונטית חצי-קלאסית כפי שהוצע על ידי קיטאייב כדי למזער את מספר הקיוביטים והשערים הפיזיים. עבודה עדכנית ביותר [7] ביצעה גם הדגמת הוכחת קונספט לפירוק גורמים של המספר השלם 21. הדגמה זו כללה גם שימוש בגרסה מקומפלת של שגרת הערכת הפאזה הקוונטית, ונבנתה על ההדגמה הקודמת על ידי [4]. המחברים עברו מעבר לעבודה זו על ידי שימוש בתצורה של שערי Toffoli משוערים עם הסטות פאזה שיוריות. האלגוריתם יושם על מעבדים קוונטיים של IBM תוך שימוש בחמישה קיוביטים בלבד, ונוכחות של אנטגלמנט בין קיוביטי הבקרה וקיוביטי הרגיסטר אומתה בהצלחה.

קנה מידה של האלגוריתם

אנחנו מציינים שהצפנת RSA כוללת בדרך כלל גדלי מפתח בסדר גודל של 2048 עד 4096 ביט. ניסיון לפרק גורמים של מספר בן 2048 ביט עם אלגוריתם שור יביא למעגל קוונטי עם מיליוני קיוביטים, כולל תקורת תיקון שגיאות ועומק מעגל בסדר גודל של מיליארד, שהוא מעבר לגבולות של החומרה הקוונטית הנוכחית לביצוע. לכן, אלגוריתם שור ידרוש או שיטות בניית מעגלים מיטביות או תיקון שגיאות קוונטי חזק כדי להיות ישים מבחינה מעשית לשבירת מערכות קריפטוגרפיות מודרניות. אנחנו מפנים אתכם ל-[9] לדיון מפורט יותר על הערכת משאבים עבור אלגוריתם שור.

אתגר

ברכותינו על השלמת המדריך! כעת זה זמן מצוין לבדוק את ההבנה שלך. האם תוכל לנסות לבנות את המעגל לפירוק גורמים של 21? אתה יכול לבחור $a$ לפי בחירתך. תצטרך להחליט על דיוק הביטים של האלגוריתם כדי לבחור את מספר הקיוביטים, ותצטרך לעצב את אופרטורי ההעלאה בחזקה המודולריים $M_a$. אנחנו מעודדים אותך לנסות זאת בעצמך, ולאחר מכן לקרוא על המתודולוגיות שמוצגות באיור 9 של [6] ובאיור 2 של [7].

def M_a_mod21():
    """
    M_a (mod 21)
    """

    # Your code here
    pass

הפניות

1. Shor, Peter W. "Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer." SIAM review 41.2 (1999): 303-332.
2. IBM Quantum Fundamentals of Quantum Algorithms course by Dr. John Watrous.
3. Vandersypen, Lieven MK, et al. "Experimental realization of Shor's quantum factoring algorithm using nuclear magnetic resonance." Nature 414.6866 (2001): 883-887.
4. Martin-Lopez, Enrique, et al. "Experimental realization of Shor's quantum factoring algorithm using qubit recycling." Nature photonics 6.11 (2012): 773-776.
5. Monz, Thomas, et al. "Realization of a scalable Shor algorithm." Science 351.6277 (2016): 1068-1070.
6. Amico, Mirko, Zain H. Saleem, and Muir Kumph. "Experimental study of Shor's factoring algorithm using the IBM Q Experience." Physical Review A 100.1 (2019): 012305.
7. Skosana, Unathi, and Mark Tame. "Demonstration of Shor's factoring algorithm for N=21 on IBM quantum processors." Scientific reports 11.1 (2021): 16599.
8. Kitaev, A. Yu. "Quantum measurements and the Abelian stabilizer problem." arXiv preprint quant-ph/9511026 (1995).
9. Gidney, Craig, and Martin Ekerå. "How to factor 2048 bit RSA integers in 8 hours using 20 million noisy qubits." Quantum 5 (2021): 433.