HUBO and QUBO models for prime factorization

The security of the RSA cryptosystem is based on the difficulty of factoring a large number N into prime numbers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p$$\end{document}p and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q$$\end{document}q satisfying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=p\times q$$\end{document}N=p×q. This paper presents a prime factorization method using a D-Wave quantum computer that could threaten the RSA cryptosystem in the future. The starting point for this method is very simple, representing two prime numbers as qubits. Then, we set the difference between the product of the two prime numbers expressed in qubits and N as a cost function, and we find the solution when the cost function is minimized. D-Wave's quantum annealer can find the minimum value of any quadratic problem. However, the cost function must be a higher-order unconstrained optimization (HUBO) model because it contains second- or higher-order terms. We used a hybrid solver accessible via Leap, D-Wave’s real-time quantum cloud service, and the dimod package provided by the D-Wave Ocean software development kit (SDK) to solve the HUBO problem. We also successfully factorized 102,454,763 with 26 logical qubits. In addition, we factorized 1,000,070,001,221 using the range-dependent Hamiltonian algorithm.


Introduction
RSA cryptosystem is a popular public-key cryptographic algorithm for secure data transmission, first introduced in 1978 [1].Public and private keys form a pair in the RSA cryptosystem.Content encrypted with the public key can only be decrypted with the private key, and content encrypted with the private key can only be decrypted with the public key.The public key is a large bi-prime that can only be decrypted through prime factors held in the private key.The RSA cryptosystem is based on prime factoriaation [2], which finds the prime factors  and  such that  =  ×  for a large bi-prime .The algorithm is an interesting mathematical problem because the algorithm relies on principles of number theory.While the RSA cryptosystem is surprisingly simple, a large bi-prime is difficult to factoriae [3].In other words, the computational complexity of decryption is much higher than that of encryption, so security is excellent.However, with the development of quantum computers, quantum algorithms that reduce the computational complexity difference between encryption and decryption are proposed, threatening the safety of the RSA cryptosystem.
Peter Shor proposed the Shor algorithm [4], a powerful quantum algorithm that can factoriae an integer  in polynomial time, so can attack the RSA cryptosystem in polynomial time in 1994 [5].There have been many simulations on quantum computers for practical use-moreover, many attempts to run the Shor algorithm on quantum computer hardware.Geller et al. applied the Shor algorithm to factoriae 51 and 85 using Fermat numbers and eight qubits [6].However, there are still limitations to applying Shor's algorithm to a large number.On the other hand, quantum annealing, which underlies a -Wave's quantum computer, allows the factoriaation of even larger numbers.-Wave's quantum annealing can find the minimum value of the quadratic unconstrained binary optimiaation (QUBO) or Ising model [7].The objective formulated as the problem Hamiltonian of a -Wave quantum computer is achieved by defining the linear and quadratic coefficients of a binary quadratic model (BQM) that maps those values to the qubits and couplers of the QPU.ridi et al. presented a new autonomous algorithm that factoriaes all bi-primes up to 200,099 by reducing Hamiltonian's degree using Grobner bases [8].Jiang et al. developed a frame that transforms the prime factoriaation problem for arbitrary integers into the Ising model using ancillary variables and factoriaed the integer 376,289 with 94 logical qubits via a -Wave 2000Q System [9].ue to the limitation of the number of usable qubits of the -Wave 2000Q System, many researchers have studied factoriaing of a larger integer using the -Wave's hybrid quantum/classical simulator qbsolv.Peng et al. improved Jiang et al.'s work and factored 1,005,973 with 89 qubits by using qbsolv [10].Wang et al. successfully factoriaed all integers within 10,000 and demonstrated 1,028,171 with 88 qubits via qbsolv [11].Also, Wang et al. recently factoriaed 1,630,729 (an 11-bit prime factor multiplied by an 11-bit prime factor) [12].Unfortunately, a generaliaed method for factoring all integers into primes has not been developed due to the limitations of current quantum computers.Jiang et al. (2018) proposed a direct method for prime factoriaation.In this paper, we propose a direct method to formulate a higher-order unconstrained optimiaation (HUBO) model for prime factoriaation.To obtain the HUBO model, it is sufficient to express  and , which are prime factors of bi-prime N, as sums of qubits, respectively.The minimum value of the HUBO model can be obtained via a Leap's quantum-classic hybrid solver provided by the -Wave Ocean software development kit (S K) [13].The hybrid solver is an updated version of qbsolv.The hybrid solver can solve arbitrary problems formulated as quadratic models.Therefore, a process of converting the HUBO model to the QUBO model is required to solve a problem using a hybrid solver.The composed sampler in the -Wave Ocean's dimod package creates a BQM from a higher order problem.Hence, if these are used, prime factoriaation of large integers is possible using the -Wave Ocean S K. We demonstrated the factoriaation of 102,454,763 with 26 logical qubits using the -Wave hybrid solver as a successful example of our method.In addition, we applied the range dependent Hamiltonian algorithm that divides the domain to the new HUBO model.This algorithm is a method of finding a solution with a small number of qubits by dividing the domain into certain subintervals.The number of 1,000,070,001,221 was prime factoriaed using the HUBO model to which this algorithm was applied.

