はじめに

量子コンピュータの話題になると、「今使われている暗号が解読されてしまうかもしれない」という話をよく耳にします。
少し調べてみると、その理由としてよく登場するのが「ショアのアルゴリズム」と「素因数分解」、そして「RSA暗号」です。

RSA暗号は、公開鍵が2つの素数の積で構成されており、その数を素因数分解して元の素数を求められてしまうと暗号が破られてしまいます。現在の古典コンピュータでは、RSAで用いられるような大きな整数の素因数分解は現実的な時間では困難であるため、この暗号方式は安全に利用されています。

一方で、量子コンピュータを用いると、ショアのアルゴリズムによってこの素因数分解を効率的に行える可能性があることが知られています。そのため、量子コンピュータが実用化されるとRSA暗号が解読されるのではないか、という議論につながっています。

ショアのアルゴリズム自体の概要を解説した記事は多くありますが、その中心となる計算である「べき乗剰余（modular exponentiation）」を量子回路としてどのように実装するのかまで踏み込んで説明した例はあまり多くありません。

この記事では、量子回路の基礎を理解している、あるいは自分で調べながら読み進められる方を対象に、ショアのアルゴリズムを完全に実装するための量子回路の構成方法を、実装レベルまで分解して説明していきます。

ショアのアルゴリズムとは

ショアのアルゴリズムの全体像としては以下となります。

\( qx, qy \) は複数の量子ビットで構成される量子レジスタを示しています。
\( qx \) を重ね合わせの状態にして、oracleという計算処理によって \( qx \) を入力とした計算結果を \( qy \) に設定します。
その後、\( qx \) に対して逆量子フーリエ変換である QFT† を適用し、\( qx \) を測定することで \( qy \) の周期を導くことができます。
例えば \( qy=1 \) となるのが \( x=0,4,8,12 \) で、\( qy=2 \) となるのが \( x=1,5,9,13 \) のような場合、同じ値を繰り返す間隔である \( 4 \) という値がここでいう周期になります。

RSA暗号解読で必要となる素因数分解を行うためのoracleはどのような計算かというと \( qy=a^{qx} \pmod{N} \) です。
任意の \( a \) という定数を用意し、\( a^{qx} \) を公開鍵 \( N \) で割った余りを求めており、この計算を「べき乗剰余」と呼びます。
このべき乗剰余の計算の周期を求めることで、その後の古典計算により \( N \) の素因数分解ができます。

oracleの計算内容をイメージするためにPythonで記述すると以下のような処理となります。

1
def oracle(qx):
2
  return (a**qx) % N

例えば \( a=2, N=15 \) だとすると、ショアのアルゴリズムの全体像をNumPyを使ってPythonで記述すると以下のイメージとなります（実装イメージなので qft_dg や measure という処理が別途あるものと思ってください）。

1
a = 2
2
N = 15
3

4
# 4bit=16個の重ね合わせ
5
qx = numpy.arange(16)
6

7
# oracle
8
qy = oracle(qx)
9

10
# QFT†
11
qx = qft_dg(qx)
12

13
# 測定
14
x = measure(qx)

Pythonで記述すると非常に簡単なプログラムとなりますが、実はこれを量子回路で実装したものはインターネット上を探してもほとんど見つけることができません。
ショアのアルゴリズムは有名であるため、量子回路の代表例として実装されることが多いですが、\( a=2, N=15 \) を量子回路で実装された例はあるものの、\( N \) がこれよりも大きくなった場合の量子回路はほとんど見当たりません。

それはなぜでしょうか？

実は、一見すると簡単に見える oracle の計算処理である \( qy=a^{qx} \pmod{N} \) を、量子回路で実装するのが非常に難しいからです。
この実装が難しい oracle について、 \( N=15 \) よりも大きい任意の \( N \) について計算するための量子回路を実現する方法を説明します。

oracleの実装はなぜ難しいのか

Pythonでは、四則演算や余り計算・べき乗計算などは、プログラム言語やNumPy内部で実現されているため、一般のプログラマーは実装する必要はありません。
しかし量子回路では、1つの量子ビットに対して単純な操作しかできない量子ゲートというものを組み合わせて、足し算や掛け算のような基本的な計算自体も自分で実装します。
量子ゲートはいくつか種類はあるものの、素因数分解で利用する oracle で必要な「べき乗剰余」という計算では、基本的に Xゲートをベースにした回路だけを使います（他の量子ゲートで実現する方法はあるものの、今回は最も分かりやすいXゲートをベースに説明します）。

