HUBO式による素因数分解

2つの素数の積の素因数分解のための HUBO

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

$\sqrt{35}=5.91$ かつ $35/2=17.5$ であるため、$p$ と $q$ の探索範囲を以下のように制限できます:

\[\begin{aligned} 2 \leq &p \leq 5 \\ 6 \leq &q \leq 17 \end{aligned}\]

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

\[\begin{aligned} f(p,q) &= (pq-35)^2 \end{aligned}\]

整数変数 $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] ...

[omitted]

... +995284220088838892027904*p[19]*p[20]*q[19]*q[20]
sol = Sol(energy=0, {p[0]: 0, p[1]: 1, p[2]: 1, p[3]: 1, p[4]: 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 など）をインポートします。