A Factorisation Algorithm in

Adiabatic Quantum Computation

Tien D. Kieu∗
Centre for Quantum and Optical Science
Swinburne University of Technology, Victoria, Australia
(Dated: February 28, 2019)

Abstract
The problem of factorising positive integer N into two integer factors x and y is first reformu-
lated as an optimisation problem over the positive integer domain of either of the Diophantine
polynomials QN (x, y) = N 2 (N − xy)2 + x(x − y)2 or RN (x, y) = N 2 (N − xy)2 + (x − y)2 + x, of
√
each of which the optimal solution is unique with x ≤ N ≤ y, and x = 1 if and only if N is
prime. An algorithm in the context of Adiabatic Quantum Computation is then proposed for the
general factorisation problem.

∗
[email protected]; [email protected]

1
Factoring an integer into its prime constituents has attracted much interest since the
advance of the RSA public-private key encryption [1]. It is suspected that factorisation is
NP-intermediate, that is, in the NP class but may be not quite NP-complete. While there
does not yet exist any polynomial-time algorithm for the problem on a classical/Turing
computer, the discovery of Shor’s quantum algorithm with quantum circuits [2] has been
one of the main motivations for research into quantum computation and building of quantum
computers.
In this paper we also consider the factorisation problem in the realm of quantum com-
putation but with adiabatic processes, in complementary addition to the computation with
quantum circuits. Also recently, the authors of [3] have considered the problem with quan-
tum annealing. We first reformulate in the next section the factorisation into two integer
factors as an optimisation of some corresponding Diophantine polynomials over the integer
domain. The optimisation could also be repeatedly applied to any integer having more than
two prime factors. Based on this reformulation, we then present an algorithm in the context
of AQC (Adiabatic Quantum Computation) for the general factorisation problem. Following
that are some numerical illustrations of the algorithm and discussion on the lower bound of
the computing time with the help of an energy-time uncertainty relation. The paper is then
concluded with some remarks.

FACTORISATION AS AN OPTIMISATION PROBLEM

We first consider the problem of factorising a natural integer N into two integer factors x
and y. We propose that this problem can be reformulated as an optimisation problem over
the integer domain of the following Diophantine polynomial

min QN (x, y) ≡ N 2 (N − xy)2 + x(x − y)2 .

(1)
x,y∈N

Without the second term in QN above, the optimal solutions contain the trivial unity factor.
We could eliminate this triviality with this second term replaced by (x − y)2 , but there still
remains a symmetry between x and y. We now show that the polynomial QN (x, y) in (1) is
minimised if and only if xy = N where x ≤ y and x is nearest to, but not exceeding, the
√
integer part of N . The optimal solution thus has no trivial factor x = 1 unless N is prime.
For xy = N the first term on the rhs of (1) vanishes and the remaining second term is
√
obviously smaller for 1 ≤ x ≤ N ≤ y. In this case, the second term is a one-variable

2
function x(x − N/x)2 , and it is not difficult to see that this term is a decreasing function for
√
1 ≤ x ≤ N . We then have

2
QN (x, y)|N =xy = QN (x, 1/x) ≤ max √ x(x − N/x) ,
1≤x≤ N

≤ x(x − N/x)2 x=1 ,

≤ (N − 1)2 . (2)

The second line results from the fact that x(x − N/x)2 obtains its maximum value at the
√
boundary where x = 1. The closer x is to the integer part of N the smaller the value of
x(x − N/x)2 .
Now if xy 6= N then obviously (N − xy)2 ≥ 1 and the first term of (1) is consequently
not less than N 2 . Thus

QN (x, y)|N 6=xy ≥ N 2 + x(x − y)2 ,

≥ N 2. (3)

Combining (2) and (3), we have for N > 1

QN (x, y)|N 6=xy > QN (x, y)|N =xy . (4)