「Xゲート」は、量子ビットの \( \ket{0} \) と \( \ket{1} \) を反転させるという効果があります。
例えば3量子ビット \( \ket{q_2 q_1 q_0}=\ket{100} \) の状態に対して、\( \ket{q_0}, \ket{q_1}, \ket{q_2} \) に対してXゲートを適用すると以下の変化が起きます。

$$ \begin{align} \ket{100} \xrightarrow{X(\ket{q_0})} \ket{101} \cr \ket{100} \xrightarrow{X(\ket{q_1})} \ket{110} \cr \ket{100} \xrightarrow{X(\ket{q_2})} \ket{000} \end{align} $$

これだけでは計算は無理ですが、量子ゲートには「制御ビット」と呼ばれる「他の量子ビットがすべて \( \ket{1} \) の場合のみ量子ゲートを適用する」という条件式のようなものを付与できるため、それらも組み合わせて実現します。

$$ \begin{array}{ll} \ket{100} \xrightarrow{ctrl(\ket{q_1}) \otimes X(\ket{q_0})} \ket{100} & (\ket{q_1}=\ket{0}なので変化なし) \cr \ket{100} \xrightarrow{ctrl(\ket{q_2}) \otimes X(\ket{q_0})} \ket{101} & (\ket{q_2}=\ket{1}なので \ket{q_0} が反転) \cr \ket{100} \xrightarrow{ctrl(\ket{q_2 q_1}) \otimes X(\ket{q_0})} \ket{100} & (\ket{q_2 q_1}=\ket{10}なので変化なし) \cr \ket{100} \xrightarrow{X(\ket{q_1})} \ket{110} \xrightarrow{ctrl(\ket{q_2 q_1}) \otimes X(\ket{q_0})} \ket{111} & (\ket{q_2 q_1}=\ket{11}なので \ket{q_0} が反転) \end{array} $$

このXゲートと制御ビットの組み合わせで実現できる「あるビットを反転させる」「他の量子ビットがすべて \( \ket{1} \) のときだけ対象ビットを反転させる」という単純な量子ゲートの組み合わせだけで、複雑な計算となる「べき乗剰余」を実装するのが難しいのです。

すべての基本はインクリメント

oracleとして用意する四則演算などの計算を行うためには、まずはインクリメント回路の実現から始めます。

$$ \begin{array}{ccccc} & q_3 & q_2 & q_1 & q_0 \cr {}+ & & & & 1 \cr \hline & q_3 & q_2 & q_1 & q_0 \end{array} $$

インクリメントを2進数の筆算で考えると上記のようになります。
\( q_3 \sim q_0 \) には \( 0 \) または \( 1 \) が入り、計算によって値が変化します。

最下位ビットに注目すると、インクリメントによって \( q_0 \) は \( 0 \to 1 \) あるいは \( 1 \to 0 \) と変わることになり、結果としては単に「\( q_0 \) のビットを反転させる」というXゲートと同じ動きになります。
次に最上位ビットの \( q_3 \) を考えると、このビットが変わるのは桁上がりがあるときだけであり、桁上がりは \( q_2 q_1 q_0 = 111 \) と3つのビットがすべて 1 のときとなります。
これは「\( q_2 q_1 q_0 \) のビットすべてが1のときに \( q_3 \) を反転する」という制御ビットを使った量子ゲートの動きそのものになります。

インクリメント回路としては、最上位ビットから制御ビットを利用した反転を繰り返し、最下位ビットは単なる反転のXゲートを適用することで実現できます。

1
from qiskit import QuantumCircuit
2
from qiskit.circuit.library import XGate
3
from functools import lru_cache
4

5
@lru_cache
6
def inc(nbits):
7
    qc = QuantumCircuit(nbits, name = f"inc_{nbits}")
8
    for i in reversed(range(1, nbits)):
9
        qc.append(XGate().control(i), range(i + 1))
10
    qc.x(0)
11
    return qc.to_gate()

次に足し算を考えてみます。
インクリメントと同じように、2進数の筆算で表現すると以下となります。

$$ \begin{array}{ccccc} & q_3 & q_2 & q_1 & q_0 \cr {}+ & r_3 & r_2 & r_1 & r_0 \cr \hline & q_3 & q_2 & q_1 & q_0 \end{array} $$

実はこれは以下のような分解のイメージで表現が可能です。

