Pell's Equation（ペル方程式）〜Quantum Zooやっていく〜 準備編【量子コンピューター Advent Calendar 2022 10 日目】

Dec. 10, 2022, 9:06 a.m. edited Jan. 2, 2023, 1:08 p.m.

#Quantum Zoo  #量子情報  #量子力学 

量子コンピューター Advent Calendar 2022 10 日目の記事です。昨年書いた記事：量子回路で色々な内積計算【量子コンピューター Advent Calendar 2021 15 日目】

有名な量子アルゴリズムである Grover の探索アルゴリズムや Shor の素因数分解アルゴリズムはググればたくさん日本語記事が出てくる一方、 Hallgren の Pell 方程式を解くアルゴリズムはほぼ存在しない。そこで、昨年 7 月に本ブログの企画である〜Quantum Zooやっていく〜に沿って挑戦してみたのだが難しすぎて挫折していた。それから完全に放置していたのだが、今年の量子アドベントカレンダーのネタとして良いだろうと思い再挑戦したというもの。したがって、投稿が当日夜ギリギリになったのも致し方ない（なお投稿当日は後述する準備 1. まで記述していた）。そして普段付属している Qiskit コードが次回の実装編に回るのも仕方ない。

背景

ある非平方数の \(d\) が与えられたとき、整数 \(x\), \(y\) についての方程式 \(x^2-dy^2=1\) をペル方程式という。この解 \((x,y)\) のうち、 \(x+y\sqrt{d}\) を最小にするものを \((x_1,y_1)\) とすると、他の解は \(x_k+dy_k=(x_1+dy_1)^k\) (\(k=1,2,3,\dots\)) を満たす。また、 \(d\) は非平方数であるため \(x+dy\) と \((x,y)\) は一対一対応にできる。ゆえに、ペル方程式を解くとは、 \(x_1+dy_1\) (これを fundamental solution という) を見つけるということになる。

古典的によく知られたペル方程式の解法として連分数展開を用いたものがある。

ここで、 \((x_1,y_1)\) を表すには一般に \(d\) と違い指数関数的に多くのビットを必要とする。これは例えば

\(d\)	\(x_1\)	\(y_1\)
2	3	1
29	9801	1820
6009	131634010632725315892594469510599473884013975	1698114661157803451688949237883146576681644
16383	128	1

のように、 \(d\) と比べて \(x_1,y_1\) が指数関数的に大きくなる場合があるためである。ゆえに、そのような \((x_1,y_1)\) を多項式時間で発見することはできない。その代わり、 \(R=\ln(x_1+y_1\sqrt{d})\) (これを regulator という) を \(n\) 桁の精度で計算するアプローチがとられる（\(R\) に最も近い整数 \(\lfloor R\rceil\) から \((x_1,y_1)\) を求めることができる¹）。この \(R\) を見つける最速の古典アルゴリズムの時間計算量は \(O(e^\sqrt{\log d}\mathrm{poly}(n))\) である一方、 Hallgren は \(O(\mathrm{poly}(\log d,n))\) の時間計算量で見つける量子アルゴリズムを開発した。

問題

背景で説明した \(R\) を求める。

方針

いつもはこの次にアルゴリズムの節が来るが、今回はかなり難解なのでまずは大きな方針から始める。本稿は原著論文 [S. Hallgren, 2002] および、それを代数学の背景含めて自己完結にまとめられた [R. Jozsa, 2003] に従っている。

今回のアルゴリズムでも Shor の素因数分解アルゴリズムや離散対数問題アルゴリズムと同様に、周期が \(R\) である周期関数を利用して量子フーリエ変換からの周期発見で \(R\) を見つけるというものである。しかし、これまでのアルゴリズムと異なり、そのような周期関数を構成するのに割と高度な前準備が必要となる。ゆえに、途中で迷子にならないように、本節で周期関数を構築するまでの大まかな方針を立てる。なお、この方針内で色々用語が出てくるが、それらは準備の節で出てくる。

準備の方針として

