第1章: アルゴリズムからパルスまでのスタック

🌐 JP | 🇬🇧 EN | Last sync: 2026-08-13

基礎数理道場 > 量子ソフトウェアスタック入門 > 第1章

アルゴリズムは数学的な対象です。機械が受け取るのはアナログの波形です。その間にあるものは何ひとつ自明ではなく、その間にあるものすべてがソフトウェアです。本章はそのソフトウェアの見取り図です。7つの層、それぞれが何のためにあるのか、そのうちどれがハードウェアについて何も知らずに書けるのか、そして — この主題を扱えるものにしている点ですが — ある層が仕事を果たしたことをどう検査するのか。

最後の点は本コース全体の土台なので、先に述べておく価値があります。ここでのコンパイラパスは回路から回路への関数であり、2つの回路が同じことを意味しているときに正しいとされます。量子回路にとって「同じことを意味する」には厳密な定義があります。大域位相を除いて同じユニタリであること、です。この定義は検査できます。小さな回路の行列を作り、書き換えた回路の行列を作り、大域位相を除去して差を取ればよいのです。第2章と第3章のすべての書き換え規則、第4章のすべてのパルス、第5章のすべての緩和式には、この種の検査が付いています。ハードウェアに触れずに本コースを誠実に書けるのはこの検査があるからであり、本コースから持ち帰る習慣として最も役に立つのもこれです。

本章は3つのものを作ります。1.2節が回路の中間表現（IR）を固定します。以降のすべての章が操作するデータ構造であり、層の境界を層の境界にしているのはゲート集合の選択であることを論じます。1.3節がコンパイルを意味を保つ書き換えとして定義し、等価性チェッカを作ります。1.4節が商用SDKがこれらの層を何と呼んでいるかの地図です。ただし一般的な言い方に限ります。フレームワーク間でも版の間でも言葉は違い、層は違わないからです。

学習目標

本章を修了すると、以下のことができるようになります。

• アルゴリズムと制御波形の間にある7つの層を挙げ、各層が何を消費して何を生産するかを述べ、各層がハードウェアのどの情報を必要とするかを正確に言える
• 中間表現がなぜ存在するのか — $N \times M$ の議論 — を説明し、層の境界を固定するのがデータ構造ではなくゲート集合の選択であることを述べられる
• 本コースの回路IR — ゲートタプル、run_circuit、circuit_depth、gate_counts — を入門コースのビッグエンディアン規約のもとで実装できる
• コンパイラパスに要求しうる3つの正しさの関係 — 大域位相を除いて厳密、量子ビットの置換を除いて厳密、$\varepsilon$ の範囲で近似 — を述べ、どの層がどれを使うかを言える
• 大域位相除去付きのユニタリ等価性チェッカを作り、Hilbert-Schmidt重なりが最適な位相を与える理由を説明できる
• 制御ブロックの内側では大域位相を捨ててはいけない理由を説明し、そこから生じるバグを数値的に示せる
• 費用の指標 — 総ゲート数、2量子ビットゲート数、深さ、実時間、誤り予算 — を区別し、それらが一致しないことを示して、目的関数を指定しない「最適化」が無意味であることを述べられる
• 任意の量子SDKのドキュメントの目次を読み、その各項目がどの層の話かを言える

引き継ぐ規約

3つの規約を量子コンピューティング入門からそのまま引き継ぎます。

ビッグエンディアンの量子ビット順序。 量子ビット0がケットの最左の記号であり、振幅の添字の最上位ビットです。したがって3量子ビットでは $\lvert 000 \rangle, \lvert 001 \rangle, \ldots, \lvert 111 \rangle$ が添字 $0$ から $7$ を占めます。SDKの文献の多くは逆の規約を用いており、この食い違いはこの分野で結果が黙って狂う最大の原因です。そしてこれは本章の主題のよい例でもあります。量子ビット順序は物理の性質ではなく層と層の界面の性質であり、事実ではなく規約であるからこそ、実際のデバッグ時間を食うのです。

シミュレータ。 すべてのコード例は同コース第2章の99行の状態ベクトルシミュレータ上で走ります。Code Example 1 は本コースが使う関数のみを逐語で再掲します。

簡約単位系。 ハミルトニアンが現れるところでは $\hbar = 1$ とします。本章では現れず、第4章では至るところに現れます。

新しい規約が1つあり、それが本コースの技術的な契約です。1.2節の回路IRです。第2章から第5章までがこれを一切変えずに使うので、ある章のために書いたパスを別の章の回路に適用できます。

1.1 なぜ層があるのか

埋めるべき距離

10量子ビットのGroverのアルゴリズムを取ります。数学としては2つの鏡映の積を $\lfloor \pi/(4\theta) \rfloor$ 回反復したものです。実験室で起きることとしては、数万個の整形されたマイクロ波あるいはレーザーパルスであり、それぞれに振幅・持続時間・周波数・位相・包絡線があり、ナノ秒精度のスケジュールで届けられ、最後に電圧をビットに弁別しなければならないアナログ測定が来ます。

後者を手で書く人はいませんし、書きたい人もいません。それは1台の機械の1日に固有のものです。$\pi$ 回転を起こすパルス振幅はドリフトし、毎日測り直されるからです。前者を直接ハードウェアに書く人もいません。前者はどの量子ビット同士が相互作用できるかについて何も言っていないからです。逃げ道は古典計算と同じで、表現を間に挟むことです。各表現は上の層が気にすべきでない詳細を捨て、下の層が要求する詳細を加えます。以下の7層が本コースの使う層です。名前はフレームワークによって違い、分業は違いません。

消費/生産の列は、第2層から第5層までがすべて同じ言語を話していることを示しています。ゲート列が入り、ゲート列が出る。これが中間表現の要点そのものです。すべてのパスが同じ型をもつので、パスは合成でき、並べ替えられ、1つずつテストできます。型が変わるのは第6層で、だからこそパルスレベルのプログラミングは回路レベルのものとまったく違う感触がします。知らなければならないことの列が、層の境界がどこにあるべきかの判定基準です。ハードウェアの情報を必要としないパスは一度書けば永久に使えます。結合グラフを必要とするパスは機械ごとに走らせ直します。今日の較正値を必要とするパスは毎朝走らせ直します。パスを入力の生鮮度で並べると、ちょうど上の層構造になります。そしてこれが、量子コンパイラのうち最適化器が最も可搬でパルスコンパイラが最も可搬でないという、一見不思議な事実の説明にもなっています。

$N \times M$ の議論

IRの経済的な理由は、古典のコンパイラIRを生んだのと同じものです。$N$ 個のアルゴリズムと $M$ 台の機械があって中間表現がなければ、$N \times M$ 個の翻訳を書きます。1つあれば、$N$ 個のフロントエンドと $M$ 個のバックエンドです。第2の、もっと目立たない利点があり、量子計算ではこちらの方が重要です。

IRは測る場所です。回路がデータ構造になった瞬間、ゲート数を数え、深さを計算し、2量子ビットゲート数を数え、$T$ ゲート数を数え、誤りを見積ることができます。何かを走らせるのは高価であり、量子アルゴリズムについての興味ある問いはたいてい「機械はどれだけ大きくなければならないか」なので、この分野の実際の仕事のほとんどはIRの層で起きています。

量子のコンパイルは何が違うか

3つあります。本コースが古典コンパイラ講義の翻訳ではない理由を説明するので、名前を付けておく価値があります。

正しさは小さいサイズでは決定可能で、大きいサイズでは高価である。 2つの古典プログラムが等価であるとは同じ関数を計算することであり、これは一般に決定不能です。2つの量子回路が等価であるとは位相を除いて同じユニタリを実装することであり、これは2つの $2^n \times 2^n$ 行列についての主張、すなわち有限の計算です。これは贈り物であり、1.3節はそれを使い切ります。落とし穴は $2^n$ です。全数検査は小さい $n$ のための道具であり、20量子ビットあたりより上は構造的あるいは統計的に検証しなければなりません。

最適化は任意選択ではない。 2倍長い古典プログラムは2倍遅く走ります。2倍長い量子回路は走らないかもしれません。出力が情報を担う限界を超えて誤りが積み上がるからです。Code Example 6 がこれを定量化します。したがって量子計算においてコンパイルは便利のための層ではなく、計算が実行可能かどうかを決める一部です。2量子ビットゲート数の2倍は、必要な機械の大きさの2倍です。

最下層は物理実験である。 第6層は表を引いてゲートをパルスに翻訳するのではありません。その機械上で最近測られた母数を使って翻訳し、その較正実験自体がスタックの中のソフトウェアです。第4章がそのループを実装します。自分の機械を測ってから自分を書き換えるソフトウェアは珍しく、この分野の本当に特徴的な点の1つです。

層が漏れるところ

上の層構造は工学上の約束事であって定理ではなく、面白い研究はその継ぎ目にあります。3つの例を挙げます。いずれも後の章で戻ってきます。

• ハードウェアを意識した最適化。 第3層は機械について何も知らないはずですが、結合グラフを知っている最適化器は、後のルーティング費用を減らす書き換えを選べます。第3層と第4層を統合すると、よい回路と、はるかに難しい問題が得られます。
• パルスレベルのコンパイル。 隣接する2つのゲートを別々にコンパイルすれば2つのパルスになりますが、一緒にコンパイルすれば1つの短いパルスになりうります。第5層と第6層を潰すと本物の忠実度が買え、可搬性は完全に失われます。
• 誤り訂正。 誤り耐性のスタックは第3層と第4層の間に丸ごと1つ層を挿し込みます。そこでは1つの論理ゲートが数千の物理操作になり、ハードウェアが何をネイティブに実装していようとネイティブゲート集合はClifford$+T$になります。第5章がそのリソース計算を扱います。

この分野の妥当な読み方は、層は別々の人が別々の問題に取り組めるようにするために存在しており、数年ごとに誰かが層の1つを無視して2倍を稼いでいる、というものです。

1.2 回路IRの設計

回路はコードではなくデータである

本コースが使う表現は意図的に最小限です。

$$ \text{回路} = \left[ g_1, g_2, \ldots, g_L \right], \qquad g_i = (\text{名前}, \text{引数} \ldots) $$

回路はタプルのPythonリストです。ゲートは、第1要素がゲート名の文字列、残りが（あれば）角度と、作用する量子ビットからなるタプルです。全体は次のとおりです。

ゲートは左から右へ作用させるので、回路の行列はゲート行列の逆順の積になります。

$$ U(\left[g_1, \ldots, g_L\right]) = U_{g_L} \cdots U_{g_2} U_{g_1} $$

これは本コースのすべての図が使う順序であり、行列を書く順序とは逆です。逆にすると意図した回路の逆順になり、多くの回路ではそれがほとんど正しいので気づきにくいのです。

3つの性質は擁護しておく価値があります。どれも本物のIRも下さなければならない設計判断だからです。

オブジェクトではなくデータである。 ゲートはタプルなので、== で比較でき、辞書のキーにでき、印字でき、パターン照合できます。するとパスは状態をもたないリストからリストへの普通の関数になり、2つのパスは関数合成で合成されます。本物のIRはもっと豊かな構造 — 可換なゲートに人工的な順序を付けないための有向非巡回グラフ — を使います。パスがリスト上で隣接するのではなくある量子ビット上で隣接するゲートを探したいとき、リスト表現が何を失っているかは第2章で論じます。

量子ビットは整数であり、整数はワイヤである。 それらはレジスタ内の位置であって、物理量子ビットの名前ではありません。プログラムの量子ビットから機械の量子ビットへの写像がレイアウトで、第4層で選ばれ、別の対象です。第3章がそれを明示します。両者を混同するのがトランスパイラの出力を読むときの混乱の最大の原因です。ルーティングの後、ワイヤ $k$ はもうプログラマが $k$ と呼んだ量子ビットを保持していないからです。

IRが意図的に省いているもの

4つあります。どれも本番のIRは含んでおり、どれも本コースの後の章で戻します。

一般的な論点は、IRは何を表現するかと同じくらい、何を表現するのを拒むかによって定義されるということです。何でも表現できる表現はプログラミング言語であり、それは最適化できません。

ゲート集合が層の境界である理由