$$ \begin{array}{cccccl} & q_3 & q_2 & q_1 & q_0 & \cr {}+ & & & & 1 & (r_0=1の場合) \cr \hline & q_3 & q_2 & q_1 & & \cr {}+ & & & 1 & & (r_1=1の場合) \cr \hline & q_3 & q_2 & & & \cr {}+ & & 1 & & & (r_2=1の場合) \cr \hline & q_3 & & & & \cr {}+ & 1 & & & & (r_3=1の場合) \cr \hline & q_3 & q_2 & q_1 & q_0 & \end{array} $$

この分解が示す意味としては、\( r_0 \sim r_3 \) の各ビットが 1 の場合に、そのビット位置から始まる上位ビットに対してインクリメント回路を実行しているのと同義になります。
そのため、以下のようにインクリメント回路に対して制御ビットを付与したものを組み合わせることで足し算回路が実現できます。

1
@lru_cache
2
def adder_reg(xbits, ybits):
3
    qc = QuantumCircuit(xbits + ybits, name = f"qy_{ybits}=qy_{ybits}+qx_{xbits}")
4
    for i in range(min(xbits, ybits)):
5
        qc.append(inc(ybits - i).control(annotated = True), [i, *range(xbits+i, xbits + ybits)])
6
    return qc.to_gate()
7

8
@lru_cache
9
def adder(nbits, a):
10
    qc = QuantumCircuit(nbits, name = f"(+{a})")
11
    for i in range(nbits):
12
        if a & (1 << i):
13
            qc.append(inc(nbits - i), range(i, nbits))
14
    return qc.to_gate()

このように、足し算回路はインクリメント回路の組み合わせで実現ができます。
また後述するように、掛け算は足し算の組み合わせで実現できることなどから、インクリメント回路がすべての計算の基本となることが分かると思います。

べき乗剰余実現へ向けて

素因数分解で利用する oracle は、\( qy=a^{qx} \pmod{N} \) という計算でした。
一見すると計算が難しそうな \( a^{qx} \) というのは、\( qx \) を \( qx_0 \sim qx_3 \) にビット分解することで以下のように展開できます。

$$ a^{qx} = a^{qx_0} \cdot a^{2 qx_1} \cdot a^{4 qx_2} \cdot a^{8 qx_3} $$

\( 1 \) から開始して、\( qx_0=1 \) の場合は \( a \) を掛け算する、次に \( qx_1=1 \) の場合は \( a^2 \) を掛け算する、続けて \( qx_2=1 \) の場合は \( a^4 \) を掛け算し、最後に \( qx_3=1 \) の場合は \( a^8 \) を掛け算します。
この「もし1の場合は掛け算する」という「1の場合は」の部分は、制御ビットで実現できると予想できます。
さらに余りの世界の特徴として、4回掛け算してから \( \bmod{N} \) するのも、1回の掛け算ごとに \( \bmod{N} \) するのも計算結果としては変わりません。

$$ a \cdot b \cdot c \cdot d \pmod{N} = ((a \cdot b \pmod{N}) \cdot c \pmod{N}) \cdot d \pmod{N} $$

要するに \( x=x \cdot a \pmod{N} \) のような「剰余乗算」を繰り返し適用することで「べき乗剰余」が実現できます。

$$ \begin{align} x &= 1 \cr x &= x \cdot a^1 \pmod{N} \quad (qx_0=1 の場合) \cr x &= x \cdot a^2 \pmod{N} \quad (qx_1=1 の場合) \cr x &= x \cdot a^4 \pmod{N} \quad (qx_2=1 の場合) \cr x &= x \cdot a^8 \pmod{N} \quad (qx_3=1 の場合) \end{align} $$

次に行うのは掛け算の実現方法の検討です。
掛け算も2進数の筆算にすることで分かりやすくなります。 ここでは \( z = a \cdot x \) の掛け算を表現しています。

$$ \begin{array}{cccccccccc} & & & & & a_3 & a_2 & a_1 & a_0 \cr \times & & & & & x_3 & x_2 & x_1 & x_0 \cr \hline & z_7 & z_6 & z_5 & z_4 & z_3 & z_2 & z_1 & z_0 \end{array} $$

これは以下のように足し算の組み合わせに分解できます。

$$ \begin{array}{ccccccccccl} & & & & & a_3 & a_2 & a_1 & a_0 \cr \times & & & & & x_3 & x_2 & x_1 & x_0 \cr \hline {}+ & & & & & a_3 & a_2 & a_1 & a_0 & (x_0=1 の場合) \cr {}+ & & & & a_3 & a_2 & a_1 & a_0 & & (x_1=1 の場合) \cr {}+ & & & a_3 & a_2 & a_1 & a_0 & & & (x_2=1 の場合) \cr {}+ & & a_3 & a_2 & a_1 & a_0 & & & & (x_3=1 の場合) \cr \hline & z_7 & z_6 & z_5 & z_4 & z_3 & z_2 & z_1 & z_0 & \end{array} $$

