5
\$\begingroup\$

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?

\$\endgroup\$
0

0