LogoQCoder Programming Contest 004
Sign in
← Previous Problem

C5: Modular Arithmetic IV

Time Limit: 10 sec

Memory Limit: 512 MiB

Score: 600

Writer: PSL

How QPC004's problems relate to Shor's algorithm is described in the follow-up at the bottom of this page.

Editorial

First, expand kk in binary representation:

k=k0+2k1+⋯+2n−1kn−1(k0,k1,…,kn−1∈{0,1})\begin{equation} k = k_0 + 2k_1 + \cdots + 2^{n-1}k_{n-1} \quad (k_0, k_1, \ldots, k_{n-1} \in \{0, 1\}) \end{equation}

Then, the desired state can be expanded as follows:

12m∑k=02n−1∣k⟩n∣ak mod L⟩2m+1=12m∑k=02n−1∣k0⟩∣k1⟩⋯∣kn−1⟩∣ak0+2k1+⋯2n−1kn−1 mod L⟩2m+1=12m∑k=02n−1∣k0⟩∣k1⟩⋯∣kn−1⟩∣ak0a2k1⋯a2n−1kn−1 mod L⟩2m+1\begin{align} &\frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^n-1} \ket{k}_{n} \ket{a^k \ \text{mod} \ L}_{2m+1} \nonumber\\ &= \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^n-1} \ket{k_0} \ket{k_1} \cdots \ket{k_{n-1}} \ket{a^{k_0 + 2k_1 + \cdots 2^{n-1}k_{n-1}} \ \text{mod} \ L}_{2m+1} \nonumber\\ &= \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^n-1} \ket{k_0} \ket{k_1} \cdots \ket{k_{n-1}} \ket{a^{k_0}a^{2k_1}\cdots a^{2^{n-1}k_{n-1}} \ \text{mod} \ L}_{2m+1} \end{align}

Therefore, by using the operation of problem C4 with ∣ki⟩\ket{k_i} as the control bit, the desired operation can be realized by using the controlled circuit CC4\mathrm{CC4} of problem C4.

∣0⟩n∣0⟩2m+1→H⊗n12m∑k=02n−1∣k0⟩∣k1⟩⋯∣kn−1⟩∣0⟩2m+1→X12m∑k=02n−1∣k0⟩∣k1⟩⋯∣kn−1⟩∣1⟩2m+1→CC4(a)12m∑k=02n−1∣k0⟩∣k1⟩⋯∣kn−1⟩∣ak0 mod L⟩2m+1→CC4(a2)12m∑k=02n−1∣k0⟩∣k1⟩⋯∣kn−1⟩∣ak0a2k1 mod L⟩2m+1⋯→CC4(a2n−1)12m∑k=02n−1∣k0⟩∣k1⟩⋯∣kn−1⟩∣ak0a2k1⋯a2n−1kn−1 mod L⟩2m+1\begin{align} &\ket{0}_n \ket{0}_{2m+1} \nonumber\\ &\xrightarrow{H^{\otimes n}} \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^n-1} \ket{k_0} \ket{k_1} \cdots \ket{k_{n-1}} \ket{0}_{2m+1} \nonumber\\ &\xrightarrow{X} \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^n-1} \ket{k_0} \ket{k_1} \cdots \ket{k_{n-1}} \ket{1}_{2m+1} \nonumber\\ &\xrightarrow{\mathrm{CC4}(a)} \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^n-1} \ket{k_0} \ket{k_1} \cdots \ket{k_{n-1}} \ket{a^{k_0} \ \text{mod} \ L}_{2m+1} \nonumber\\ &\xrightarrow{\mathrm{CC4}(a^2)} \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^n-1} \ket{k_0} \ket{k_1} \cdots \ket{k_{n-1}} \ket{a^{k_0}a^{2k_1} \ \text{mod} \ L}_{2m+1} \nonumber\\ &\cdots \nonumber \\ &\xrightarrow{\mathrm{CC4}(a^{2^{n-1}})} \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^n-1} \ket{k_0} \ket{k_1} \cdots \ket{k_{n-1}} \ket{a^{k_0}a^{2k_1} \cdots a^{2^{n-1}k_{n-1}} \ \text{mod} \ L}_{2m+1} \end{align}

The circuit depth is O(n3)O(n^3).

Sample Code