本節の主張はこうです。1.1節の表の層の境界は、ちょうどゲート集合が変わる点です。本コースには3つのゲート集合が現れ、その間の各境界は費用を伴う合成問題です。

集合を変えるたびに書き換え問題が生じ、それぞれに特徴的な費用があります。記述用集合を $\lbrace R_z, R_y, \mathrm{CX}\rbrace$ の形のネイティブ集合に変換する費用はゲートあたり定数倍です。Code Example 7 は小さな回路でおよそ $1.2$ と測ります。ほとんどのゲートが1個か2個のネイティブに展開されるからです。Clifford$+T$への変換は種類が違います。一般の回転 $R_z(\theta)$ は精度 $\varepsilon$ に対して長さが $\log(1/\varepsilon)$ で伸びるClifford$+T$ゲートの列を要するので、費用は定数倍ではなく要求する精度の関数です。この非対称性は誤り耐性のリソース見積りについて最も重要な事実であり、第5章がそれを使います。

系として、実際に気にする理由が出てきます。SDKのドキュメントがあるバックエンドの「基底ゲート集合」を述べているとき、それはその機械についてこの境界の1つがどこにあるかを教えており、したがってそのコンパイラがどの合成問題を解かなければならなかったかを教えています。そのインタフェースについて、これほど重要なことは他にありません。

深さと、それが仮定していること

契約の最後の部分は費用モデルです。2つの関数が回路を要約します。

$$ \mathrm{depth}(C) = \text{ゲートを貪欲に詰めたときの層数}, \qquad \mathrm{counts}(C) = \lbrace \text{名前} \mapsto \text{個数} \rbrace $$

貪欲な層分けは自明なアルゴリズムです。各量子ビットについて最初に空く層の添字を保ち、各ゲートを、それが触れるすべての量子ビットが空いている最も早い層に置き、それらの量子ビットを進める。結果は量子ビットを共有するゲートの最長の鎖の長さ、すなわち回路のクリティカルパスです。

深さは誤りが待ち時間から来るときに正しい指標です。早く終わる回路はデコヒーレンスする時間が短いからです。そして誤りがゲートから来るときには誤った指標です。10個のゲートからなる層が1個の層と同じに数えられるからです。さらに、単純に偽である仮定に乗っています。すべてのゲートが同じ時間を要するという仮定です。実際のハードウェアでは2量子ビットゲートは1量子ビットゲートより典型的に1桁遅いので、回路の実時間は重み付きの深さです。Code Example 6 が両方を計算して、どれだけ離れるかを示します。第4章がその重みをパルスの持続時間で置き換えます。

1.3 コンパイルとは意味を保つ書き換えである

3つの正しさの関係

コンパイラパスは回路 $C$ を回路 $C' = P(C)$ に写す関数 $P$ です。$C'$ が $C$ の意味していたことを意味するとき正しいとされます。この主張には3つの版が使われており、本コースのすべてのパスはどれを満たすかを宣言します。

大域位相を除いて厳密。 最も強い関係であり、既定です。

