Quantum Annealing for Prime Factorization.

We have developed a framework to convert an arbitrary integer factorization problem to an executable Ising model by first writing it as an optimization function then transforming the k-bit coupling (k ≥ 3) terms to quadratic terms using ancillary variables. Our resource-efficient method uses [Formula: see text] binary variables (qubits) for finding the factors of an integer N. We present how to factorize 15, 143, 59989, and 376289 using 4, 12, 59, and 94 logical qubits, respectively. This method was tested using the D-Wave 2000Q for finding an embedding and determining the prime factors for a given composite number. The method is general and could be used to factor larger integers as the number of available qubits increases, or combined with other ad hoc methods to achieve better performances for specific numbers.

Introduction

Integer factorization reduces an integer N to its factors p and q such that pq = N. While this fundamental problem in number theory is computationally hard in practice, integer factorization is not believed to belong to the class of NP-hard problems. However, all known classical factoring algorithms which are deterministic and don’t have unproven hypotheses require time exponential in logN. For example, the fastest, known classical algorithm for integer factorization is the general number field sieve method[1], which scales exponentially in the number of operations required with respect to the integer N. Thus, the integer factorization problem has been used as a basic hardness assumption for many encryption methods including the widely deployed RSA cryptographic system. With broad applications in cryptographic data storage and communications[2], identifying new methods for integer factorization plays an important role in modern information security.

Quantum computing theory has the potential to reduce the number of operations required for solving the integer factorization problem. Within the circuit model of quantum computation, Shor’s algorithm is perhaps the most well-known method for integer factorization, in which the number of operations to factorize an integer N is polynomial in the size logN[3]. This exponential speedup over the general number field sieve is achieved by reducing factorization to the order-finding problem. Several experimental demonstrations using quantum computing hardware have validated the correctness of Shor’s algorithm for small integer values, including early work by Vandersypen et al.[4] to factorize N = 15 using seven spin-1/2 nuclei in a molecule as qubits. Subsequent experiments by Lanyon et al.[5], Lu et al.[6], and Politi et al.[7] implemented compiled versions of Shor’s algorithm using photonic systems for factoring 15. Martín-López et al.[8] factored 21 using qubit recycling, and Lucero et al.[9] used superconducting qubits to factor 15. Geller et al.[10] used a simplified version of Shor’s algorithm for factoring products of the Fermat primes 3, 5, 17, 257, and 65537. Recent work from Grosshans et al. has shown how factoring safe semi-primes using the quantum order-finding algorithm can reduce the failure probability[11].

An equally powerful model of quantum computing is the adiabatic quantum computing (AQC) model[12,13], which can also solve the integer factorization problem. Peng et al. first developed integer factorization within AQC by reducing it to unconstrained optimization and solving this problem using adiabatic quantum dynamics. They further validated these ideas experimentally using a three-qubit NMR quantum processor for the case of N = 21[14], while Xu et al. subsequently factored 143 using similar NMR technology[15]. Schaller et al. developed a novel approach based on multiplication tables that can be cast as an optimization, which they have demonstrated for biprimes up to N = 217[16]. Dridi et al. furthered these ideas by using Gröbner bases to reduce the number of auxiliary variables required and simplify equations, thereby enabling a demonstration of factorization up to 223357[17].

In this contribution, we introduce a new procedure for solving the integer factorization problem using quantum annealing[18,19] which utilizes adiabatic quantum computation. Differ from the recent work[20] which sketched the hardware design of reversible multiplier to achieve factorization, we provide specific mathematical derivations to be tested on the existing hardware. We begin by describing a direct method for integer factorization that reduces the problem to unconstrained optimization. We review how this optimization problem can be reduced to a quadratic form Ising Hamiltonian and be solved using quantum annealing. We then describe a modified multiplication table method that reduces the overall resource requirements on the optimization problem and permits methods to account for constraints that may appear in quantum annealing hardware, such as connectivity and number of qubits. Our modification also reduces the range of coefficients in the underlying cost function without increasing the number of qubits required. This method avoided the time-consuming preprocessing steps[17] to achieve comparable results, and could be combined with other ad-hoc function simplification methods to further reduce the number of qubits or other aspects. We finally tested both methods with results from experimental demonstrations using quantum annealing hardware.

