QAOA を用いた素因数分解法について (Part 2: 平方差法による素因数分解)

2022 年 12 月に Yan 等の中国の研究者グループが発表した "[Factoring integers with sublinear resources on a superconducting quantum processor](https://arxiv.org/abs/2212.12372)" という論文 (プレプリント) を 6 回に分けて解説しています. この論文は QAOA という量子計算を使った素因数分解アルゴリズム SQIF (Sublinear Quantum Integer Factorization) を提案しています. 本記事では SQIF の背景となる, 平方差法による素因数分解法の原理を紹介します. この原理は[二次篩法](https://en.wikipedia.org/wiki/Quadratic_sieve)や[数体篩法](https://ja.wikipedia.org/wiki/%E4%B8%80%E8%88%AC%E6%95%B0%E4%BD%93%E7%AF%A9%E6%B3%95)でも使用されているので, おなじみの方も多いと思います.

[Part 1: 論文の概要](https://qiita.com/izutetsuya/items/f40854aeb63d623410b8)
Part 2: 平方差法による素因数分解  (本記事)
[Part 3: Schnorr の素因数分解法](https://qiita.com/izutetsuya/items/a82aa95f11cf0de6ee6b)
[Part 4: SQIF の詳細](https://qiita.com/izutetsuya/items/2ceeb8c8dd581799b4da) 
[Part 5: SQIF の数値実験の分析](https://qiita.com/izutetsuya/items/15a537c6573153c200c3)
[Part 6: まとめ](https://qiita.com/izutetsuya/items/3b54b70337a529b1d8d8)

https://arxiv.org/abs/2212.12372

# 更新履歴
- 2023年5月29日：Typo を修正するとともに、わかりにくい部分を変更しました。
- 2023年5月9日：一連の記事へのリンクを追加しました。また各記事の用語やトーンを揃えるため、全体的に加筆変更しました。

# 平方差法の原理

## アイディア

平方差法は, 合成数 $N$ を素因数分解する際に, 平方を $N$ で割った余りが等しくなるような異なる自然数 $x, y$ (つまり $x^2 \equiv y^2 \pmod{N}$ を満たす自然数 $x, y$) を見つけることを目標とします. このような自然数 $x,y$ が見つかったならば $x^2-y^2$ は $N$ の倍数となるので, $x^2-y^2$ の多項式としての因数 $x+y$ が自然数としての $N$ の因数に一致することに期待して両者の最大公約数 $\gcd(x+y,N)$ を計算し, $N$ の因数を得ようとします. 平方の差を利用することから平方差法 (Difference of Square Method) と呼ばれます. 平方差法では, $\gcd$ の計算結果が $1$ または $N$ 自身になってしまい素因数分解に失敗することもあるのですが, その場合には $x,y$ を取り替えていくことで素因数分解が可能になります.

## 平方差法の数値例

例として $N=140086253=10559 \times 13267$ の場合で考えてみましょう. $x=23778411$, $y=83479911$ のとき $x^2$, $y^2$ を $N$ で割った余りはともに $57396393$ となる (つまり $x^2 \bmod{N}=y^2 \bmod{N}=57396393$ となる)ので,
```math
  \begin{eqnarray}
    \gcd(x+y,N) &=& \gcd(23778411+83479911,140086253) \\ 
                &=& \gcd(107258322,140086253) \\
                &=& 10559
  \end{eqnarray}
```
となり, 素因数 $10559$ を見つけることに成功します. しかし $y'=116307842$ のときは確かに $y'^2 \bmod{N}=57396393$ となりますが,
```math
  \begin{eqnarray}
    \gcd(x+y',N) &=& \gcd(23778411+116307842,140086253) \\
                 &=& \gcd(140086253,140086253) \\
                 &=& 140086253
  \end{eqnarray}
```
となり, 素因数分解に失敗します. この場合は $x+y'=N$ になっていることが失敗の原因です.

## 自然数の見つけ方

ではどうやって平方を $N$ で割った余りが等しい自然数 $x,y$ を求めれば良いでしょうか. 

量子素因数分解アルゴリズムとして有名な [Shor アルゴリズム](https://en.wikipedia.org/wiki/Shor%27s_algorithm) は入力された自然数 $a,N$ に対し $a^r \equiv 1 \pmod{N}$ となる自然数 $r$ を求めるアルゴリズムです (このような $r$ を位数と呼びます) ので, $r$ が偶数のときに $x=a^{\frac{r}{2}}$, $y=1$ とおくことで平方差法を適用できます. しかし一般にはこのような $x,y$ を直接求めるのは困難です. そこで「べき積」を利用して上記のような $x,y$ を構成するアプローチを用いるのが一般的です.

# 関係式

素数のべき乗を使った積をべき積と呼びます. 例えば $2^3 \times 3^2~(=72)$ はべき積ですし, $3 \times 5~(=15)$ のようにべき指数が 1 の場合もべき積と呼びます. ここで, 平方差を求める代わりに, 以下のように平方が $\bmod N$ で素数のべき積となる自然数 $x$ を求める問題を考えます. 

```math
 x^2 \equiv q_1^{f_1} \times q_2^{f_2} \times \dots \times q_t^{f_t} \bmod{N}
```

上記の式を「関係式 (relation)」と呼びます. 平方差法の式 (合同式) と比べると, 左辺は同じですが右辺はべき積になっており, 平方差を求めるよりも優しくなっています. 

## 関係式の数値例

関係式を何個か得られた場合, それらから目標の合同式を構成することができます. 例えば $N=140086253=10559 \times 1326$ に対し, 以下の3つの関係式が得られたとします.
```math
 28367^2 \equiv 2^{6} \times 3^{4} \times 7^{1} \times 13^{2} \times 17^{1} \bmod{N} \\
 28549^2 \equiv 2^{3} \times 3^{3} \times 7^{4} \times 13^{1} \times 17^{1} \bmod{N}  \\
 29448^2 \equiv 2^{1} \times 3^{1} \times 7^{1} \times 13^{3} \times 17^{2} \bmod{N}
 ```
このとき左辺同士, 右辺同士の積を考えることで
```math
  \begin{eqnarray}
 (28367 \times 28549 \times 29448)^2 &\equiv& 
 2^{10} \times 3^{8} \times 7^{6} \times 13^{6} \times 17^{4} \bmod{N} \\
 &\equiv& 
 (2^{5} \times 3^{4} \times 7^{3} \times 13^{3} \times 17^{2})^2 \bmod{N}
  \end{eqnarray}
```
という式を得ます. よって
```math
  x=28367 \times 28549 \times 29448 \bmod N=23778411, \\
  y=2^5 \times 3^4 \times 7^3 \times 13^3 \times 17^2 \bmod{N}=83479911, 
```
とおくことで, 合同式 $x^2 \equiv y^2 \pmod{N}$ を得ることができました.

## 関係式から合同式を求める計算

さきほどの数値例では右辺の積がたまたま平方数になりましたが, 一般には平方数にならないでしょうから, その場合の計算方法が必要になります. そこで以下のように関係式が $t+1$ 個得られた場合を考えます:

```math
 x_0^2 \equiv q_1^{f_{0,1}} \times q_2^{f_{0,2}} \times \dots \times q_t^{f_{0,t}} \bmod{N} \\
 x_1^2 \equiv q_1^{f_{1,1}} \times q_2^{f_{1,2}} \times \dots \times q_t^{f_{1,t}} \bmod{N} \\
 x_2^2 \equiv q_1^{f_{2,1}} \times q_2^{f_{2,2}} \times \dots \times q_t^{f_{2,t}} \bmod{N} \\
 \vdots \\
 x_t^2 \equiv q_1^{f_{t,1}} \times q_2^{f_{t,2}} \times \dots \times q_t^{f_{t,t}} \bmod{N} \\
```
われわれの目標は, 適切な式を掛け合わせることで右辺が平方数になるようにすることです. そこで右辺のべき指数から生成できるベクトルを $f_i=(f_{i,1},f_{i,2},\dots,f_{i,t}) \bmod{2}$ と定めると, 右辺の積が平方数になるような関係式を求めるという問題は $u_0 f_0 + u_1 f_1 + \dots + u_t f_t = (0,0,\dots,0) \bmod{2}$ となる係数 $u_0, u_1, \dots, u_t ~(u_i=0 ~\textrm{or}~ 1)$ を求める問題に言い換えることができます. このとき 
```math
  \begin{eqnarray}
    x &=& x_0^{u_0} \times x_1^{u_1} \times x_2^{u_2} \times \dots \times x_t^{u_t} \bmod{N} \\
    y &=& q_1^{v_1} \times q_1^{v_2} \times \dots \times q_t^{v_t} \bmod{N}
  \end{eqnarray}
```
とおくことで, 欲しかった合同式 $x^2 \equiv y^2 \pmod{N}$ を得ることができます. ただし $(v_1,v_2,\dots,v_t)$ $=\frac{1}{2}(u_0 f_0 + u_1 f_1 + \dots + u_t f_t)$ です.

## 掃き出し法

係数 $u_0, u_1, \dots, u_t ~(u_i=0 ~\textrm{or}~ 1)$ を求めるには, ベクトル $f_i$ から生成される以下の $(t+1)\times t$ 行列に対して掃き出し法を適用します:
```math
  \begin{pmatrix}
    f_{0,1} & f_{0,2} & \cdots & f_{0,t} \\
    f_{1,1} & f_{1,2} & \cdots & f_{1,t} \\
    f_{2,1} & f_{2,2} & \cdots & f_{2,t} \\
    \vdots & \vdots & \ddots & \vdots \\
    f_{t,1} & f_{t,2} & \cdots & f_{t,t} \\
  \end{pmatrix}
```
この行列は $t$ 列なので, 全ての行が線型独立ならば, 係数 $u_0, u_1, \dots, u_t ~(u_i=0 ~\textrm{or}~ 1)$ を求めることができます.

## 掃き出し法の数値例

$N=140086253$ のときに, 以下の7個の関係式が得られたとしましょう
```math
  16891^2 \equiv 2^0 \times 3^5 \times 5^3 \times 7^0 \times 13^2 \times 17^0 \pmod{N} \\
  27665^2 \equiv 2^7 \times 3^3 \times 5^1 \times 7^0 \times 13^1 \times 17^2 \pmod{N} \\
  28367^2 \equiv 2^6 \times 3^4 \times 5^0 \times 7^1 \times 13^2 \times 17^1 \pmod{N} \\
  28549^2 \equiv 2^3 \times 3^3 \times 5^0 \times 7^4 \times 13^1 \times 17^1 \pmod{N} \\
  29448^2 \equiv 2^1 \times 3^1 \times 5^0 \times 7^1 \times 13^3 \times 17^2 \pmod{N} \\
  44299^2 \equiv 2^0 \times 3^5 \times 5^0 \times 7^0 \times 13^0 \times 17^3 \pmod{N} \\
  52940^2 \equiv 2^2 \times 3^8 \times 5^1 \times 7^1 \times 13^0 \times 17^0 \pmod{N} \\
```
このとき以下のような行列が得られます (どの行とどの行を演算したかを記録するために, 最後の列にベクトル名を追加しました):
```math
  \begin{pmatrix}
    0 & 1 & 1 & 0 & 0 & 0 & f_0 \\
    \fbox{1} & 1 & 1 & 0 & 1 & 0 & f_1 \\
    0 & 0 & 0 & 1 & 0 & 1 & f_2 \\
    1 & 1 & 0 & 0 & 1 & 1 & f_3 \\
    1 & 1 & 0 & 1 & 1 & 0 & f_4 \\
    0 & 1 & 0 & 0 & 0 & 1 & f_5 \\
    0 & 0 & 1 & 1 & 0 & 0 & f_6 \\
  \end{pmatrix}
```
以下, この行列に掃き出し法を適用します. まず1列目に着目すると, 非ゼロの成分は2行目にあります(四角で囲まれた値) ので, 4行目と5行目に2行目を加えることで, 行列は以下のようになります. 演算は $\bmod{2}$ で考えていることに注意してください. 2 行目は以後使用しないので, わかるように色をつけています:
```math
  \begin{pmatrix}
    0 & \fbox{1} & 1 & 0 & 0 & 0 & f_0 \\
    \color{red}{1} & \color{red}{1} & \color{red}{1} & \color{red}{0} & \color{red}{1} & \color{red}{0} & \color{red}{f_1} \\
    0 & 0 & 0 & 1 & 0 & 1 & f_2 \\
    0 & 0 & 1 & 0 & 0 & 1 & f_1+f_3 \\
    0 & 0 & 1 & 1 & 0 & 0 & f_1+f_4 \\
    0 & 1 & 0 & 0 & 0 & 1 & f_5 \\
    0 & 0 & 1 & 1 & 0 & 0 & f_6 \\
  \end{pmatrix}
```
次に2列目に着目すると, 非ゼロの成分は1行目にある (四角で囲まれた値) ので, 6行目に1行目を加えることで, 行列は以下のようになります (2行目は使用しない行なので手を加える必要はありません).
```math
  \begin{pmatrix}
    \color{red}{0} & \color{red}{1} & \color{red}{1} & \color{red}{0} & \color{red}{0} & \color{red}{0} & \color{red}{f_0} \\
    \color{red}{1} & \color{red}{1} & \color{red}{1} & \color{red}{0} & \color{red}{1} & \color{red}{0} & \color{red}{f_1} \\
    0 & 0 & 0 & 1 & 0 & 1 & f_2 \\
    0 & 0 & \fbox{1} & 0 & 0 & 1 & f_1+f_3 \\
    0 & 0 & 1 & 1 & 0 & 0 & f_3+f_4 \\
    0 & 0 & 1 & 0 & 0 & 1 & f_0+f_5 \\
    0 & 0 & 1 & 1 & 0 & 0 & f_6 \\
  \end{pmatrix}
```
さらに3列目に着目すると, 非ゼロ成分は4行目にある (四角で囲まれた値) ので, 5行目, 6行目, 7行目に4行目を加えることで以下の行列を得ます. 
```math
  \begin{pmatrix}
    \color{red}{0} & \color{red}{1} & \color{red}{1} & \color{red}{0} & \color{red}{0} & \color{red}{0} & \color{red}{f_0} \\
    \color{red}{1} & \color{red}{1} & \color{red}{1} & \color{red}{0} & \color{red}{1} & \color{red}{0} & \color{red}{f_1} \\
    0 & 0 & 0 & \fbox{1} & 0 & 1 & f_2 \\
    \color{red}{0} & \color{red}{0} & \color{red}{1} & \color{red}{0} & \color{red}{0} & \color{red}{1} & \color{red}{f_1+f_3} \\
    0 & 0 & 0 & 1 & 0 & 1 & f_3+f_4 \\
    0 & 0 & 0 & 0 & 0 & 0 & f_0+f_1+f_3+f_5 \\
    0 & 0 & 0 & 1 & 0 & 1 & f_6 \\
  \end{pmatrix}
```
このとき6行目は全ての成分がゼロになっているので, $f_0+f_1+f_3+f_5=0 \bmod{2}$ という式が得られます. 実際
```math
  f_0+f_1+f_3+f_5=(10,16,4,4,4,6)
```
なので, 
```math
  \begin{eqnarray}
    x &=& 16891 \times 27665 \times 28549 \times 44299 \bmod{140086253}=20820499 \\
    y &=& 2^5 \times 3^8 \times 5^2 \times 7^2 \times 13^2 \times 17^3 \bmod{N} = 119265754
  \end{eqnarray}
```
となりますが, 
```math
  \begin{eqnarray}
    \gcd(x+y,N) &=& (20820499+119265754,140086253) \\
                &=& \gcd(140086253,140086253) \\
                &=& 140086253
  \end{eqnarray}
```
となり, 素因数分解に失敗します. そこで行列の掃き出しに戻り, 3行目に着目することで以下の行列を得ます.
```math
  \begin{pmatrix}
    \color{red}{0} & \color{red}{1} & \color{red}{1} & \color{red}{0} & \color{red}{0} & \color{red}{0} & \color{red}{f_0} \\
    \color{red}{1} & \color{red}{1} & \color{red}{1} & \color{red}{0} & \color{red}{1} & \color{red}{0} & \color{red}{f_1} \\
    \color{red}{0} & \color{red}{0} & \color{red}{0} & \color{red}{1} & \color{red}{0} & \color{red}{1} & \color{red}{f_2} \\
    \color{red}{0} & \color{red}{0} & \color{red}{1} & \color{red}{0} & \color{red}{0} & \color{red}{1} & \color{red}{f_1+f_3} \\
    0 & 0 & 0 & 0 & 0 & 0 & f_2+f_3+f_4 \\
    0 & 0 & 0 & 0 & 0 & 0 & f_0+f_1+f_3+f_5 \\
    0 & 0 & 0 & 0 & 0 & 0 & f_2f_6 \\
  \end{pmatrix}
```
すると5行目より新たに $f_2+f_3+f_4=(10,8,0,6,6,4)$ が得られるので
```math
  \begin{eqnarray}
    x &=& 28367 \times 28549 \times 29448 \bmod{140086253} = 23778411 \\
    y &=& 2^5 \times 3^4 \times 7^3 \times 13^3 \times 17^2 \bmod{140086253} = 83479911
  \end{eqnarray}
```
となり,
```math
  \begin{eqnarray}
    \gcd(x+y,N) &=& \gcd(23778411+83479911,140086253) \\
                &=& \gcd(107258322,140086253) \\
                &=& 10559
  \end{eqnarray}
```
より素因数 $10559$ を得ることができました.

## 関係式の収集方法

これまでの説明により, 十分な個数の関係式を収集できれば, 行列計算によって平方差で必要な合同式を求められることがわかりました. このため次は関係式の集め方が課題となります. Dixon法, 連分数法, 二次篩法などの素因数分解法はすべて平方差法を利用しており, 関係式の生成法の違いが方式の違いになっています (が, 具体的な生成法は別の文献を参照ください).

# 拡張関係式

平方差法によって巨大な合成数 $N$ を素因数分解する場合, 関係式の左辺が平方数に限定されている点が足かせになってきます. そこで, 合同式の右辺を平方数からべき積に条件緩和したのと同様に, 左辺も平方数からべき積に条件緩和した
```math
  p_1^{e_1} \times p_2^{e_s} \times \dots \times p_s^{e_s} \equiv q_1^{f_1} \times q_2^{f_2} \times \dots \times q_t^{f_t} \pmod{N}
```
という形の関係式を集めることで, より大きな合成数の素因数分解が可能となります. このような関係式を拡張関係式と呼ぶことにします (が, これはこの一連の記事だけの造語であり, 通常は単に関係式と呼ばれています).

十分な個数の拡張関係式が得られた場合, べき指数から生成されるベクトル $(e_1,e_2,\dots,e_s,f_1,f_2,\dots,f_t) \bmod{2}$ を考えて, そのベクトルたちが生成する行列に対して掃き出し法を適用することで, 関係式の場合と同様に合同式を得ることができます.

## 拡張関係式の収集法

最後は拡張関係式の収集法です. 数体篩法や Schnorr 法は拡張関係式を利用した素因数分解法です. また, 新しい量子素因数分解アルゴリズムである SQIF も拡張関係式を利用しています. 収集法の違いが方式の違いになっています. Schnorr 法と SQIF の拡張関係式の収集法は次回以降の記事で紹介します.

# まとめ

本記事では平方差法による素因数分解法の原理を紹介し, 素因数分解問題が関係式/拡張関係式を大量に求める問題に変換できることを説明しました. また, 関係式/拡張関係式から平方差に必要な合同式を具体的に求める手順も紹介しました.