Discrete Log(離散対数問題) 〜Quantum Zooやっていく〜

July 19, 2021, 3:14 p.m. edited Dec. 4, 2021, 2:42 p.m.

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

$$ \def\bra#1{\mathinner{\left\langle{#1}\right|}} \def\ket#1{\mathinner{\left|{#1}\right\rangle}} \def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}} $$

まさかの第二回。

なお、いきなり結論に向かいたい読者は アルゴリズム(測定を最後にした場合) および 測定で得られた整数値から \(d\) を計算(おそらく完全版)、実装も見る場合は Qiskit による実装 を見れば良い。

では、なぜこのような構成になっているかというと、本記事で最初に紹介するアルゴリズムは色々なところで Shor の離散対数問題を解くアルゴリズムとして載っているため、最初に参考にして書いたものである。にもかかわらず、途中に測定があるために Qiskit での実装が難しく、そのために測定を最後におこなうバージョンとして書いたものが次のアルゴリズム。しかしながら、ここまでの方法では実は求まらないという事実が判明し、結局 Shor の原論文に立ち向かうのが一番の近道だったと判明するのが後半である。要は、私の辿った道をともに歩んでほしいという文学的(?)なノリである。

背景

Shor のアルゴリズム における周期発見アルゴリズムの応用先の一つ。最速の古典アルゴリズムでは \(n\) (定義は次節の問題にて)に関する超多項式時間1 (Shor の論文内には古典では \(\exp(O(\log p)^{1/3}(\log\log p)^{2/3})\) の時間計算量がかかるとある) かかるのに対し、 Peter Shor により \(\mathrm{poly}(n)\) 時間の計算量で解ける量子アルゴリズムが得られた。

問題

\(p\) を \(n\) ビットの素数とし2、 \(Z_p^\ast=\{1,\dots,p-1\}\) とする。 \(\alpha\) を \(Z_p^\ast\) の生成子とし3、 \(\beta\in Z_p^\ast\) とする。 \(p\), \(\alpha\), \(\beta\) が与えられたとき、 \(\alpha^d=\beta\) を満たす \(d\in[0, p-1]\) を求める4

アルゴリズム

内容は基本的に Shor による周期発見アルゴリズムと同じである。ただし、二変数関数でそれぞれに周期が存在するというところが異なる。今回こそ [P. Shor, 1996] (p.20) に従いたかったが、よくわからなかったので、楕円曲線への適用をした [J. Proos and C. Zalka, 2004] のレビューにあった内容 (2章) に従う。このアルゴリズムは 2 つの量子アルゴリズムから構成される。

1 つ目の量子アルゴリズム: \(\alpha^q=1\) を満たす \(q\) を探す

本題の \(d\) を見つける前に、まずはじめに \(\alpha^q=1\) を満たす \(q\) を見つける必要がある。すなわち、 \(\alpha^x=\alpha^{x+q}\) という周期関数であるため、前回のショアの Factoring アルゴリズムにより得ることができる。

2 つ目の量子アルゴリズム: \(\alpha^d=\beta\) を満たす \(d\) を探す

1 つ目の量子アルゴリズムにより得られた \(q\) を考えると、 \(0\leq d\leq q\) となる。そのうえで \(f(x,y)=\alpha^x\beta^y\) を整数 \(x\), \(y\) について定義する。この関数 \(f\) は 2 つの独立な周期をもつ5

$$f(x+q,y)=f(x,y),\ f(x+d,y-1)=f(x,y)$$

この関数 \(f\) を用いて、重ね合わせ状態

$$\frac{1}{q}\sum_{x=0}^{q-1}\sum_{y=0}^{q-1}\ket{x}\ket{y}\ket{0}$$

から

$$\frac{1}{q}\sum_{x=0}^{q-1}\sum_{y=0}^{q-1}\ket{x}\ket{y}\ket{\alpha^x\beta^y}\tag{1}$$

を得る。ここで、最後のレジスタ \(\ket{\alpha^x\beta^y}\) に測定をし、 \(\alpha^{x_0}\) が得られたとする。すると、

$$\alpha^x\beta^y=\alpha^x(\alpha^d)^y=\alpha^{x_0}$$