Below is a sample program:

import math
 
from qiskit import QuantumCircuit, QuantumRegister
 
 
def inv_mod(a: int, L: int) -> int:
    def extended_euclidean(a: int, b: int) -> tuple[int, int, int]:
        # when b == 0:
        if b == 0:
            return a, 1, 0
 
        # recursive step
        gcd, x1, y1 = extended_euclidean(b, a % b)
 
        # update x and y
        x = y1
        y = x1 - (a // b) * y1
 
        return gcd, x, y
 
    _, x, _ = extended_euclidean(a, L)
    return x % L
 
 
def beki(n: int, a: int, L: int) -> list[int]:
    result = [a]
    for _ in range(n - 1):
        result.append(result[-1] ** 2 % L)
 
    return result
 
 
def qft(n: int) -> QuantumCircuit:
    qc = QuantumCircuit(n)
 
    for i in reversed(range(n)):
        qc.h(i)
        for j in reversed(range(i)):
            qc.cp(math.pi / 2 ** (i - j), j, i)
 
    for i in range(n // 2):
        qc.swap(i, n - i - 1)
 
    return qc
 
 
def mcadd(n: int, a: int, c: int) -> QuantumCircuit:
    qc = QuantumCircuit(n + c)
 
    qc.compose(qft(n), qubits=range(n), inplace=True)
    for i in range(n):
        theta = 2 * math.pi * a * 2**i / 2**n
        if c == 0:
            qc.p(theta, i)
        elif c == 1:
            qc.cp(theta, n, i)
        else:
            qc.mcp(theta, list(range(n, n + c)), i)
    qc.compose(qft(n).inverse(), qubits=range(n), inplace=True)
 
    return qc
 
 
def mcadd_mod(n: int, a: int, L: int, c: int) -> QuantumCircuit:
    qc = QuantumCircuit(n + 1 + c)
 
    qc.compose(mcadd(n + 1, 2**n - L + a, c), qubits=range(n + 1 + c), inplace=True)
 
    qc.x(n)
    qc.compose(mcadd(n, -(2**n - L + a), c + 1), qubits=range(n + 1 + c), inplace=True)
    qc.x(n)
 
    qc.compose(mcadd(n, 2**n - a, c + 1), qubits=range(n + 1 + c), inplace=True)
 
    qc.compose(mcadd(n + 1, a, c), qubits=range(n + 1 + c), inplace=True)
 
    return qc
 
 
def cswap(n: int) -> QuantumCircuit:
    qc = QuantumCircuit(2 * n + 1)
 
    for i in range(n):
        qc.cswap(2 * n, i, i + n)
 
    return qc
 
 
def solve(n: int, m: int, a: int, L: int) -> QuantumCircuit:
    x, y = QuantumRegister(n), QuantumRegister(2 * m + 1)
    qc = QuantumCircuit(x, y)
 
    qc.h(x)
    qc.x(y[0])
 
    lst = beki(n, a, L)
    a_inv = inv_mod(a, L)
    lst_inv = beki(n, a_inv, L)
    lst_minus_inv = [(-v) % L for v in lst_inv]
    for i in range(n):
        b = lst[i]
        for j in range(m):
            qc.compose(mcadd_mod(m, b, L, 2), qubits=[*y[m:], x[i], y[j]], inplace=True)
            b = 2 * b % L
 
        qc.compose(cswap(m), qubits=[*y[: 2 * m], x[i]], inplace=True)
 
        b = lst_minus_inv[i]
        for j in range(m):
            qc.compose(mcadd_mod(m, b, L, 2), qubits=[*y[m:], x[i], y[j]], inplace=True)
            b = 2 * b % L
 
    return qc

Follow-up

In the original draft of the contest problems, there were the following problems, D1 and D2, at the end. Although they were not used in the contest because of limitations in the judging system, we include the problems and their explanations here to deepen understanding of Shor’s algorithm.

Problem D1

You are given integers nn, mm, aa, and LL. The inputs satisfy 1≤L≤2n1 \leq L \leq 2^n and 0≤a<L0 \leq a < L, and aa and LL are coprime.

Implement, on a quantum circuit qc\mathrm{qc} with n+2m+1{n + 2m + 1} qubits, an operation that creates the superposition state ∣B⟩\ket{B} from the zero state. Then, by measuring the leftmost nn qubits, estimate the period rr of ak mod La^k \bmod L.

The superposition state ∣B⟩\ket{B} is defined as follows:

∣B⟩n+2m+1=(IQFT⊗I)12m∑k=02m−1∣k⟩n∣ak mod L⟩2m+1\ket{B}_{n + 2m +1} = \left( \mathrm{IQFT} \otimes I \right) \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} \ket{k}_{n} \ket{a^k \bmod L}_{2m + 1}

Editorial for D1

Open

We begin with

∣C5⟩=12m∑k=02m−1∣k⟩∣gk mod L⟩.\begin{equation} \ket{\mathrm{C5}} = \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} \ket{k} \ket{g^k \bmod L}. \end{equation}