これから分かるように、\( z=a \cdot x \) というのは、\( z = a \cdot x_0 + 2a \cdot x_1 + 4a \cdot x_2 + 8a \cdot x_3 \) のように足し算の繰り返しとなります。
この値を \( N \) で割った余りを計算することになりますが、これも余りの世界の特徴として、すべてを足してから \( \bmod{N} \) するのも、1回の足し算ごとに \( \bmod{N} \) するのも同じ結果となるため、以下のように繰り返し計算することで実現できます。

$$ \begin{align} z &= 0 \cr z &= z + a \pmod{N} \quad (x_0=1 の場合) \cr z &= z + 2a \pmod{N} \quad (x_1=1 の場合) \cr z &= z + 4a \pmod{N} \quad (x_2=1 の場合) \cr z &= z + 8a \pmod{N} \quad (x_3=1 の場合) \end{align} $$

このように、足し算して余りを求める「剰余加算」を実現することで、最終的には「べき乗剰余」も実現できるようになります。

量子回路の可逆性

べき乗剰余の実現方法は分かってきましたが、量子回路にはもう1つ重要なルールがあります。
可逆性という言葉で表現されている「逆操作で元に戻すことができる」「計算後に情報が損失されてはならない」というものです。

どういうことかというと、例えば掛け算で \( x \times y \to z \) と量子回路を実現する際は先に説明したような足し算等を組み合わせて行いますが、計算前後で \( z=0 \) の状態から \( z=x \cdot y \) に変わります。
\( x,y \) は計算前と同じなので情報が失われていないのは明白ですが、\( z \) は値が変わっています。
「逆操作で元に戻すことができる」というのは \( x,y \) を使って \( z \) を計算した逆の操作を行うことで元の \( z=0 \) に戻すことができることです。 また、「計算後に情報が損失されてはならない」というのは \( z=x \cdot y \) という計算後の値から計算前は \( z=0 \) であったことが求められる状態を指します。
この例の \( x \times y \to z \) は可逆となっています。

次に \( x \times y \to y \) というのを考えてみます。
\( x \) の値は計算前後で変わらないため問題ないとして、可逆であるためには計算後に変わってしまった \( y \) について計算前の値が何であったか求められ、元に戻せる必要があります。
一見すると、計算後の \( y \) を \( x \) で割れば計算前の値に戻るように思えますが、\( x=0 \) の場合に困ります。
\( x=0 \) だと計算後は \( y=0 \) となりますが、これだと計算前の \( y \) はどんな値でも取り得ることになり、元の値が分からないという情報が失われた状態になっています。
そのため \( x \times y \to y \) というのは可逆ではなくなり量子回路として実現できません。

今回実現したい剰余加算は、\( \bmod N \) という \( N \) で割った余りを使います。
ある値 \( x \) を \( N \) で割った余りを求めるには、\( x \pmod{N} \to x \) という量子計算ができればよいですが、これは可逆ではありません。
なぜかというと、例えば計算後の値が \( x=1 \) だったとしたら、計算前は \( x=1 \) あるいは \( x=N+1 \) という2つの可能性があり、そのどちらであるかを判定できる情報が失われているためです。
では余りの計算ができないかというとそうではありません。
\( x+y \pmod{N} \to x \) で計算前の 「\( x \) は必ず \( N \) 未満である」 と分かっている場合は可逆となります。
\( y \) は計算前後で変わらないため情報は保持されており、その \( y \) を使って計算後の \( x \) から引き算することで、\( x-y \lt 0 \) ならば \( x-y+N \to x \) で復元でき、\( x-y \ge 0 \) ならば \( x-y \to x \) となります。

上記の例の \( x+y \pmod N \to x \) はまさに剰余加算のことであり、量子回路で実現可能な計算です。

剰余加算の実現

実現したい剰余加算は \( x + a \pmod{N} \to x \) であり、\( x \) は量子レジスタで \( a \) は定数となります。
計算で使う量子ビットは、桁上がり用の \( \ket{c} \) を1ビット、\( N \) と同じビット数を持つ量子レジスタ \( \ket{x} \) を用意します。

