HUBO式による因数分解¶

2つの素数の積を因数分解するためのHUBO¶

ここでは、2つの素数の積である整数の因数分解を考えます。
例えば、積 \(pq=35\) が与えられたとき、その2つの素因数 \(p=5\) と \(q=7\) を求めることが目標になります。

\(\sqrt{35}=5.91\) かつ \(35/2=17.5\) であることから、探索範囲を次のように制限できます。

\[\begin{split} 2 &≤ p ≤ 5,\\ 6 &≤ q ≤ 17. \end{split}\]

このような整数変数に対して、35 の因数分解問題は次のペナルティ式を用いて定式化できます。

\[ f(p, q)=(pq − 35)² \]

整数変数 \(p\) および \(q\) は2値変数の線形結合として実装されるため、その積 \(pq\) は2次式になります。したがって \(f(p, q)\) は4次式になります。\(f(p, q)\) は、\(p\) と \(q\) が 35 の正しい因数であるときにちょうど最小値 0 をとります。

因数分解のためのQUBO++プログラム

以下のQUBO++プログラムは、HUBO式 \(f(p, q)\) を構築し、Easy Solver を用いて最適化問題を解きます。

hubo-factorization-program1.cpp¶

#include <qbpp/qbpp.hpp>
#include <qbpp/easy_solver.hpp>

int main() {
  auto p = 2 <= qbpp::var_int("p") <= 5;
  auto q = 6 <= qbpp::var_int("q") <= 17;

  auto f = p * q == 35;
  f.simplify_as_binary();
  std::cout << "f = " << f << std::endl;
  std::cout << "f.body() = " << f.body() << std::endl;

  auto solver = qbpp::EasySolver(f);
  auto sol = solver.search({{"target_energy", 0}});

  std::cout << "sol = " << sol << std::endl;
  std::cout << "p = " << sol(p) << std::endl;
  std::cout << "q = " << sol(q) << std::endl;
  std::cout << "f(sol) = " << f(sol) << std::endl;
  std::cout << "f.body(sol) = " << f.body(sol) << std::endl;
}

このプログラムでは、式 p * q == 35 は自動的に qbpp::sqr(p * q - 35) に変換されます。そして、等式が満たされたときにエネルギー値 0 を達成します。

このプログラムの出力は次のとおりです。

f = 529 -240*p[0] -408*p[1] -88*q[0] -168*q[1] -304*q[2] -304*q[3] +144*p[0]*p[1] -5*p[0]*q[0] +40*p[0]*q[2] +40*p[0]*q[3] +16*p[1]*q[0] +56*p[1]*q[1] +208*p[1]*q[2] +208*p[1]*q[3] +16*q[0]*q[1] +32*q[0]*q[2] +32*q[0]*q[3] +64*q[1]*q[2] +64*q[1]*q[3] +128*q[2]*q[3] +52*p[0]*p[1]*q[0] +112*p[0]*p[1]*q[1] +256*p[0]*p[1]*q[2] +256*p[0]*p[1]*q[3] +20*p[0]*q[0]*q[1] +40*p[0]*q[0]*q[2] +40*p[0]*q[0]*q[3] +80*p[0]*q[1]*q[2] +80*p[0]*q[1]*q[3] +160*p[0]*q[2]*q[3] +48*p[1]*q[0]*q[1] +96*p[1]*q[0]*q[2] +96*p[1]*q[0]*q[3] +192*p[1]*q[1]*q[2] +192*p[1]*q[1]*q[3] +384*p[1]*q[2]*q[3] +16*p[0]*p[1]*q[0]*q[1] +32*p[0]*p[1]*q[0]*q[2] +32*p[0]*p[1]*q[0]*q[3] +64*p[0]*p[1]*q[1]*q[2] +64*p[0]*p[1]*q[1]*q[3] +128*p[0]*p[1]*q[2]*q[3]
sol = 0:{{p[0],1},{p[1],1},{q[0],1},{q[1],0},{q[2],0},{q[3],0}}
p = 5
q = 7

出力から、式 f には4次の項が含まれていることが確認できます。これにより、f がHUBO式であることが分かります。また、ソルバーは正しく素因数 \(p=5\) および \(q=7\) を見つけています。

大きな数の素因数分解における無制限の大きな係数

デフォルトでは、QUBO++ における式の係数およびエネルギー値のデータ型は、それぞれ int32_t および int64_t です。これらの型は、マクロ COEFF_TYPE および ENERGY_TYPE を定義することで変更できます。