1. ペル方程式の fundamental solution を求めることと、 quadratic number field \(\mathbb{Q}[\sqrt{d}]\) の代数的整数の fundamental unit \(\epsilon_0\) を見つけることは等価である
2. \(\mathbb{Q}[\sqrt{d}]\) の principal reduced ideals である \(\{J_i\}_{i\in\mathbb{Z}}\) は周期的であり、その周期に \(\epsilon_0\) が関係
3. 距離および長距離を移動するための大ジャンプを導入し、そして周期関数 \(h\) を構成

という流れでいく。今回のアドベントカレンダーである本稿は 1., 2. をカバーし、次回の実装編でそれ以降を書く予定。なお、各種証明は先述の [R. Jozsa, 2003] を参照してほしい。

準備 1. ペル方程式を解くことと quadratic number field \(\mathbb{Q}[\sqrt{d}]\) の代数的整数の fundamental unit \(\epsilon_0\) を求めることは等価

ペル方程式の解よりも広い集合として quadratic number field

$$\mathbb{Q}[\sqrt{d}]=\{r_1+r_2\sqrt{d} \mid r_1,r_2\in\mathbb{Q}\}$$

を定義する。広い、といったのは、ペル方程式の解は整数 (\(\mathbb{Z}\)) のみに対し、これは有理数 (\(\mathbb{Q}\)) にて定義されているためである。この部分集合として

$$\mathcal{O}=\{\xi \mid \xi^n+a_{n-1}\xi^{n-1}+\cdots+a_1\xi+a_0=0,a_0,a_1,\dots,a_{n-1}\in\mathbb{Z},\xi\in\mathbb{Q}[\sqrt{d}]\}$$

と定義する。この要素 \(\xi\) を代数的整数 (algebraic integer) という（なお、\(\xi=x+y\sqrt{d}\) のとき \(\overline{\xi}=x-y\sqrt{d}\)）。すると、驚くべきことに、

$$\mathcal{O}=\{m+n\omega \mid m,n\in\mathbb{Z}\}$$

ただし、

$$ \omega=\begin{cases} \frac{-1+\sqrt{d}}{2} & (d\equiv 1 \pmod 4) \\ \sqrt{d} & (d\equiv 2,3 \pmod 4) \end{cases} $$

が成り立つ。ここで、自身の逆数もまた代数的整数となる代数的整数を unit という。すると、

\(\xi=x+y\sqrt{d}\in\mathcal{O}\) は unit \(\Leftrightarrow\) \(2x\in\mathbb{Z}\) かつ \(x^2-dy^2=\pm 1\)

が成り立つ。さらに、 \(\mathcal{O}\) の 1 より大きい最小の unit を \(\epsilon_0\) (これを fundamental unit という) とすると、すべてのその unit は \(\{\pm\epsilon_0^k \mid k\in\mathbb{Z}\}\) と表される。

ここまでの準備を踏まえると、次のことがいえる。

\(d\equiv 2,3\pmod 4\) つまり \(\mathcal{O}=\mathbb{Z}[\sqrt{d}]\) のとき
\(\epsilon_0=m+n\sqrt{d}\) について、
\(m^2-dn^2=1\) ならば、それすなわちペル方程式の解。
\(m^2-dn^2=-1\) ならば、 \(\epsilon_1=\epsilon_0^2\) について \((2m^2+1)^2-d(2mn)^2=\epsilon_1\overline{\epsilon}_1=\epsilon_0^2\overline{\epsilon}_0^2=(-1)^2=1\) より、ペル方程式の解が得られる。
\(d\equiv 1 \pmod 4\) つまり \(\mathcal{O}=\mathbb{Z}\left[\frac{-1+\sqrt{d}}{2}\right]\) のとき
\(\epsilon_0=m+n\sqrt{d}\) について、
\(m,n\) が \(\frac{1}{2}\) の倍数になることがあるが、何乗かすればペル方程式の解が得られる。

したがって、ペル方程式と解くことと、 \(\mathbb{Q}[\sqrt{d}]\) の代数的整数の fundamental unit \(\epsilon_0\) を見つけることは等価である。