Background

Quantum Annealing was introduced[18] to solve optimization problems using quantum fluctuations to transit to the ground state, compared to simulated annealing which uses thermal fluctuations to get to the global minimum. Quantum fluctuations such as quantum tunneling[21] provide ways of transitions between states. The transverse field controls the rate of the transition, as the role of temperature played in simulated annealing.

Farhi et al.[12] remodeled the procedure as Adiabatic Quantum Computing(AQC), which finds the energetic ground state of a problem Hamiltonian by adiabatically evolving the quantum state. If the system begins at the ground state, after the adiabatic evolution, the system will remain at the ground state. By combining the initial Hamiltonian H and the problem Hamiltonian H linearly, the system could be defined as a time-dependent Hamiltonian , where the duration T defines the time-scale for evolution and controls the rate at which the time-dependent Hamiltonian changes, the initial Hamiltonian H is of the form over L qubits with the Pauli operator defining the x-basis of the i-th qubit. The problem Hamiltonian H is in the form of Ising model over L qubits aswhere defines the z-basis for the i-th qubit and the local fields h and couplings J define the problem instance.

Computation within the AQC model evolves the L-qubit quantum state under the time-dependent Hamiltonian H(t) according to the Schrödinger equation , where |ψ(t)〉 is the state of the system at time t ∈ [0, T] and we will set ℏ = 1. Let |ϕ(t)〉 be the i-th instantaneous eigenstate, such that H(t)|ϕ(t)〉 = E(t) |ϕ(t)〉, and let the initial state of the system be the ground state at time t = 0, such that |ψ(0)〉 = |ϕ0(0)〉. According to the adiabatic theorem[22], the system state will remain in the instantaneous ground state of the time-dependent Hamiltonian provided the evolution is sufficiently slow to prevent excitations to higher-lying states. Under these idealized adiabatic conditions, the system will evolve into the energetic ground state of the problem Hamiltonian as |ψ(T)〉 = |ϕ0(T)〉. This prepared quantum state of L qubits is then measured to generate a classical string of L bits that represents the solution to the encoded factorization problem.

Several practical considerations limit the applicability of the AQC model for solving optimization problems. Foremost is the requirement that changes in the quantum state must be adiabatic, i.e, slow, relative to the internal timescales of the instantaneous Hamiltonian H(t). Theoretical analyses of this requirement provide a best lower bound on the time T as O(Δ−2), where Δ is the minimum energy gap within the instantaneous eigen spectrum of H(t)[13]. However, the minimum spectral gap is dependent on the specific instances of the initial and final Hamiltonians and the interpolation between them. A priori knowledge of spectral gap information would provide a significant insight into the underlying optimization problem, if not the solution directly, and therefore is an impractical expectation for a computational method. In addition, the pure-state model for AQC fails to account for finite-temperature effects observed in actual hardware as well as unexpected environmental coupling, unpredictable control noise, unwanted crosstalk, and other imperfections.

