HUBO式による因数分解¶
2つの素数の積を因数分解するためのHUBO
2つの素数の積である整数の因数分解を考えます. 例えば,積 \(pq = 35\) が与えられたとき,2つの素因数 \(p=5\) と \(q=7\) を求めることが目標です.
\(\sqrt{35}=5.91\) かつ \(35/2=17.5\) であるため,\(p\) と \(q\) の探索範囲を以下のように制限できます:
このような整数変数に対して,\(35\) の因数分解問題は以下のペナルティ式を用いて定式化できます:
整数変数 \(p\) と \(q\) はバイナリ変数の線形式として実装されるため,その積 \(pq\) は2次式となり,したがって \(f(p,q)\) は4次式になります. 明らかに,\(f(p,q)\) は \(p\) と \(q\) が35の正しい因数であるとき,かつそのときに限り最小値0を達成します.
因数分解のためのHi-QUBOプログラム
以下のHi-QUBOプログラムはHUBO式 \(f(p,q)\) を構築し,Easy Solverを使って最適化問題を解きます:
#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 は制約式の qbpp::Expr で,展開後のペナルティ qbpp::sqr(p * q - 35) と元の式 p * q の両方を保持します.f 自身は式 qbpp::sqr(p * q - 35) を表し,f.body() は元の式 p * q を返します.f.simplify_as_binary() は f 自身（ペナルティ）と f.body()（元の式）の両方を同時に簡約します.
このプログラムの出力は以下の通りです:
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]
f.body() = 12 +6*p[0] +12*p[1] +2*q[0] +4*q[1] +8*q[2] +8*q[3] +p[0]*q[0] +2*p[0]*q[1] +4*p[0]*q[2] +4*p[0]*q[3] +2*p[1]*q[0] +4*p[1]*q[1] +8*p[1]*q[2] +8*p[1]*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(sol) = 0
f.body(sol) = 35
出力から,式 f が4次の項を含んでおり,HUBO式であることが確認できます.
ソルバーは素因数 \(p=5\) と \(q=7\) を正しく求めています.
大きな数の素因数分解のための無制限の大きな係数
Hi-QUBOでは,式の係数とエネルギー値のデータ型はデフォルトでそれぞれ int32_t と int64_t です.
これらの型は INTEGER_TYPE_* マクロ（INTEGER_TYPE_C32E32,INTEGER_TYPE_C32E64（デフォルト）,INTEGER_TYPE_C64E64,INTEGER_TYPE_C64E128,INTEGER_TYPE_C128E128,INTEGER_TYPE_CPP_INT）のいずれかを定義することで切り替えられます — 詳細は VAREXPR 参照.
さらに,Hi-QUBOは任意の大きさの係数とエネルギー値を持つ式を INTEGER_TYPE_CPP_INT でサポートします.このマクロを定義すると coeff_t と energy_t の両方が qbpp::cpp_int になります.
以下のHi-QUBOプログラムは,2つの大きな素数の積を因数分解します:
#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/qbpp.hpp をインクルードする前に INTEGER_TYPE_CPP_INT を定義し,coeff_t と energy_t を任意精度整数 cpp_int に設定します.
定数整数は qbpp::integer() に10進文字列を渡して指定します（int * int の中間オーバーフローを避けるため,積は qbpp::integer("1000039") * qbpp::integer("1000079") と書いています）.
このプログラムは以下の結果を出力します:
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 任意精度整数を使用するには,ヘッダのインクルード前に
INTEGER_TYPE_CPP_INTを定義します（あるいは-DINTEGER_TYPE_CPP_INTをコンパイラフラグで渡します）.他の型はINTEGER_TYPE_C32E32,INTEGER_TYPE_C32E64（デフォルト）,INTEGER_TYPE_C64E64,INTEGER_TYPE_C64E128,INTEGER_TYPE_C128E128で選択します — 全対応表は VAREXPR 参照.
2つの素数の積の素因数分解のための HUBO
2つの素数の積である整数の素因数分解を考えます. 例えば,積 \(pq = 35\) が与えられたとき,2つの素因数 \(p=5\) と \(q=7\) を求めることが目標です.
\(\sqrt{35}=5.91\) かつ \(35/2=17.5\) であるため,\(p\) と \(q\) の探索範囲を以下のように制限できます:
このような整数変数に対して,\(35\) の素因数分解問題はペナルティ式を用いて以下のように定式化できます:
整数変数 \(p\) と \(q\) はバイナリ変数の線形式として実装されるため,その積 \(pq\) は2次式となり,したがって \(f(p,q)\) は4次式になります. 明らかに,\(f(p,q)\) は \(p\) と \(q\) が 35 の正しい因数であるときに限り最小値 0 を達成します.
素因数分解の PyQBPP プログラム
以下のプログラムは HUBO 式 \(f(p,q)\) を構築し,Easy Solver を用いて最適化問題を解きます:
import pyqbpp as qbpp
p = qbpp.var("p", between=(2, 5))
q = qbpp.var("q", between=(6, 17))
f = (p * q == 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))
このプログラムでは,qbpp.var("p", between=(2, 5)) で値域 \(\{2, 3, 4, 5\}\) をとる整数変数 \(p\) を生成します（内部的にはバイナリ変数 p[0], p[1] の線形結合として表現されます）.同様に,値域 \(\{6, 7, \dots, 17\}\) をとる \(q\) はバイナリ変数 q[0], q[1], q[2], q[3] に展開されます.式 (p * q == 35) は自動的に sqr(p * q - 35) に変換され,等式が満たされたときにエネルギー値 0 を達成します.f は制約式で,展開後のペナルティ sqr(p * q - 35) と元の式 p * q の両方を保持します.f 自身は式 sqr(p * q - 35) を表し,f.body は元の式 p * q を返します.f.simplify_as_binary() は f 自身（ペナルティ）と f.body（元の式）の両方を同時に簡約します.整数変数は内部のバイナリ変数に対して線形なので,積 \(pq\) は2次式となり,その2乗である \(f(p,q)\) は4次式となります.
このプログラムの出力は以下の通りです:
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]
f.body = 12 +6*p[0] +12*p[1] +2*q[0] +4*q[1] +8*q[2] +8*q[3] +p[0]*q[0] +2*p[0]*q[1] +4*p[0]*q[2] +4*p[0]*q[3] +2*p[1]*q[0] +4*p[1]*q[1] +8*p[1]*q[2] +8*p[1]*q[3]
sol = Sol(energy=0, {p[0]: 1, p[1]: 1, q[0]: 1, q[1]: 0, q[2]: 0, q[3]: 0})
p = 5
q = 7
f(sol) = 0
f.body(sol) = 35
出力から,式 f が4次の項を含むことが確認でき,これが HUBO 式であることがわかります.
ソルバーは素因数 \(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.EasySolver などの呼び出しは,選択した係数・エネルギー型を一貫して使用します.
以下のプログラムは2つの大きな素数の積を因数分解します:
import pyqbpp.cppint as qbpp
p = qbpp.var("p", between=(2, 2000000))
q = qbpp.var("q", between=(2, 2000000))
f = (p * q == 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など）をインポートします.