import numpy as np
from numpy.random import randint
from fractions import Fraction
from qiskit import assemble, Aer, transpile, QuantumCircuit, ClassicalRegister
from diskit import *

Distribute Shor’s Algorithm

The code for building the monolithic circuit for factoring 15 is taken from Qiskit Library. https://qiskit.org/textbook/ch-algorithms/shor.html

In this example we will solve the period finding problem for \(a=7\) and \(N=15\). We provide the circuits for \(U\) where:

\[U|y\rangle = |ay\bmod 15\rangle\]

without explanation. To create \(U^x\), we will simply repeat the circuit \(x\) times. In the next section we will discuss a general method for creating these circuits efficiently. The function c_amod15 returns the controlled-U gate for a, repeated power times.

def c_amod15(a, power):
    """Controlled multiplication by a mod 15"""
    if a not in [2,4,7,8,11,13]:
        raise ValueError("'a' must be 2,4,7,8,11 or 13")
    U = QuantumCircuit(4)
    for iteration in range(power):
        if a in [2,13]:
            U.swap(2,3)
            U.swap(1,2)
            U.swap(0,1)
        if a in [7,8]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
        if a in [4, 11]:
            U.swap(1,3)
            U.swap(0,2)
        if a in [7,11,13]:
            for q in range(4):
                U.x(q)
    U = U.to_gate()
    U.name = "%i^%i mod 15" % (a, power)
    c_U = U.control()
    return c_U

We also import the circuit for the QFT (you can read more about the QFT in the quantum Fourier transform chapter):