Practically quantum annealing relaxes the AQC with guarantee that the observed final state corresponds to the energetic ground state. However, the probability to observe the energetic ground state may be reduced due to physical noise and non-adiabatic dynamics. The resulting error rate p(L) characterizes the quantum annealing dynamics, which is most accurately modeled as an open quantum system of L-qubits evolving the presence of an uncontrolled environment. The quantum annealing model is therefore more robust to the above practical considerations but it is necessarily a probabilistic computational model. Statistical sampling of a quantum annealing computation is always necessary to gather confidence in the observed result. Quantum annealing may also be interpreted as a meta-heuristic for managing noisy AQC computation, whereby the aggregate likelihood of success p is determined by the number of samples S as p = 1−(p). The number of samples necessary to achieve a desired probability of success is therefore S ≥ log(1−p)/log(p). For a fixed annealing duration T and probability of success p, we may expect the probability of error p to increase as the size of the system increases, i.e., as L increases. The rate at which the sample number S increases with system size plays an important role in determining the computational complexity of using the quantum annealing model. For example, an error rate that increases exponentially with system size, i.e., p∝ exp(L), yields a sampling rate that increases linearly.

A related practical consideration is the resource efficiency with which quantum annealing can be implemented. Specifically, the number of qubits necessary to realize the problem Hamiltonian H influences not only the number of samples required but also the feasibility of demonstrating the method on available hardware. The most general case of L-qubit Hamiltonian may include all-to-all connectivity, whereby each qubit must interact with every other qubits. However, most of the existing hardware does not permit such connectivity directly, and methods for realizing implicit connections have been developed[23,24]. In our implementation of integer factorization using the Ising Hamiltonian, it is necessary to compose a problem Hamiltonian in terms of pairwise interactions, and we develop an efficient transformation of the factoring problem Hamiltonian into pair-wise coupling.

Methods

We describe two methods for implementing integer factorization within the quantum annealing model. We found these two corresponding Hamiltonians H to encode the factors of an input integer N, such that the energetic ground state corresponds to factorization of the input. The first is a direct method to compute the factors of N = pq by constructing the associated optimization problem as an Ising Hamiltonian. The second method is based on the modified multiplication tables to translate the problem into the Ising Hamiltonian. We test our methods using experimental quantum computing hardware appropriate for quantum annealing. The D-Wave 2000Q processor natively implements an Ising model Hamiltonian and provides programmable control over the parameters h and J as well as the annealing duration T.

Direct Method

Our direct method factors N = pq, where p and q are prime numbers. We set and . Because p and q are prime numbers, we use the binary representations and , where l1 > l2 and x ∈ {0, 1} for i = 1 to l1 + l2−2. We define the cost function , and explicit multiplication of the binary representations for p and q yields a sum of binary products. We reduce the resulting 3-local terms to 2-local terms using the following equivalence[24]: for x, y, z ∈ {0, 1}, xy = z iff xy−2xz−2yz + 3z = 0, and xy ≠ z iff xy−2xz−2yz + 3z > 0. Therefore,andthus, the x1x2x3 term may be transformed to quadratic form by replacing x1x2 with x4 plus a constraint in the form of a penalty term:

By introducing a new variable and adding the penalty term, we are able to transform 3-local terms to 2-local terms.

For integer factorization, we require auxiliary variables to form a quadratic cost function, and when l1 = l2 = l this number is l × (l−1). Counting the variables to denote the factors themselves, the quadratic function requires L = 2 × (l−1) + l × (l−1) = (l + 2) × (l−1) binary variables in total. Since , . We could also let when lengths of p and q are prefixed.

We illustrate this direct method of factorization for the case of N = 15. Because , we define p = (x11)2 and q = (x2x31)2. The objective function f(x1, x2, x3) = (N−pq)2 may then be reduced to the 3-local form:

We reduce the 3-local terms to 2-local terms using the method described above to obtainwhere

This result is a quadratic unconstrained binary optimization (QUBO) problem over L = 4 variables that may be transformed into an equivalent Ising Hamiltonian as defined in Eq. 1 by identifying the binary variable x with the i-th spin state s = 2x − 1. For N = 15, the local fields h and couplings J of the Ising problem Hamiltonian are then determined to beand

It is notable that the L × L coupling matrix J is generally dense on the upper triangle, indicating that L(L−1)/2 couplings are necessary. Similarly Ising parameters may be generated for other integers N by the appropriate quadritization into a QUBO problem and then reduced to Ising form.