そして、 regulator \(R\) は \(R=\ln\epsilon_0\) と定義される。

準備 2. イデアルおよび principle cycle の導入

ここからはイデアルというものを導入して、周期関数の構築に迫っていく。

まず、\(A,B\subseteq\mathcal{O}\) について

$$A\cdot B=\{a_1b_1+\cdots a_nb_n \mid a_1,\dots,a_n\in A,b_1,\dots,b_n\in B,n\in\mathbb{N}\}$$

と定義する（例えば、仮に \(A=\{p, q\}, B=\{r,s\}\) とすると、 \(A\cdot B=\{pr,ps,qr,qs,pr+qs,ps+qr\}\) ということ）。すると、

\(I\cdot\mathcal{O}=I\) かつ \(\alpha,\beta\in I\) ならばすべての \(m,n\in\mathbb{Z}\) について \(m\alpha+n\beta\in I\) を満たす

• \(I\subseteq\mathcal{O}\) を integral ideal
• \(I\subseteq\mathbb{Q}[\sqrt{d}]\) を fractional ideal

という。このとき、

• \(\gamma\in\mathcal{O}\) ならば \(\gamma\mathcal{O}=\{\gamma\xi \mid \xi\in\mathcal{O}\}\) は integral ideal となる
• \(\gamma\in\mathbb{Q}[\sqrt{d}]\) ならば \(\gamma\mathcal{O}=\{\gamma\xi \mid \xi\in\mathcal{O}\}\) は fractional ideal となる

また、このような \(\gamma\mathcal{O}\) を principal ideal という。

ここまででだいぶ謎概念が出てきていてお腹いっぱい感があるが、馴染みある形になるまでもう少し。

\(I\) を fractional ideal とすると \(\alpha\in I,\alpha>0\) について、 \(|\beta|<|\alpha|\) かつ \(|\overline{\beta}|<|\overline{\alpha}|\) となる \(0\neq\beta\in I\) が存在しないとき、 \(\alpha\) を \(I\) の minimum という。また、 \(I\) の minimum に 1 があるとき、 \(I\) は reduced であるという。

すると、 \(I\) が reduced であるとき、その標準形は

$$I=\mathbb{Z}+\frac{b+\sqrt{D}}{2a}\mathbb{Z}$$

ただし、 \(a,b\in\mathbb{Z}\) かつ \(0<a<\sqrt{D}\) かつ \(b<\sqrt{D}\) を満たす（最後の \(b\) についてはもう少し条件があるが割愛）。それらにより、 reduced ideal の個数は有限である。

こうして、 reduced fractional ideal とペル方程式の解に近い形があることが見えてきた感じがする。

次に、 \(\gamma(I)=\frac{b+\sqrt{D}}{2a}\) とし、 \(\rho(I)=\frac{1}{\gamma(I)}I\) とする。すると、 \(4a\) は \(b^2-D\) を割り切るなどの条件を利用してうまく式変形し、新たに \(a',b'\) を定義することで

$$\rho(I)=\mathbb{Z}+\frac{b'+\sqrt{D}}{2a'}\mathbb{Z}$$

と再び標準形を得ることができる。こうして、 principal reduced ideal である \(J_i,i\in\mathbb{Z}\) について \(J_{i+1}=\rho(J_i)\) としてやると、 \( \{J_i\}_{i\in\mathbb{Z}}\) は周期的となる。つまり、 \(J_i= J_j \Leftrightarrow i\equiv j\pmod {k_0}\) をすべての \(i,j\in\mathbb{Z}\) について満たす最小の \(k_0\) が存在するということ。この \(\{J_0=\mathcal{O}, J_1, \dots, J_{k_0-1}\}\) を principal cycle という。さらに、 \(J_i=\alpha_i\mathcal{O}=\mathbb{Z}+\gamma_i\mathbb{Z}\) （ただし、 \(\alpha_0=1\in\mathcal{O}\)）を満たす \(\alpha_i\) を考えることができ、 \(\alpha_{k_0}=\epsilon_0\) を満たす（\(\mathcal{O}=J_{k_0}=\alpha_{k_0}\mathcal{O}\)）。また、 \(I\) を reduced principal fractional ideal とすると、任意の \(I\) は principal cycle に含まれる（つまり、どれかの \(k\) について \(I=J_k\) となる）。

