最大公約数 (GCD)¶

\(P\) と \(Q\) を2つの正の整数とします. 最大公約数 (GCD) の計算はHUBO問題として定式化できます.

\(p\),\(q\),\(r\) を以下の制約を満たす正の整数とします:

\[\begin{split} \begin{aligned} p\cdot r &= P \\ q\cdot r &=Q \end{aligned} \end{split}\]

明らかに,\(r\) は \(P\) と \(Q\) の公約数です. したがって,これらの制約を満たす \(r\) の最大値が \(P\) と \(Q\) の最大公約数となります. そのような \(r\) を求めるために,HUBO定式化において \(-r\) を目的関数として使用します.

Hi-QUBO プログラム

上記のアイデアに基づき,以下のHi-QUBOプログラムは2つの整数 P = 858 と Q = 693 の最大公約数を計算します:

greatest-common-divisor-program1.cpp¶

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

int main() {
  const int P = 858;
  const int Q = 693;
  auto p = 1 <= qbpp::var_int("p") <= 1000;
  auto q = 1 <= qbpp::var_int("q") <= 1000;
  auto r = 1 <= qbpp::var_int("r") <= 1000;

  auto constraint = (p * r == P) + (q * r == Q);
  auto f = -r + constraint * 1000;

  f.simplify_as_binary();

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

  std::cout << "GCD = " << sol(r) << std::endl;
  std::cout << sol(p) << " * " << sol(r) << " = " << P << std::endl;
  std::cout << sol(q) << " * " << sol(r) << " = " << Q << std::endl;
}

このプログラムでは,p,q,r は \([1,1000]\) の範囲の整数変数として定義されています. 二乗ペナルティ項を展開すると係数が \(10^{13}\) 程度に達し,デフォルトの32ビット coeff_t を超えるため,ヘッダのインクルード前に INTEGER_TYPE_C64E64 を定義して64ビットの係数・エネルギーを使用します（整数型マクロの一覧は FACTORIZATION 参照）. 式 constraint は,両方の制約が満たされたときにゼロに評価されるように構成されています.

目的関数 -r はペナルティ係数 1000 を掛けた制約項と組み合わされ,結果の式が f に格納されます.

EasySolverは f を最小化する解を探索します. 得られた p,q,r の値は以下のように出力されます:

この出力から,858と693の最大公約数が正しく33として得られたことが確認できます.

\(P\) と \(Q\) を2つの正の整数とします. 最大公約数 (GCD) の計算は HUBO 問題として定式化できます.

\(p\),\(q\),\(r\) を以下の制約を満たす正の整数とします:

\[\begin{split} \begin{aligned} p\cdot r &= P \\ q\cdot r &=Q \end{aligned} \end{split}\]

明らかに,\(r\) は \(P\) と \(Q\) の公約数です. したがって,これらの制約を満たす \(r\) の最大値が \(P\) と \(Q\) の GCD です. そのような \(r\) を求めるために,HUBO 定式化において \(-r\) を目的関数として使用します.

PyQBPP プログラム

上記の考え方に基づき,以下の PyQBPP プログラムは2つの整数 P = 858 と Q = 693 の GCD を計算します:

greatest-common-divisor-program1.py¶

import pyqbpp.c64e64 as qbpp

P = 858
Q = 693
p = qbpp.var("p", between=(1, 1000))
q = qbpp.var("q", between=(1, 1000))
r = qbpp.var("r", between=(1, 1000))

constraint = (p * r == P) + (q * r == Q)
f = -r + constraint * 1000

f.simplify_as_binary()

solver = qbpp.EasySolver(f)
sol = solver.search(time_limit=1.0)

print(f"GCD = {sol(r)}")
print(f"{sol(p)} * {sol(r)} = {P}")
print(f"{sol(q)} * {sol(r)} = {Q}")

このプログラムでは,p,q,r は範囲 \([1,1000]\) の整数変数として定義されています. 二乗ペナルティ項を展開すると係数が \(10^{13}\) 程度に達し,デフォルトの pyqbpp モジュールの32ビット係数を超えるため,pyqbpp.c64e64 をインポートして64ビットの係数・エネルギーを使用します（利用可能な型バリアントは FACTORIZATION 参照）. 式 constraint は,両方の制約が満たされたときにゼロと評価されるように構築されています.

目的関数 -r はペナルティ係数 1000 を掛けた制約項と組み合わされ,結果の式は f に格納されます.

EasySolver は f を最小化する解を探索します. 得られた p,q,r の値は以下のように出力されます:

GCD = 33
26 * 33 = 858
21 * 33 = 693

この出力から,858 と 693 の GCD が 33 として正しく求められたことが確認できます.