The least-squares problem for prime factorizations
The Ising model is a mathematical model for ferromagnetism in statistical mechanics.The energy Hamiltonian (the cost function) is formulated as follows: where  ⃗ = ( 1 , ⋯ ,   )  and   ∈ {+1, −1}.
QUBO is a combinatorial optimiaation problem in computer science.In this problem, a cost function f is defined on an  −dimensional binary vector space   onto ℝ.

𝑓(𝑞
where Q is an upper diagonal matrix,  ⃗ = ( 1 , ⋯ ,   )  , and   is a binary element of  ⃗.In this paper, the matrix Q is referred to as the QUBO matrix.The problem is to find  * ⃗⃗⃗⃗ which minimiaes the cost function  among vectors  ⃗.Since we have   2 =   , the cost function is reformulated as follows: In Q, the diagonal terms  , and the off-diagonal terms  , represent the linear terms and the quadratic terms, respectively.The unknowns of the Ising model σ and the unknowns of the QUBO model q have the linear relation Assume that the integer  is the product of two prime numbers  and  .To calculate  and  , consider the least square problem as below: Equation ( 1) satisfies the minimum value 0 when  = .To compute Eq. ( 5) conveniently, we apply 2-norm square to it.

HUBO model
While solving the binary least-squares problem,  and  are represented by combinations of qubits   ∈ {0,1}.The radix 2 representation of the positive integer values,  and , given by where the integer,  + 1, denotes the number of binary digits of  and  [14].We use qubits from  = 0 to use the same equation in the range dependent Hamiltonian algorithm in Chapter 2.4.
To derive a HUBO model, we insert Eq. ( 7) into Eq.( 6).This yields the summation terms of the first term in Eq. ( 6), as indicated below: In Eq. ( 8), since (  ) 2 =   , equation ( 9) can be obtained.The second term of Eq. ( 6) is calculated as follows: The HUBO model consists of the sum of Eqs.(10) and (12).And the global minimum energy we need to get is − 2 .

QUBO model
The HUBO model for prime factoriaations consists of quadratic, cubic, and quartic terms.To reformulate a non-quadratic (higher degree) polynomial into QUBO form, terms of the form, , where c is a real number, are substituted with one of the following quadratic terms [9]: For all , ,  ∈ {0,1},  can be transformed into a combination of linear and quadratic terms by adding a new qubit w to every cubic term.Similarly, Eq. ( 13) can be applied twice to convert quartic terms into QUBO formulations.The quartic terms with positive coefficients of this HUBO model are calculated as follows: For all  1 ,  2 ,  1 ,  2 , ∈ {0,1},  1  2  1  2 can be transformed into a combination of linear and quadratic terms by adding new seven qubits,  1 ,  2 , ⋯ ,  7 , for each quartic term.

HUBO model with the range dependent Hamiltonian algorithm
Recently, the range dependent Hamiltonian algorithm was proposed [15].This algorithm divides the domain into subregions that can be represented by the desired number of qubits.Applying this algorithm,  and  can be expressed as follows: and  ≈ ∑ 2   + −1 =0 +   (16) where   = 2  and  is an integer.
To derive a HUBO model, we insert Eq. ( 16) into Eq.( 6).This yields the summation terms of the first term in Eq. ( 6), as indicated below: In Eq. ( 8), since (  ) 2 =   , equation ( 9) can be obtained.The second term of equation ( 6) is calculated as follows: The HUBO model consists of the sum of Eqs. ( 19) and ( 21).And the global minimum energy we need to get is − 2 −   2   2 + 2    .

HUBO model with the qubit decomposition algorithm
Jun proposed the qubit decomposition algorithm [16].This algorithm is a method of finding Eq. ( 7) by continuously using  and  for a specific number of qubits.The advantage of this algorithm is that it can find a given equation with one qubit.We show here how  and  each find a solution with one qubit.This algorithm first finds the minimum energy from the highest order qubit for an exponent of 2. The next step is to find the minimum energy at each step by lowering the order one by one.The larger the energy difference for each step of a given problem, the easier it is to find a solution.

Experiments
We use a solver, dimod.ExactPolySolver().sample_hubo(), in the -Wave Ocean S K to test the HUBO model for the linear systems.This solver represents the energy for every possible number of cases each qubit can have.We tested it on 100 different numbers belonging to the first thousand prime numbers.In all cases, we were able to find pairs of  and  using the new HUBO model.The largest number we've tested is 102,454,763.The solver factored the number into two prime numbers, 10,111 and 10,133.The pseudo-code used in the test is shown in Algorithm 1.

Algorithm 1
The HUBO model for prime factoriaation etermine the number of qubits to be used for each prime number: NQ Represent  and  as the sum of qubits: Global minimum energy for solution: − 2 ▷ Quadratic Terms in Eqs.(10) and ( 12) ▷Cubic Terms in Eq.( 10) for  3 = 0:  − 1 do for  1 = 0:  − 2 do for  2 =  1 + 1:  − 1 do ▷Quartic Terms in Eq. ( 10) To calculate the first quadratic terms, we multiply 2  1 and 2  2 and then square them for  2 loops.The total number of flops is 2( + 1) 2 .In a similar way, the number of 3( + 1) 2 flops can be obtained for the second quadratic terms.The total number of flops for quadratic terms is 5( + 1) 2 .We can consider the change in terms of the first quadratic terms. 2 2( 1 + 2 ) varies from 4 0 to 4 2−2 .We can represent all of these terms as a sum of 2 + 1.Rather than directly calculating each coefficient, if we find the amount of change and put it into the coefficient appropriately, we can reduce the amount of calculation by about the square root.In a similar way, calculating the number of flops as a change amount for the cubic term and the quartic term requires 3 + 1 and 4 + 1 , respectively.We can reduce the amount of calculation once more here.The coefficients of the first quadratic terms and cubic terms are included in the coefficients of quartic terms.Therefore, the total optimiaed number of flops for the HUBO model is (4 + 3) + (2 + 3).
Both cubic and quartic terms of this HUBO model have positive coefficients.Therefore, the QUBO model for the factoriaation can be formulated using Eqs.(13) and (15).The pseudo-code is shown in Algorithm 2. In Algorithm 2, the cubic terms are divided into linear terms, quadratic terms, and constant terms by Eq. ( 13).The constant term generated here is not included in the QUBO model and appears in the reduced form of the global minimum energy.Reduction algorithm: Two cubic terms to the combination of linear and quadratic terms Two cubic terms with positive coefficients can be formulated into 10 linear and quadratic terms using Eq. ( 13).In the process of converting the two cubic terms into the form of the QUBO model, two new qubits are added.Eq. ( 15) is used to formulate each quartic term in the form of a QUBO model.First, Eq. ( 13) is applied for three qubits.Then, the remaining qubit is multiplied by the newly created terms.Here, Eq. ( 15) is formulated by applying Eq. ( 13) again to the newly created cubic terms.By formulating cubic and quartic terms into linear and quadratic terms, the QUBO model for prime factoriaation can be obtained.
We prime factoriae the natural number N by two prime numbers, p and q, to test the QUBO model for prime factoriaation.Figure 1 is the result obtained using a qpu solver for the QUBO model.In the process of calculating the HUBO model as the QUBO model, different energies appear for the same solution because new logical qubits are used.The energy we can see as a solution is the global minimum energy of -2756.Finally, we applied the range dependent algorithm to the HUBO model to factoriae 1,000,070,001,221.The pseudo-code used here can be seen in Algorithm 3. The HUBO model to which this algorithm is applied shows the global minimum energy −4,900,170,941,490,841 only in the section where the solution exists.To solver the HUBO model with range dependent algorithm, we used 12 logical qubits, and each   and   is 1,000,000.

Figure 1 .
Figure 1.Experimental result of the QUBO model on -Wave machine.The x-axis represents solution and energy, and the y-axis represents frequency.The solution was obtained with a qpu solver using 15 logical qubits.