さらに、QUBO++ は任意に大きな係数およびエネルギー値を持つ式をサポートしています。このオプションを有効にするには、両方のマクロを qbpp::cpp_int に設定します。以下の QUBO++ プログラムは、2つの大きな素数の積を因数分解します。

hubo-factorization-program2.cpp¶

#define INTEGER_TYPE_CPP_INT

#include <qbpp/qbpp.hpp>
#include <qbpp/easy_solver.hpp>

int main() {
  auto p = 2 <= qbpp::var_int("p") <= qbpp::integer("2000000");
  auto q = 2 <= qbpp::var_int("q") <= qbpp::integer("2000000");

  auto f = p * q == qbpp::integer("1000039") * qbpp::integer("1000079");
  f.simplify_as_binary();
  std::cout << "f = " << f << std::endl;

  auto solver = qbpp::EasySolver(f);
  auto sol = solver.search({{"target_energy", 0}});

  std::cout << "sol = " << sol << std::endl;
  std::cout << "p = " << sol(p) << std::endl;
  std::cout << "q = " << sol(q) << std::endl;
}

qbpp.hpp をインクルードする前に、マクロ COEFF_TYPE と ENERGY_TYPE を qbpp::cpp_int に設定しています。整数定数は、文字列リテラルを引数とする qbpp::cpp_int() を用いて指定します。

このプログラムは、次のような出力を生成します。

f = 1000236020078726181467929 -4000472012304*p[0] -8000944024600*p[1] -16001888049168*p[2] -32003776098208*p[3] -64007552195904*p[4] -128015104389760*p[5] -256030208771328*p[6] -512060417509888*p[7] -1024120834888704*p[8] -2048241669253120*p[9] -4096483336409088*p[10] -8192966664429568*p[11] -16385933295304704*p[12] -32771866456391680*p[13] -65543732375912448*p[14] -131087462604341248*p[15] -262174916618747904*p[16] -524349798877757440*p[17] -1048699460316561408*p[18] -2097398370877308928*p[19] -3806137462543214568*p[20] -4000472012304*q[0] -8000944024600*q[1] -16001888049168*q[2] -32003776098208*q[3] -64007552195904*q[4] -128015104389760*q[5] -256030208771328*q[6] -512060417509888*q[7] -1024120834888704*q[8] -2048241669253120*q[9] -4096483336409088*q[10] -8192966664429568*q[11] -16385933295304704*q[12] -32771866456391680*q[13] -65543732375912448*q[14] -131087462604341248*q[15] -262174916618747904*q[16] -524349798877757440*q[17] -1048699460316561408*q[18] -2097398370877308928*q[19] -3806137462543214568*q[20] 

[omitted]

+17139313073086877663232*p[19]*p[20]*q[14]*q[19] +31102631877776215375872*p[19]*p[20]*q[14]*q[20] +4284828268271719415808*p[19]*p[20]*q[15]*q[16] +8569656536543438831616*p[19]*p[20]*q[15]*q[17] +17139313073086877663232*p[19]*p[20]*q[15]*q[18] +34278626146173755326464*p[19]*p[20]*q[15]*q[19] +62205263755552430751744*p[19]*p[20]*q[15]*q[20] +17139313073086877663232*p[19]*p[20]*q[16]*q[17] +34278626146173755326464*p[19]*p[20]*q[16]*q[18] +68557252292347510652928*p[19]*p[20]*q[16]*q[19] +124410527511104861503488*p[19]*p[20]*q[16]*q[20] +68557252292347510652928*p[19]*p[20]*q[17]*q[18] +137114504584695021305856*p[19]*p[20]*q[17]*q[19] +248821055022209723006976*p[19]*p[20]*q[17]*q[20] +274229009169390042611712*p[19]*p[20]*q[18]*q[19] +497642110044419446013952*p[19]*p[20]*q[18]*q[20] +995284220088838892027904*p[19]*p[20]*q[19]*q[20]
sol = 0:{{p[0],0},{p[1],1},{p[2],1},{p[3],1},{p[4],0},{p[5],0},{p[6],0},{p[7],0},{p[8],0},{p[9],1},{p[10],1},{p[11],1},{p[12],1},{p[13],1},{p[14],0},{p[15],1},{p[16],0},{p[17],0},{p[18],0},{p[19],0},{p[20],1},{q[0],0},{q[1],1},{q[2],1},{q[3],0},{q[4],0},{q[5],1},{q[6],1},{q[7],1},{q[8],1},{q[9],0},{q[10],1},{q[11],1},{q[12],1},{q[13],1},{q[14],0},{q[15],1},{q[16],0},{q[17],0},{q[18],0},{q[19],0},{q[20],1}}
p = 1000079
q = 1000039