一気に多くの概念が押し寄せた感があるが、こうして周期性が得られ、かつそこに準備 1. で出てきた fundamental unit \(\epsilon_0\) が関係しているのでとても重要そうにみえる。

準備 3. 距離の導入、そして \(\ast\) 演算により多項式時間で計算できる周期 \(R\) の関数 \(h\) の構築

周期関数を与える前に、それぞれの fractional principal ideal の間の距離を考える。 \(I_1,I_2\) が fractional principal ideal で \(I_1=\gamma I_2\) であるとき、その間の距離 \(\delta(I_1,I_2)\) を

$$\delta(I_1,I_2)=\ln|\gamma| \mod R$$

と定義する。特に、 unit ideal \(\mathcal{O}\) との距離は \(\delta(I)=\delta(\mathcal{O},I)\) と略記する。すると、 principle cycle に含まれる \(J_i=\alpha_i\mathcal{O}=\mathbb{Z}+\gamma_i\mathbb{Z}\) について

$$\delta(J_i)=\ln |\alpha_i| \mod R$$

が成り立ち、そして

$$\frac{3}{32D}\leq \delta(J_i,\rho(J_i))=\ln\gamma_i\leq\frac{1}{2}\ln D$$

が成り立つ。

ここで、 Hallgren のアルゴリズムで用いる周期関数として、

$$h(x)=(I_x, \tilde{x}-\delta(I_x))$$

と定義する。ただし、 \(\tilde{x}\equiv x\mod R\), \(0\leq\tilde{x}<R\)、 そして、 \(I_x\) は \(\delta(I_x)<\tilde{x}\) を満たしつつ距離 \(\delta(I_x)\) が最も大きくなる reduced principal fractional ideal である。

この周期関数は何者なのかと思うが、第一成分の \(I_x\) は \(I_0=\mathcal{O}\) から始まり、 \(I_R=I_0=\mathcal{O}\) で戻ってくるので、確かに我々が欲しい周期関数のように感じられる（それなら第一成分のみで十分に感じるが、第二成分は \(h\) が各周期で一対一となるために必要らしい）。

したがって、あとは与えられた \(x\)（\(\tilde{x}\)）に対応する \(I_x=J_i\) が得られれば良いのだが、実は \(\delta(J_i)=\ln |\alpha_i|\mod R\) の \(\alpha_i\) は多項式時間で計算できない。一方、 \(\delta(J_i, \rho(J_i))=\ln \gamma_i\) は多項式時間で計算できるので、これを \(i=1\) から繰り返していく方法も考えられるが、添字の最後である \(k_0\) は指数関数的に大きいので、結局多項式時間で計算できなくなってしまう。

そこで、そのような大ジャンプを実現するために \(\ast\) 演算を導入する。 \(I=\mathbb{Z}+\frac{b+\sqrt{D}}{2a}\mathbb{Z}\) を principle cycle 内の reduced ideal とし、これを用いて