Modified Multiplication Table Method

The second method is based on modified multiplication table. It reduces the range of Ising parameter values used as coefficients for the local fields and couplings. At the meantime, it considers to use a smaller number of carry variables without complicated preprocessing.

The modified multiplication table method uses local minimizations over the products of individual binary substring bits representing the integers p and q. It divides the multiplication table into several blocks, and considers each block individually. We could also choose the size of each block to get the desired range of parameters and the number of variables, or make a balance between them. A detailed analysis of the range of coefficients is shown in the last section of the supplemental material. Note that the modified multiplication table method does not eliminate the need for quadratization of 4-body and 3-body product terms and auxiliary variables are required. However, this approach does reduce the number of these higher-order terms compared to the direct method, which makes it possible to embed larger problem sizes on currently available quantum hardware.

We describe this method by an illustrative example of N = 143, for which p = 13 and q = 11. Past approaches for integer factorization constructed a system of equations from each column or part of each column in the multiplication table[17,25,26], which accounts for a carry bit or several carry bits for each part. In our approach, we divide the multiplication table into blocks so that it is only needed to use carries between blocks. This greatly reduced the number of carries, thus the total number of variables reduced too.

As shown in Table 1 for N = 143, we introduce two sets of carry bits. We denote them using c ∈ {0, 1}, and the two-bit numbers (c2c1)2 = c2 × 2 + c1 and (c4c3)2 = c4 × 2 + c3 represent the carry bits for each of the divided columns in the table. Note that the columns are composed along two-bit domains, so that addition within each block is over two-bit numbers. But the sums are over four-bit numbers. The resulting block system of equations derived from Table 1 is

Table 1

Multiplication table for 13 × 11 or 11 × 13 = 143 in binary.

	27	26	25	24	23	22	21	20
p					1	p 2	p 1	1
q					1	q 2	q 1	1
					1	p 2	p 1	1
				q 1	p 2 q 1	p 1 q 1	q 1
			q 2	p 2 q 2	p 1 q 2	q 2
		1	p 2	p 1	1
carries		c 4	c 3	c 2	c 1
p × q = 143	1	0	0	0	1	1	1	1

Because this modified multiplication table method calculates carries only for each block, it avoids requiring carry bits for each column in the multiplication table. This reduces the overall complexity of the computation by reducing the number of carry bits as well as the number of couplings between bits. In the limit of a single column per block, the conventional multiplication table is recovered, while in the limit of a single equation the direct method is recovered. Instead of making the sum of each column equal to every each bit of the number to be factored as in a conventional multiplication table, we make each block of the multiplication table equal to the corresponding block of the number to be factored. As shown in the appendix material, the equations for these blocks may be reduced to the non-negative cost function

This form may be expanded and further simplified using the property , while the remaining cubic and higher-order terms like c1p1q1 and p1p2q1q2 can be reduced to quadratic form by introducing auxiliary variables. In particular, we note that the quadratization of the negative term is similar to the position term, e.g.,

as detailed in the appendix material, the conversion to QUBO form leads to the parameters for the Ising Hamiltonian. For N = 143, this yields the local fieldsand the upper triangular coupling matrix

Our approach requires a decision to partition the columns of the multiplication table into blocks, and this choice must balance the number of unknown variables (carries) against the range of coefficients in the problem Hamiltonian. We illustrate this choice for the factorization of biprimes 59989 and 376289 in Tables 2 and 3. Our approach is to set the bit-length of the carry variable for each block based on the largest possible number of that block(the right neighboring columns). For example, the maximum carry for the right-most block in the multiplication table of N = 59989 is 3 which requires two bits to represent. Thus, the bit-length of the carry variable for this block is 2, i.e., (c2c1)2. Similarly, for N = 376289, the bit-length of the carry variables for the right-most block is 3, while the bit-length of the carry variable for the third block is 4. Because this bit-length is larger than the size of the fourth block, which has a bit-length of 3, the most significant bit of the carry is included in the neighboring block, i.e., the fifth block in this example.

