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?