Project Q: Shor's Algorithm

1from __future__ import print_function
2
3import strangeworks
4
5import math
6import random
7import sys
8from fractions import Fraction
9try:
10    from math import gcd
11except ImportError:
12    from fractions import gcd
13
14from builtins import input
15
16import projectq.libs.math
17import projectq.setups.decompositions
18from projectq.backends import Simulator, ResourceCounter
19from projectq.cengines import (AutoReplacer, DecompositionRuleSet,
20                               InstructionFilter, LocalOptimizer,
21                               MainEngine, TagRemover)
22from projectq.libs.math import (AddConstant, AddConstantModN,
23                                MultiplyByConstantModN)
24from projectq.meta import Control
25from projectq.ops import (All, BasicMathGate, get_inverse, H, Measure, QFT, R,
26                          Swap, X)
27
28
29def run_shor(eng, N, a, verbose=False):
30    """
31    Runs the quantum subroutine of Shor's algorithm for factoring.
32    Args:
33        eng (MainEngine): Main compiler engine to use.
34        N (int): Number to factor.
35        a (int): Relative prime to use as a base for a^x mod N.
36        verbose (bool): If True, display intermediate measurement results.
37    Returns:
38        r (float): Potential period of a.
39    """
40    n = int(math.ceil(math.log(N, 2)))
41
42    x = eng.allocate_qureg(n)
43
44    X | x[0]
45
46    measurements = [0] * (2 * n)  # will hold the 2n measurement results
47
48    ctrl_qubit = eng.allocate_qubit()
49
50    for k in range(2 * n):
51        current_a = pow(a, 1 << (2 * n - 1 - k), N)
52        # one iteration of 1-qubit QPE
53        H | ctrl_qubit
54        with Control(eng, ctrl_qubit):
55            MultiplyByConstantModN(current_a, N) | x
56
57        # perform inverse QFT --> Rotations conditioned on previous outcomes
58        for i in range(k):
59            if measurements[i]:
60                R(-math.pi/(1 << (k - i))) | ctrl_qubit
61        H | ctrl_qubit
62
63        # and measure
64        Measure | ctrl_qubit
65        eng.flush()
66        measurements[k] = int(ctrl_qubit)
67        if measurements[k]:
68            X | ctrl_qubit
69
70        if verbose:
71            print("{}".format(measurements[k]), end="")
72            sys.stdout.flush()
73
74    All(Measure) | x
75    # turn the measured values into a number in [0,1)
76    y = sum([(measurements[2 * n - 1 - i]*1. / (1 << (i + 1)))
77             for i in range(2 * n)])
78
79    # continued fraction expansion to get denominator (the period?)
80    r = Fraction(y).limit_denominator(N-1).denominator
81
82    # return the (potential) period
83    return r
84
85
86# Filter function, which defines the gate set for the first optimization
87# (don't decompose QFTs and iQFTs to make cancellation easier)
88def high_level_gates(eng, cmd):
89    g = cmd.gate
90    if g == QFT or get_inverse(g) == QFT or g == Swap:
91        return True
92    if isinstance(g, BasicMathGate):
93        return False
94        if isinstance(g, AddConstant):
95            return True
96        elif isinstance(g, AddConstantModN):
97            return True
98        return False
99    return eng.next_engine.is_available(cmd)
100
101
102if __name__ == "__main__":
103    # build compilation engine list
104    resource_counter = ResourceCounter()
105    rule_set = DecompositionRuleSet(modules=[projectq.libs.math,
106                                             projectq.setups.decompositions])
107    compilerengines = [AutoReplacer(rule_set),
108                       InstructionFilter(high_level_gates),
109                       TagRemover(),
110                       LocalOptimizer(3),
111                       AutoReplacer(rule_set),
112                       TagRemover(),
113                       LocalOptimizer(3),
114                       resource_counter]
115
116    # make the compiler and run the circuit on the simulator backend
117    eng = MainEngine(Simulator(), compilerengines)
118    
119    # what number do you want to factor?
120    N=35
121
122    # choose a base at random:
123    a = int(random.random()*N)
124    if not gcd(a, N) == 1:
125        print("\n\n\tOoops, we were lucky: Chose non relative prime"
126              " by accident :)")
127        print("\tFactor: {}".format(gcd(a, N)))
128    else:
129        # run the quantum subroutine
130        r = run_shor(eng, N, a, True)
131
132        # try to determine the factors
133        if r % 2 != 0:
134            r *= 2
135        apowrhalf = pow(a, r >> 1, N)
136        f1 = gcd(apowrhalf + 1, N)
137        f2 = gcd(apowrhalf - 1, N)
138        if ((not f1 * f2 == N) and f1 * f2 > 1 and
139                int(1. * N / (f1 * f2)) * f1 * f2 == N):
140            f1, f2 = f1*f2, int(N/(f1*f2))
141        if f1 * f2 == N and f1 > 1 and f2 > 1:
142            print("\n\n\tFactors found :-) : {} * {} = {}"
143                  .format(f1, f2, N))
144        else:
145            print("\n\n\tBad luck: Found {} and {}".format(f1,
146                                                                          f2))
147
148        print(resource_counter)  # print resource usage