Table 2

Multiplication table for 251 × 239 = 59989 in binary.

	215	214	213	212	211	210	29	28	27	26	25	24	23	22	21	20
p									1	p6	p5	p4	p3	p2	p1	1
q									1	q6	q5	q4	q3	q2	q1	1
									1	p6	p5	p4	p3	p2	p1	1
								q1	p6q1	p5q1	p4q1	p3q1	p2q1	p1q1	q1
							q2	p6q2	p5q2	p4q2	p3q2	p2q2	p1q2	q2
						q3	p6q3	p5q3	p4q3	p3q3	p2q3	p1q3	q3
					q4	p6q4	p5q4	p4q4	p3q4	p2q4	p1q4	q4
				q5	p6q5	p5q5	p4q5	p3q5	p2q5	p1q5	q5
			q6	p6q6	p5q6	p4q6	p3q6	p2q6	p1q6	q6
		1	p6	p5	p4	p3	p2	p1	1
	c11	c10	c9	c8	c7	c6	c5	c4	c3		c2	c1
	1	1	1	0	1	0	1	0	0	1	0	1	0	1	0	1

Table 3

Multiplication table for 659 × 571 = 376289 in binary.

218	217	216	215	214	213	212	211	210	29	28	27	26	25	24	23	22	21	20
									1	p8	p7	p6	p5	p4	p3	p2	p1	1
									1	q8	q7	q6	q5	q4	q3	q2	q1	1
									1	p8	p7	p6	p5	p4	p3	p2	p1	1
								q1	p8q1	p7q1	p6q1	p5q1	p4q1	p3q1	p2q1	p1q1	q1
							q2	p8q2	p7q2	p6q2	p5q2	p4q2	p3q2	p2q2	p1q2	q2
						q3	p8q3	p7q3	p6q3	p5q3	p4q3	p3q3	p2q3	p1q3	q3
					q4	p8q4	p7q4	p6q4	p5q4	p4q4	p3q4	p2q4	p1q4	q4
				q5	p8q5	p7q5	p6q5	p5q5	p4q5	p3q5	p2q5	p1q5	q5
			q6	p8q6	p7q6	p6q6	p5q6	p4q6	p3q6	p2q6	p1q6	q6
		q7	p8q7	p7q7	p6q7	p5q7	p4q7	p3q7	p2q7	p1q7	q7
	q8	p8q8	p7q8	p6q8	p5q8	p4q8	p3q8	p2q8	p1q8	q8
1	p8	p7	p6	p5	p4	p3	p2	p1	1
	c14			c10	c9	c8	c7	c6	c5	c4	c3	c2	c1
		c13	c12	c11
1	0	1	1	0	1	1	1	1	0	1	1	1	1	0	0	0	0	1

We estimate the number of variables needed to construct the Ising Hamiltonian that encodes the factorization problem for N = 59989 and 376289. This requires quadratization of the resulting systems of factoring equations followed by reduction to the Ising form, exactly as discussed explicitly above. For the case of N = 59989, we have l1 = l2 = 6, and therefore 12 variables are required to represent the factors themselves plus 11 variables to denote the carries, while 36 auxiliary variables are required for quadratization of the factoring equations. The total number variables is 59. For N = 376289, we have l1 = l2 = 8 with 14 carries and 64 auxiliary variables. As noted above, sometimes vary bits in the multiplication table will overlap, as is the case for column 214 of 376289 shown in Table 3. In such circumstances, we just simply add these carries in the table and then use the same method as before to find the corresponding Ising Hamiltonian. Thus, this problem Hamiltonian requires 94 qubits.