式 f には非常に大きな係数が含まれていることが分かります。また、大きな合成数の因数分解が正しく求められていることも確認できます。

TIP: COEFF_TYPE および ENERGY_TYPE には、int8_t, int16_t, int32_t, int64_t, qbpp::int128_t, qbpp::int256_t, qbpp::int512_t, qbpp::int1024_t, qbpp::cpp_int を設定できます。

因数分解のためのPyQBPPプログラム

以下のPyQBPPプログラムは、HUBO式 \(f(p, q)\) を構築し、Easy Solver を用いて最適化問題を解きます。

hubo-factorization-program1.py¶

import pyqbpp as qbpp

p = qbpp.var("p", between=(2, 5))
q = qbpp.var("q", between=(6, 17))

f = qbpp.constrain(p * q, equal=35)
f.simplify_as_binary()
print("f =", f)
print("f.body =", f.body)

solver = qbpp.EasySolver(f)
sol = solver.search(target_energy=0)

print("sol =", sol)
print("p =", sol(p))
print("q =", sol(q))
print("f(sol) =", f(sol))
print("f.body(sol) =", f.body(sol))

このプログラムでは、式 p * q == 35 は自動的に sqr(p * q - 35) に変換されます。そして、等式が満たされたときにエネルギー値 0 を達成します。

このプログラムの出力は次のとおりです。

f = 529 -240*p[0] -408*p[1] -88*q[0] -168*q[1] -304*q[2] -304*q[3] +144*p[0]*p[1] -5*p[0]*q[0] +40*p[0]*q[2] +40*p[0]*q[3] +16*p[1]*q[0] +56*p[1]*q[1] +208*p[1]*q[2] +208*p[1]*q[3] +16*q[0]*q[1] +32*q[0]*q[2] +32*q[0]*q[3] +64*q[1]*q[2] +64*q[1]*q[3] +128*q[2]*q[3] +52*p[0]*p[1]*q[0] +112*p[0]*p[1]*q[1] +256*p[0]*p[1]*q[2] +256*p[0]*p[1]*q[3] +20*p[0]*q[0]*q[1] +40*p[0]*q[0]*q[2] +40*p[0]*q[0]*q[3] +80*p[0]*q[1]*q[2] +80*p[0]*q[1]*q[3] +160*p[0]*q[2]*q[3] +48*p[1]*q[0]*q[1] +96*p[1]*q[0]*q[2] +96*p[1]*q[0]*q[3] +192*p[1]*q[1]*q[2] +192*p[1]*q[1]*q[3] +384*p[1]*q[2]*q[3] +16*p[0]*p[1]*q[0]*q[1] +32*p[0]*p[1]*q[0]*q[2] +32*p[0]*p[1]*q[0]*q[3] +64*p[0]*p[1]*q[1]*q[2] +64*p[0]*p[1]*q[1]*q[3] +128*p[0]*p[1]*q[2]*q[3]
sol = 0:{{p[0],1},{p[1],1},{q[0],1},{q[1],0},{q[2],0},{q[3],0}}
p = 5
q = 7

大きな数の素因数分解のための無制限の大きな係数

PyQBPP のデフォルトモジュール pyqbpp（pyqbpp.c32e64 の別名）は、32 ビット係数と 64 ビットエネルギー値を使用し、最も高速に動作します。大きな合成数の素因数分解では、ペナルティ式の係数が 32 ビット範囲を超える場合があります。そのような場合は pyqbpp.cppint をインポートすることで、任意精度整数演算（coeff_t と energy_t の両方を cpp_int に設定）へ切り替えることができます。

PyQBPP では、係数型・エネルギー型の選択をインポートするサブモジュール名（pyqbpp.c32e64、pyqbpp.cppint など）によって指定します。

以降の qbpp.var、qbpp.constrain、qbpp.EasySolver などの呼び出しは、選択した係数型・エネルギー型を一貫して使用します。

以下のプログラムは、2つの大きな素数の積を因数分解する例です。

hubo-factorization-program2.py¶

import pyqbpp.cppint as qbpp

p = qbpp.var("p", between=(2, 2000000))
q = qbpp.var("q", between=(2, 2000000))

f = qbpp.constrain(p * q, equal=1000039 * 1000079)
f.simplify_as_binary()
print("f =", f)