The second term x(x − y)2 of QN (x, y) in (1) is thus designed to introduce an asymmetry
between x and y, enforcing x ≤ y at the optimum value, and to eliminate the trivial unity
factor of N unless N is prime. The optimisation problem (1) thus has a unique solution
√
(x, y) for xy = N with x ≤ N ≤ y. In general, x is the integer factor nearest to, but not
√
exceeding, N . Consequently, x = 1 if and only if N is prime.
The optimisation of QN (x, y) can thus also determine the primality of N .
For N having more than one pair of two prime factors, the optimisation of QN (x, y) still
admits only one unique solution to yield a non-trivial factoring pair of N . For N having
more than two prime factors we can apply the optimisation repeatedly for the successive
factors to obtain the constituent prime factors.
It can be shown that the objective function in (1) is non-convex as depicted in FIG. 1
for a slice of the polynomial for N = 6 along the plane y + x = 5. (In general, determining
convexity for multivariate quartic polynomials is NP-hard.)
It is also known that the problem of minimising a general non-convex degree four poly-
nomial over a two-dimensional convex polygon like that in (1) is classically NP-hard – even

3
FIG. 1. Illustration of non-convexity of QN (x, y) along the plane x + y = 5 for N = 6.

though there exist polynomial-time approximations for some lower and upper bounds of the
optimal values; see, for example, [4] and references therein. The condition of convexity or
lack of it is an important element here; were QN (x, y) convex then we would already have a
polynomial time algorithm for the factorisation!
In the next section we propose an AQC algorthim for optimising this type of polynomials
over the integer domain.

AN AQC ALGORITHM FOR THE FACTORISATION PROBLEM

AQC [5, 6] makes use of some appropriate time-varying Hamiltonians and an initial
state. The computation starts with the readily constructible ground state |gI i of an initial
Hamiltonian HI which is then adiabatically extrapolated to the final target Hamiltonian
HP whose ground state |gi encodes the desirable solution of the problem and could be

4
then obtained with reasonably high probability. The interpolation between HI and HP is
facilitated by a time-dependent Hamiltonian in the time interval 0 ≤ t ≤ T ,

H(t) = f (t/T )HI + g(t/T )HP , (5)

either in a temporally linear manner (that is, f (t/T ) = (1 − t/T ) and g(t/T ) = t/T ); or
otherwise with f (0) = 1 = g(1) and f (1) = 0 = g(0). We also assume that both f and g
are continuous, and g is semi-positive for all t ∈ [0, T ]. Such a time evolution is captured
by the following Schrödinger equation in the time interval [0, T ]

i∂t |ψ(t)i = H(t) |ψ(t)i, (6)

|ψ(0)i = |gI i.

As the rate of the evolution of the Hamiltonian approaches zero, the end state |ψ(T )i asymp-
totically converges to the target state |gi as asserted by the quantum adiabatic theorem [7].

The initial Hamiltonian HI and its ground state |gI i

Corresponding to the variable x, we introduce the operators a†x and ax , which respectively
create and annihilate the number states |nx i, for nx = 0, 1, . . .:

ax , a†x = 1,

(7)
√
a†x |nx i = nx + 1|nx + 1i,
√
ax |nx i = nx |nx − 1i, nx ≥ 1,
ax |0x i = 0,

and the number operators n̂x are constructed in the usual way:

n̂x = a†x ax ,

n̂x |nx i = nx |nx i. (8)

Similarly for the variable y with a†y , ay and n̂y .

We propose to start the AQC with the following initial Hamiltonian HI :

HI = (a†x − θx∗ )(ax − θx ) + (a†y − θy∗ )(ay − θy ), (9)

5
for some c-numbers θx and θy . The states in the total Hilbert space can be now decomposed in
terms of |nx i⊗|my i, the direct products of the two bases of number states. This Hamiltonian
admits the readily constructible direct product of coherent states as the ground state |gI i,

|gI i = |θx i ⊗ |θy i, shortened as |θx θy i, (10)

with |θi the canonical quantum coherent state

∞ n
|θ|2 X θ
− 2
|θi = e √ |ni,
n=0 n!
where:

a |θi = θ |θi.

The target Hamiltonian HP

Building on (1), we introduce the target Hamiltonian HP ,

HP = N 2 (N − n̂x n̂y )2 + n̂x (n̂x − n̂y )2 , (11)