より、 \(x+dy\equiv x_0\pmod q\) 、つまり、 \(x=(x_0-dy)\pmod q\) となる。したがって、状態は

$$\frac{1}{\sqrt{q}}\sum_{y=0}^{q-1}\ket{x_0-dy}\ket{y}$$

となる。ここに量子フーリエ変換を適用し、

$$ \begin{align} &\frac{1}{\sqrt{q}}\frac{1}{q}\sum_{x',y'=0}^{q-1}\sum_{y=0}^{q-1}e^{2\pi i((x_0-dy)x'+yy')/q}\ket{x'}\ket{y'}\tag{2}\\ =&\frac{1}{\sqrt{q}}\sum_{x'=0}^{q-1}e^{2\pi ix_0 x'/q}\ket{x'}\ket{y'\equiv dx'\pmod q}\tag{3} \end{align} $$

を得る。ここで、式(2)から式(3)への変形では

$$\sum_{a=0}^{n-1}e^{\frac{2\pi i}{n}ab}=n\delta_{b,0}$$

を用いている(\(\delta\) はデルタ関数、\(a\), \(b\), \(n\) は整数)。つまり、 \(dx'-y'\) が \(q\) で割り切れる場合にのみ \(y\) についての和の結果は \(q\) となり、そのときの \(y'\) の値は \(y'\equiv dx'\pmod q\) を満たす値でなくてはならない。そしてその値が 2 番目のレジスタに入っていることになる。

これを測定し、得られた \(x'\neq 0\), \(y'\) を用いて

$$d=y'(x')^{-1}\bmod q$$

と求まる。・・・のだが、次節で測定を最後に移動した場合6を書く。

アルゴリズム(測定を最後にした場合)

次に実装例を書こうと思ったら Qiskit で途中で測定する方法がよくわからなかったので、測定を最後に持っていった場合の 2 つ目(というか 1 つ目はもはや不要)の量子アルゴリズムを記す7。要は Factoring のときと同じ感じになる。

式(1)から始めるが、どうせ量子回路組むときはビット数のべき乗 \(N\) を分母にすることになるので、もうそのように書くと

$$\frac{1}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1}\ket{x}\ket{y}\ket{\alpha^x\beta^y}$$

量子フーリエ変換を適用し、

$$\frac{1}{N^2}\sum_{x'=0}^{N-1}\sum_{y'=0}^{N-1}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1}e^{2\pi i(xx'+yy')/N}\ket{x'}\ket{y'}\ket{\alpha^x\beta^y}$$

3 つのレジスタ \(\ket{x'}\), \(\ket{y'}\), \(\ket{\alpha^{x_0}\beta^{y_0}}\) が測定される確率は

$$ \begin{align} &\frac{1}{N^4}\left|\sum_{x,y\ \text{where}\ \alpha^x\beta^y=\alpha^{x_0}\beta^{y_0}}e^{2\pi i(xx'+yy')/N}\right|^2\tag{4}\\ =&\frac{1}{N^4}\left|\sum_k\sum_\ell e^{2\pi i\{(x+kq+\ell d)x'+(y-\ell)y'\}/N}\right|^2\tag{5} \end{align} $$