We wish to transform ∣C5⟩\ket{\mathrm{C5}}; to do so, first we prepare the following equations, labeled A, B, C, and D.

We define the operator GG as follows:

For any xx,

G∣x mod L⟩=∣gx mod L⟩⋯A\begin{equation} G \ket{x \bmod L} = \ket{gx \bmod L} \cdots A \end{equation}

Now, for 0≤s<r0 \leq s < r, define

∣us⟩=1r∑k=0r−1exp⁡(−2πiksr)∣gk mod L⟩⋯B\begin{equation} \ket{u_{s}} = \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} \exp\left(\frac{-2\pi i ks}{r}\right) \ket{g^k \bmod L} \cdots B \end{equation}

Then the following properties C and D hold for ∣us⟩\ket{u_{s}}:

1r∑s=0r−1∣us⟩=∣1⟩⋯C\begin{equation} \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \ket{u_{s}} = \ket{1} \cdots C \end{equation}G∣us⟩=exp⁡(2πisr)∣us⟩⋯D\begin{equation} G \ket{u_{s}} = \exp\left(\frac{2\pi i s}{r}\right) \ket{u_{s}} \cdots D \end{equation}

Proof of C:

1r∑s=0r−1∣us⟩=1r∑s=0r−11r∑k=0r−1exp⁡(−2πiksr)∣gk mod L⟩=∑k=0r−1(1r∑s=0r−1exp⁡(−2πiksr))∣gk mod L⟩=∣1⟩\begin{align} \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \ket{u_{s}} &= \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} \exp\left(\frac{-2\pi i ks}{r}\right) \ket{g^k \bmod L} \nonumber\\ &= \sum_{k=0}^{r-1} \left( \frac{1}{r} \sum_{s=0}^{r-1} \exp\left(\frac{-2\pi i ks}{r}\right) \right) \ket{g^k \bmod L} \nonumber\\ &= \ket{1} \end{align}

Note: The final equality holds because:

  • For k=0k = 0
1r∑s=0r−1exp⁡(−2πiksr)=1r∑s=0r−11=1\frac{1}{r} \sum_{s=0}^{r-1} \exp\left(\frac{-2\pi i ks}{r}\right) = \frac{1}{r} \sum_{s=0}^{r-1} 1 = 1
  • For k≠0k \neq 0
1r∑s=0r−1exp⁡(−2πiksr)=1r1−exp⁡(−2πikrr)1−exp⁡(−2πikr)=1r1−11−exp⁡(−2πikr)=0\frac{1}{r} \sum_{s=0}^{r-1} \exp\left(\frac{-2\pi i ks}{r}\right) = \frac{1}{r} \frac{1 - \exp\left(\frac{-2\pi i kr}{r}\right)}{1 - \exp\left(\frac{-2\pi i k}{r}\right)} = \frac{1}{r} \frac{1 - 1}{1 - \exp\left(\frac{-2\pi i k}{r}\right)} = 0

Proof of D:

