再帰による素因数分解

HUBO 式による因数分解では、2 つの素数の積を 2 つの因数に分けました。このページでは、これを繰り返して、整数を素因数に分解します。 1 回の分解では、ネイティブ整数変数（integer=）と qbpp.cons() を使い、$xy = n$ となる $x$、$y$ を Easy Solver で探します。分解できたら、$x$ と $y$ のそれぞれを同じように分解します。

1 回の分解

合成数 $n$ は、$2 \le x \le y$ を満たす 2 つの整数の積 $n = xy$ で表せます。 $x \le y$ なら $x^2 \le xy = n \le y^2$ なので、$x \le \sqrt{n} \le y$ です。また $x \ge 2$ なので $y \le n/2$ です。そこで、次の範囲で

\[2 \le x \le \lfloor\sqrt{n}\rfloor, \qquad \lfloor\sqrt{n}\rfloor \le y \le \lfloor n/2 \rfloor\]

制約

\[xy = n\]

を満たす $x$、$y$ を探します。 $x \le y$ は、$x$ と $y$ を入れ替えただけの解を除くための条件です。

素数判定

$n$ が素数のときは、この制約を満たす $x$、$y$ はありません。ところが、ソルバーは「解がない」ことを示せません。解が見つからないのが、$n$ が素数だからなのか、探しきれなかったからなのかを区別できず、制限時間いっぱいまで探し続けます。そこで、素数かどうかは PyQBPP を使わずに sympy の isprime で判定し、合成数だけをソルバーで分解します。 sympy は pip install sympy でインストールできます。

プログラム

import math
from sympy import isprime
import pyqbpp.c64e128 as qbpp


def split(n):
    r = math.isqrt(n)
    x = qbpp.var("x", integer=(2, r))
    y = qbpp.var("y", integer=(r, n // 2))
    f = qbpp.cons(x * y == n)
    f.simplify_as_binary()
    sol = qbpp.EasySolver(f).search(time_limit=10, target_energy=0)
    return sol(x), sol(y)


def factorize(n):
    if isprime(n):
        return [n]
    if n == 4:  # x and y would both be the constant 2
        return [2, 2]
    x, y = split(n)
    if x * y != n:
        raise RuntimeError(f"no factor of {n} found within the time limit")
    print(f"{n} = {x} * {y}")
    return factorize(x) + factorize(y)


for n in [4294967297, 999999999999]:
    factors = sorted(factorize(n))
    print(n, "=", " * ".join(str(p) for p in factors))

split(n) は、上の範囲のネイティブ整数変数 x、y を作り、qbpp.cons(x * y == n) を Easy Solver で解きます。 target_energy を 0 にしているので、制約を満たす解が見つかった時点で探索を止めます。
factorize(n) は、n が素数ならそのまま返し、そうでなければ split(n) で 2 つに分けて、それぞれを再帰的に分解します。 10 秒以内に分解が見つからなければ、例外を投げて止まります。
n が 4 のときは、x と y の範囲がどちらも 2 だけになり、式に変数が残らないので、ソルバーに渡す前に扱います。
pyqbpp.c64e128 は、係数を 64 ビット、エネルギーを 128 ビットの整数にするモジュールです（理由は後で説明します）。

このプログラムの出力は、たとえば次のとおりです:

4294967297 = 641 * 6700417
4294967297 = 641 * 6700417
999999999999 = 3 * 333333333333
333333333333 = 3 * 111111111111
111111111111 = 3 * 37037037037
37037037037 = 7 * 5291005291
5291005291 = 11 * 481000481
481000481 = 13 * 37000037
37000037 = 37 * 1000001
1000001 = 101 * 9901
999999999999 = 3 * 3 * 3 * 7 * 11 * 13 * 37 * 101 * 9901

$4294967297 = 2^{32}+1$ は、フェルマーが素数だと予想した数で、オイラーが 1732 年に $641$ で割り切れることを示しました。 Easy Solver は乱数を使うので、どの約数の組が見つかるかは実行するたびに変わり、途中の分け方は変わることがあります。最後の素因数の並びは変わりません。

係数とエネルギーの型

$n$ は 64 ビットの整数に収まりますが、探索の途中ではもっと大きな値が現れます。 $x$ と $y$ がそれぞれの範囲の上限にあるとき、$xy$ は約 $n^{1.5}/2$ になり、制約の値（違反量の 2 乗）$(xy-n)^2$ は約 $n^3/4$ になります。 $n < 10^{12}$ ならこの値は約 $2.5\times10^{35}$ 以下で、128 ビットの整数（約 $1.7\times10^{38}$ まで）に収まります。そこで、係数を 64 ビット、エネルギーを 128 ビットとする pyqbpp.c64e128 を使います。既定の pyqbpp（係数 32 ビット、エネルギー 64 ビット）では足りません。あふれても警告は出ないので、型は探索の途中に現れる最大の値で選びます。 $n$ が約 $9\times10^{12}$ を超えると 128 ビットにも収まらなくなるので、桁数に上限のない pyqbpp.cppint が必要です（データ型を参照）。

時間がかかる数

この方法で速く分解できるかどうかは、小さいほうの因数 $x$ の位置で決まります。

因数が小さい（2 に近い）とき、または 2 つの因数が近い（$x$ が $\sqrt{n}$ に近い）ときは、範囲の端に答えがあるので、すぐに見つかります。 12 桁のランダムな合成数 30 個は、どれも 0.15 秒以内に最初の分解が見つかりました。
小さいほうの因数が範囲の中ほどにあると、時間がかかります。 4〜6 桁の素数ともう 1 つの素数の積で 12 桁になる数 30 個では、20 秒の制限で分解できたのは 20 個（中央値 6.6 秒）で、10 個は見つかりませんでした。たとえば $508069567969 = 7793 \times 65195633$ は、20 秒では分解できませんでした（20 コアの CPU で計測）。

$(xy-n)^2$ は、$x$ が因数に近づいても小さくなるとは限らず、どこに因数があるかの手がかりになりません。そのため、ソルバーはほぼ手探りで $x$ を探すことになります。これは QUBO で素因数分解するときの本質的な限界で、このページの方法は実用的な素因数分解の方法ではありません（試し割りや楕円曲線法などのほうがはるかに速く分解できます）。