Generally, for factoring a biprime number, we use approximately log(N) binary variables to encode the integer factors and about log(N) binary variables to denote the carries where N is the number to be factored. An additional log2(N)/4 auxiliary binary variables are required for quadratization. Therefore, a total of approximately log2(N)/4 binary variables are required to represent the problem Hamiltonian and, consequently, a similar number of qubits must be available within hardware. As a point of reference, applying this method to the current factoring record for RSA-768 would require approximately 147456 qubits.

Experiments

We tested both the direct and modified multiplication table methods using quantum annealing hardware from D-Wave Systems. The D-Wave System hardware consists of a programmable platform of integrated superconducting flux qubits designed to operate within the quantum annealing model. In particular, the hardware system accepts as input a problem Hamiltonian H presented in Ising form and the parameters h and J. The hardware also enables interpolation from the starting Hamiltonian H and the final Hamiltonian at a rate controlled by the annealing duration T. Measurements of the resulting quantum state are performed in σ basis for each qubits.

While the latest 2000Q system contains up to 2048 qubits arranged in a connectivity pattern expressed as a 16-by-16 Chimera graph, this sparse connectivity pattern requires additional resources to ensure the required interactions between the logical variables defining the problem Hamiltonian. This is accomplished by embedding the problem Hamiltonian into the hardware graph while maintaining the logical form of the cost function[27-29]. From the coupling matrix J, we define a graph G that represents the variables as vertices and non-zero coupling as edges. A minor embedding of G(V, E) into a hardware graph G′(V′, E′) is defined by a mapping such that each vertex v ∈ G is mapped to a connected subtree T of G′ and if (u, v) ∈ E then there exist i, i ∈ G′ such that i ∈ T, i ∈ T and (i, i) ∈ E′. If such a mapping ϕ exists between G and G′, we say G is a minor of G′ and we use G ≤ G′ to denote such relationship.

The logical parameters for local fields and couplers should also be considered. In the parameter setting problem[30], we assign each node and each edge in the minor embedding graph such that: (1) for each node in the tree T expanded by the same vertex i, its value satisfies , (2) for each edge in the tree T expanded by the same vertex i, the value needs to be large enough to make sure all physical qubits that correspond to the same logical qubit to be of the same value and (3) for each edge in the minor embedding graph which is in the original graph, we could use the same J value.

Our programmed implementation of these methods were written in C/C++ using the XACC programming framework[31]. XACC enables integration of the D-Wave solver application programming interface (SAPI) using a directive-based programming model. Pre-processing of the input N generated the Ising parameters for a logical Hamiltonian that was then embedded into the hardware graph structure. For minor embedding, we use the sapi_findembedding method included in the D-Wave 2000Q control software, SAPI version 3.0. This embedding methods is based on a randomized algorithm from Cai, Macready and Roy[32], while access to these methods were managed using the XACC dwsapi-embedding plugin[31]. The corresponding biases and couplings for the embedded problem were generated using the logical Ising parameters. The output of the embedding was a program implementation of the physical Ising model that was submitted for execution on the D-Wave processor. Additional parameters for the execution included the number of samples S and the annealing duration T. The default annealing schedule for the 2000Q was used for all executions. The output from each of the S executions was a measured binary string designating ±1 values for each spin variable. The number of samples was S = 10,000. Each returned string was then classified according to the corresponding energy for the physical Ising model and subsequently decoded into the factors p and q. A histogram of all solutions returned for a specific annealing time was recorded.

Figure 1(a,b) shows the frequency of each decoded solutions to the factorization problem for N = 15 and 21 using the direct method. These observed solutions are decoded using the inverse of the embedding with majority vote used to resolve any ambiguity in results. The plot presents the decoded results in order of lowest energy to highest energy (left to right). For these two examples, the lowest energy solution corresponds to the correct factors. In addition, several other computed solutions decoded into the correct factors as the associated errors were resolved by the decoding method. For example, several solutions are labeled as (3,7) because the observed bit strings corresponded to high-energy states before decoding. Only the first (leftmost) corresponds to the lowest energy state. The others were higher energy solutions, thus can’t be counted as the correct solution.