G∣us⟩=1r∑k=0r−1exp⁡(−2πiksr)G∣gk mod L⟩=1r∑k=0r−1exp⁡(−2πiksr)∣gk+1 mod L⟩=exp⁡(2πisr)1r∑k=0r−1exp⁡(−2πi(k+1)sr)∣gk+1 mod L⟩=exp⁡(2πisr)1r∑k′=0r−1exp⁡(−2πik′sr)∣gk′ mod L⟩=exp⁡(2πisr)∣us⟩\begin{align} G \ket{u_{s}} &= \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} \exp\left(\frac{-2\pi i ks}{r}\right) G \ket{g^k \bmod L} \nonumber\\ &= \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} \exp\left(\frac{-2\pi i ks}{r}\right) \ket{g^{k+1} \bmod L} \nonumber\\ &= \exp\left(\frac{2\pi i s}{r}\right) \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} \exp\left(\frac{-2\pi i (k+1)s}{r}\right) \ket{g^{k+1} \bmod L} \nonumber\\ &= \exp\left(\frac{2\pi i s}{r}\right) \frac{1}{\sqrt{r}} \sum_{k'=0}^{r-1} \exp\left(\frac{-2\pi i k's}{r}\right) \ket{g^{k'} \bmod L} \nonumber\\ &= \exp\left(\frac{2\pi i s}{r}\right) \ket{u_{s}} \end{align}

Now that we have prepared (A), (B), (C), and (D), we use them to transform ∣C5⟩\ket{\mathrm{C5}}:

∣C5⟩=12m∑k=02m−1∣k⟩∣gk mod L⟩=12m∑k=02m−1∣k⟩Gk∣1⟩=12m∑k=02m−1∣k⟩Gk1r∑s=0r−1∣us⟩=12m∑k=02m−1∣k⟩exp⁡(2πiskr)1r∑s=0r−1∣us⟩=12m∑k=02m−1exp⁡(2πiskr)∣k⟩1r∑s=0r−1∣us⟩\begin{align} \ket{\mathrm{C5}} &= \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} \ket{k} \ket{g^k \bmod L} \nonumber\\ &= \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} \ket{k} G^{k} \ket{1} \nonumber\\ &= \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} \ket{k} G^{k} \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \ket{u_{s}} \nonumber\\ &= \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} \ket{k} \exp\left(\frac{2\pi i s k}{r}\right) \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \ket{u_{s}} \nonumber\\ &= \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} \exp\left(\frac{2\pi i s k}{r}\right) \ket{k} \frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \ket{u_{s}} \end{align}

Although it is not completely rigorous, if we ignore the factor 1r∑s=0r−1∣us⟩\frac{1}{\sqrt{r}} \sum_{s=0}^{r-1} \ket{u_{s}}, we obtain

12m∑k=02m−1exp⁡(2πi(2ms/r)k2m)∣k⟩\begin{equation} \frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m-1} \exp\left(\frac{2\pi i (2^{m}s/r) k}{2^{m}}\right) \ket{k} \end{equation}

Now, if we apply the IQFT and obtain a measurement result λ\lambda, then λ≈2msr\lambda \approx \frac{2^{m}s}{r}, so that λ2m≈sr\frac{\lambda}{2^{m}} \approx \frac{s}{r}.

Note that both rr and ss are unknown. However, even this information is sufficient to narrow down the candidates for rr. For example, suppose λ2m=34138192\frac{\lambda}{2^{m}} = \frac{3413}{8192}, then, by using a continued fraction expansion, we can approximate as follows:

34138192=0+12+12+12+1170+14≈0+12+12+12+0=512\begin{align} \frac{3413}{8192} &= 0 + \cfrac{1}{2 + \cfrac{1}{2 + \cfrac{1}{2 + \cfrac{1}{170 + \frac{1}{4}}}}} \nonumber\\ &\approx 0 + \cfrac{1}{2 + \cfrac{1}{2 + \cfrac{1}{2 + 0}}} \nonumber\\ &= \frac{5}{12} \end{align}

Thus, we estimate that r=12r = 12.

Problem D1

Factorize 1515 into its prime factors.

Editorial for D2

Open
  1. Choose gg in the range 2≤g<L2 \leq g < L.
  2. If gg is a factor of LL, the factorization is complete.
  3. Find rr such that gr=1 mod Lg^r = 1 \bmod L using the method from D1.
  4. If rr is odd, start over.
  5. If rr is even, it is likely that (gr/2+1)(gr/2−1)=0 mod L(g^{r/2}+1)(g^{r/2}-1) = 0 \bmod L.
  6. Either gcd(gr/2+1,L)gcd(g^{r/2}+1, L) or gcd(gr/2−1,L)gcd(g^{r/2}-1, L) is likely to be a factor of LL. If the factorization is successful, the process is complete.
  7. If the factorization fails, start over.