$$ \begin{align} \ket{c}\ket{x} + a &\to \ket{c}\ket{x} \tag{1}\cr \ket{c}\ket{x} - N &\to \ket{c}\ket{x} \tag{2}\cr ctrl(\ket{c}) \otimes (\ket{x}+N) &\to \ket{x} \tag{3}\cr \ket{c}\ket{x} - a &\to \ket{c}\ket{x} \tag{4}\cr \ket{x} + a &\to \ket{x} \tag{5}\cr X(\ket{c}) &\to \ket{c} \tag{6} \end{align} $$

上記の (3) までは何となく意味は読み取れるかと思います。
(1)は \( x+a \) を計算し、(2)で一度 \( x+a-N \) となるように \( N \) を引き算します。
もし \( x + a \lt N \) であればマイナスになり \( \ket{c} \) が \( 1 \) となります。
これを利用して、(3)では \( x+a-N+N (x + a \lt N の場合) \) を行っています。

この時点で \( \ket{x} \) には正しい答えが入っていますが、問題なのは \( \ket{c} \) の状態です。
\( \ket{c} \) は場合によって \( 0 \) や \( 1 \) の可能性があり、これを残したままにすると剰余加算を繰り返して剰余乗算を実現する際に、続きの計算が期待したものになりません。
そのため、\( \ket{c} \) は確実に \( 0 \) に戻しておきたいのです。

この \( \ket{c} \to 0 \) にするための計算が (4) 以降なのです。
いろいろなパターンで検証すると分かりますが、(4)によって \( \ket{c} \) は必ず \( 1 \) となります。
(5)で一度引き算した \( a \) を再び足して、最後に \( \ket{c} \) をXゲートで \( 0 \) に戻します。
剰余加算の量子回路を図で示すと以下となります。

1
@lru_cache
2
def adder_mod(a, N):
3
    nbits = N.bit_length()
4
    qc = QuantumCircuit(nbits + 1, name = f"(+{a} mod {N})")
5
    qc.append(adder(nbits + 1, a), range(nbits + 1))
6
    qc.append(adder(nbits + 1, N).inverse(), range(nbits + 1))
7
    qc.append(adder(nbits, N).control(annotated = True), [nbits, *range(nbits)])
8
    qc.append(adder(nbits + 1, a).inverse(), range(nbits + 1))
9
    qc.append(adder(nbits, a), range(nbits))
10
    qc.x(nbits)
11
    return qc.to_gate()

剰余乗算の実装

剰余乗算 \( z=a \cdot x \pmod{N} \) は、以下のように剰余加算を繰り返すことで実現できました。

$$ \begin{align} z &= 0 \cr z &= z + a \pmod{N} \quad (x_0=1 の場合) \cr z &= z + 2a \pmod{N} \quad (x_1=1 の場合) \cr z &= z + 4a \pmod{N} \quad (x_2=1 の場合) \cr z &= z + 8a \pmod{N} \quad (x_3=1 の場合) \end{align} $$

これを量子レジスタ \( \ket{x} \ket{z} \) と桁上がり用の量子ビット \( \ket{c} \) を用意して同じような順番の回路を作成します（\( \ket{z} \) の初期値は \( \ket{0} \)）。

$$ \begin{align} ctrl(\ket{x_0}) \otimes (\ket{c}\ket{z} + a \pmod{N}) &\to \ket{c}\ket{z} \cr ctrl(\ket{x_1}) \otimes (\ket{c}\ket{z} + 2a \pmod{N}) &\to \ket{c}\ket{z} \cr ctrl(\ket{x_2}) \otimes (\ket{c}\ket{z} + 4a \pmod{N}) &\to \ket{c}\ket{z} \cr ctrl(\ket{x_3}) \otimes (\ket{c}\ket{z} + 8a \pmod{N}) &\to \ket{c}\ket{z} \end{align} $$

量子回路で示すと以下となります。

1
@lru_cache
2
def multi_mod(a, N):
3
    nbits = N.bit_length()
4
    qc = QuantumCircuit(nbits * 2 + 1, name = f"qy=(qx * {a}) mod {N}")
5
    for i in range(nbits):
6
        qc.append(adder_mod(a, N).control(annotated = True), [i, *range(nbits, nbits * 2 + 1)])
7
        a = (a * 2) % N
8
    return qc.to_gate()

べき乗剰余の実装

ついにべき乗剰余の実装までたどり着きました。