def qft_dagger(n):
    """n-qubit QFTdagger the first n qubits in circ"""
    qc = QuantumCircuit(n)
    # Don't forget the Swaps!
    for qubit in range(n//2):
        qc.swap(qubit, n-qubit-1)
    for j in range(n):
        for m in range(j):
            qc.cp(-np.pi/float(2**(j-m)), m, j)
        qc.h(j)
    qc.name = "QFT†"
    return qc

Factoring from Period Finding

Not all factoring problems are difficult; we can spot an even number instantly and know that one of its factors is 2. In fact, there are specific criteria for choosing numbers that are difficult to factor, but the basic idea is to choose the product of two large prime numbers.

A general factoring algorithm will first check to see if there is a shortcut to factoring the integer (is the number even? Is the number of the form \(N = a^b\)?), before using Shor’s period finding for the worst-case scenario. Since we aim to focus on the quantum part of the algorithm, we will jump straight to the case in which N is the product of two primes.

To see an example of factoring on a small number of qubits, we will factor 15, which we all know is the product of the not-so-large prime numbers 3 and 5.

We will use 8 counting qubits, so in the topology the number of qubits must be \(\geq\) 8 + 4

N = 15

The first step is to choose a random number, \(a\), between \(1\) and \(N-1\):

np.random.seed(1) # This is to make sure we get reproduceable results
a = randint(2, 15)
print(a)

7

Next we quickly check it isn’t already a non-trivial factor of \(N\):

from math import gcd # greatest common divisor
gcd(a, N)

1

At First we see the monolithic implementation.

def qpe_amod15(a):
    n_count = 8
    qc = QuantumCircuit(4+n_count, n_count)
    for q in range(n_count):
        qc.h(q)     # Initialize counting qubits in state |+>
    qc.x(3+n_count) # And auxiliary register in state |1>
    for q in range(n_count): # Do controlled-U operations
        qc.append(c_amod15(a, 2**q),
                 [q] + [i+n_count for i in range(4)])
    qc.append(qft_dagger(n_count), range(n_count)) # Do inverse-QFT
    qc.measure(range(n_count), range(n_count))
    # Simulate Results
    aer_sim = Aer.get_backend('aer_simulator')
    # Setting memory=True below allows us to see a list of each sequential reading
    t_qc = transpile(qc, aer_sim)
    qobj = assemble(t_qc, shots=1)
    result = aer_sim.run(qobj, memory=True).result()
    readings = result.get_memory()
    print("Register Reading: " + readings[0])
    phase = int(readings[0],2)/(2**n_count)
    print("Corresponding Phase: %f" % phase)
    return phase, qc

We can modify the monolithic version to a distributed version of the function conveniently with minimal changes in structure and data

def qpe_amod15_dist(a, topology):
    qregs = topology.get_regs()
    n_count = 8
    qc = QuantumCircuit(*qregs, ClassicalRegister(n_count))
    for q in range(n_count):
        qc.h(q)     # Initialize counting qubits in state |+>
    qc.x(3+n_count) # And auxiliary register in state |1>
    for q in range(n_count): # Do controlled-U operations
        qc.append(c_amod15(a, 2**q),
                 [q] + [i+n_count for i in range(4)])
    qc.append(qft_dagger(n_count), range(n_count)) # Do inverse-QFT
    qc.measure(range(n_count), range(n_count))
    remapper = CircuitRemapper(circuit_topo)
    dist_circ = remapper.remap_circuit(qc, decompose=True)

    # Simulate Results
    aer_sim = Aer.get_backend('aer_simulator')

    # Setting memory=True below allows us to see a list of each sequential reading
    t_qc = transpile(dist_circ, aer_sim)
    qobj = assemble(t_qc, shots=1)
    result = aer_sim.run(qobj, memory=True).result()
    readings = result.get_memory()

    ## Need to manipulate readouts to exclude cat_measurement registers
    final_readings = readings[0].split()
    for reading in final_readings:
        if len(reading) == n_count:
            final_reading = reading

    print("Register Reading: " + final_reading)
    phase = int(final_reading,2)/(2**n_count)
    print("Corresponding Phase: %f" % phase)
    return phase,  dist_circ

Again, Topology qubits must be :math:`geq` 8 + 4, we have taken 13 qubits with distribution over 2 qpus

circuit_topo = Topology()
circuit_topo.create_qmap(2, [6, 7],"sys_shor")
circuit_topo.qmap, circuit_topo.emap

({'sys_shor0': [Qubit(QuantumRegister(6, 'sys_shor0'), 0),
   Qubit(QuantumRegister(6, 'sys_shor0'), 1),
   Qubit(QuantumRegister(6, 'sys_shor0'), 2),
   Qubit(QuantumRegister(6, 'sys_shor0'), 3),
   Qubit(QuantumRegister(6, 'sys_shor0'), 4),
   Qubit(QuantumRegister(6, 'sys_shor0'), 5)],
  'sys_shor1': [Qubit(QuantumRegister(7, 'sys_shor1'), 0),
   Qubit(QuantumRegister(7, 'sys_shor1'), 1),
   Qubit(QuantumRegister(7, 'sys_shor1'), 2),
   Qubit(QuantumRegister(7, 'sys_shor1'), 3),
   Qubit(QuantumRegister(7, 'sys_shor1'), 4),
   Qubit(QuantumRegister(7, 'sys_shor1'), 5),
   Qubit(QuantumRegister(7, 'sys_shor1'), 6)]},
 {'sys_shor0': Qubit(QuantumRegister(1, 'com_sys_shor0'), 0),
  'sys_shor1': Qubit(QuantumRegister(1, 'com_sys_shor1'), 0)})

From this phase, we can easily find a guess for \(r\):

phase, circ = qpe_amod15_dist(a, circuit_topo)
Fraction(phase).limit_denominator(15) # Denominator should (hopefully!) tell us r

Single level decomposition cutoff of 1 minute reached. Performing transpilation with basis gates: cx, u1, u2, u3, id
Register Reading: 11000000
Corresponding Phase: 0.750000

Fraction(3, 4)

frac = Fraction(phase).limit_denominator(15)
s, r = frac.numerator, frac.denominator
print(r)

4

Now we have \(r\), we might be able to use this to find a factor of \(N\). Since:

\[a^r \bmod N = 1\]

then:

\[(a^r - 1) \bmod N = 0\]

which means \(N\) must divide \(a^r-1\). And if \(r\) is also even, then we can write:

\[a^r -1 = (a^{r/2}-1)(a^{r/2}+1)\]

(if \(r\) is not even, we cannot go further and must try again with a different value for \(a\)). There is then a high probability that the greatest common divisor of \(N\) and either \(a^{r/2}-1\), or \(a^{r/2}+1\) is a proper factor of \(N\) [2]:

guesses = [gcd(a**(r//2)-1, N), gcd(a**(r//2)+1, N)]
print(guesses)

[3, 5]

The cell below repeats the algorithm until at least one factor of 15 is found.

a = 7
factor_found = False
attempt = 0
while not factor_found:
    attempt += 1
    print("\nAttempt %i:" % attempt)
    phase,_ = qpe_amod15_dist(a, circuit_topo) # Phase = s/r
    frac = Fraction(phase).limit_denominator(N) # Denominator should (hopefully!) tell us r
    r = frac.denominator
    print("Result: r = %i" % r)
    if phase != 0:
        # Guesses for factors are gcd(x^{r/2} ±1 , 15)
        guesses = [gcd(a**(r//2)-1, N), gcd(a**(r//2)+1, N)]
        print("Guessed Factors: %i and %i" % (guesses[0], guesses[1]))
        for guess in guesses:
            if guess not in [1,N] and (N % guess) == 0: # Check to see if guess is a factor
                print("*** Non-trivial factor found: %i ***" % guess)
                factor_found = True

Attempt 1:
Single level decomposition cutoff of 1 minute reached. Performing transpilation with basis gates: cx, u1, u2, u3, id
Register Reading: 01000000
Corresponding Phase: 0.250000
Result: r = 4
Guessed Factors: 3 and 5
*** Non-trivial factor found: 3 ***
*** Non-trivial factor found: 5 ***

References

1. Stephane Beauregard, Circuit for Shor’s algorithm using 2n+3 qubits, arXiv:quant-ph/0205095

2. 13. Nielsen and I. Chuang, Quantum Computation and Quantum Information, Cambridge Series on Information and the Natural Sciences (Cambridge University Press, Cambridge, 2000). (Page 633)