Figure 1

Experimental results on D-Wave machine: rates of getting different solutions. For example, the (3, 5) in the x-axis denotes the factorization of 15 is 3 multiplied by 5, the number in y-axis denotes the rate to get this factorization.

Using the modified multiplication table method for factoring 143, we embed the problem Hamiltonian to D-Wave machine using the following method. Suppose n qubits are needed in the Hamiltonian, we divide n into groups. For each group, we use 4 copies of the nodes with each . We assign each edge in the tree T the negative number with largest absolute value to make it a penalty term. This method guarantees the nodes correspond to the same original qubit have the same value. We assign each edge corresponding to the original edge in the problem graph the same J value. The embedded graph to D-Wave machine is in Fig. 2.

Figure 2

Embedding the factoring instance N = 143 to Chimera graph. The nodes with the same color denote the same original qubit, with their connected lines corresponding to strong couplings. The left footnotes refer to which spin the node was embedded.

The results graph are shown in Fig. 3. The final state of the part of the system which represents the problem solution will be |1 −1 −1 1〉 or |−1 1 1−1〉 with relatively high probability, which corresponds to solutions p = (1p2p11) = (1101)2 = 13, q = (1q2q11) = (1011)2 = 11 or p = 11, q = 13.

Figure 3

Experimental results on D-Wave machine: final ground state of factoring 143. Nodes colored red denote +1, nodes colored blue denote −1. (a) This graph shows s1 = 1, s2 = −1, s3 = −1, s4 = 1 which means p = 1101, q = 1011. (b) This graph shows s1 = −1, s2 = 1, s3 = 1, s4 = −1 which means p = 1011, q = 1101.

For factoring larger numbers like 376289 using D-Wave, we embed the Hamiltonians into the Chimera hardware using the predefined function find_embedding and embed_problem available in the vendor’s software developer packages, because the problem graph can’t be embedded directly like the case 143 shown above which has 12 qubits, due to the size limitation of the current Chimera hardware. For now, the largest number experimented for finding an embedding in D-Wave 2000Q is 249919 which equals to 509 × 491. It uses 74 qubits in the final Ising Hamiltonian, and embeds to 1803 physical qubits in Chimera graph.

Conclusions

We have presented two general methods for factoring integers using quantum annealing for optimizing a cost function that is reduced to an Ising Hamiltonian. Both methods requires qubits in total, where N is the number to be factored. The novelty of our demonstration of quantum annealing for prime factorization is based on the reduction in quantum resources required to execute factoring and the experimental verification of the algorithmic accuracy using currently available hardware. As a proof-of-concept, we have demonstrated these methods by factoring integers using the D-Wave 2000Q quantum annealing hardware, but these methods may be used on any other quantum annealing system with a similar number of qubits, qubit degree of connectivity, and hardware parameter precision. Assuming that quantum annealing hardware systems will continue to grow both in the number of qubits and bits of precision capabilities, our methods offer a promising path toward factor much larger numbers in the future. It is also good to combine our method with other ad hoc methods to achieve significantly better performances for specific numbers.

Finally, we note that while our demonstrations of factoring have made use of currently available quantum annealers, there is an outstanding question regarding the asymptotic complexity for this approach. It is well known that algorithmic complexity within the AQC model depends on the minimum spectral gap between the ground and first-excited states of the underlying time-dependent Hamiltonian. Attempts to classify the complexity of the spectral gap with respect to system size have not yet succeed and, indeed, Cubitt, Perez-Garcia, and Wolf have proven that the problem of claiming a Hamiltonian has a gap is undecidable in general[33]. Nonetheless, there is hope that our resource-efficient algorithms may find use in pre-processing potential factors for noisy factorization algorithms, e.g., as suggested by Patterson et al. within the context of RSA[34].

Supplemental material