ここで、式(4)から式(5)への変形には上で述べた 2 つの独立な周期を用いている。 Factoring のときと同様に測定で得られるのは \(\frac{qx'}{N}\) が整数かつ \(\frac{dx'-y'}{N}\) が整数を満たすものがほとんどなので、特に後者の条件より

$$dx'\equiv y'\pmod N$$

\(\gcd(x',N)=1\) のとき

$$d=y'(x')^{-1}\pmod N\tag{6}$$

より \(d\) は求まる。式(6)の解き方としては、拡張ユークリッド互除法より

$$x'd_0+Nd_1=1$$

を満たす整数 \(d_0,d_1\) を求め、その両辺に \(y'\) を掛けることで

$$x'd_0y'+Nd_1y'=y'$$

より

$$x'd_0y'\equiv y'\pmod N$$

であるから、 \(d=d_0y'\) と求まる。

これでは不完全なことがわかる例

\(p=7\) として \(Z_7^\ast=\{1,2,3,4,5,6\}\) とする。生成子 \(\alpha=3\), \(\beta=6\) とすると、答えは \(\alpha^3\equiv\beta\pmod 7\) より、 \(d=3\) となる。これを shor のアルゴリズムにより求める。

今回のアルゴリズム(Qiskit プログラムは後述)を動かすことにより以下のプロットを得る。


(やはり若干読みづらいが、例えば右端 100111 は \(\ket{x'}\ket{y'}=\ket{111}\ket{100}=\ket{7}\ket{4}\) を表す)

最も大きいピークが 2 本、次点で大きいピークが 4 本存在する。そう、これも前回と同じである。ゆえに、やはり \(u,v\) について表を書くと次のようになる。この表は簡単に

def f(u,v):
    x = u * 8 / 6
    y = 3 * x - 8 * v
    return x,y

for u in range(20):
    for v in range(20):
        x,y = f(u,v)
        if 0<=x<8 and 0<=y<8:
            print(f'(u,v)=({u},{v}): (x,y)=({x},{y})')

で求めた。

\((u,v)\) \(x'\) \(y'\)
(0,0) 0 0
(1,0) 1.333 4
(2,1) 2.666 0
(3,1) 4 4
(4,2) 5.333 0
(5,2) 6.666 4

より、やはり最も大きいピークは \(x',y'\) が整数のところであり、次点はそれに近い小数値となっている。ゆえに、測定で得られた \(x'\) を用いて拡張ユークリッド互除法をおこなっても正しい \(d\) を得ることは難しい。現に、今回のピークの値を用いてもうまく \(d\) を求めることはできない。

では、 \(N=8^2=64\) でおこなうとどうだろうか。再び後述する Qiskit プログラムでピークを描いてみると以下のようになる。

ここで得られた大きなピークを用いて上で述べたように拡張ユークリッド互除法をおこなっても、やはり残念ながら正しく \(d\) を得ることはできない。

測定で得られた整数値から \(d\) を計算(おそらく完全版)

仕方がないので、結局 Shor の原論文を参考にするのがこういうときは良いだろう。論文の p.21 を見ると

$$T=dx'+y'-\frac{d}{p-1}\{x'(p-1)\}_N$$

という式が出てくる(本記事に各記号は合わせた)。ここで、 \(\{z\}_N\) は \(z\pmod N\) のことである。

これは次の不等式を満たす。

$$|\{T\}_N|=|dx'+y'-\frac{d}{p-1}\{x'(p-1)\}_N-jN|\leq\frac{1}{2}$$

ここで、 \(j\) は \(T/N\) に最も近い整数である。ゆえに、

$$ \begin{align} -\frac{1}{2}&\leq dx'+y'-\frac{d}{p-1}\{x'(p-1)\}_N-jN \leq\frac{1}{2}\\ -\frac{1}{2N}&\leq \frac{y'}{N} + d\left(\frac{x'(p-1)-\{x'(p-1)\}_N}{(p-1)N}\right)-j\leq\frac{1}{2N}\\ -\frac{1}{2N}&\leq \frac{y'}{N} + d\left(\frac{x'(p-1)-\{x'(p-1)\}_N}{(p-1)N}\right)\leq\frac{1}{2N}\pmod 1 \end{align} $$

最後の変形では \(\pmod 1\) を用いた。これにより \(j\in\mathbb{Z}\) は \(1\) で割った余りが \(0\) なので消える。ここで、 \(\{x'(p-1)\}_N=x'(p-1)-c'N\) (\(c'\in\mathbb{Z}\)) なので

$$-\frac{1}{2N}\leq \frac{y'}{N} + \frac{dc'}{p-1}\leq\frac{1}{2N}\pmod 1$$

それから \(\frac{y'}{N}\) を連分数展開(など)により分母を \(p-1\) にする(このとき得られる分子を \(y''\in\mathbb{Z}\) とする)

$$-\frac{1}{2N}\leq \frac{y''+dc'}{p-1}\leq\frac{1}{2N}\pmod 1$$

ここから原論文では \(c'\) で割る \(\pmod {p-1}\) とあったのだが、(ちょっとなにいってるかわからなかったので)ここでは \(N\gg 1\) つまり、左辺と右辺がほぼ 0 であることと \(\pmod 1\) を活かして $$\frac{y''+dc'}{p-1}=m\in\mathbb{Z}$$ となるような \(d\) を求めるアプローチをとる。式変形して $$ \begin{align} y''+dc'&=m(p-1)\\ dc'-m(p-1)&=-y'' \end{align} $$ 最後の不定方程式を拡張ユークリッド互除法により解くことで \(d\) が求まる。

Qiskit による実装

ここまでの内容を Qiskit を用いて実装したのが以下になる。コードの全体は GitHub にある

from typing import Optional
import math

import numpy as np
import matplotlib.pyplot as plt
from sympy import gcdex
from qiskit import QuantumCircuit, Aer, transpile, assemble
from qiskit.visualization import plot_histogram
from sympy import Rational
from sympy.ntheory.continued_fraction import continued_fraction, continued_fraction_convergents

from qft import qft
from elementary import ax_modM

def discrete_log(alpha: int, beta: int, p: int, N_len: Optional[int] = None, show_hist: Optional[bool] = True) -> int:
    # find d s.t. alpha^d = beta
    if N_len is None:
        N_len = int(np.ceil(np.log2(p ** 2)))
    N = 2 ** N_len
    qc = QuantumCircuit(11 * N_len - 2, N_len * 2)

    qc.h(range(N_len * 2))

    qc.append(ax_modM(a=alpha, M=p, N_len=N_len), list(range(N_len)) + list(range(N_len * 2, 11 * N_len - 2)))
    qc.append(ax_modM(a=beta, M=p, N_len=N_len, x_0_at_first=False), range(N_len, 11 * N_len - 2))

    qc.append(qft(n=N_len), range(N_len))
    qc.append(qft(n=N_len), range(N_len, N_len * 2))

    qc.measure(range(N_len * 2), range(N_len * 2))

    backend = Aer.get_backend('aer_simulator_matrix_product_state')
    qc = transpile(qc, backend)
    job = backend.run(qc, shots=10000)
    hist = job.result().get_counts()

    if show_hist:
        plot_histogram(hist)
        plt.show()

    for measured_key, _ in sorted(hist.items(), key=lambda x: x[1], reverse=True):
        x = int(measured_key[-N_len:], 2)
        y = int(measured_key[-N_len * 2:-N_len], 2)

        for fraction in continued_fraction_convergents(continued_fraction(Rational(y, N))):
            if fraction.denominator == p - 1:
                yy = fraction.numerator
            elif (p - 1) % fraction.denominator == 0:
                yy = fraction.numerator * ((p - 1) // fraction.denominator)
            else:
                continue

            c = Rational(x * (p - 1) - (x * (p - 1) % N)) / N
            if c == 0:
                continue

            d, m, _ = gcdex(c, 1 - p)
            d *= -yy
            m *= -yy
            if alpha ** d % p == beta:
                return d
            d += 1 - p
            if alpha ** d % p == beta:
                return d

    print('d was not found?!')

if __name__ == '__main__':
    print(discrete_log(alpha=3, beta=6, p=7, N_len=3))

これを実行すると、確かに \(d=3\) が得られる。


  1. 多項式時間ではできないという計算量 

  2. 基本的には \(p\) が素数でない(素数に近くないと)と古典でも Pohlig–Hellman algorithm により簡単に解かれてしまう ので、ここでは素数に限定している 

  3. つまり、 \(\alpha^0, \alpha^1, \dots, \alpha^{p-1}\in Z_p^\ast\) 

  4. これと \(\alpha^d\equiv\beta\pmod p\) は等価。また、より一般的な離散対数問題は代数学でもっと一般化されて表されるらしい 

  5. 1 つ目の周期は \(\alpha^{x+q}=\alpha^{x}\)、 2 つ目の周期は \(\alpha^{x+d}\beta^{y-1}=\alpha^x\beta\beta^{y-1}=\alpha^x\beta^y\) 

  6. in principle, measurements can always be deferred to the end of a quantum algorithm (p.4)

  7. あとで確認したら IPA の文書に載ってるのと大体同じものとなった