$$I_2=I\cdot I=\frac{1}{k'}\left(\mathbb{Z}+\frac{b'+\sqrt{D}}{2a'}\mathbb{Z}\right)$$

とする（\(k',a',b'\) は適切に定義する）。これは reduced ideal の標準形とは違うので、括弧の中身を抜き出すと

$$I_2'=\mathbb{Z}+\frac{b'+\sqrt{D}}{2a'}\mathbb{Z}$$

が得られる。実はこの \(I_2'\) は一般に reduced ではないのだが、 \(\rho\) を \(\left\lceil\log_2\frac{a}{\sqrt{D}}\right\rceil+1\)（つまり \(a<\sqrt{D}\) より定数のオーダー）の回数作用させると reduced となる。その reduced ideal を \(I_2''\) とすると、 principal cycle に含まれるようになる。そして、この \(I_2''\) に \(\rho\) を \(2n\) 回（\(n<\frac{3\ln D}{2\ln 2}\)）作用させると \(\delta(J_k)>2\delta(I)\) を満たす最初の \(J_k\) が得られる。ここまでの演算を \(J_k=I\ast I\) と表すこととする。この \(I\ast I\) にかかる時間計算量は \(O(\mathit{poly}\log D)\) と表される。

さらに、 reduced principle ideal \(I\) について \(\ast-\)squaring を

$$I\to I^{(2)}=I\ast I\to I^{(4)}=I^{(2)}\ast I^{(2)}\to\cdots\to I^{(2^n)}=I^{(2^{n-1})}\ast I^{(2^{n-1})}$$

とする。こうして得られる \(I^{(2^n)}\) の距離は \(\delta(I^{(2^n)})>2^n\delta(I)\) と指数関数的に大きくなり、そのうえ \(I^{(2^n)}\) は \(O(\mathit{poly}\log D,n)\) で得られる。この大ジャンプを利用することで \(h\) を多項式時間で計算可能となる。より詳しくいうと、 \(x\) について \(10^{-n}\) の精度で与えると、 \(\tilde{x}-\delta(I_x)\) は \(10^{-n}\) の精度で \(\mathit{poly}(\log D,\log x,n)\) の時間計算量で計算できる。

その方法としては、まず \(\mathcal{O}\) として \(\mathcal{O}=\mathbb{Z}+\frac{\tau(D,2)+\sqrt{D}}{2}\mathbb{Z}\) を用いる。ただし、 \(\tau(D,2)\) とは \(\tau\equiv D \mod 2\cdot 2\) および \(-2<\tau\leq 2\) \(\text{if}\) \(2>\sqrt{D}\) \(\text{otherwise}\) \(\sqrt{D}-2\cdot 2<\tau\leq\sqrt{D}\) を満たす整数 \(\tau\) のことである。この \(\mathcal{O}\) に \(\rho\) を 2 回作用させて \(I_0=\rho^2(\mathcal{O})\) を作ると、 \(\ln D>\delta(I_0)>\ln 2\) を満たす。ここから、 \(\ast-\)squaring により \(I_{k+1}=I_k\ast I_k\) としていくと、

$$2\delta(I_{k-1})\leq\delta(I_k)\leq 2\delta(I_{k-1})+\frac{1}{2}\ln D$$

を満たす。ゆえに、 \(\delta(I_k)\geq\delta(I_0)2^k\) となる。この繰り返しは \(\delta(I_{N+1})>x\) と初めてなるような \(k=N\) で終了する。この \(N\) は \(N<\left\lceil\log_2\frac{x}{\delta(I_0)}\right\rceil\) を満たす。ここからさらにギリギリまで \(x\) に近づけるため、まず \(J=I_N\) とする。そして \(\delta(J\ast I_{k-1})< x\) ならば \(J:=J\ast I_{k-1}\) で更新し \(k:=k-1\) と更新する。 \(\delta(J\ast I_{k-1})>x\) であったならば \(k:=k-1\) と更新するのみ。こうして \(k=1\) となるまで繰り返し（繰り返し回数は \(N\) 回）、 \(x-\delta(J)\leq\delta(I_0)+\frac{1}{2}\ln D\leq \frac{3}{2}\ln D\) を満たす \(J=J_\ast\) を得られる。それから最後に \(\rho^2\) を何回か（多くとも \(\frac{3\ln D}{2\ln 2}\) 回） \(J_\ast\) に作用させ、 \(x\) を超える直前の \(J\) が求める \(I_x\) となる。

ここまでをまとめると、以下の図のようになるはず：

次回予告

次回は Qiskit で実際に動く Python コードとして実装してやりたいところ。・・・そんなこと本当にできるのか？？？（実は \(h\) をそのまま実装するわけではないのである→擬周期的な関数の離散化問題）

次回更新は、未定（2023 年中にはしたい）！