which has a unique ground state, i.e. no degeneracy. This ground state would give us the
solution for the desirable optimisation problem (1) over the integer domain, thus solving the
factorisation for the integer N .
Our proposed AQC may be physically realised by quantum optics or other means, with
some suitably chosen end time T . Then the solution factors for N would be obtained from
the final state |ψ(T )i by its measurements in the number state basis. Even though such
measurements only yield probabilities, any factor candidates so obtained could be verified
straightforwardly and efficiently by a multiplication of nx and ny , or by a division of N by
nx or by ny .
After the release of a draft of this present work, the algorithm presented herein has been
physically implemented with qubits on D-Wave Systems [8].

NUMERICAL SIMULATIONS

For the purpose of illustration, we simulate the AQC algorithm for N = 6 by numeri-
cally solving the Schrödinger equation (6) in a suitably truncated Hilbert space of dimen-
√
sions 23 × 23, with linear extrapolation in (5) and parameters θx = θy = 4 N (so that
hθx θy |n̂x n̂y |θx θy i = N ).

6
FIG. 2. Spectral flow of H(t) with linear extrapolation in (5) for N = 6. The two lowest branches
respectively correspond at t = T to |2x 3y i and |3x 2y i, both of which are non-trivial factor pairs.

FIG. 2 depicts the spectral flow of H(t). The two lowest branches respectively correspond
at t = T to |2x 3y i and |3x 2y i, both of which are non-trivial factor pairs of N = 6. Note the
relatively wide energy gaps between these solution states and the non-solution nx my 6= N .
This would be helpful in shortening the running time with an appropriately chosen g(t/T )
for the AQC algorithm, in accordance with (13) in the next Section.

FIG. 3 shows the probabilities |hψ(T )|nx my i|2 versus T for the two largest probabili-
ties. In all instances of T there, the number state |2x 3y i has the maximum probability
among all the number states. This probability increases with T , the inverse of the evolution
rate, in agreement with the quantum adiabatic theorem [7]. The next highest probability
corresponds to the other solution |3x 2y i.

We note that the increased probability for |3x 2y i at T = 100 compared to that at T = 30

7
FIG. 3. Probability |hψ(T )|nx my i|2 versus T in the simulated AQC for the factorisation of N = 6
with a suitably truncated Hilbert space of dimensions 23 × 23. Only the two largest probabilities
are shown.

is an artifact of our truncation of the Hilbert space. With a truncation at |nmax i such that
a† |nmax i = 0, the commutation relation (7) is violated

hnmax |aa† |nmax i = 0,

hnmax |a† a + 1|nmax i = nmax + 1. (12)

The effect of such a truncation artifact becomes more prominent for larger T , unless we
increase nmax accordingly.

8
COMPUTATIONAL COMPLEXITY

We have presented elsewhere [9] a necessary condition for a lower time limit required in
a general AQC for an initial state to evolve into an orthogonal state under the dynamics of
H. Namely, it is necessary that the evolution time cannot be less than T⊥ ,
!
1
T⊥ ∼ O R1 , (13)
∆I EP 0 g(τ )dτ
where ∆I EP is the energy spread of the initial state |gI i in terms of the target Hamiltonian
HP ,
q
∆I EP ≡ hgI |HP2 |gI i − hgI |HP |gI i2 . (14)

It is important to note that only the initial eigenstate |gI i (which is of course time-
independent), and neither the instantaneous eigenstates nor the full time-dependent wave
function at any other times, is required for the time condition (13). This hallmark of our
results in [9] enables their wider applicability and usefulness.
We note that T⊥ is particularly a function of the parameters |θ|,

T⊥ = T⊥ (|θ|). (15)

The parameters θ’s give our proposed AQC algorithm some advantage that is denied or
not evident elsewhere, and this could be exploited in order to reduce the lower time limit
T⊥ for the computation to be O(1) or even less! But this reduction would at the same
time incur an appropriate increase in the energy cost Ecost = max0≤t≤T hgI |H(t)|gI i. For
instance, for θ ∼ O(1) the cost Ecost would be, as evident from (11), of order O(N 4 ) – that
is, exponentially in the number of bits in the binary representation of N .
It is not that surprising that we could reduce the computation time with more energy
resources. Other examples of this delicate balance between the energy required and the
lower time limit for a well known quantum algorithm are given elsewhere [9, 10].

