「第2回 Microsoft Q# Coding Contest 2019冬」 に参加しました。今のところ唯一の量子プログラミングでの競技プログラミングコンテストで、使用言語はMicrosoft社開発のQ#です。
上位50位に入るとMicrosoft Quantumの限定Tシャツが貰えます。今回は13問中12問正解、総合48位という結果でギリギリ入賞出来ました。半年に1回開催されそう(公式アナウンス一覧)なので、半年後のために備忘録としてQ#の使い方・コンテスト問題の考え方解き方などのまとめを残しておきます。
・Microsoft Q# Coding Contest 2019冬 (問題 / 解答 / 順位表)
・Microsoft Q# Coding Contest 公式アナウンス一覧
・Q# 公式ドキュメント
・Quantum Ketas
・Q# Language Quick Reference
・Quantum logic gate - wikipedia
Q#備忘録
開発環境はMicrosoft公式HPを参考に、MacOS向けの「Visual Studio Code」等をインストールした。コマンドラインでdotnet new console -lang Q# --output Solutionと入力するとプロジェクトフォルダが作成される。フォルダ内にはC#で書かれた「Driver.cs」とQ#で書かれた「Operations.qs」が入っており、この2つのファイルを編集する。
using System;
using Microsoft.Quantum.Simulation.Core;
using Microsoft.Quantum.Simulation.Simulators;
namespace Solution{
class Driver{
static void Main(string[] args){
using (var qsim = new QuantumSimulator()) {
var (a,b) = Output.Run(qsim).Result; // Operations.qs内のOutput()関数を実行
}
}
}
}
Driver.csの方はこれ以上弄る必要は無い。
namespace Solution {
open Microsoft.Quantum.Canon;
open Microsoft.Quantum.Primitive;
open Microsoft.Quantum.Extensions.Diagnostics; // DumpMachine()関数に必要
open Microsoft.Quantum.Extensions.Math; // Sqrt,Sin関数などが使える
open Microsoft.Quantum.Extensions.Convert; // 型変換に必要 (int->doubleなど)
// 1つのキュービットを"desired"の状態に変化させる(※測定を使用する)
operation Set (desired: Result, q1: Qubit) : Unit {
let current = M(q1);
if (desired != current) {X(q1);}
}
//
operation Solve (qs : Qubit[]) : Unit {
// your code here
let N = Length(qs); // Qubit数取得
let theta = ArcSin(1.0/Sqrt(ToDouble(3))); // 算術演算
let array = IntArrayFromRange(1..1..N-1); // 配列生成 1~N-1まで、第二引数はstep幅
X(qs[0]); // Xゲート |0> <-> |1> のFlip
H(qs[0]); // Hゲート
Ry(theta, qs[0]); // Rゲート 回転行列 |0> -> cosθ/2|0> + sinθ/2|1> 等
let M0 = M(qs[0]); // 測定Measurement
RAll1(ArcSin(0.5), qs); // e^(iθ)|111><111| 全1の状態だけに位相因子を付ける
Controlled X ([qs[1],qs[2]], qs[0]); // Controlledゲート ([qubit配列],qubit) []が全てOneの場合にのみ実行
Controlled X (qs[1 .. N-1], qs[0]); // ↑と同じ
Controlled X (Subarray(array,qs), qs[0]); // ↑と同じ
Message($"Value = {M0}"); // コンソールに文字列・値を出力
DumpMachine("dump1.txt"); // 量子状態をテキストファイルに出力
}
// Driver.csから呼び出される関数、Solve(x)関数のqubit引数xを生成・代入するためだけのもの
operation Output () : (Int,Int) { // 返り値が無い場合は"Unit"と書く
body (...) {
let N = 3; // 不変変数の定義
mutable a = 1; // 可変変数の定義
set a = 2; // 可変変数に代入
using (x = Qubit[N]) { // 新たにqubitを生成するにはusing()関数を使う 生成時には全て|0>
Set(One, x[0]); // |0>|0> -> |0>|1> // qubitの0,1は「Zero」,「One」で表す
Set(Zero, x[0]); // |0>|1> -> |0>|0>
Solve(x); // Solve()関数を実行
ResetAll(x); // using()関数終了時には新たに定義したqubitを全て|0>にして開放しないとエラーが出る
}
return (1,1);
}
}
}
Operations.qsを主に編集する。こちらはQ#の文法で書かれている。今回主に使用した文法・関数を上のコード内に列挙している。コンテスト提出の際は、ローカル実行用に置いているだけのOutput()関数部分を除いたコードを提出する。
コマンドラインからdotnet runで計算を実行できる。
Value = Zero
$$
\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}}
$$
問題A
与えられるinput状態から指定された重ね合わせ状態を作成する問題。やることはほぼ決まっているので簡単。
A1. Generate state |00⟩ + |01⟩ + |10⟩
次の変換を実行せよ。
\begin{eqnarray}
\ket{00} &\rightarrow& \frac{1}{\sqrt3} \left( \ket{00} + \ket{01} + \ket{10} \right)
\end{eqnarray}
(解答1)
Hゲートで4状態の重ね合わせが作れる。後は$\ket{11}$状態を消去すればよい。新たなqubitを追加し、Controlled Xゲートにより$\ket{11}$状態の時のみそれをFlipさせてエンタングルメント状態を作る。
\begin{eqnarray}
\ket{00} &\rightarrow& \frac{1}{\sqrt4} \left( \ket{00} + \ket{01} + \ket{10} + \ket{11} \right) \\
&\rightarrow& \frac{1}{\sqrt4} \left( \ket{00} + \ket{01} + \ket{10} + \ket{11} \right) \otimes \ket{0} \\
&\rightarrow& \frac{1}{\sqrt4} \left( \ket{00} + \ket{01} + \ket{10} \right) \otimes \ket{0} + \frac{1}{\sqrt4} \ket{11} \otimes \ket{1} \\
\end{eqnarray}
第2qubitを測定して0が観測されれば終了(確率3/4)、1が観測された場合は初期状態$\ket{00}$にリセットして繰り返す。
namespace Solution {
open Microsoft.Quantum.Primitive;
open Microsoft.Quantum.Canon;
operation Set (desired: Result, q1: Qubit) : Unit {
let current = M(q1);
if (desired != current) {X(q1);}
}
operation Solve (qs : Qubit[]) : Unit {
// your code here
mutable flag = 1;
repeat { // repeat関数 until(条件)で終了 fixup{}が必要
using (y = Qubit[1]) {
H(qs[0]);
H(qs[1]);
Controlled X (qs, y[0]);
let measurement = M(y[0]);
if (measurement==One) {
set flag = 0;
X(y[0]);
Set(Zero, qs[0]);
Set(Zero, qs[1]);
}
}
}
until (flag == 1)
fixup {}
}
}
(解答2)
測定やリセットを行わなくても可能。第2qubitをRyゲートで回転。次に$\ket {00}$の第1qubitをControlled Hゲートで2分割する。
\begin{eqnarray}
\ket{00} &\rightarrow& \ket{0} \otimes \left( \sqrt{\frac{2}{3}} \ket{0} + \frac{1}{\sqrt3} \ket{1} \right) \\
&\rightarrow& \sqrt{\frac{2}{3}} \ket{00} + \frac{1}{\sqrt3} \ket{01} \\
&\rightarrow& \frac{1}{\sqrt3} \ket{00} + \frac{1}{\sqrt3} \ket{10} + \frac{1}{\sqrt3} \ket{01}
\end{eqnarray}
A2. Generate equal superposition of four basis states
N個qubitの$\ket{00...0}$状態から、指定される4状態$\ket{\psi_0},\ket{\psi_1},\ket{\psi_2},\ket{\psi_3}$の線形結合状態を作成せよ。
\begin{eqnarray}
\ket{00...0} &\rightarrow& \frac{1}{\sqrt4} \left( \ket{\psi_0} + \ket{\psi_1} + \ket{\psi_2} + \ket{\psi_3} \right)
\end{eqnarray}
A1の解答1手法を使用すると、正解確率は$\frac{4}{2^N}$になる。制約が$N\le16$かつ実行時間1sであり、$N>10$程度でTLEを起こす。4状態しかないので、解答2のように地道に作っていけばよい。
namespace Solution {
open Microsoft.Quantum.Primitive;
open Microsoft.Quantum.Canon;
operation Set (desired: Result, q1: Qubit) : Unit {
let current = M(q1);
if (desired != current) {X(q1);}
}
operation Solve (qs : Qubit[], bits : Bool[][]) : Unit {
let length = Length(qs);
using (y = Qubit[2]) {
H(y[0]); H(y[1]);
for (i in 0..length-1) {
// 0番目
if (bits[0][i]==true) {
X(y[0]); X(y[1]);
Controlled X (y, qs[i]);
X(y[0]); X(y[1]);
}
// 1番目
if (bits[1][i]==true) {
X(y[0]);
Controlled X (y, qs[i]);
X(y[0]);
}
// 2番目
if (bits[2][i]==true) {
X(y[1]);
Controlled X (y, qs[i]);
X(y[1]);
}
// 3番目
if (bits[3][i]==true) {
Controlled X (y, qs[i]);
}
}
// yを戻す
for (j in 0..length-1) {if (bits[0][j]==false) {X(qs[j]);} } // qsの一部を反転
Controlled X (qs, y[0]); Controlled X (qs, y[1]);
for (j in 0..length-1) {if (bits[0][j]==false) {X(qs[j]);} } // qsの一部を再反転して戻す
for (j in 0..length-1) {if (bits[1][j]==false) {X(qs[j]);} } // qsの一部を反転
Controlled X (qs, y[0]);
for (j in 0..length-1) {if (bits[1][j]==false) {X(qs[j]);} } // qsの一部を再反転して戻す
for (j in 0..length-1) {if (bits[2][j]==false) {X(qs[j]);} } // qsの一部を反転
Controlled X (qs, y[1]);
for (j in 0..length-1) {if (bits[2][j]==false) {X(qs[j]);} } // qsの一部を再反転して戻す
X(y[0]); X(y[1]);
}
}
}
B問題
上位勢の正答率から見ても、ここが一番難しい。正答数は「B2 << D6 < B1 << その他」という感じだった。
B1. Distinguish three-qubit states
次の2状態$\ket{\psi_0},\ket{\psi_1}$の識別を行う。
\begin{eqnarray}
\ket{\psi_0} &=& \frac{1}{\sqrt3} ( \ket{100} + \omega \ket{010} + \omega^2 \ket{001} ) \\
\ket{\psi_1} &=& \frac{1}{\sqrt3} ( \ket{100} + \omega^2 \ket{010} + \omega \ket{001} ) \\
\omega &=& e^{\frac{2\pi i}{3}}
\end{eqnarray}
係数の並びが{$1,\omega,\omega^2$}と{$1,\omega^2,\omega$}なので、古典的には3状態に対する離散フーリエ変換を行えればそれぞれ一つの状態に収束する。しかし、2^3基底の3状態であり、そのまま量子離散フーリエ変換に乗せることは出来なかった。
\begin{eqnarray}
\begin{pmatrix}
1 & 1 & 1 \\
1 & \omega & \omega^2 \\
1 & \omega^2 & \omega \\
\end{pmatrix}
\begin{pmatrix}
\ket{100} \\
\ket{010} \\
\ket{001} \\
\end{pmatrix}
\end{eqnarray}
しかし、同じような離散フーリエ変換っぽい操作の逆変換を行えばいけそうであることがわかった。順変換として
\begin{eqnarray}
\ket{100} &\rightarrow& \frac{1}{\sqrt3} \ket{100} + \frac{\omega}{\sqrt3} \ket{010} + \frac{\omega^2}{\sqrt3} \ket{001}
\end{eqnarray}
を行うような操作を作り、これを逆に辿れば逆変換となる。この変換は問題A1の解法2とほぼ同じことをやれば良い。
具体的には、
\begin{eqnarray}
\ket{100} &\rightarrow& \frac{1}{\sqrt3} \ket{100} + \frac{1}{\sqrt3} \ket{010} + \frac{1}{\sqrt3} \ket{001}
\end{eqnarray}
のW状態を作り、位相部分をRAll1オペレータで選択的に調整する。W状態の作成はQuantumKatas/Superposition/ReferenceImplementation.qsのTask13のコードを参考にした。状態$\ket{\psi_1}$に対して完全な逆フーリエ変換を実現しているわけではないが、テスト段階で2状態の識別が成功していたのでこれで提出した。後から理由を考えると、直交する2状態$\ket{\psi_0},\ket{\psi_1}$に対してユニタリ変換を行った場合その直交性は保たれる
\begin{eqnarray}
\ket{\psi_0}, \ket{\psi_1} &\rightarrow& U \ket{\psi_0}, U \ket{\psi_1} \\
\braket{\psi_0}{\psi_1}=0 &\rightarrow& \bra{\psi_0} U^\dagger U \ket{\psi_1}=\braket{\psi_0}{\psi_1}=0
\end{eqnarray}
ので、$\ket{\psi_0}$の方さえある1状態に収束させてしまえば、$\ket{\psi_1}$は自動的にその状態以外の重ね合わせとなっているのでOKみたいな感じでした。離散フーリエ変換の話はあまり本質的では無かったです。
operation Solve (qs : Qubit[]) : Int {
body (...) {
let N = Length(qs);
// |W_N⟩ = |0⟩|W_(N-1)⟩ + |1⟩|0...0⟩
// do a rotation on the first qubit to split it into |0⟩ and |1⟩ with proper weights
// |0⟩ -> sqrt((N-1)/N) |0⟩ + 1/sqrt(N) |1⟩
let theta = ArcSin(1.0 / Sqrt(ToDouble(3)));
Ry(2.0 * theta, qs[0]);
// do a zero-controlled W-state generation for qubits 1..N-1
X(qs[0]);
let theta2 = ArcSin(1.0 / Sqrt(ToDouble(2)));
Controlled Ry ([qs[0]], (2.0*theta2, qs[1]) );
X(qs[0]);
// 000 -> 100
X(qs[0]); X(qs[1]);
Controlled X ([qs[1],qs[0]], qs[2]);
X(qs[0]); X(qs[1]);
// 010 -> w 010
X(qs[2]); X(qs[0]);
RAll1(4.0*ArcSin(0.5), qs);
X(qs[2]); X(qs[0]);
// 001 -> w^2 001
X(qs[2]); X(qs[1]);
RAll1(8.0*ArcSin(0.5), qs);
X(qs[2]); X(qs[1]);
DumpMachine("dump1.txt");
// w^2 001 -> 001
X(qs[2]); X(qs[1]);
RAll1(-8.0*ArcSin(0.5), qs);
X(qs[2]); X(qs[1]);
// w 010 -> 010
X(qs[2]); X(qs[0]);
RAll1(-4.0*ArcSin(0.5), qs);
X(qs[2]); X(qs[0]);
// 100 -> 000
X(qs[0]); X(qs[1]);
Controlled X ([qs[1],qs[0]], qs[2]);
X(qs[0]); X(qs[1]);
X(qs[0]);
Controlled Ry ([qs[0]], (-2.0*theta2, qs[1]) );
X(qs[0]);
Ry(-2.0*theta, qs[0]);
let M0 = M(qs[2]);
let M1 = M(qs[1]);
let M2 = M(qs[0]);
DumpMachine("dump2.txt");
if (M0==Zero && M1==Zero && M2==Zero) {
return 1;
}
else {return 0;}
}
//adjoint invert;
//controlled distribute;
// controlled adjoint distribute;
}
前半が順変換、後半が逆変換で元に戻ることを確認している。解答としては後半部分を提出すれば良い。(逆変換はadjoint autoの機能を使えば自動で複素共役転置行列を作ってくれるっぽい:参考 (問題文も参照))
B2. Not A, not B or not C?
次の3状態$\ket A, \ket B, \ket C$の識別を行う。完全に特定する必要はなく、「どれかでない」ことを識別できれば良い。
\begin{eqnarray}
\ket{A} &=& \frac{1}{\sqrt2} ( \ket0 + \ket1 ) \\
\ket{B} &=& \frac{1}{\sqrt2} ( \ket0 + \omega \ket1 ) \\
\ket{C} &=& \frac{1}{\sqrt2} ( \ket0 + \omega^2 \ket1 ) \\
\omega &=& e^{\frac{2\pi i}{3}}
\end{eqnarray}
解けず。
C問題
与えられる量子状態$\ket{x}$がbit列としての条件(交互にならんでいるかなど)を満たす場合にTarget register$\ket{y}$をFlipさせる。
\begin{eqnarray}
\ket{x} \otimes \ket{y} &\rightarrow& \ket{x} \otimes \ket{f(x) \oplus y}
\end{eqnarray}
関数$f(x)$は$x$が条件を満たすとき1、満たさない時0になる関数で、$\oplus$はXOR演算子である。
やることは簡単。実装が大変がち。
C1. Alternating bits oracle
$\ket{x}$が交互に並んでいるか判定する。パターンとして$\ket{101010...}$か$\ket{010101...}$のどちらか。1,3,5...番目をFlipさせると$\ket{00000...}$か$\ket{11111...}$になる。後はControlled X$(x,y)$とControlled X$(\bar x,y)$を作用させればよい。最後にここまで$\ket{x}$に行った操作を逆に辿って初期状態の$\ket{x}$に戻す(測定を用いないユニタリ変換しか行っていないので必ず逆変換が可能)。
C2. "Is the bit string periodic?" oracle
$\ket{x}$が周期的に並んでいるか判定する。bit数($2\le N \le 7$)に対して有効な周期は限られるので、全列挙して場合分け・判定をした。
C3. "Is the number of ones divisible by 3?" oracle
$\ket{x}$に含まれる$\ket 1$の個数が3の倍数であるかどうかを判定する。counter用レジスターを用意して、$x[0]$から順に$\ket 1$が含まれている場合+1をする。倍数判定をして$\ket y$をFlipさせる。
関連する話題として、$\ket{\psi_0}= F \ket{00...0}$が$\ket 1$をいくつ含んでいるかを数えるアルゴリズムがある。$\ket 1$が0または1つしか無いと分かっている場合、古典的なデータ列であればbit数Nに対して線形$O(N)$のデータアクセス回数が必要であるが、量子的なinputオラクル関数Fとして与えられている場合$O(\sqrt N)$回のクエリ回数で判定が出来る(Groverアルゴリズム)。より一般に$n$個の$\ket 1$が含まれている場合に、その個数$n$を確率的に数える(古典より)効率的な量子アルゴリズムも存在する(「Quantum Counting」,「Groverの量子探索アルゴリズムの応用」)。
D問題
指定されたユニタリ行列の形を満たすような量子操作を見つける。例えば2qubitに作用する$4\times4$ユニタリ行列が次のように指定される。ここで、$X$はnon-zeroの任意の値であることを表す($X$は別に同じ値で無くて良い)。
\begin{eqnarray}
\begin{pmatrix}
X & X & 0 & 0 \\
X & X & 0 & 0 \\
0 & 0 & X & 0 \\
0 & 0 & 0 & X \\
\end{pmatrix}
\begin{pmatrix}
\ket{00} \\
\ket{01} \\
\ket{10} \\
\ket{11} \\
\end{pmatrix}
\end{eqnarray}
これは要素ごとに分解すると、
\begin{eqnarray}
\ket{00} &\rightarrow& \ket{00} + \ket{01} \\
\ket{01} &\rightarrow& \ket{00} + \ket{01} \\
\ket{10} &\rightarrow& \ket{10} \\
\ket{11} &\rightarrow& \ket{11} \\
\end{eqnarray}
よって、第1qubitが$\ket 0$の場合のみ第2qubitにHゲートを作用させる操作で実現出来る。より具体的な例はWarmup-Roundの問題と解答を参照するとよい。
ポイントとして、測定が禁止されているという制約上、指定される行列は全てユニタリである(あらゆるユニタリ変換は量子操作で実現可能(EMAN物理学 2入力量子ゲート)であるので、ユニタリ行列ならば必ず解が存在する)。ユニタリ行列の特徴として、
・列ベクトルは正規直交基底である
・行ベクトルは正規直交基底である
という性質があるので、$X$の値をほぼ推測することが可能である。例えば、上の例であれば右下$2\times2$ブロックにある2つの$X$は正規直交基底となるためには必ず1となっているはずである。
問題文で指定されている通り、実装・デバッグは「DumpUnitary tool」を使用して行った。DumpUnitaryフォルダをコピーして来て、dotnet runでそのまま動かせる。「Operations.qs」ファイルの// Apply the operationの行の下に望みの量子操作を書き込めんで実行すれば、「DumpUnitaryPattern.txt」ファイルに行列の形が出力される。コンテスト提出も書き込んだ量子操作部分をほぼそのまま提出する形でいける。
D1. Block diagonal matrix
省略
D2. Pattern of increasing blocks
\begin{eqnarray}
\begin{pmatrix}
X & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & X & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & X & X & 0 & 0 & 0 & 0 \\
0 & 0 & X & X & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & X & X & X & X \\
0 & 0 & 0 & 0 & X & X & X & X \\
0 & 0 & 0 & 0 & X & X & X & X \\
0 & 0 & 0 & 0 & X & X & X & X \\
\end{pmatrix}
\begin{pmatrix}
\ket{000} \\
\ket{001} \\
\ket{010} \\
\ket{011} \\
\ket{100} \\
\ket{101} \\
\ket{110} \\
\ket{111} \\
\end{pmatrix}
\end{eqnarray}
$\ket{00?}$の場合何もしない(恒等演算,$I$)、$\ket{01?}$の場合$\ket?$にControlled Hゲート、$\ket{1??}$の場合2つの$\ket?$qubitにControlled Hゲートでおしまい。
// Apply the operation
for (i in 1..N-1) {
for (j in i+1..N-1) {X(qs[j]);}
let array = IntArrayFromRange(i..1..N-1);
for (j in 0..i-1) {Controlled H (Subarray(array, qs), qs[j]);}
for (j in i+1..N-1) {X(qs[j]);}
}
D3. X-wing fighter
\begin{eqnarray}
\begin{pmatrix}
X & 0 & 0 & 0 & 0 & 0 & 0 & X \\
0 & X & 0 & 0 & 0 & 0 & X & 0 \\
0 & 0 & X & 0 & 0 & X & 0 & 0 \\
0 & 0 & 0 & X & X & 0 & 0 & 0 \\
0 & 0 & 0 & X & X & 0 & 0 & 0 \\
0 & 0 & X & 0 & 0 & X & 0 & 0 \\
0 & X & 0 & 0 & 0 & 0 & X & 0 \\
X & 0 & 0 & 0 & 0 & 0 & 0 & X \\
\end{pmatrix}
\begin{pmatrix}
\ket{000} \\
\ket{001} \\
\ket{010} \\
\ket{011} \\
\ket{100} \\
\ket{101} \\
\ket{110} \\
\ket{111} \\
\end{pmatrix}
\end{eqnarray}
2つの対角線はそれぞれ恒等操作とSWAP(全Flip)操作に対応している。よって次のように、自身と全Flipしたものに分裂させればよい。
\begin{eqnarray}
\ket{000} &\rightarrow& \frac{1}{\sqrt2} \left( \ket{000} + \ket{111} \right)
\end{eqnarray}
具体的な手順としては、Hゲートで2つに分割、$\ket{1}$の方の状態をControlled Xで残りqubitをFlipさせる。
\begin{eqnarray}
\ket{000} &\rightarrow& \frac{1}{\sqrt2} \left( \ket{000} + \ket{001} \right) \\
&\rightarrow& \frac{1}{\sqrt2} \left( \ket{000} + \ket{111} \right) \\
\end{eqnarray}
これだと$\ket{001}$の状態と混成してうまくいかない。第三qubitが$\ket1$のものは予めControlled Xで他qubitを全Flipさせておけばうまくいく。
\begin{eqnarray}
\ket{001} &\rightarrow& \ket{111} \\
&\rightarrow& \frac{1}{\sqrt2} \left( \ket{110} + \ket{111} \right) \\
&\rightarrow& \frac{1}{\sqrt2} \left( \ket{110} + \ket{001} \right) \\
\end{eqnarray}
// Apply the operation
for (i in 1..N-1) {Controlled X([qs[0]], qs[i]);}
H(qs[0]);
for (i in 1..N-1) {Controlled X([qs[0]], qs[i]);}
D4. TIE fighter
省略
D5. Creeper
省略
D6. Hessenberg matrix
対角線のもう一本下(subdiagonal)まで要素が存在し、それより下が0となるような行列。$N=2,3,4$の場合についてやる。かなり難しい上に実装も一番大変だった。もう二度とやりたくないレベル(効率的な実装方法が他にありそう)。
\begin{eqnarray}
\begin{pmatrix}
X & X & X & X \\
X & X & X & X \\
0 & X & X & X \\
0 & 0 & X & X \\
\end{pmatrix}
\begin{pmatrix}
\ket{00} \\
\ket{01} \\
\ket{10} \\
\ket{11} \\
\end{pmatrix}
\end{eqnarray}
一番簡単な$N=2$の場合の実装を考える。普通に試行錯誤しても一向に解決方法が見つからなかった。趣向を変えて作成するべきユニタリ行列の成分が具体的にどうなるか考える。すると、列ベクトルが正規直交基底になるのは次のような場合があることがわかる。
\begin{eqnarray}
\begin{pmatrix}
\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{4}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} \\
-\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{4}} & \frac{1}{\sqrt{8}} & \frac{1}{\sqrt{8}} \\
0 & -\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{4}} & \frac{1}{\sqrt{4}} \\
0 & 0 & -\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\
\end{pmatrix}
\end{eqnarray}
これを実装すればよい。ユニタリ行列の下の要素から順に実現していく。まず4行目の変換
\begin{eqnarray}
\ket{11} &\rightarrow& \frac{1}{\sqrt2} \ket{10} + \frac{1}{\sqrt2} \ket{11} \\
\ket{10} &\rightarrow& \frac{1}{\sqrt2} \ket{10} - \frac{1}{\sqrt2} \ket{11}
\end{eqnarray}
はこれだけ。次に
\begin{eqnarray}
\ket{10} &\rightarrow& \frac{1}{\sqrt2} \ket{01} + \frac{1}{\sqrt2} \ket{10} \\
\ket{01} &\rightarrow& \frac{1}{\sqrt2} \ket{01} - \frac{1}{\sqrt2} \ket{10}
\end{eqnarray}
を行うと、全体としては、
\begin{eqnarray}
\ket{11} &\rightarrow& \frac{1}{\sqrt4} \ket{01} + \frac{1}{\sqrt4} \ket{10} + \frac{1}{\sqrt2} \ket{11} \\
\ket{10} &\rightarrow& \frac{1}{\sqrt4} \ket{01} + \frac{1}{\sqrt4} \ket{10} - \frac{1}{\sqrt2} \ket{11} \\
\ket{01} &\rightarrow& \frac{1}{\sqrt2} \ket{01} - \frac{1}{\sqrt2} \ket{10} \\
\ket{00} &\rightarrow& \ket{00}
\end{eqnarray}
となっており、これで3,4行目が完成した。最後に残った$\ket{01},\ket{00}$を分解する操作
\begin{eqnarray}
\ket{01} &\rightarrow& \frac{1}{\sqrt2} \ket{00} + \frac{1}{\sqrt2} \ket{01} \\
\ket{00} &\rightarrow& \frac{1}{\sqrt2} \ket{00} - \frac{1}{\sqrt2} \ket{01}
\end{eqnarray}
を加えれば、
\begin{eqnarray}
\ket{11} &\rightarrow& \frac{1}{\sqrt8} \ket{00} + \frac{1}{\sqrt8} \ket{01} + \frac{1}{\sqrt4} \ket{10} + \frac{1}{\sqrt2} \ket{11} \\
\ket{10} &\rightarrow& \frac{1}{\sqrt8} \ket{00} + \frac{1}{\sqrt4} \ket{01} + \frac{1}{\sqrt4} \ket{10} - \frac{1}{\sqrt2} \ket{11} \\
\ket{01} &\rightarrow& \frac{1}{\sqrt4} \ket{00} + \frac{1}{\sqrt2} \ket{01} - \frac{1}{\sqrt2} \ket{10} \\
\ket{00} &\rightarrow& \frac{1}{\sqrt2} \ket{00} - \frac{1}{\sqrt2} \ket{01}
\end{eqnarray}
となり問題のユニタリ変換が実現している。$N=3,4$の場合も同様に出来るが、変換を辿って実装するのがクソ大変になっていく。
// Apply the operation
//let N = Length(qs);
if (N==2) {
Controlled H ([qs[1]], qs[0]);
Controlled X ([qs[1]], qs[0]);
Controlled H ([qs[0]], qs[1]);
X(qs[1]);
Controlled H ([qs[1]], qs[0]);
X(qs[1]);
X(qs[0]);
}
elif (N==3) {
// 1
Controlled H ([qs[2],qs[1]], qs[0]);
// 2
Controlled X ([qs[2],qs[1]], qs[0]);
Controlled H ([qs[2],qs[0]], qs[1]);
Controlled X ([qs[2],qs[1]], qs[0]);
// 3
X(qs[1]);
Controlled H ([qs[2],qs[1]], qs[0]);
X(qs[1]);
// 4
X(qs[1]); X(qs[0]);
Controlled X ([qs[1],qs[0]], qs[2]);
X(qs[1]); X(qs[0]);
X(qs[2]); X(qs[0]);
Controlled X ([qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
X(qs[2]);
Controlled H ([qs[2],qs[1]], qs[0]);
X(qs[2]);
X(qs[2]); X(qs[0]);
Controlled X ([qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
X(qs[1]); X(qs[0]);
Controlled X ([qs[1],qs[0]], qs[2]);
X(qs[1]); X(qs[0]);
// 5
X(qs[2]);
Controlled H ([qs[2],qs[1]], qs[0]);
X(qs[2]);
// 6
X(qs[2]); X(qs[0]);
Controlled X ([qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
X(qs[2]); X(qs[1]);
Controlled H ([qs[2],qs[1]], qs[0]);
X(qs[2]); X(qs[1]);
X(qs[2]); X(qs[0]);
Controlled X ([qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
// 7
X(qs[2]); X(qs[1]);
Controlled H ([qs[2],qs[1]], qs[0]);
X(qs[2]); X(qs[1]);
}
elif (N==4) {
// 1
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
// 2
Controlled X ([qs[3],qs[2],qs[1]], qs[0]);
Controlled H ([qs[3],qs[2],qs[0]], qs[1]);
Controlled X ([qs[3],qs[2],qs[1]], qs[0]);
// 3
X(qs[1]);
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
X(qs[1]);
// 4
X(qs[1]); X(qs[0]);
Controlled X ([qs[3],qs[1],qs[0]], qs[2]);
X(qs[1]); X(qs[0]);
X(qs[2]); X(qs[0]);
Controlled X ([qs[3],qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
X(qs[2]);
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
X(qs[2]);
X(qs[2]); X(qs[0]);
Controlled X ([qs[3],qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
X(qs[1]); X(qs[0]);
Controlled X ([qs[3],qs[1],qs[0]], qs[2]);
X(qs[1]); X(qs[0]);
// 5
X(qs[2]);
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
X(qs[2]);
// 6
X(qs[2]); X(qs[0]);
Controlled X ([qs[3],qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
X(qs[2]); X(qs[1]);
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
X(qs[2]); X(qs[1]);
X(qs[2]); X(qs[0]);
Controlled X ([qs[3],qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
// 7
X(qs[2]); X(qs[1]);
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
X(qs[2]); X(qs[1]);
// 8 !!
X(qs[2]); X(qs[1]);
Controlled X ([qs[3],qs[2],qs[1]],qs[0]);
X(qs[2]); X(qs[1]);
X(qs[2]);
Controlled X ([qs[3],qs[2],qs[0]],qs[1]);
X(qs[2]);
Controlled X ([qs[3],qs[1],qs[0]],qs[2]);
Controlled H ([qs[2],qs[1],qs[0]],qs[3]);
Controlled X ([qs[3],qs[1],qs[0]],qs[2]);
X(qs[2]);
Controlled X ([qs[3],qs[2],qs[0]],qs[1]);
X(qs[2]);
X(qs[2]); X(qs[1]);
Controlled X ([qs[3],qs[2],qs[1]],qs[0]);
X(qs[2]); X(qs[1]);
X(qs[3]);
// 1
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
// 2
Controlled X ([qs[3],qs[2],qs[1]], qs[0]);
Controlled H ([qs[3],qs[2],qs[0]], qs[1]);
Controlled X ([qs[3],qs[2],qs[1]], qs[0]);
// 3
X(qs[1]);
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
X(qs[1]);
// 4
X(qs[1]); X(qs[0]);
Controlled X ([qs[3],qs[1],qs[0]], qs[2]);
X(qs[1]); X(qs[0]);
X(qs[2]); X(qs[0]);
Controlled X ([qs[3],qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
X(qs[2]);
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
X(qs[2]);
X(qs[2]); X(qs[0]);
Controlled X ([qs[3],qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
X(qs[1]); X(qs[0]);
Controlled X ([qs[3],qs[1],qs[0]], qs[2]);
X(qs[1]); X(qs[0]);
// 5
X(qs[2]);
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
X(qs[2]);
// 6
X(qs[2]); X(qs[0]);
Controlled X ([qs[3],qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
X(qs[2]); X(qs[1]);
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
X(qs[2]); X(qs[1]);
X(qs[2]); X(qs[0]);
Controlled X ([qs[3],qs[2],qs[0]], qs[1]);
X(qs[2]); X(qs[0]);
// 7
X(qs[2]); X(qs[1]);
Controlled H ([qs[3],qs[2],qs[1]], qs[0]);
X(qs[2]); X(qs[1]);
X(qs[3]);
}