$$ U(C') = e^{i\varphi}\, U(C) \quad \text{ある } \varphi \in \mathbb{R} \text{ について} $$

位相を割り出すのは、それが観測できないからです。測定確率は $\lvert \langle x \rvert \psi \rangle \rvert^2$ に、期待値は $\langle \psi \rvert A \lvert \psi \rangle$ に依存し、どちらも全体の $e^{i\varphi}$ を見ません。最適化パスと基底変換はこの関係を満たし、以下のチェッカが検査するのもこれです。

量子ビットの置換を除いて厳密。 SWAPを挿入して戻さないルータは、プログラムの量子ビットが別のワイヤに載った状態で終わります。

$$ U(C') = e^{i\varphi}\, \Pi\, U(C)\, \Pi'^{\dagger} $$

ここで $\Pi, \Pi'$ は置換行列であり、コンパイラはこれを記録して測定結果を解釈する層に渡さなければなりません。記録を忘れるのは典型的なバグで、その症状は、まったくもっともらしい出力分布のビットが入れ替わっていることです。第3章は置換を明示的に持ち回ります。Code Example 7 のルータはSWAPを戻すことでこの問題を回避します。正しく、そして無駄です。

$\varepsilon$ の範囲で近似。 離散ゲート集合への合成は厳密にはできないので、保証は距離になります。

$$ \min_{\varphi} \left\lVert U(C') - e^{i\varphi} U(C) \right\rVert \le \varepsilon $$

作用素ノルムで、回路が測定やノイズを含むならダイヤモンドノルムで測ります。これが第2章のClifford$+T$合成と第4章の有限時間パルスが必要とする関係です。構造上重要な事実は $\varepsilon$ が足し算されることです。それぞれ $\varepsilon$ で合成された $L$ 個のゲートは、ノルムの劣加法性により $L\varepsilon$ の精度の回路を与えるので、長い回路が要求するゲートごとの精度はその長さとともに厳しくなります。

大域位相を正しく除去する

大きさ $d = 2^n$ の2つの行列 $U$ と $V$ が与えられたとき、位相ができる範囲で両者を最も近づける位相と、残差の距離が欲しいわけです。Hilbert-Schmidt内積を書きます。

$$ \langle V, U \rangle = \mathrm{tr}\left( V^{\dagger} U \right) $$

最適な位相は $e^{i\varphi} = \langle V, U \rangle / \lvert \langle V, U \rangle \rvert$ であり、これは標準的な最小二乗の答えです。$\lVert U - e^{i\varphi}V\rVert_F^2$ を $\varphi$ について最小化することは $e^{i\varphi}$ を重なりに揃えることを意味します。ユニタリな $U, V$ に対してCauchy-Schwarzは

$$ \left\lvert \mathrm{tr}\left(V^{\dagger}U\right) \right\rvert \le \sqrt{\mathrm{tr}\left(U^{\dagger}U\right)\,\mathrm{tr}\left(V^{\dagger}V\right)} = d $$

を与え、等号が成り立つのはちょうど $U = e^{i\varphi} V$ のときです。したがって正規化した重なり $\lvert \mathrm{tr}(V^{\dagger}U)\rvert / d$ 自身が完全な判定になります。等価な回路では $1$ に等しく、そうでなければ真に小さいのです。これは退化する場合も片付けます。重なりがゼロなら位相は定義されませんが、重なりがゼロの2つの回路はこれ以上ないほど非等価なので、どんな位相でも構いません。

大域位相はいつでも大域的とは限らないという規則

位相を割り出す便宜がバグに変わるのはここだけなので、独立した規則にしておきます。

$$ R_z(\pi) = \begin{pmatrix} e^{-i\pi/2} & 0 \cr 0 & e^{i\pi/2} \end{pmatrix} = -i \begin{pmatrix} 1 & 0 \cr 0 & -1 \end{pmatrix} = -i\, Z $$

$R_z(\pi)$ と $Z$ は等価で、両者を区別する実験はありません。ところが両方を別の量子ビットで制御してみます。制御版は制御が $\lvert 1 \rangle$ のときにだけ因子 $-i$ を掛けるので、条件付きで掛かる位相は相対位相であり、これは完全に観測可能です。

$$ \mathrm{C}\text{-}R_z(\pi) = \mathrm{diag}\left(1, 1, -i, i\right) \ne \mathrm{diag}\left(1,1,1,-1\right) = \mathrm{CZ} $$

Code Example 4 はこの差を $0.765$ と測ります。丸め誤差ではなく、別のゲートです。したがって規則です。どの量子コンパイラもどこかで強制しており、どの実装者も一度は間違えます。

大域位相を除いてのみ正しい書き換えを、制御ゲートの中身に適用してはいけません。位相をデータとして持ち回るか、制御版を最初から合成するかのどちらかです。

演習3が上の例を修復する補正を求め、捨てた位相が制御量子ビット上の正真正銘のゲートとして再登場することを見ます。

本コースの検証の方法論

以上はすべて1つの再利用可能なルーチンになり、Code Example 4 以降はすべてのパスがそれを通して呼ばれます。

1. $2^n$ 個の基底状態それぞれで回路を走らせて $U(C)$ を作る。費用はシミュレーション $2^n$ 回。
2. 同じ方法で $U(C')$ を作る。
3. トレースの式で大域位相を除去し、$\max_{jk} \lvert U_{jk} - e^{i\varphi} V_{jk} \rvert$ を取る。
4. 許容値 $10^{-10}$ と比べる。これは浮動小数点の床 $10^{-16}$ の8桁上、誤りと呼ぶに値するどんな誤りよりも8桁下です。

誠実な限界は手順1の $2^n$ です。$n = 8$ ならシミュレーション256回でミリ秒です。$n = 20$ なら百万振幅の状態を百万回シミュレートすることになり、手が届きません。全数検査の範囲より上では2つの選択肢があり、本コースは両方を使います。ランダムな入力状態: 両方の回路を1つのHaarランダム状態に走らせて重なりを比べます。ユニタリを変えるパスは確率1でランダム状態の出力を変えるので、ランダム状態1個は驚くほど強いテストです。演習4は、試したすべての回路で全数検査と一致することを、$2^n$ 回ではなく1回のシミュレーションで測ります。構造的な議論: 書き換え規則を、それが触れる2つか3つの量子ビット上で一度証明し、そのうえでパスがその前提が成り立つ場所でしか適用しないことを論じます。本物のコンパイラは後者に頼り前者でテストしており、だからこの分野のコンパイラのバグはたいてい誤った規則ではなく、誤った場所で適用された規則です。Code Example 5 はまさにその種のバグを作り、捕まえます。

「よい」とは何か

コードの前にもう1つ語彙が必要です。意味を保つパスはそれだけでは改善ではありません。改善には目的関数が要り、使われているものが少なくとも5つあります。

これらは一致せず、Code Example 6 がどれだけ違うかを示します。同じユニタリを2通りにコンパイルしたものは、2量子ビットゲートが速い機械では時間で $60\%$ 安く、遅い機械では $1.8\%$ 安いのです。したがって最適化器はいつでも何かのための最適化器であり、何のためかを指定させないコンパイラは誤ったものを最適化します。

1.4 SDKの地図

本コースはSDKを教えません。それでもSDKを使う人のために書かれているので、本節では層をフレームワークの部品に対応させます。どのフレームワークの今後数版にも耐えるだけ一般的な言い方でです。

層の対応

どの回路モデルのフレームワークも、名前は違っても同じ部品をもっています。

この表について、表そのものより価値のある観察が3つあります。

機械の記述こそが面白い対象である。 第4層・第5層・第6層はすべて同じもので母数化されています。結合グラフ、ネイティブゲート集合、較正データを含むターゲットの記述です。フレームワークがハードウェアをまたいで可搬なのは、その3つの情報が1つのオブジェクトに括り出されているからであり、そのオブジェクトの定義を読むのが、そのフレームワークが実際に何を狙えるのかを知る最短経路です。

フレームワークの違いは、ほとんどがどの層を前面に出すかの違いである。 回路を中心に作られたフレームワークは第2層を前面に出し、勾配をサービスとして扱います。微分可能プログラムを中心に作られたものはアルゴリズムの層を前面に出し、回路を実装の詳細として扱い、パラメータシフト則と逆伝播を軸に構成されます。これは名前の違いではなく本当に別の設計です。パルスを中心に作られたものは第6層を前面に出します。3つともすべての層を含んでいます。

最適化レベルはパス順序のメニューである。 ほとんどの利用者が触る唯一のつまみ — 「optimization level」と名付けられた整数 — はパスの列を選びます。だからこそ第2章はそのつまみではなく個々の書き換え規則に時間を使います。

SDKが与えてくれて本コースが与えないもの

これについて公正であることは本節の要点の一部です。本番のフレームワークはハードウェアへのアクセスとキューを提供し、ここで書くものよりはるかに優れた、実機に対して開発されたパスを提供し、較正が変わるたびに更新される機械の記述を対応バックエンドごとに提供し、そしてあなたがこれから踏む地雷を既に踏んだ人々の共同体を提供します。どれも代替できませんし、どれも本コースの目的ではありません。

本コースが提供するのはもう半分、パスが何をしたのかを知っていることです。その知識があれば、トランスパイル後の回路を読んでなぜ膨れたのかが見え、最適化レベルが効いているかを判断でき、結果が誤ったときにルーティングの問題か合成の問題かを見分けられ、コードを1行も書く前にアルゴリズムの費用を見積れます。具体的には、この5章を終えたら、任意の量子フレームワークのドキュメントを開き、その目次を読んで、各項目を上の表の行に割り当てられるはずです。それが成果物です。

1.5 層を作る

本章の残りは実装です。Code Example 1 と 2 がIRを作り、3がそれを動かし、4がチェッカを作り、5が本物のパスを書いて本物のバグを捕まえ、6が費用モデルを付け、7が第2章から第5章が埋めていくパイプラインを組み立てます。

Code Example 1: シミュレータの再掲

これは量子コンピューティング入門第2章の状態ベクトルシミュレータであり、本コースが使う関数のみを逐語で再掲したものです。qcsim.py として保存してください。以下のIRモジュールは from qcsim import * から始まります。ここに新しいものは何もなく、変更も一切ありません。入門コースのファイルをすでにお持ちなら、そちらを使ってください。

"""Minimal state-vector simulator (big-endian: qubit 0 = leftmost = most significant).

Save this file as qcsim.py; every later example does `from qcsim import *`.
"""
import numpy as np

# ---- 1量子ビットゲート --------------------------------------------------
I2 = np.eye(2, dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
S = np.array([[1, 0], [0, 1j]], dtype=complex)
T = np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]], dtype=complex)


def rx(theta):
    c, s = np.cos(theta / 2), np.sin(theta / 2)
    return np.array([[c, -1j * s], [-1j * s, c]], dtype=complex)


def ry(theta):
    c, s = np.cos(theta / 2), np.sin(theta / 2)
    return np.array([[c, -s], [s, c]], dtype=complex)


def rz(theta):
    e = np.exp(-1j * theta / 2)
    return np.array([[e, 0], [0, np.conj(e)]], dtype=complex)


# ---- 状態 ---------------------------------------------------------------
def ket(bits: str) -> np.ndarray:
    """'01' -> 4次元の基底状態 |01>（ビッグエンディアン）"""
    n = len(bits)
    psi = np.zeros(2 ** n, dtype=complex)
    psi[int(bits, 2)] = 1.0
    return psi


def apply_gate(state, U, targets, n):
    """n量子ビット状態の targets に 2^k x 2^k ユニタリ U を作用させる"""
    k = len(targets)
    psi = state.reshape([2] * n)          # 1. n添字テンソルとして見る
    psi = np.moveaxis(psi, targets, range(k))   # 2. 標的軸を先頭へ
    rest = psi.shape[k:]
    psi = psi.reshape(2 ** k, -1)         # 3. 平坦化して行列積
    psi = U @ psi
    psi = psi.reshape(list((2,) * k) + list(rest))
    psi = np.moveaxis(psi, range(k), targets)   # 4. 軸を元に戻す
    return psi.reshape(-1)


CNOT4 = np.array([[1, 0, 0, 0],
                  [0, 1, 0, 0],
                  [0, 0, 0, 1],
                  [0, 0, 1, 0]], dtype=complex)


def cnot(state, control, target, n):
    """任意の量子ビット対・任意の向きのCNOT"""
    return apply_gate(state, CNOT4, [control, target], n)


def probs(state):
    """Born則による全 2^n 通りの確率"""
    return np.abs(state) ** 2


def sample(state, shots, seed=None):
    """測定のシミュレーション: {ビット列: 回数}"""
    n = int(np.log2(state.size))
    rng = np.random.default_rng(seed)
    idx = rng.choice(state.size, size=shots, p=probs(state))
    out = {}
    for i in idx:
        b = format(i, f'0{n}b')
        out[b] = out.get(b, 0) + 1
    return dict(sorted(out.items()))

実際の仕事をしている関数は apply_gate だけであり、本コース全体がこれに乗っているので一文だけ触れておきます。状態を $n$ 添字テンソルに変形し、標的の軸を先頭に移し、$2^k \times 2^k$ 行列を掛け、軸を戻します。したがって量子ビットの任意の部分集合に対するゲートは $O(2^k 2^n)$ 演算で済み、追加のメモリは不要です。だから $n = 8$ は瞬時で、$n = 20$ でも可能なのです。

Code Example 2: 回路IR

これが本コースの契約です。qir.py として保存してください。以下のすべてのコード例、そして以降のすべての章が from qir import * から始まります。

"""Chapter 1, Example 2: 本コースの回路IRです。

回路はゲートタプルのリストです。ゲート名は文字列、量子ビットは整数
（ビッグエンディアン、量子ビット0が最左）です。このファイルを qir.py として
保存してください。以下のすべてのコード例は `from qir import *` から始まり、
以降のすべての章がこれを再掲します。

    ("h", q)   ("x", q)   ("z", q)   ("s", q)   ("t", q)
    ("rx", theta, q)      ("ry", theta, q)      ("rz", theta, q)
    ("cx", control, target)                     ("cz", q1, q2)
"""
import numpy as np
from qcsim import *

CZ4 = np.diag([1.0, 1.0, 1.0, -1.0]).astype(complex)

FIXED_1Q = {"h": H, "x": X, "z": Z, "s": S, "t": T}
ROT_1Q = {"rx": rx, "ry": ry, "rz": rz}
TWO_Q = ("cx", "cz")


def gate_qubits(g):
    """ゲートタプル1個が触れる量子ビットを、書かれた順に返します"""
    if g[0] in ROT_1Q:
        return (g[2],)
    if g[0] in TWO_Q:
        return (g[1], g[2])
    if g[0] in FIXED_1Q:
        return (g[1],)
    raise ValueError(f"unknown gate name {g[0]!r}")


def apply_ir_gate(state, g, n):
    """ゲートタプル1個を n量子ビット状態ベクトルに作用させます"""
    name = g[0]
    if name in FIXED_1Q:
        return apply_gate(state, FIXED_1Q[name], [g[1]], n)
    if name in ROT_1Q:
        return apply_gate(state, ROT_1Q[name](g[1]), [g[2]], n)
    if name == "cx":
        return cnot(state, g[1], g[2], n)
    if name == "cz":
        return apply_gate(state, CZ4, [g[1], g[2]], n)
    raise ValueError(f"unknown gate name {name!r}")


def run_circuit(circ, n, psi0=None):
    """ゲートタプルのリストを状態ベクトルシミュレータで実行し、最終状態を返します。

    psi0 の既定値は |00...0> です。ゲートは左から右へ作用させるので、回路の
    行列はゲート行列の逆順の積になります。
    """
    state = ket("0" * n) if psi0 is None else np.asarray(psi0, dtype=complex)
    for g in circ:
        state = apply_ir_gate(state, g, n)
    return state


def circuit_depth(circ, n):
    """量子ビットの排他性による貪欲な層分け: 回路が必要とする層数です。

    すべてのゲートが1単位時間を要すると仮定しています。これはハードウェア上では
    誤りであり、第4章で修正します。
    """
    ready = [0] * n              # 各量子ビットが最初に空く層
    for g in circ:
        qs = gate_qubits(g)
        layer = max(ready[q] for q in qs)
        for q in qs:
            ready[q] = layer + 1
    return max(ready) if n else 0


def gate_counts(circ):
    """ゲート名 -> 個数。加えてキー "2q" に2量子ビットゲートの総数を入れます"""
    counts = {}
    for g in circ:
        counts[g[0]] = counts.get(g[0], 0) + 1
    counts["2q"] = sum(counts.get(name, 0) for name in TWO_Q)
    return counts

いかに小さいかに注目してください。公開関数4つ、ディスパッチ表1つ、クラスは0。本コースの残りの内容はすべて、このデータ構造の上のパスです。あとで効いてくる細部が2つあります。run_circuit は初期状態を省略可能な引数として受け取ります。これが、各基底ベクトルで順に走らせて回路の行列を取り出すことを可能にします。Code Example 4 がまさにそれをします。そして gate_counts は名前ごとの個数と集計値 "2q" の両方を報告します。2量子ビットの個数が誤りを予測する数であり、呼び出し側が毎回計算し直すべきものではないからです。

Code Example 3: IRを通したBellとGHZ

新しい表現でまずすべきことは、答えの分かっているものを走らせることです。

"""Chapter 1, Example 3: IRを通したBell状態とGHZ状態。
Code Example 2 の続き（同一セッション）。"""
import numpy as np
from qir import *

# ---- どの量子コンピューティング入門も最初に書く2行の回路 -----------------
bell = [("h", 0), ("cx", 0, 1)]
psi = run_circuit(bell, 2)

print("Bell回路:", bell)
print(f"  振幅            : {np.round(psi.real, 6)}")
print(f"  深さ            : {circuit_depth(bell, 2)}")
print(f"  ゲート数        : {gate_counts(bell)}")
print(f"  ノルム = 1      : "
      f"{np.isclose(probs(psi).sum(), 1.0, rtol=1e-12, atol=0.0)}")
print(f"  振幅の誤差      : 1/sqrt(2) に対して "
      f"{abs(psi[0] - 2 ** -0.5):.2e}, {abs(psi[3] - 2 ** -0.5):.2e}")


# ---- 同じ状態を n量子ビット上で2通りに作ります --------------------------
def ghz_chain(n):
    """CNOTの鎖によるGHZ: 2量子ビットゲートは n - 1 個、深さは n です"""
    return [("h", 0)] + [("cx", q, q + 1) for q in range(n - 1)]


def ghz_log(n):
    """倍化によるGHZ: 2量子ビットゲートは同じ n - 1 個、深さは約 log2(n) です"""
    circ, span = [("h", 0)], 1
    while span < n:
        for q in range(min(span, n - span)):
            circ.append(("cx", q, q + span))
        span *= 2
    return circ


print("\n同一の n量子ビットGHZ状態に対する2つの回路")
head = (f"{'n':>3}{'chain depth':>13}{'log depth':>11}{'cx (chain)':>12}"
        f"{'cx (log)':>10}{'max |dpsi|':>13}")
print(head)
print("-" * len(head))
for n in range(2, 9):
    a, b = ghz_chain(n), ghz_log(n)
    pa, pb = run_circuit(a, n), run_circuit(b, n)
    err = np.max(np.abs(pa - pb))
    print(f"{n:>3}{circuit_depth(a, n):>13}{circuit_depth(b, n):>11}"
          f"{gate_counts(a)['2q']:>12}{gate_counts(b)['2q']:>10}{err:>13.2e}")

# ---- 4量子ビットGHZ状態を測定すると何が返るか --------------------------
n = 4
psi4 = run_circuit(ghz_log(n), n)
print(f"\n{n}量子ビットのGHZ、2000ショット、seed 7:")
print(f"  {sample(psi4, 2000, seed=7)}")
print(f"  振幅がゼロでない基底: "
      f"{[format(i, f'0{n}b') for i in np.flatnonzero(np.abs(psi4) > 1e-12)]}")
print(f"  深さ chain / log    : {circuit_depth(ghz_chain(n), n)} / "
      f"{circuit_depth(ghz_log(n), n)}")

# ---- 回路の深さはゲート数を n で割ったものではありません ----------------
print("\n1量子ビットゲート4個、2通りの配置:")
for label, circ in [("すべて別の量子ビット上", [("h", q) for q in range(4)]),
                    ("すべて量子ビット0上  ",
                     [("h", 0), ("x", 0), ("z", 0), ("s", 0)])]:
    print(f"  {label} ゲート数 = {len(circ)}  深さ = "
          f"{circuit_depth(circ, 4)}")

Bell回路: [('h', 0), ('cx', 0, 1)]
  振幅            : [0.707107 0.       0.       0.707107]
  深さ            : 2
  ゲート数        : {'h': 1, 'cx': 1, '2q': 1}
  ノルム = 1      : True
  振幅の誤差      : 1/sqrt(2) に対して 1.11e-16, 1.11e-16

同一の n量子ビットGHZ状態に対する2つの回路
  n  chain depth  log depth  cx (chain)  cx (log)   max |dpsi|
--------------------------------------------------------------
  2            2          2           1         1     0.00e+00
  3            3          3           2         2     0.00e+00
  4            4          3           3         3     0.00e+00
  5            5          4           4         4     0.00e+00
  6            6          4           5         5     0.00e+00
  7            7          4           6         6     0.00e+00
  8            8          4           7         7     0.00e+00

4量子ビットのGHZ、2000ショット、seed 7:
  {'0000': 1013, '1111': 987}
  振幅がゼロでない基底: ['0000', '1111']
  深さ chain / log    : 4 / 3

1量子ビットゲート4個、2通りの配置:
  すべて別の量子ビット上 ゲート数 = 4  深さ = 1
  すべて量子ビット0上   ゲート数 = 4  深さ = 4

注目すべき点。 Bell状態はちょうど正しく出てきます。そしてノルムの検査は浮動小数点数を $1.0$ と比べるのではなく相対許容を使っています。いま身につけておく価値のある習慣です。第5章は10桁にわたる量を比べ、そこでは単純な等値検査は不作法どころか誤りになります。GHZの表は本コースの中心主題の初登場です。

どちらの構成も2量子ビットゲートをちょうど $n-1$ 個使い、どちらも最後のビットまで同じ状態を作ります。違うのは深さだけです。鎖は $n$、倍化構成は $\lceil \log_2 n \rceil + 1$ で、$n = 8$ では8対4です。デコヒーレンスで制限される機械では後者が2倍よく、2量子ビットゲート誤りで制限される機械では両者は同一です。コンパイラを1つも書く前に、同じ回路を2通りにコンパイルしたものが、ある指標では違い別の指標では違わない — これが1.3節で「最適化」には目的関数が要ると言い張った理由です。最後のブロックは深さの定義を具体的にします。4量子ビット上の4ゲートは深さ1、1量子ビット上の4ゲートは深さ4であり、ゲート数では区別できません。

Code Example 4: 等価性チェッカと、それが検証する規則

ここで本コースの残りが乗る道具を作ります。

"""Chapter 1, Example 4: 等価性チェッカと、それが検証する書き換え規則。
Code Example 3 の続き（同一セッション）。"""
import numpy as np
from qir import *


def unitary_of(circ, n):
    """回路の 2^n x 2^n 行列: 各基底状態に対して1回ずつ実行します"""
    dim = 2 ** n
    U = np.empty((dim, dim), dtype=complex)
    for j in range(dim):
        e = np.zeros(dim, dtype=complex)
        e[j] = 1.0
        U[:, j] = run_circuit(circ, n, psi0=e)
    return U


def best_global_phase(U, V):
    """e^{i phi} V を U に位相だけで最も近づける位相を返します。

    これはHilbert-Schmidt重なり tr(V^dagger U) から得られます。Cauchy-Schwarzに
    より、その絶対値が 2^n に達するのは U = e^{i phi} V のときに限るので、
    重なりがほぼゼロであること自体が両者が非等価であることの証明になります。
    """
    tr = np.trace(V.conj().T @ U)
    return 1.0 + 0.0j if abs(tr) < 1e-12 else tr / abs(tr)


def phase_free_error(U, V):
    """最良の大域位相を除去したあとの max |U - e^{i phi} V|"""
    return float(np.max(np.abs(U - best_global_phase(U, V) * V)))


def assert_equivalent(a, b, n, label="", atol=1e-10):
    """本コースのすべての書き換えパスを守るテストです"""
    err = phase_free_error(unitary_of(a, n), unitary_of(b, n))
    if err > atol:
        raise AssertionError(f"{label}: circuits differ, max error {err:.3e}")
    return err


def verdict(U, V, atol=1e-10):
    """比較が返しうる3つの答えのうちどれか"""
    if np.max(np.abs(U - V)) <= atol:
        return "identical"
    if phase_free_error(U, V) <= atol:
        return "up to phase"
    return "DIFFERENT"


PI = np.pi
RULES = [
    (1, "H H = I", [("h", 0), ("h", 0)], []),
    (1, "X X = I", [("x", 0), ("x", 0)], []),
    (1, "S S = Z", [("s", 0), ("s", 0)], [("z", 0)]),
    (1, "T T = S", [("t", 0), ("t", 0)], [("s", 0)]),
    (1, "H Z H = X", [("h", 0), ("z", 0), ("h", 0)], [("x", 0)]),
    (1, "H X H = Z", [("h", 0), ("x", 0), ("h", 0)], [("z", 0)]),
    (1, "Rz(a) Rz(b) = Rz(a+b)",
     [("rz", 0.3, 0), ("rz", 0.7, 0)], [("rz", 1.0, 0)]),
    (1, "Rz(pi) = Z", [("rz", PI, 0)], [("z", 0)]),
    (1, "Rz(pi/2) = S", [("rz", PI / 2, 0)], [("s", 0)]),
    (2, "CX CX = I", [("cx", 0, 1), ("cx", 0, 1)], []),
    (2, "CZ = H(1) CX H(1)",
     [("cz", 0, 1)], [("h", 1), ("cx", 0, 1), ("h", 1)]),
    (2, "CZ(0,1) = CZ(1,0)", [("cz", 0, 1)], [("cz", 1, 0)]),
    (2, "CX(0,1) = CX(1,0)", [("cx", 0, 1)], [("cx", 1, 0)]),
    (2, "H,H CX(0,1) H,H = CX(1,0)",
     [("h", 0), ("h", 1), ("cx", 0, 1), ("h", 0), ("h", 1)], [("cx", 1, 0)]),
    (4, "GHZ chain = GHZ doubling", ghz_chain(4), ghz_log(4)),
]

head = (f"{'rule':<28}{'n':>3}{'max |U-V|':>12}{'phase-free':>12}"
        f"{'phase/pi':>10}  verdict")
print(head)
print("-" * len(head))
for n, label, left, right in RULES:
    U, V = unitary_of(left, n), unitary_of(right, n)
    phi = np.angle(best_global_phase(U, V)) / PI
    phi = phi if abs(phi) > 1e-9 else 0.0
    print(f"{label:<28}{n:>3}{np.max(np.abs(U - V)):>12.2e}"
          f"{phase_free_error(U, V):>12.2e}{phi:>10.3f}  {verdict(U, V)}")

print("\nこの表のうち興味深いのは2行です。")
print("  Rz(pi) と Z は観測できない因子 exp(-i pi/2) = -i だけ異なります。")
print("  GHZ鎖とGHZ倍化は異なるユニタリですが、回路が実際に与えられる唯一の")
print("  入力状態の上では一致します:")
for n in (4, 6, 8):
    a, b = ghz_chain(n), ghz_log(n)
    on_zero = np.max(np.abs(run_circuit(a, n) - run_circuit(b, n)))
    as_unitary = phase_free_error(unitary_of(a, n), unitary_of(b, n))
    print(f"    n = {n}: |0...0> 上の誤差 = {on_zero:.2e},  "
          f"ユニタリとして = {as_unitary:.2e}")


# ---- 罠: 制御ブロックの内側では大域位相は大域的ではありません -------------
def controlled(U):
    """量子ビット0が |1> のとき量子ビット1に U を作用させる2量子ビットゲート"""
    C = np.eye(4, dtype=complex)
    C[2:, 2:] = U
    return C


print("\nゲートが制御された瞬間に、大域位相は大域的でなくなります:")
print(f"  位相除去後の誤差  Rz(pi) vs Z    : "
      f"{phase_free_error(rz(PI), Z):.2e}")
print(f"  位相除去後の誤差  C-Rz(pi) vs CZ : "
      f"{phase_free_error(controlled(rz(PI)), CZ4):.2e}")
print("  C-Rz(pi) の対角成分: "
      + " ".join(f"{v:+.3f}" for v in np.diag(controlled(rz(PI)))))
print("  CZ の対角成分      : "
      + " ".join(f"{v:+.3f}" for v in np.diag(CZ4)))

rule                          n   max |U-V|  phase-free  phase/pi  verdict
--------------------------------------------------------------------------
H H = I                       1    2.22e-16    2.22e-16     0.000  identical
X X = I                       1    0.00e+00    0.00e+00     0.000  identical
S S = Z                       1    0.00e+00    0.00e+00     0.000  identical
T T = S                       1    2.22e-16    1.11e-16     0.000  identical
H Z H = X                     1    2.22e-16    2.22e-16     0.000  identical
H X H = Z                     1    2.22e-16    2.22e-16     0.000  identical
Rz(a) Rz(b) = Rz(a+b)         1    1.11e-16    1.11e-16     0.000  identical
Rz(pi) = Z                    1    1.41e+00    6.12e-17    -0.500  up to phase
Rz(pi/2) = S                  1    7.65e-01    1.11e-16    -0.250  up to phase
CX CX = I                     2    0.00e+00    0.00e+00     0.000  identical
CZ = H(1) CX H(1)             2    2.22e-16    2.22e-16     0.000  identical
CZ(0,1) = CZ(1,0)             2    0.00e+00    0.00e+00     0.000  identical
CX(0,1) = CX(1,0)             2    1.00e+00    1.00e+00     0.000  DIFFERENT
H,H CX(0,1) H,H = CX(1,0)     2    3.33e-16    3.33e-16     0.000  identical
GHZ chain = GHZ doubling      4    7.07e-01    7.07e-01     0.000  DIFFERENT

この表のうち興味深いのは2行です。
  Rz(pi) と Z は観測できない因子 exp(-i pi/2) = -i だけ異なります。
  GHZ鎖とGHZ倍化は異なるユニタリですが、回路が実際に与えられる唯一の
  入力状態の上では一致します:
    n = 4: |0...0> 上の誤差 = 0.00e+00,  ユニタリとして = 7.07e-01
    n = 6: |0...0> 上の誤差 = 0.00e+00,  ユニタリとして = 7.07e-01
    n = 8: |0...0> 上の誤差 = 0.00e+00,  ユニタリとして = 7.07e-01

ゲートが制御された瞬間に、大域位相は大域的でなくなります:
  位相除去後の誤差  Rz(pi) vs Z    : 6.12e-17
  位相除去後の誤差  C-Rz(pi) vs CZ : 7.65e-01
  C-Rz(pi) の対角成分: +1.000+0.000j +1.000+0.000j +0.000-1.000j +0.000+1.000j
  CZ の対角成分      : +1.000+0.000j +1.000+0.000j +1.000+0.000j -1.000+0.000j

注目すべき点。 15の規則のうち11が $10^{-16}$ の床でidenticalと出ます。教科書的な回路恒等式であり、これを信用で済ませるのはこれが最後になるべきです。残る4行が内容を担います。Rz(pi) = Z と Rz(pi/2) = S は位相を除いて等価であり、同一ではありません。素朴な比較は $1.41$ と $0.765$ という誤りを報告します。要素の絶対値が最大1のユニタリ行列の尺度では巨大です。一方で位相除去後の比較は $10^{-17}$ を報告します。大域位相除去なしのチェッカはここで正しい書き換えを棄却し、コンパイラはそのような書き換えを多数含むので、ほとんどすべてを棄却することになります。復元された位相 $-0.5\pi$ と $-0.25\pi$ はちょうど $-\theta/2$ です。回転ゲートは $R_z(\theta) = e^{-i\theta Z/2}$ という位相規約を担い、位相ゲートは担わないからです。

CX(0,1) = CX(1,0) は誤り $1.0$ で正しく棄却されます。制御と標的を入れ替えるとゲートが変わり、チェッカはそう言います。それを成り立たせる方法が次の行です。両方の量子ビットをHadamardで共役すると役割が入れ替わります。CNOTが非対称なのは計算基底においてのみだからです。

GHZの行は本章で最も教育的な失敗です。Code Example 3 の2つの構成は $\lvert 0\ldots0\rangle$ 上で最後のビットまで一致し、ユニタリとしては $0.707$ 違います。どちらも真であり、別の問いに答えています。状態等価性は回路が固定の入力をもつ状態準備であるときに欲しいものです。ユニタリ等価性は、回路が他の何かに適用される、あるいは制御される、あるいは逆にされるサブルーチンであるときに必ず必要なものです。状態等価性で正当化したパスをより大きな回路の内側で適用するのはバグであり、しかも回路が最上位で使われる限り動いてしまうので厄介です。本コースが一貫してユニタリ等価性を使うのはそのためです。

制御ゲートのブロックは1.3節の規則を数値にしたものです。$\mathrm{C}$-$R_z(\pi)$ の対角は $(1, 1, -i, i)$、CZの対角は $(1,1,1,-1)$ であり、捨てた因子 $-i$ が制御の2つの枝の間の相対位相になっています。

Code Example 5: 3つのパスと、チェッカが捕まえるバグ

最適化器を1つ、そして意図的に壊した最適化器を1つ。

"""Chapter 1, Example 5: 3つの書き換えパスと、それを守るテスト。
Code Example 4 の続き（同一セッション）。"""
import numpy as np
from qir import *

SELF_INVERSE = {"h", "x", "z", "cx", "cz"}
NAMES = ["h", "x", "z", "s", "t", "rx", "ry", "rz", "cx", "cz"]
ANGLES = [k * np.pi / 4 for k in (-3, -2, -1, 1, 2, 3, 4)]


def next_touching(circ, i, qs):
    """i より後で qs と量子ビットを共有する最初のゲートの添字、なければ len(circ)"""
    j = i + 1
    while j < len(circ) and not (qs & set(gate_qubits(circ[j]))):
        j += 1
    return j


def cancel_inverses(circ):
    """パス1: 何も間に挟まない、同一の自己逆ゲートの対を削除します"""
    out = list(circ)
    i = 0
    while i < len(out):
        g = out[i]
        if g[0] in SELF_INVERSE:
            j = next_touching(out, i, set(gate_qubits(g)))
            if j < len(out) and out[j] == g:
                del out[j], out[i]
                i = max(i - 1, 0)
                continue
        i += 1
    return out


def fuse_rotations(circ):
    """パス2: 同じ量子ビット上の同じ軸まわりの隣接回転を融合します"""
    out = list(circ)
    i = 0
    while i < len(out):
        g = out[i]
        if g[0] in ROT_1Q:
            j = next_touching(out, i, {g[2]})
            if j < len(out) and out[j][0] == g[0]:
                out[i] = (g[0], g[1] + out[j][1], g[2])
                del out[j]
                continue
        i += 1
    return out


def drop_null_rotations(circ, atol=1e-12):
    """パス3: 厳密に恒等となる 4 pi の整数倍の回転を削除します"""
    period = 4 * np.pi
    keep = []
    for g in circ:
        if g[0] in ROT_1Q:
            residual = abs((g[1] + period / 2) % period - period / 2)
            if residual < atol:
                continue
        keep.append(g)
    return keep


PASSES = [cancel_inverses, fuse_rotations, drop_null_rotations]


def optimize(circ, n, passes=PASSES, max_rounds=50):
    """パスを不動点まで走らせ、各パスの後で等価性を検査します"""
    current = list(circ)
    for _ in range(max_rounds):
        before = current
        for p in passes:
            candidate = p(current)
            assert_equivalent(current, candidate, n, label=p.__name__)
            current = candidate
        if current == before:
            break
    return current


def random_circuit(n, length, rng):
    """IRのゲート集合上のランダム回路。角度は pi/4 の整数倍です"""
    circ = []
    for _ in range(length):
        name = NAMES[int(rng.integers(len(NAMES)))]
        if name in TWO_Q:
            a, b = (int(v) for v in rng.choice(n, size=2, replace=False))
            circ.append((name, a, b))
        elif name in ROT_1Q:
            circ.append((name, ANGLES[int(rng.integers(len(ANGLES)))],
                         int(rng.integers(n))))
        else:
            circ.append((name, int(rng.integers(n))))
    return circ


# ---- 各サイズ200本のランダム回路を最適化し、検査します -------------------
head = (f"{'n':>3}{'length':>8}{'gates out':>11}{'2q out':>9}{'depth in':>10}"
        f"{'depth out':>11}{'worst error':>13}{'failures':>10}")
print(head)
print("-" * len(head))
for n, length in [(2, 20), (3, 30), (4, 40), (5, 50)]:
    rng = np.random.default_rng(1000 + n)
    tally = np.zeros(6)                  # ゲート数 in/out、2q in/out、深さ in/out
    worst, failures = 0.0, 0
    for _ in range(200):
        c = random_circuit(n, length, rng)
        o = optimize(c, n)
        err = phase_free_error(unitary_of(c, n), unitary_of(o, n))
        worst = max(worst, err)
        failures += int(err > 1e-10)
        tally += [len(c), len(o), gate_counts(c)["2q"], gate_counts(o)["2q"],
                  circuit_depth(c, n), circuit_depth(o, n)]
    print(f"{n:>3}{length:>8}{tally[1] / tally[0]:>10.1%}"
          f"{tally[3] / max(tally[2], 1):>9.1%}{tally[4] / 200:>10.2f}"
          f"{tally[5] / 200:>11.2f}{worst:>13.2e}{failures:>10}")


# ---- そしてバグ入りのパス -------------------------------------------------
def buggy_commute(circ):
    """誤り: Z が CX の標的を通り抜けて可換であるかのように扱っています。
    Z が可換なのは制御側であって、標的側ではありません。"""
    out = list(circ)
    for i in range(len(out) - 1):
        g, h = out[i], out[i + 1]
        if g[0] == "cx" and h[0] == "z" and h[1] == g[2]:
            out[i], out[i + 1] = h, g
    return out


print("\nもっともらしいバグを含むパスと、チェッカの反応:")
trap = [("h", 0), ("cx", 0, 1), ("z", 1), ("h", 1)]
print(f"  回路          : {trap}")
print(f"  パス適用後    : {buggy_commute(trap)}")
print(f"  位相除去後の誤差: "
      f"{phase_free_error(unitary_of(trap, 2), unitary_of(buggy_commute(trap), 2)):.3f}")

rng = np.random.default_rng(5)
caught = tried = 0
for _ in range(500):
    c = random_circuit(3, 12, rng)
    b = buggy_commute(c)
    if b != c:
        tried += 1
        if phase_free_error(unitary_of(c, 3), unitary_of(b, 3)) > 1e-10:
            caught += 1
print(f"  ランダム性質テスト: 500本のうち {tried} 本でパスが発火し、"
      f"{caught} 本で捕まりました")
print(f"  Z を*制御*側で通した場合: "
      f"{phase_free_error(unitary_of([('cx', 0, 1), ('z', 0)], 2), unitary_of([('z', 0), ('cx', 0, 1)], 2)):.2e}")


# ---- 効果が出るのはランダムさではなく構造です ----------------------------
def cz_star(n):
    """他のすべての量子ビットから量子ビット n-1 へのCZを、2周分並べます"""
    ctrl = list(range(n - 1))
    return [("cz", q, n - 1) for q in ctrl + ctrl]


def expand_cz(circ):
    """基底変換のパス: すべてのCZを第2量子ビット上の H CX H に置き換えます"""
    out = []
    for g in circ:
        if g[0] == "cz":
            out += [("h", g[2]), ("cx", g[1], g[2]), ("h", g[2])]
        else:
            out.append(g)
    return out


print("\n構造をもつ回路: CZの星をCXに変換し、そのあと最適化します。")
head = (f"{'n':>3}{'cz':>5}{'gates after expand':>20}{'after optimize':>16}"
        f"{'h out':>7}{'depth in':>10}{'depth out':>11}{'error':>10}")
print(head)
print("-" * len(head))
for n in (4, 5, 6, 7):
    logical = cz_star(n)
    expanded = expand_cz(logical)
    assert_equivalent(logical, expanded, n, label="expand_cz")
    tuned = optimize(expanded, n)
    assert_equivalent(logical, tuned, n, label="expand+optimize")
    err = phase_free_error(unitary_of(logical, n), unitary_of(tuned, n))
    print(f"{n:>3}{gate_counts(logical)['cz']:>5}{len(expanded):>20}"
          f"{len(tuned):>16}{gate_counts(tuned).get('h', 0):>7}"
          f"{circuit_depth(expanded, n):>10}{circuit_depth(tuned, n):>11}"
          f"{err:>10.1e}")
print("CZを1個ずつ変換して生じた内側のHの対は、次のCZが生む対と打ち消し合い、")
print("両端の2個だけが残ります。")

  n  length  gates out   2q out  depth in  depth out  worst error  failures
---------------------------------------------------------------------------
  2      20     91.1%    93.9%     14.94      13.76     4.71e-16         0
  3      30     91.8%    96.6%     16.84      15.89     5.24e-16         0
  4      40     92.2%    96.9%     18.25      17.27     4.78e-16         0
  5      50     92.7%    98.1%     19.18      18.25     5.90e-16         0

もっともらしいバグを含むパスと、チェッカの反応:
  回路          : [('h', 0), ('cx', 0, 1), ('z', 1), ('h', 1)]
  パス適用後    : [('h', 0), ('z', 1), ('cx', 0, 1), ('h', 1)]
  位相除去後の誤差: 1.000
  ランダム性質テスト: 500本のうち 19 本でパスが発火し、18 本で捕まりました
  Z を*制御*側で通した場合: 0.00e+00

構造をもつ回路: CZの星をCXに変換し、そのあと最適化します。
  n   cz  gates after expand  after optimize  h out  depth in  depth out     error
----------------------------------------------------------------------------------
  4    6                  18               8      2        18          8   2.2e-16
  5    8                  24              10      2        24         10   2.2e-16
  6   10                  30              12      2        30         12   2.2e-16
  7   12                  36              14      2        36         14   2.2e-16
CZを1個ずつ変換して生じた内側のHの対は、次のCZが生む対と打ち消し合い、
両端の2個だけが残ります。

注目すべき点。 ランダム回路の表が誠実な測定であり、誠実な答えは、ピープホール最適化器はランダムな入力に対してはほとんど効かない、というものです。ゲート数の約 $8\%$、深さの約 $5\%$ です。これはパスの欠陥ではありません。ランダム回路には利用できる構造がなく、打ち消せる対もほとんどありません。本物の回路は構造に満ちています。人間かテンプレートが生成したものだからです。そして最後のブロックがそのときに何が起きるかを示します。800本の回路が最適化され、800本が検証され、最悪の食い違いは $5.9 \times 10^{-16}$、失敗0本。この最後の列が表の要点です。

バグはこの例の最も価値ある部分です。「$Z$ はCXを通り抜けて可換」は、パスに書き込まれてしまう類の半分覚えの規則です。制御側については真であり — 最後の行でちょうど $0.00$ と検証されています — 標的側については偽です。CXが掛ける $X$ と $Z$ は反可換だからです。チェッカは4ゲートの回路で誤り $1.0$ として捕まえます。そのあとの性質テストがもっと面白いことを言います。500本のランダム回路のうち、このバグ入りパスが書き換える対象を見つけたのは19本だけで、その19本のうち1本では、偶然まだ等価な回路が出てきました。ランダムテストはバグを見つけますが、発火させるだけの標本が必要であり、1回通っただけのテストは見た目ほどの意味をもちません。

構造をもつ例が輪を閉じます。CZを1個ずつ独立に $H\,\mathrm{CX}\,H$ に変換するのは自明な基底変換であり、そして無駄です。$n = 6$ では12個で足りるところに30個のゲートを生みます。共有量子ビット上の内側のHの対がすべて冗長だからです。それらを打ち消すとゲート18個と層18個が消えます。だからこそ本物のコンパイラは合成と最適化を1回ずつ走らせるのではなく交互に走らせるのであり、それが第2章の出発点です。

Code Example 6: 時間の予算と誤りの予算

同じユニタリを実装する2つの回路の費用は同じではありません。機械の関数としてどれだけ違うかを見ます。

"""Chapter 1, Example 6: 指標が買うもの — 時間の予算と誤りの予算。
Code Example 5 の続き（同一セッション）。"""
import math
import numpy as np
from qir import *


def circuit_layers(circ, n):
    """circuit_depth が数えている層分けを、ゲートのリストとして返します"""
    ready, layers = [0] * n, []
    for g in circ:
        qs = gate_qubits(g)
        k = max(ready[q] for q in qs)
        while len(layers) <= k:
            layers.append([])
        layers[k].append(g)
        for q in qs:
            ready[q] = k + 1
    return layers


def wall_clock(circ, n, t_2q):
    """1量子ビットゲート時間を単位とした所要時間を、層ごとに足します。

    t_2q は同じ単位で表した2量子ビットゲートの所要時間です。1つの層はその中で
    最も遅いゲートと同じ時間を要するので、深さだけでは時間は決まりません。
    """
    total = 0.0
    for layer in circuit_layers(circ, n):
        total += max(t_2q if g[0] in TWO_Q else 1.0 for g in layer)
    return total


def gate_error(circ, eps_1q, eps_2q):
    """1 からゲートごとの成功確率の積を引いた値"""
    n2 = gate_counts(circ)["2q"]
    n1 = len(circ) - n2
    return 1.0 - (1.0 - eps_1q) ** n1 * (1.0 - eps_2q) ** n2


def idle_error(circ, n, t_2q, t_coh):
    """回路の所要時間だけ保持された n量子ビットのデコヒーレンス"""
    return 1.0 - math.exp(-n * wall_clock(circ, n, t_2q) / t_coh)


# ---- 層分けは circuit_depth が数えているものと同一の対象です --------------
n = 6
logical = cz_star(n)
naive, tuned = expand_cz(logical), optimize(expand_cz(logical), n)
for label, circ in [("naive", naive), ("optimized", tuned)]:
    print(f"{label:>10}: {len(circuit_layers(circ, n))} 層, circuit_depth = "
          f"{circuit_depth(circ, n)}, ゲート数 = {len(circ)}")

# ---- 同じユニタリ、2つの費用。機械が非対称になるにつれて ------------------
print("\n2つの回路は同じユニタリを実装します。1量子ビットゲート時間を単位とした")
print("費用を、2量子ビットゲートが遅くなるにつれて示します:")
head = (f"{'t_2q':>6}{'naive time':>12}{'tuned time':>12}{'saved':>8}"
        f"{'naive/depth':>13}{'tuned/depth':>13}")
print(head)
print("-" * len(head))
for t_2q in (1.0, 3.0, 10.0, 30.0, 100.0):
    a, b = wall_clock(naive, n, t_2q), wall_clock(tuned, n, t_2q)
    print(f"{t_2q:>6.0f}{a:>12.0f}{b:>12.0f}{1 - b / a:>8.1%}"
          f"{a / circuit_depth(naive, n):>13.2f}"
          f"{b / circuit_depth(tuned, n):>13.2f}")

# ---- どちらの個数が誤りを支配するか --------------------------------------
print("\n最適化後の回路の誤り予算。1量子ビット誤り率を 1e-4 に固定し、")
print("2量子ビット誤り率を数桁にわたってスイープします:")
eps_1q = 1e-4
head = (f"{'eps_2q':>9}{'ratio':>8}{'from 1q':>11}{'from 2q':>11}"
        f"{'total':>11}{'2q share':>10}")
print(head)
print("-" * len(head))
n2 = gate_counts(tuned)["2q"]
n1 = len(tuned) - n2
for eps_2q in (1e-4, 3e-4, 1e-3, 3e-3, 1e-2):
    e1 = 1.0 - (1.0 - eps_1q) ** n1
    e2 = 1.0 - (1.0 - eps_2q) ** n2
    tot = gate_error(tuned, eps_1q, eps_2q)
    print(f"{eps_2q:>9.0e}{eps_2q / eps_1q:>8.0f}{e1:>11.2e}{e2:>11.2e}"
          f"{tot:>11.2e}{e2 / (e1 + e2):>10.1%}")
print(f"  この回路は1量子ビットゲート {n1} 個、2量子ビットゲート {n2} 個です")

# ---- 時間とコヒーレンスを合わせて ----------------------------------------
print("\nデコヒーレンスを加えます。コヒーレンス時間も同じ単位です:")
head = (f"{'t_coh':>9}{'naive gate':>12}{'naive idle':>12}{'naive tot':>11}"
        f"{'tuned tot':>11}{'improvement':>13}")
print(head)
print("-" * len(head))
for t_coh in (1e3, 1e4, 1e5, 1e6):
    ga = gate_error(naive, eps_1q, 1e-3)
    ia = idle_error(naive, n, 10.0, t_coh)
    gb = gate_error(tuned, eps_1q, 1e-3)
    ib = idle_error(tuned, n, 10.0, t_coh)
    ta, tb = 1 - (1 - ga) * (1 - ia), 1 - (1 - gb) * (1 - ib)
    print(f"{t_coh:>9.0e}{ga:>12.2e}{ia:>12.2e}{ta:>11.2e}{tb:>11.2e}"
          f"{ta / tb:>12.2f}x")


# ---- 浮動小数点数を10の冪と比べるとき ------------------------------------
def relative_verdict(value, reference, rel_tol=1e-9):
    """10の冪との比較を、決して等号を使わずに行います"""
    if math.isclose(value, reference, rel_tol=rel_tol):
        return "at"
    return "below" if value < reference else "above"


print("\nしきい値の比較に相対許容を付ける理由:")
budget = 0.0
for _ in range(100):
    budget += 1e-5           # ゲート1個ずつ積み上げた誤り予算
cubed = (1e-1) ** 3          # 同じ数に別の経路で到達したもの
print(f"  リテラルの 1e-3         : {1e-3!r}")
print(f"  1e-5 を100回足した値    : {budget!r}   == 1e-3 か: "
      f"{budget == 1e-3}")
print(f"  (1e-1) ** 3             : {cubed!r}   == 1e-3 か: "
      f"{cubed == 1e-3}")
print(f"  それぞれの相対判定      : "
      f"{relative_verdict(budget, 1e-3)}, {relative_verdict(cubed, 1e-3)}")
tot = gate_error(tuned, eps_1q / 10, 1e-4)
print(f"  eps_2q = 1e-4 のこの回路: {tot:.6e}  -> 1e-3 より "
      f"{relative_verdict(tot, 1e-3)}")

     naive: 30 層, circuit_depth = 30, ゲート数 = 30
 optimized: 12 層, circuit_depth = 12, ゲート数 = 12

2つの回路は同じユニタリを実装します。1量子ビットゲート時間を単位とした
費用を、2量子ビットゲートが遅くなるにつれて示します:
  t_2q  naive time  tuned time   saved  naive/depth  tuned/depth
----------------------------------------------------------------
     1          30          12   60.0%         1.00         1.00
     3          50          32   36.0%         1.67         2.67
    10         120         102   15.0%         4.00         8.50
    30         320         302    5.6%        10.67        25.17
   100        1020        1002    1.8%        34.00        83.50

最適化後の回路の誤り予算。1量子ビット誤り率を 1e-4 に固定し、
2量子ビット誤り率を数桁にわたってスイープします:
   eps_2q   ratio    from 1q    from 2q      total  2q share
------------------------------------------------------------
    1e-04       1   2.00e-04   1.00e-03   1.20e-03     83.3%
    3e-04       3   2.00e-04   3.00e-03   3.20e-03     93.7%
    1e-03      10   2.00e-04   9.96e-03   1.02e-02     98.0%
    3e-03      30   2.00e-04   2.96e-02   2.98e-02     99.3%
    1e-02     100   2.00e-04   9.56e-02   9.58e-02     99.8%
  この回路は1量子ビットゲート 2 個、2量子ビットゲート 10 個です

デコヒーレンスを加えます。コヒーレンス時間も同じ単位です:
    t_coh  naive gate  naive idle  naive tot  tuned tot  improvement
--------------------------------------------------------------------
    1e+03    1.19e-02    5.13e-01   5.19e-01   4.63e-01        1.12x
    1e+04    1.19e-02    6.95e-02   8.06e-02   6.89e-02        1.17x
    1e+05    1.19e-02    7.17e-03   1.90e-02   1.62e-02        1.17x
    1e+06    1.19e-02    7.20e-04   1.26e-02   1.08e-02        1.18x

しきい値の比較に相対許容を付ける理由:
  リテラルの 1e-3         : 0.001
  1e-5 を100回足した値    : 0.001000000000000002   == 1e-3 か: False
  (1e-1) ** 3             : 0.0010000000000000002   == 1e-3 か: False
  それぞれの相対判定      : at, at
  eps_2q = 1e-4 のこの回路: 1.019530e-03  -> 1e-3 より above

注目すべき点。 最初の表が1.3節の問いに答えます。最適化後の回路は、2量子ビットゲートが1量子ビットゲートと同じ速さの機械では素朴版より $60\%$ 速く、100倍遅い機械では $1.8\%$ 速いのです。最適化器が取り除いたのがすべて1量子ビットゲートであり、2量子ビットゲートがクロックを支配すると1量子ビットゲートは効かなくなるからです。同じ書き換え、同じ2つの回路で、その価値が30倍違います。最適化器に目的関数が要るとはこういうことです。目的はターゲット機械の性質であり、コンパイラには教えてやらなければなりません。

naive/depth と tuned/depth の列は層1つの平均所要時間を与えており、2量子ビットゲートが遅くなるにつれて1から34と84へ登ります。深さが時間の代理になるのは層が似ているときだけで、最適化の後では似ていません。最適化後の回路は2量子ビットゲートについて密なので、層1つあたりの費用が個別に高いのです。安いゲートを取り除いて深さを下げたパスは、誤った数を改善したことになります。誤りの表は一方向の結果です。

2量子ビット誤りの分担は、2つの誤り率が等しいときすでに $83\%$ です。最適化後の回路が2量子ビットゲート10個と1量子ビットゲート2個からなるからです。そして2量子ビット率が100倍悪いときには $99.8\%$ に達します。だからこそ本物のコンパイラの費用関数は2量子ビット数に支配され、第2章は総ゲート数ではなくCNOT数に労力を使うのです。デコヒーレンスの表がクロスオーバーを示します。コヒーレンス時間がゲート時間の1000倍しかないとき、待ちが誤りを完全に支配し — ゲート由来の $0.0119$ に対して $0.513$ — 最適化器による回路の短縮は全体で $1.12$ 倍の価値です。コヒーレンス時間 $10^6$ では待ちの項は無視できるようになり、改善は $1.18$ に落ち着きます。この2つの領域は違う目的関数を欲しています。待ちが支配するなら深さ、ゲートが支配するなら2量子ビット数です。

最後のブロックは本コースが至るところで適用する浮動小数点の規則です。$10^{-5}$ を100回足して積み上げた誤り予算は $0.001000000000000002$ になり、$(10^{-1})^3$ は $0.0010000000000000002$ になります。どちらもリテラルの 1e-3 と等しくなく、== で書いたしきい値検査、あるいは演算がちょうど乗るはずの境界に対する不等号は狂います。本コースの10の冪に対する比較はすべて相対許容を通します。10の冪を掛け合わせて物理量子ビット数に到達する第5章が、これが細かい話でなくなる場所です。

Code Example 7: スタックをパイプラインとして

最後の例は、以降の章が埋めていく形にすべてを組み立てます。

"""Chapter 1, Example 7: スタック全体を、守られたパスのパイプラインとして。
Code Example 6 の続き（同一セッション）。"""
import numpy as np
from qir import *


def swap_gates(a, b):
    """SWAPをCX 3個で。IRのゲート集合で量子ビットを動かす唯一の方法です"""
    return [("cx", a, b), ("cx", b, a), ("cx", a, b)]


def route_on_line(circ, n):
    """第3層のスタブ: すべての2量子ビットゲートを直線上の隣接対に載せます。

    離れたゲートはSWAPで寄せてから、そのSWAPを戻します。したがって量子ビットと
    ワイヤの対応は最初と最後で同じであり、回路は厳密に等価なままです。第3章では
    より安価なこと、すなわち置換をそのまま残してコンパイラの状態に持つことを
    行います。
    """
    out = []
    for g in circ:
        if g[0] in TWO_Q and abs(g[1] - g[2]) > 1:
            c, t = g[1], g[2]
            d = 1 if c > t else -1
            pairs, p = [], t
            while p + d != c:
                pairs.append((p, p + d))
                p += d
            forward = [s for pair in pairs for s in swap_gates(*pair)]
            out += forward + [(g[0], c, p)] + forward[::-1]
        else:
            out.append(g)
    return out


# ---- 第5層のスタブ: ネイティブ基底 {rz, ry, cx} への変換 ------------------
PI = np.pi
NATIVE_SET = {"rz", "ry", "cx"}
NATIVE_RULES = {
    "z": lambda q: [("rz", PI, q)],
    "s": lambda q: [("rz", PI / 2, q)],
    "t": lambda q: [("rz", PI / 4, q)],
    "x": lambda q: [("rz", PI, q), ("ry", PI, q)],
    "h": lambda q: [("rz", PI, q), ("ry", PI / 2, q)],
}


def translate_to_native(circ):
    """すべてのゲートを {rz, ry, cx} に書き換えます。各規則は位相を除いて成立します"""
    out = []
    for g in circ:
        name = g[0]
        if name in NATIVE_SET:
            out.append(g)
        elif name in NATIVE_RULES:
            out += NATIVE_RULES[name](g[1])
        elif name == "rx":
            theta, q = g[1], g[2]
            out += [("rz", PI / 2, q), ("ry", theta, q), ("rz", -PI / 2, q)]
        elif name == "cz":
            out += translate_to_native([("h", g[2])])
            out.append(("cx", g[1], g[2]))
            out += translate_to_native([("h", g[2])])
        else:
            raise ValueError(f"no native rule for {name!r}")
    return out


print("すべての変換規則を、使う前に検査します:")
head = f"{'rule':<34}{'phase-free error':>18}"
print(head)
print("-" * len(head))
for name, circ in [("z  -> rz(pi)", [("z", 0)]),
                   ("s  -> rz(pi/2)", [("s", 0)]),
                   ("t  -> rz(pi/4)", [("t", 0)]),
                   ("x  -> rz(pi) ry(pi)", [("x", 0)]),
                   ("h  -> rz(pi) ry(pi/2)", [("h", 0)]),
                   ("rx(0.7) -> rz ry rz", [("rx", 0.7, 0)])]:
    err = phase_free_error(unitary_of(circ, 1),
                           unitary_of(translate_to_native(circ), 1))
    print(f"{name:<34}{err:>18.2e}")
err = phase_free_error(unitary_of([("cz", 0, 1)], 2),
                       unitary_of(translate_to_native([("cz", 0, 1)]), 2))
print(f"{'cz -> h cx h, then h expanded':<34}{err:>18.2e}")
print("どの規則も大域位相を除いてのみ厳密です。だからこそ、どれ1つとして")
print("制御ゲートの中身に適用してはいけません。")

# ---- パイプライン --------------------------------------------------------
STACK = [("optimize", lambda c, n: optimize(c, n)),
         ("route on a line", route_on_line),
         ("translate to native", lambda c, n: translate_to_native(c)),
         ("optimize again", lambda c, n: optimize(c, n))]


def compile_stack(circ, n, verbose=True):
    """スタックを走らせ、各段の後で入力との等価性を検査します"""
    rows = [("input", circ)]
    current = circ
    for label, stage in STACK:
        current = stage(current, n)
        assert_equivalent(circ, current, n, label=label)
        rows.append((label, current))
    if verbose:
        head = (f"{'stage':<22}{'gates':>7}{'2q':>5}{'depth':>7}"
                f"{'time@10':>9}{'error':>12}{'vs input':>11}")
        print(head)
        print("-" * len(head))
        for label, c in rows:
            print(f"{label:<22}{len(c):>7}{gate_counts(c)['2q']:>5}"
                  f"{circuit_depth(c, n):>7}{wall_clock(c, n, 10.0):>9.0f}"
                  f"{gate_error(c, 1e-4, 1e-3):>12.2e}"
                  f"{phase_free_error(unitary_of(circ, n), unitary_of(c, n)):>11.1e}")
    return current


# 2量子ビットゲートを意図的に非局所にした5量子ビット回路です。
n = 5
algorithm = [("h", 0), ("h", 1), ("t", 2),
             ("cz", 0, 4), ("cx", 0, 3), ("s", 4),
             ("rx", PI / 3, 2), ("cx", 4, 1), ("h", 4), ("h", 4),
             ("cz", 1, 3), ("rz", PI / 8, 0), ("rz", -PI / 8, 0)]
print(f"\n5量子ビット回路1本をスタック全体に通します"
      f"（入力は {len(algorithm)} ゲート）:")
final = compile_stack(algorithm, n)
print(f"\n出力に現れるネイティブゲート名: "
      f"{sorted(set(g[0] for g in final))}")
print(f"2量子ビットゲート: 論理 {gate_counts(algorithm)['2q']} 個 -> "
      f"物理 {gate_counts(final)['2q']} 個、"
      f"倍率 {gate_counts(final)['2q'] / gate_counts(algorithm)['2q']:.1f}")

# ---- 接続性が違えば同じ回路がどうなるか ----------------------------------
print("\n機械の接続性が違った場合の同じ回路:")
head = (f"{'connectivity':<20}{'2q out':>8}{'depth out':>11}{'error':>11}"
        f"{'check':>10}")
print(head)
print("-" * len(head))
for label, routed in [("all-to-all", optimize(algorithm, n)),
                      ("line 0-1-2-3-4", route_on_line(optimize(algorithm, n), n))]:
    tuned = optimize(translate_to_native(routed), n)
    print(f"{label:<20}{gate_counts(tuned)['2q']:>8}{circuit_depth(tuned, n):>11}"
          f"{gate_error(tuned, 1e-4, 1e-3):>11.2e}"
          f"{assert_equivalent(algorithm, tuned, n, label=label):>10.1e}")
print("この回路が要求する対は "
      f"{sorted((min(g[1], g[2]), max(g[1], g[2])) for g in algorithm if g[0] in TWO_Q)}")

すべての変換規則を、使う前に検査します:
rule                                phase-free error
----------------------------------------------------
z  -> rz(pi)                                6.12e-17
s  -> rz(pi/2)                              1.21e-16
t  -> rz(pi/4)                              2.48e-16
x  -> rz(pi) ry(pi)                         6.12e-17
h  -> rz(pi) ry(pi/2)                       1.19e-16
rx(0.7) -> rz ry rz                         1.24e-16
cz -> h cx h, then h expanded               2.09e-16
どの規則も大域位相を除いてのみ厳密です。だからこそ、どれ1つとして
制御ゲートの中身に適用してはいけません。

5量子ビット回路1本をスタック全体に通します（入力は 13 ゲート）:
stage                   gates   2q  depth  time@10       error   vs input
-------------------------------------------------------------------------
input                      13    4      6       42    4.89e-03    0.0e+00
optimize                    9    4      5       41    4.49e-03    1.2e-16
route on a line            57   52     53      521    5.12e-02    1.2e-16
translate to native        69   52     63      531    5.23e-02    3.6e-16
optimize again             69   52     63      531    5.23e-02    3.6e-16

出力に現れるネイティブゲート名: ['cx', 'ry', 'rz']
2量子ビットゲート: 論理 4 個 -> 物理 52 個、倍率 13.0

機械の接続性が違った場合の同じ回路:
connectivity          2q out  depth out      error     check
------------------------------------------------------------
all-to-all                 4         10   5.59e-03   4.3e-16
line 0-1-2-3-4            52         63   5.23e-02   3.6e-16
この回路が要求する対は [(0, 3), (0, 4), (1, 3), (1, 4)]

注目すべき点。 まず変換の表です。7つの規則、それぞれ使う前に検証され、大域位相を除去したあとではすべて $10^{-16}$ で厳密であり、除去する前はどれ1つも厳密ではありません。$\lbrace R_z, R_y, \mathrm{CX}\rbrace$ というネイティブ基底はIRのすべてのゲートを表現するのに十分です。これが、記述用集合とネイティブ集合が同じものに対する別の語彙であるという主張の内容です。ただし一方から他方への変換はどの行でも位相について雑であり、だからこそ制御ブロックについての規則が規則でなければならないのです。

パイプラインの表が本章から覚えておくべき数です。2量子ビットゲート4個を含む論理ゲート13個が入ります。最適化後は9ゲート。直線上へのルーティング後は2量子ビットゲート52個を含む57ゲート — 誤りを支配する数について13倍、実時間について10倍で、接続性以外に原因はありません。誤りの見積りもそれに従います。ルーティング前が $4.9 \times 10^{-3}$、後が $5.1 \times 10^{-2}$ です。そしてすべての段が入力に対して $10^{-16}$ で検査されています。

誠実な観察が3つ。第1に、このルータは意図的に素朴です。量子ビットとワイヤの対応を変えず等価性検査を厳密にできるように、自分のSWAPを戻すので、SWAP数が2倍になります。第3章は置換をそのまま残して半分にし、さらによい初期レイアウトを選んでずっとよくします。第2に、最後の optimize 段は何も見つけません。69ゲート入り、69ゲート出です。これは水増しではなく、Code Example 5 のピープホール規則がルータの後片付けには弱すぎるという測定です。ルータが生む打ち消せる対は、リスト上で隣接するのではなく他のゲートで隔てられているからです。第2章が可換性を意識した打ち消しでこれを直します。パスの順序とパスの強さが別の問題であることの好例です。第3に、all-to-all の行が反実仮想です。4つのゲートを直接実行できる機械上の同じ回路は2量子ビットゲート4個と深さ10で済み、52個と63に対する数字です。要求された対 — $(0,3), (0,4), (1,3), (1,4)$ — はまさに直線がもたない対です。

このパイプラインのすべての段はスタブであり、本コースの残りがそれらを1つずつ置き換えていきます。その間、防護柵である assert_equivalent はちょうどそのままの位置にとどまります。

演習

演習1: IRを拡張する

IRに2つのゲートを追加してください。Pauli $Y$ のための ("y", q) と、SWAPゲートのための ("swap", a, b) です。

1. 2つそれぞれについて、qir.py のどの部分を変えなければなりませんか。gate_qubits、apply_ir_gate、circuit_depth、gate_counts それぞれについて答えてください。
2. 両方を実装し、("swap", 0, 1) がCX 3個に等しいこと、("y", 0) が位相を除いて $Z$ のあとの $X$ に等しいことをチェッカで検証してください。
3. [("swap", 0, 1), ("y", 0), ("cx", 0, 1)] に対して gate_counts は何を報告しますか。ネイティブなSWAPをもたない機械にとって、それは正しい答えですか。
4. swap を使って書いた回路と、SWAPをCX 3個に展開した同じ回路の深さを比べてください。スケジューラはどちらの数を信じるべきですか。

import qir                                  # 自分の名前空間ではなくモジュールを書き換えます

qir.FIXED_1Q["y"] = Y                       # 1. 固定1量子ビットゲートを1つ追加
qir.TWO_Q = ("cx", "cz", "swap")            # 2. gate_qubits と gate_counts が追従
SWAP4 = np.eye(4, dtype=complex)[[0, 2, 1, 3]]
_base_apply = qir.apply_ir_gate


def apply_ir_gate(state, g, n):             # 3. ディスパッチに分岐を1つ追加
    if g[0] == "swap":
        return apply_gate(state, SWAP4, [g[1], g[2]], n)
    return _base_apply(state, g, n)


qir.apply_ir_gate = apply_ir_gate           # run_circuit は呼び出し時にこれを解決します

swap3 = [("cx", 0, 1), ("cx", 1, 0), ("cx", 0, 1)]
print(f"{phase_free_error(unitary_of([('swap', 0, 1)], 2), unitary_of(swap3, 2)):.1e}")
print(f"{phase_free_error(unitary_of([('y', 0)], 1), unitary_of([('z', 0), ('x', 0)], 1)):.1e}")
print(gate_counts([("swap", 0, 1), ("y", 0), ("cx", 0, 1)]))
print(circuit_depth([("swap", 0, 1), ("y", 0)], 2),
      circuit_depth(swap3 + [("y", 0)], 2))
# 0.0e+00
# 0.0e+00
# {'swap': 1, 'y': 1, 'cx': 1, '2q': 2}
# 2 4

演習2: 深さ、時間、そして貪欲な層分けが仮定していること

4量子ビット上の回路 $C = [(h,0), (cx,0,1), (h,2), (cx,2,3), (z,1), (z,3)]$ を取ります。

1. circuit_depth(C, 4) を手で計算し、各層の中身を挙げてください。そのあと検算してください。
2. 4量子ビット上に、深さが同じで実時間が $t_{2q}$ 倍まるごと違う2つの回路を作ってください。深さを目的関数にすることについて、これは何を言っていますか。
3. 貪欲な層分けは実際の実行時間の下限を与えます。それが省いているものを3つ挙げてください。
4. Code Example 6 の2つの誤りモデルのうち、深さが正しい目的関数になるのはどちらで、2量子ビットゲート数が正しいのはどちらですか。機械がどちらの領域にいるかを決めるのは何ですか。

c = [("h", 0), ("cx", 0, 1), ("h", 2), ("cx", 2, 3), ("z", 1), ("z", 3)]
print(circuit_depth(c, 4), [len(v) for v in circuit_layers(c, 4)])
heavy = [("cz", 0, 1), ("cz", 2, 3), ("cz", 0, 1), ("cz", 2, 3)]
light = [("h", q) for q in range(4)] * 2
for label, circ in [("2q only", heavy), ("1q only", light)]:
    print(f"{label}: depth {circuit_depth(circ, 4)}, "
          f"time {wall_clock(circ, 4, 10.0):.0f}")
# 3 [2, 2, 2]
# 2q only: depth 2, time 20
# 1q only: depth 2, time 2

演習3: 大域位相の規則と、その修復

Code Example 4 は $R_z(\pi)$ と $Z$ が等価であり、$\mathrm{C}$-$R_z(\pi)$ と CZ はそうでないことを示しました。

1. $R_z(\pi) = e^{i\varphi} Z$ となる $\varphi$ を、手計算と best_global_phase の両方で求めてください。
2. $\mathrm{C}$-$R_z(\pi)$ と CZ の対角成分を書き下し、余分な因子を制御量子ビットだけに作用するゲートとして同定してください。
3. IRのゲートだけを使って、$\mathrm{C}$-$R_z(\pi)$ に厳密に等しい回路を作ってください。位相を除いてではありません。$S^{\dagger}$ が必要になりますが、これはゲート集合にないので、あるゲートで表してください。
4. コンパイラが大域位相を捨ててよいのはいつかを一般に述べ、したがって位相をデータとして持ち回らなければならないパスの例を1つ挙げてください。

print(f"{np.angle(best_global_phase(rz(PI), Z)) / PI:+.3f}")
print(f"{phase_free_error(controlled(rz(PI)), CZ4):.3f}")
fix = [("s", 0), ("s", 0), ("s", 0), ("cz", 0, 1)]     # S^3 = S^dagger
print(f"{np.max(np.abs(controlled(rz(PI)) - unitary_of(fix, 2))):.1e}")
# -0.500
# 0.765
# 6.1e-17

演習4: 検証の費用と、ランダムテストが買うもの

全数検査は回路を $2^n$ 回走らせます。これをランダムな入力状態での検査で置き換えることを考えます。

1. Haarランダム状態を引き、両方の回路を走らせ、重なりの絶対値 $\lvert \langle \psi_b \vert \psi_a \rangle \rvert$ を1と比べる equivalent_random(a, b, n, trials) を書いてください。なぜ状態の差ではなく重なりの絶対値なのですか。
2. $U(a) \ne e^{i\varphi} U(b)$ ならHaarランダム状態1個が確率1で両者を区別することを論じてください。浮動小数点演算でその議論が崩れるのはどこですか。
3. Code Example 5 のバグ入りパスに対して走らせ、数百本の回路で全数検査の判定と比べてください。
4. どの $n$ で切り替えますか。そしてどちらの方法も使えない $n = 30$ より上ではどうしますか。

def equivalent_random(a, b, n, trials=1, rng=None, atol=1e-10):
    """2^n 個すべてではなくランダムな入力状態で2つの回路を比べます"""
    rng = np.random.default_rng() if rng is None else rng
    for _ in range(trials):
        psi = rng.normal(size=2 ** n) + 1j * rng.normal(size=2 ** n)
        psi /= np.linalg.norm(psi)
        overlap = abs(np.vdot(run_circuit(b, n, psi0=psi),
                             run_circuit(a, n, psi0=psi)))
        if abs(overlap - 1.0) > atol:
            return False
    return True


rng = np.random.default_rng(11)
fired = agree = missed = 0
for _ in range(600):
    c = random_circuit(3, 12, rng)
    b = buggy_commute(c)
    if b == c:
        continue
    fired += 1
    exact = phase_free_error(unitary_of(c, 3), unitary_of(b, 3)) <= 1e-10
    single = equivalent_random(c, b, 3, trials=1, rng=rng)
    agree += int(exact == single)
    missed += int(single and not exact)
print(f"発火 {fired} 本、ランダム状態1個は全数検査と {agree} 本で一致、"
      f"見逃し {missed} 本")
# 発火 26 本、ランダム状態1個は全数検査と 26 本で一致、見逃し 0 本

演習5: スタックを読む

以下にコンパイラパスの説明が6つあり、フレームワーク名は付いていません。それぞれについて、1.1節のどの層に属するか、どのハードウェア情報を必要とするか、そしてその出力を別の機械で再利用できるかを述べてください。

1. 「同じ量子ビット上に連続する1量子ビットゲートを1つの $U(2)$ にまとめ、それを再合成する。」
2. 「SWAPの推定個数を最小化するように、プログラムの量子ビットが最初にどの物理量子ビットに載るかを選ぶ。」
3. 「ターゲットの命令集合にない2量子ビットゲートを、集合内の等価な列で置き換える。」
4. 「ゲートの開始時刻が波形のサンプル周期の整数倍になるよう、遅延命令を挿入する。」
5. 「各基底状態を準備して測定して得た行列の逆行列を、測定されたビット列の分布に掛ける。」
6. 「各2量子ビットゲートを自分自身の3つのコピーで置き換え、コピー数ゼロへ外挿する。」

まとめ

要点

1. 7つの層、入力の生鮮度で並べたもの

• アルゴリズム、回路IR、最適化、配置とルーティング、ゲート合成、パルスとスケジューリング、読み出しと誤り緩和。境界を固定する基準は、各パスが機械について知らなければならないことです。何も、結合グラフ、ネイティブゲート集合、あるいは今日の較正値。
• 第2層から第5層は同じ対象 — ゲート列 — を消費し生産するので、パスは合成でき1つずつテストできます。型が変わるのは第6層で、だからパルスのプログラミングは感触が違います。

2. IRはデータであり、そのゲート集合が層の境界である

• 回路はタプルのリスト、ゲートは名前と省略可能な角度と整数のワイヤ、ビッグエンディアン。run_circuit、circuit_depth、gate_counts、そしてチェッカの4つの関数がインタフェースの全体です。
• IRは何を表現するのを拒むかによって定義されます。測定、時間、割り当て、注釈はここではすべて欠けており、それぞれを後の章が戻します。
• ゲート集合は3つ現れます。便宜で選ばれた記述用集合、物理で選ばれたネイティブ集合、誤り訂正符号で選ばれたClifford$+T$。前2つの間の変換は定数倍、3つ目への変換は $\log(1/\varepsilon)$ で伸びる倍率です。

3. コンパイルは意味を保つ書き換えであり、その意味は検査できる

• 3つの正しさの関係。大域位相を除いて厳密（最適化、合成）、量子ビットの置換を除いて厳密（ルーティング）、$\varepsilon$ の範囲で近似（離散合成、パルス）。$\varepsilon$ は回路にわたって足し算されます。
• 最適な位相はHilbert-Schmidt重なりから来て、$\lvert\mathrm{tr}(V^{\dagger}U)\rvert = 2^n$ 自身が等価性の完全な判定になります。
• 状態等価性はユニタリ等価性より弱いです。2つのGHZ構成は $\lvert 0\ldots 0\rangle$ 上で厳密に一致し、ユニタリとしては $0.707$ 違います。弱い関係で正当化したパスを強い関係のところで使うのは、隠れるバグです。

4. 大域位相は制御ブロックの内側では大域的でない

• $R_z(\pi) = -iZ$ は等価で、$\mathrm{C}$-$R_z(\pi) \ne \mathrm{CZ}$ が $0.765$。観測できない因子が観測できる相対位相になります。Code Example 7 のすべての変換規則が位相について雑であり、だからこの規則は原理ではなく実務で効きます。
• 演習3で求めた修復は、捨てた位相が制御量子ビット上の本物のゲートとして再登場することです。$\mathrm{C}$-$R_z(\pi) = \mathrm{CZ}\cdot(S^{\dagger}\otimes I)$。

5. 検証は小さい $n$ では安価で、それより上では性格を変えなければならない

• 全数検査の費用はシミュレーション $2^n$ 回。$n \approx 12$ までは快適で、$n = 20$ では手が届きません。Haarランダム入力状態1個は非等価な回路を確率1で区別し、バグ入りパスが発火した26本すべてで全数検査の判定と一致しました。
• ランダムテストは証明ではありません。バグ入りパスが発火したのは500本のうち19本だけで、そのうち1本では偶然等価な回路が出ました。シミュレートできる範囲より上では検証は構造的になり、コンパイラのバグは誤った文脈で適用された正しい規則になります。

6. 目的関数のない「最適化」は無意味である

• 総ゲート数、2量子ビットゲート数、深さ、実時間、$T$ 数は一致しません。同じ書き換えが、2量子ビットゲートが速い機械では実行時間の $60\%$、遅い機械では $1.8\%$ の価値でした。
• 誤りに占める2量子ビットの分担は、誤り率が等しいとき $83\%$、100倍差で $99.8\%$ でした。だからコンパイラの費用関数は2量子ビット数に支配されます。
• 深さは待ちが支配するときに正しい目的関数で、ゲートが支配するときは誤り。例の回路のクロスオーバーはコヒーレンス時間がゲート時間の $10^4$ 倍から $10^5$ 倍の間でした。

実践上の含意

• パスを書く前に等価性の検査を書いてください。本コースの残りのすべての書き換え規則には検査が付いており、費用は数行です。
• 計算した浮動小数点数を10の冪と == や素の不等号で比べないでください。ループで積み上げた $100 \times 10^{-5}$ は 1e-3 ではありません。相対許容を使ってください。
• 2量子ビットゲート数は、どのゲート集合とどの接続性を仮定しているかを述べてから引用してください。5量子ビット回路を直線にルーティングしたら13倍になりました。トランスパイル後の回路が予想外に大きいときは、どの層が膨らませたのかを問うてください。ルーティングと合成は違う理由で回路を膨らませ、違う手段で直ります。

この先へ

Code Example 7 のパイプラインが本コースの縮図であり、そのすべての段がスタブです。第2章が最適化器を置き換えます。ここのピープホール規則は隣接する同一ゲートしか打ち消せず、本物には可換性の規則、Euler分解による1量子ビット合成、正しいCNOT数での2量子ビット合成、そして $T$ ゲートが高価なゲートである理由の誠実な説明が必要です。第2章はまた、本章で最も目立った非結果 — パイプラインの最後の optimize 段がルーティング後に取り除くものを何も見つけなかったこと、そしてそれを見つけるにはパスが何を知らなければならないか — も説明します。

← シリーズトップ 第2章: 回路最適化とゲート合成 →

免責事項

• 本コンテンツは教育・研究・情報提供のみを目的としており、専門的な助言(法律・会計・技術的保証など)を提供するものではありません。
• 本コンテンツおよび付随するCode examplesは「現状有姿(AS IS)」で提供され、明示または黙示を問わず、商品性、特定目的適合性、権利非侵害、正確性・完全性、動作・安全性等いかなる保証もしません。
• 本コースに現れるゲートの所要時間・誤り率・コヒーレンス時間は、スケーリングと定数倍を露わにするために数桁にわたってスイープした無次元の母数です。装置仕様・測定値ではなく、特定の機械についての予測でもありません。また1.4節の層の対応はフレームワークの構造を一般的に述べたものであり、特定の製品や版について述べたものではありません。
• 外部リンク、第三者が提供するデータ・ツール・ライブラリ等の内容・可用性・安全性について、作成者および東北大学は一切の責任を負いません。
• 本コンテンツの利用・実行・解釈により直接的・間接的・付随的・特別・結果的・懲罰的損害が生じた場合でも、適用法で許容される最大限の範囲で、作成者および東北大学は責任を負いません。
• 本コンテンツの内容は、予告なく変更・更新・提供停止されることがあります。
• 本コンテンツの著作権・ライセンスは明記された条件(例: CC BY 4.0)に従います。当該ライセンスは通常、無保証条項を含みます。