ANOTHER DIOPHANTINE POLYNOMIAL

We could also employ another Diophantine polynomial RN (x, y), in addition to QN (x, y)
of (1), to construct HP in the search to find a pair of factors for N ,

min RN (x, y) ≡ N 2 (N − xy)2 + x(x − y) + x .

(16)
x,y∈N

9
The proof that the optimisation of RN (x, y) could yield a non-trivial pair of factors is similar
to that for (1).
For xy 6= N then obviously (N − xy)2 ≥ 1 and the first term of (16) is consequently not
less than N 2 ,

RN (x, y)|N 6=xy ≥ N 2 . (17)

For xy = N the first term on the rhs of (16) vanishes and the remaining terms, for
√
x ≤ N ≤ y,
2
N
RN (x, y)|N =xy = RN (x, 1/x) = x − +x (18)
x
√
form a convex function (i.e. having non-negative second derivative) for x ∈ [1, N ]. It then
follows that such a function is bounded by its maximum value obtained at one of the end
points

RN (x, y)|N =xy = RN (x, 1/x) ≤ max RN (x, 1/x)|x=1 , RN (x, 1/x)|x=√N ,
n √ o
≤ max (1 − N )2 + 1, N , (19)

≤ (1 − N )2 + 1,
≤ N 2,
≤ RN (x, y)|N 6=xy .

A pair of factors of N so obtained is not only unique but also non-trivial, because from (19)
it follows that x > 1 when N is not prime. Only when N is prime then we would have x = 1.

CONCLUDING REMARKS

We have mapped the factorisation problem into the optimisation over the integer domain
of some corresponding Diophantine polynomials. The optimal solutions are unique, and
contain the unity factor if and only if the integers to be factored are prime. Based on this, we
then propose an algorithm in the context of Adiabatic Quantum Computation for the general
factorisation. We also indicate that the lower bound on the running time could be reduced
at the expense of more energy input in order to carry out the physical computation. Some
numerical simulations of the algorithm for a simple case are also presented for illustration
purpose.

10
After the release of a draft of this present work, the algorithm presented herein has been
physically implemented with qubits on D-Wave Systems [8].
We are grateful for discussions with Peter Hannaford, Adolfo del Campo and Richard
Warren, who also brought the work [3] to our attention after the completion of this paper.
The publication fee for this work is covered by Swinburne University.

[1] R. L. Rivest, A. Shamir, and L. Adleman, Commun. ACM 21, 120 (1978).
[2] P. W. Shor, “Algorithms for quantum computation: Discrete logarithms and factoring,” Pro-
ceedings of the 35th Annual Symposium on Foundations of Computer Science, IEEE Computer
Society Press, Los Alamitos, CA (1994).
[3] S. Jiang, K. A. Britt, A. J. McCaskey, T. S. Humble, and S. Kais, “Quantum annealing for
prime factorization,” arXiv:1804.02733 [quant-ph] (2018).
[4] J. A. De Loera, R. Hemmecke, M. Köppe, and R. Weismantel, Mathematics of Operations
Research 31, 147 (2006), https://doi.org/10.1287/moor.1050.0169.
[5] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, Science 292, 472
(2001), arXiv:quant-ph/0104129.
[6] E. Farhi, J. Goldstone, and S. Gutmann, “A numerical study of the performance of a quantum
adiabatic evolution algorithm for satisfiability,” arXiv:quant-ph/0007071 (2000).
[7] A. Messiah, Quantum Mechanics (John Wiley & Sons New York, 1966).
[8] R. Warren and D. Dahl, private communication (2018).
[9] T. D. Kieu, “A new class of time-energy uncertainty relations for time-dependent hamiltoni-
ans,” arXiv:1702.00603 [quant-ph] (2017).
[10] T. D. Kieu, “The travelling salesman problem and adiabatic quantum computation: An algo-
rithm,” arXiv:1801.07859 [quant-ph] (2018).

11