solver = qbpp.EasySolver(f)
sol = solver.search(target_energy=0)

print("sol =", sol)
print("p =", sol(p))
print("q =", sol(q))

Python は任意精度整数をネイティブにサポートしているため、目標値 1000039 * 1000079 は Python の int としてそのまま記述できます。pyqbpp.cppint をインポートすれば、PyQBPP は係数とエネルギー値をすべて Python の int として保持・操作するため、上限値 2000000 や目標の積を Python 整数として自然に記述できます。

このプログラムは以下の結果を出力します:

f = 1000236020078726181467929 -4000472012304*p[0] -8000944024600*p[1] -16001888049168*p[2] -32003776098208*p[3] -64007552195904*p[4] -128015104389760*p[5] -256030208771328*p[6] -512060417509888*p[7] -1024120834888704*p[8] -2048241669253120*p[9] -4096483336409088*p[10] -8192966664429568*p[11] -16385933295304704*p[12] -32771866456391680*p[13] -65543732375912448*p[14] -131087462604341248*p[15] -262174916618747904*p[16] -524349798877757440*p[17] -1048699460316561408*p[18] -2097398370877308928*p[19] -3806137462543214568*p[20] -4000472012304*q[0] -8000944024600*q[1] -16001888049168*q[2] -32003776098208*q[3] -64007552195904*q[4] -128015104389760*q[5] -256030208771328*q[6] -512060417509888*q[7] -1024120834888704*q[8] -2048241669253120*q[9] -4096483336409088*q[10] -8192966664429568*q[11] -16385933295304704*q[12] -32771866456391680*q[13] -65543732375912448*q[14] -131087462604341248*q[15] -262174916618747904*q[16] -524349798877757440*q[17] -1048699460316561408*q[18] -2097398370877308928*q[19] -3806137462543214568*q[20] 

[omitted]

+17139313073086877663232*p[19]*p[20]*q[14]*q[19] +31102631877776215375872*p[19]*p[20]*q[14]*q[20] +4284828268271719415808*p[19]*p[20]*q[15]*q[16] +8569656536543438831616*p[19]*p[20]*q[15]*q[17] +17139313073086877663232*p[19]*p[20]*q[15]*q[18] +34278626146173755326464*p[19]*p[20]*q[15]*q[19] +62205263755552430751744*p[19]*p[20]*q[15]*q[20] +17139313073086877663232*p[19]*p[20]*q[16]*q[17] +34278626146173755326464*p[19]*p[20]*q[16]*q[18] +68557252292347510652928*p[19]*p[20]*q[16]*q[19] +124410527511104861503488*p[19]*p[20]*q[16]*q[20] +68557252292347510652928*p[19]*p[20]*q[17]*q[18] +137114504584695021305856*p[19]*p[20]*q[17]*q[19] +248821055022209723006976*p[19]*p[20]*q[17]*q[20] +274229009169390042611712*p[19]*p[20]*q[18]*q[19] +497642110044419446013952*p[19]*p[20]*q[18]*q[20] +995284220088838892027904*p[19]*p[20]*q[19]*q[20]
sol = 0:{{p[0],0},{p[1],1},{p[2],1},{p[3],1},{p[4],0},{p[5],0},{p[6],0},{p[7],0},{p[8],0},{p[9],1},{p[10],1},{p[11],1},{p[12],1},{p[13],1},{p[14],0},{p[15],1},{p[16],0},{p[17],0},{p[18],0},{p[19],0},{p[20],1},{q[0],0},{q[1],1},{q[2],1},{q[3],0},{q[4],0},{q[5],1},{q[6],1},{q[7],1},{q[8],1},{q[9],0},{q[10],1},{q[11],1},{q[12],1},{q[13],1},{q[14],0},{q[15],1},{q[16],0},{q[17],0},{q[18],0},{q[19],0},{q[20],1}}
p = 1000079
q = 1000039

式 f が非常に大きな係数（64ビットの範囲を大きく超える）を含み、大きな合成数の因数分解が \(1000079 \times 1000039\) として正しく得られていることがわかります。整数変数 \(p\) と \(q\) はそれぞれ範囲 \([2, 2000000]\) をカバーするために21個のバイナリ変数（p[0]-p[20]、q[0]-q[20]）に展開されています。

TIP: 任意精度整数を使用するには import pyqbpp.cppint as qbpp とします。他の整数型バリアントを使用する場合は、対応するサブモジュール（pyqbpp.c64e64、pyqbpp.c128e128 など）をインポートします。