$$ \begin{align} x &= 1 \cr x &= x \cdot a^1 \pmod{N} \quad (qx_0=1 の場合) \cr x &= x \cdot a^2 \pmod{N} \quad (qx_1=1 の場合) \cr x &= x \cdot a^4 \pmod{N} \quad (qx_2=1 の場合) \cr x &= x \cdot a^8 \pmod{N} \quad (qx_3=1 の場合) \end{align} $$

上記のように繰り返して掛け算をしたいため \( \ket{x} \times a \pmod{N} \to \ket{x} \) を実現したいのですが、掛け算は可逆性の制限から \( \ket{x} \times a \pmod{N} \to \ket{z} \) のように別の量子レジスタに結果を設定する必要があります。
しかし、\( \ket{x} \) を上書きしながら更新していきたいので、以下のような工夫を行います。

$$ \begin{align} \ket{x} \times a \pmod{N} &\to \ket{z} \tag{1} \cr \ket{x} &\leftrightarrow \ket{z} \tag{2} \cr (\ket{x} \times a^{-1} \pmod{N} )^{\dagger} &\to \ket{z} \tag{3} \end{align} $$

(1)で計算した後に (2) で \( \ket{x} \) と \( \ket{z} \) を入れ替えています。
この時点で \( \ket{x} = \ket{z} \times a \pmod{N} \) であり \( \ket{z} = \ket{x} \times a^{-1} \pmod{N} \) になっています。
\( \ket{z} = \ket{x} \times a^{-1} \pmod{N} \) になっているということは、\( \ket{x} \times a^{-1} \pmod{N} \to \ket{z} \) という計算を逆回路にすることで \( \ket{z}=\ket{0} \) に戻ることを意味します。
この逆回路を適用しているのが (3) になります。
量子回路で示すと以下となります。

1
@lru_cache
2
def multi_mod_inplace(a, N):
3
    nbits = N.bit_length()
4
    qc = QuantumCircuit(nbits * 2 + 1, name = f"qx=(qx * {a}) mod {N}")
5
    qc.append(multi_mod(a, N), range(nbits * 2 + 1))
6
    qc.swap(range(nbits), range(nbits, nbits * 2))
7
    qc.append(multi_mod(pow(a, -1, N), N).inverse(), range(nbits * 2 + 1))
8
    return qc.to_gate()

\( \ket{x} \times a \pmod{N} \to \ket{x} \) が実現できたので、あとはこれを繰り返して \( qy=a^{qx} \pmod{N} \) を実現するだけです。

1
@lru_cache
2
def pow_mod(xbits, a, N):
3
    nbits = N.bit_length()
4
    qc = QuantumCircuit(xbits + nbits * 2 + 1, name = f"qy={a}^qx mod {N}")
5
    qc.x(xbits)
6
    for i in range(xbits):
7
        qc.append(multi_mod_inplace(a, N).control(annotated = True), [i, *range(xbits, xbits + nbits * 2 + 1)])
8
        a = (a * a) % N
9
    return qc.to_gate()

\( qy_0 \) の最初にXゲートがあるのは、\( \ket{qy}=\ket{1} \) に初期化するためです。

終わりに

Xゲートという \( \ket{0} \) と \( \ket{1} \) を反転させる単純な量子ゲートと、制御ビットの組み合わせだけで \( qy=a^{qx} \pmod{N} \) という計算を実現するのがどれほど難しいか伝わったでしょうか？
これを実現するためには、インクリメント→足し算→剰余加算→剰余乗算→べき乗剰余というように、計算の要素を分解して組み立てるスキルが必要となります。
これはプログラミングの本質的なスキルに通じるものであり、このような考えができるかどうかがプログラミングスキルの高さを示すことにもつながります。

今回紹介した手法以外にも効率の良い量子回路の実装方法はあるので、興味のある方は他の実装方法にチャレンジしてみてはいかがでしょうか。

最後に、べき乗剰余を利用したショアのアルゴリズムのプログラム（Qiskitを利用した量子回路の構築部分）を掲載いたします。

1
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
2
from qiskit.circuit.library import QFT
3

4
a = 2
5
N = 15
6
nbits = N.bit_length()
7

8
qx = QuantumRegister(nbits * 2, "qx")
9
qy = QuantumRegister(nbits, "qy")
10
qz = QuantumRegister(nbits + 1, "qz")
11
cx = ClassicalRegister(len(qx), "cx")
12
qc = QuantumCircuit(qx, qy, qz, cx)
13
qc.h(qx)
14
qc.append(pow_mod(len(qx), a, N), [*qx, *qy, *qz])
15
qc.append(QFT(len(qx)).inverse(), qx)
16
qc.measure(qx, cx)