# Optimizing SymPy Implementation of prime factorization in form of QUBO

I'm trying to reproduce a paper on Prime Factorization. This paper converts the problem into a QUBO form, which then we can map it to an Ising Minimizer. I've basically done everything, and I've obtained the final QUBO expression of page 3. My only problem is that my implementation gets slow for larger numbers really fast. The bottleneck of my code is in the reduced_f function. I also suspect that all of these operations are possible without symbolic calculations and only with matrix multiplication, but I don't know how to implement it like that.

from sympy import IndexedBase, expand, Indexed, Mul
import numpy as np
from tqdm import tqdm

x = IndexedBase('x')

def high_degree_f(N, l1 = 2, l2 = 3):
# p and q are binary numbers of length l1 and l2
p = 1
q = 1

for i in range(1, l1):
p += x[i] * 2**i

for idx, l in enumerate(range(l1, l1+l2 - 1)):
q += x[l] * 2**(idx+1)

f = (N - p * q ) **2
f = expand(f)
r = len(f.free_symbols) + 1

# replace x[i]**k by x[i] as x[i] is binary
for i in range(1, r):
for k in range(2, 4):
f = f.subs(x[i]**k, x[i])
return (f, p, q)

def out_degrees_more_than_2(f):
# Iterate over the terms in the expanded polynomial
for term in f.as_ordered_terms():
# Check if the term is a product
if isinstance(term, Mul):
# Count the number of Indexed instances in the term
variable_count = sum(isinstance(factor, Indexed) for factor in term.args)
# Check if there are more than two variables in the product
if variable_count > 2:
# Print the term
yield (variable_count, term)

def max_degree(f):
max_degree = 0
for (variable_count, term) in out_degrees_more_than_2(f):
max_degree = max(max_degree, variable_count)
return max_degree

def reduced_f(N, l1 = 2, l2 = 3):

(f, p, q) = high_degree_f(N, l1 = l1, l2 = l2)
number_of_variables = len(f.free_symbols) - 1

while max_degree(f) > 2:
for (variable_count, term) in tqdm(out_degrees_more_than_2(f)):
if variable_count == 3:

# extract the numerical coefficient of the term:
coefficient, v = term.as_coeff_Mul()
variables = v.args

new_term = coefficient* (variables[2] * x[number_of_variables+1] + 2 * ( variables[0] * variables[1] - 2 * variables[0] * x[number_of_variables+1] - 2 * variables[1] * x[number_of_variables+1] + 3 * x[number_of_variables+1]))
number_of_variables += 1

# substitute the term with the new term
f = f.subs(term, new_term)
elif variable_count == 4:
pass
else:
raise ValueError("Unexpected number of variables in term")
f = expand(f)

return (f, p, q)

N = 15
l1 = 2
l2 = 3

(f, p, q) = reduced_f(N, l1 = l1, l2 = l2)

print("N:\t\t", N)
print("p:\t\t", p)
print("q:\t\t", q)
print("Reduced f:\t", f)

N:               15
p:               2*x[1] + 1
q:               2*x[2] + 4*x[3] + 1
Reduced f:       200*x[1]*x[2] - 48*x[1]*x[3] - 512*x[1]*x[4] - 52*x[1] + 16*x[2]*x[3] - 512*x[2]*x[4] - 52*x[2] + 128*x[3]*x[4] - 96*x[3] + 768*x[4] + 196


expression from paper:

$$\begin{array}{rcl}f^{\prime} (x) & = & 200{x}_{1}{x}_{2}-48{x}_{1}{x}_{3}-512{x}_{1}{x}_{4}+16{x}_{2}{x}_{3}-512{x}_{2}{x}_{4}+128{x}_{3}{x}_{4}\\ & & -52{x}_{1}-52{x}_{2}-96{x}_{3}+768{x}_{4}+\mathrm{196,}\end{array}$$

But for N = 10403, l1 = 7, l2 = 7, it takes minutes to get the final expression. Can anyone suggest any improvement to my implementation?

## 1 Answer

To make it fast, I would guess a good option is to ditch sympy and simply use a dictionary and strings to create the QUBO, probably something like dic[(x_1,x_2,x_3)] = alpha_(x_1,x_2,x_3) which takes the (sorted) tuple of variables to corresponding coeffecient in the QUBO and then worry about making all tuples of length 2.

From a pragmatic point of view, I would define all the 'product variables' from $$p\times q$$ first and sub them into $$(N-pq)^2$$. This way, $$N-p\times q$$ is linear and and no further reductions will be needed. How these variables are chosen may change the final QUBO though which might not be what you are after.

• The whole point of using SymPy is to do the reduction. otherwise, the coefficients are pretty easy to obtain even in a matrix multiplication form (which is even faster than dictionary). I don't understand how your suggestion can skip the reduction. Please read the paper and if you think you can somehow improve the embedding, don't hesitate to write it as a detailed answer. Commented Jun 11 at 11:21
• My point is you are creating a 4th degree polynomial and then reducing it to a quadratic. The alternative I'm suggesting is to make the quadratic expression linear first and then square it, this should result in less terms to deal with. I'm not sure how sympy stores equations but if the long times are being caused by list manipulations. Moving to something like dictionaries and doing the reductions yourself could be a lot better.
– john
Commented Jun 11 at 11:32
• please write your approach as an answer. Commented Jun 11 at 16:16