GCC Code Coverage Report

Directory: ./
File: Circuit/ThreeQubitConversion.cpp
Date: 2022-10-15 05:10:18
Warnings: 1 unchecked decisions!
Exec Total Coverage
Lines: 340 350 97.1%
Functions: 16 16 100.0%
Branches: 596 1155 51.6%
Decisions: 15 18 83.3%
Line Branch Decision Exec Source
1 // Copyright 2019-2022 Cambridge Quantum Computing
2 //
3 // Licensed under the Apache License, Version 2.0 (the "License");
4 // you may not use this file except in compliance with the License.
5 // You may obtain a copy of the License at
6 //
7 // http://www.apache.org/licenses/LICENSE-2.0
8 //
9 // Unless required by applicable law or agreed to in writing, software
10 // distributed under the License is distributed on an "AS IS" BASIS,
11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 // See the License for the specific language governing permissions and
13 // limitations under the License.
14
15 #include "ThreeQubitConversion.hpp"
16
17 #include <array>
18 #include <cmath>
19 #include <complex>
20 #include <optional>
21 #include <stdexcept>
22 #include <tkassert/Assert.hpp>
23
24 #include "CircPool.hpp"
25 #include "CircUtils.hpp"
26 #include "Circuit.hpp"
27 #include "Circuit/Command.hpp"
28 #include "Gate/GatePtr.hpp"
29 #include "Gate/Rotation.hpp"
30 #include "OpType/OpType.hpp"
31 #include "Utils/Constants.hpp"
32 #include "Utils/CosSinDecomposition.hpp"
33 #include "Utils/EigenConfig.hpp"
34 #include "Utils/MatrixAnalysis.hpp"
35 #include "Utils/UnitID.hpp"
36
37 namespace tket {
38
39 // Return a 3-qubit circuit implementing the unitary
40 // [ D ]
41 // [ D* ]
42 // using 4 Rz and 4 CX operations, where D is a 4x4 diagonal unitary matrix.
43 // The circuit consists of Rz operations on qubit 0 and CX operations with
44 // target qubit 0.
45 732 static Circuit two_qubit_diag_adjoint_plex(const Eigen::Matrix4cd &D) {
46 // Compute angles a_i such that D_ii = e^{-i pi/2 a_i}.
47 static const double f = -2 / PI;
48 732 double a0 = f * std::arg(D(0, 0));
49 732 double a1 = f * std::arg(D(1, 1));
50 732 double a2 = f * std::arg(D(2, 2));
51 732 double a3 = f * std::arg(D(3, 3));
52 732 double t0 = (a0 + a1 + a2 + a3) / 4;
53 732 double t1 = (a0 + a1 - a2 - a3) / 4;
54 732 double t2 = (a0 - a1 - a2 + a3) / 4;
55 732 double t3 = (a0 - a1 + a2 - a3) / 4;
56
1/2
✓ Branch 2 taken 732 times.
✗ Branch 3 not taken.
 732 Circuit circ(3);
57
3/6
✓ Branch 3 taken 732 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 732 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 732 times.
✗ Branch 10 not taken.
 732 circ.add_op<unsigned>(OpType::Rz, t0, {0});
58
2/4
✓ Branch 3 taken 732 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 732 times.
✗ Branch 7 not taken.
 732 circ.add_op<unsigned>(OpType::CX, {1, 0});
59
3/6
✓ Branch 3 taken 732 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 732 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 732 times.
✗ Branch 10 not taken.
 732 circ.add_op<unsigned>(OpType::Rz, t1, {0});
60
2/4
✓ Branch 3 taken 732 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 732 times.
✗ Branch 7 not taken.
 732 circ.add_op<unsigned>(OpType::CX, {2, 0});
61
3/6
✓ Branch 3 taken 732 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 732 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 732 times.
✗ Branch 10 not taken.
 732 circ.add_op<unsigned>(OpType::Rz, t2, {0});
62
2/4
✓ Branch 3 taken 732 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 732 times.
✗ Branch 7 not taken.
 732 circ.add_op<unsigned>(OpType::CX, {1, 0});
63
3/6
✓ Branch 3 taken 732 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 732 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 732 times.
✗ Branch 10 not taken.
 732 circ.add_op<unsigned>(OpType::Rz, t3, {0});
64
2/4
✓ Branch 3 taken 732 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 732 times.
✗ Branch 7 not taken.
 732 circ.add_op<unsigned>(OpType::CX, {2, 0});
65 732 return circ;
66 }
67
68 // Return a 3-qubit circuit implementing the unitary
69 // [ D ]
70 // [ D* ]
71 // using 1-qubit and 4 TK2 operations, where D is a 4x4 diagonal unitary matrix.
72 // The circuit consists of Rz operations on qubit 0 and TK2 operations.
73 86 static Circuit two_qubit_diag_adjoint_plex_tk(const Eigen::Matrix4cd &D) {
74 // Compute angles a_i such that D_ii = e^{-i pi/2 a_i}.
75 static const double f = -2 / PI;
76 86 double a0 = f * std::arg(D(0, 0));
77 86 double a1 = f * std::arg(D(1, 1));
78 86 double a2 = f * std::arg(D(2, 2));
79 86 double a3 = f * std::arg(D(3, 3));
80 86 double t0 = (a0 + a1 + a2 + a3) / 4;
81 86 double t1 = (a0 + a1 - a2 - a3) / 4;
82 86 double t2 = (a0 - a1 - a2 + a3) / 4;
83 86 double t3 = (a0 - a1 + a2 - a3) / 4;
84
1/2
✓ Branch 2 taken 86 times.
✗ Branch 3 not taken.
 86 Circuit circ(3);
85
3/6
✓ Branch 3 taken 86 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 86 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 86 times.
✗ Branch 10 not taken.
 86 circ.add_op<unsigned>(OpType::Rz, t0, {0});
86
5/10
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 86 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 172 times.
✓ Branch 15 taken 86 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 344 circ.append_with_map(
87
4/8
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 86 times.
✗ Branch 11 not taken.
 258 CircPool::CX_using_TK2(), {{Qubit(0), Qubit(1)}, {Qubit(1), Qubit(0)}});
88
3/6
✓ Branch 3 taken 86 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 86 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 86 times.
✗ Branch 10 not taken.
 86 circ.add_op<unsigned>(OpType::Rz, t1, {0});
89
5/10
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 86 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 172 times.
✓ Branch 15 taken 86 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 344 circ.append_with_map(
90
4/8
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 86 times.
✗ Branch 11 not taken.
 258 CircPool::CX_using_TK2(), {{Qubit(0), Qubit(2)}, {Qubit(1), Qubit(0)}});
91
3/6
✓ Branch 3 taken 86 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 86 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 86 times.
✗ Branch 10 not taken.
 86 circ.add_op<unsigned>(OpType::Rz, t2, {0});
92
5/10
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 86 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 172 times.
✓ Branch 15 taken 86 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 344 circ.append_with_map(
93
4/8
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 86 times.
✗ Branch 11 not taken.
 258 CircPool::CX_using_TK2(), {{Qubit(0), Qubit(1)}, {Qubit(1), Qubit(0)}});
94
3/6
✓ Branch 3 taken 86 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 86 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 86 times.
✗ Branch 10 not taken.
 86 circ.add_op<unsigned>(OpType::Rz, t3, {0});
95
5/10
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 86 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 172 times.
✓ Branch 15 taken 86 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 344 circ.append_with_map(
96
4/8
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 86 times.
✗ Branch 11 not taken.
 258 CircPool::CX_using_TK2(), {{Qubit(0), Qubit(2)}, {Qubit(1), Qubit(0)}});
97 86 return circ;
98 }
99
100 122 static const std::vector<Eigen::Matrix4cd> &get_conj_unitaries() {
101 1 static const std::vector<Eigen::Matrix4cd> vec = ([]() {
102
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Circuit c1(2);
103
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Circuit c2(2);
104
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1 c2.add_op<unsigned>(OpType::SWAP, {0, 1});
105
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Circuit c3(2);
106
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1 c3.add_op<unsigned>(OpType::CX, {0, 1});
107
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Circuit c4(2);
108
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1 c4.add_op<unsigned>(OpType::CX, {1, 0});
109
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Circuit c5(2);
110
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1 c5.add_op<unsigned>(OpType::CX, {0, 1});
111
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1 c5.add_op<unsigned>(OpType::CX, {1, 0});
112
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 Circuit c6(2);
113
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1 c6.add_op<unsigned>(OpType::CX, {1, 0});
114
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1 c6.add_op<unsigned>(OpType::CX, {0, 1});
115 return std::vector<Eigen::Matrix4cd>{
116 get_matrix_from_2qb_circ(c1), get_matrix_from_2qb_circ(c2),
117 get_matrix_from_2qb_circ(c3), get_matrix_from_2qb_circ(c4),
118
7/14
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
 2 get_matrix_from_2qb_circ(c5), get_matrix_from_2qb_circ(c6)};
119
4/8
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 121 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 123 })();
120 122 return vec;
121 }
122
123 // Return a 3-qubit circuit and a unit complex number z which together implement
124 // the unitary
125 // [ U0 ]
126 // [ U1 ]
127 // where U0 and U1 are 4x4 unitaries.
128 //
129 // The unitary is implemented by the circuit followed by the diagonal operator
130 // diag(z, z*, z*, z) on qubits 1 and 2.
131 //
132 // If `extract_final_diagonal` is false, then z=1 and the circuit implements the
133 // unitary exactly, using 9 CX gates.
134 //
135 // If `extract_final_diagonal` is true, then the circuit uses 8 CX gates.
136 122 static std::pair<Circuit, Complex> two_qubit_plex(
137 const Eigen::Matrix4cd &U0, const Eigen::Matrix4cd &U1,
138 bool extract_final_diagonal) {
139 // 1. Decompose U0 U1* as L T L* where L and T are unitary and T is diagonal.
140
3/6
✓ Branch 1 taken 122 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 122 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 122 times.
✗ Branch 8 not taken.
 122 Eigen::Matrix4cd U = U0 * U1.adjoint();
141
1/2
✓ Branch 1 taken 122 times.
✗ Branch 2 not taken.
 122 Eigen::ComplexSchur<Eigen::Matrix4cd> schur(U);
142
1/2
✓ Branch 2 taken 122 times.
✗ Branch 3 not taken.
 122 Eigen::Matrix4cd L = schur.matrixU();
143
1/2
✓ Branch 2 taken 122 times.
✗ Branch 3 not taken.
 122 Eigen::Matrix4cd T = schur.matrixT();
144 // By construction T is unitary and upper-triangular, hence diagonal.
145 TKET_ASSERT(T.isDiagonal());
146 // 2. Let D = sqrt(T)
147
2/4
✓ Branch 1 taken 122 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 122 times.
✗ Branch 5 not taken.
 122 Eigen::Matrix4cd D = Eigen::Matrix4cd::Zero();
148
2/2
✓ Branch 0 taken 488 times.
✓ Branch 1 taken 122 times.
2/2
✓ Decision 'true' taken 488 times.
✓ Decision 'false' taken 122 times.
 610 for (unsigned i = 0; i < 4; i++) {
149
2/4
✓ Branch 1 taken 488 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 488 times.
✗ Branch 6 not taken.
 488 D(i, i) = sqrt(T(i, i));
150 }
151 // 3. Compute R such that U0 = L D R and U1 = L D* R.
152
4/8
✓ Branch 1 taken 122 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 122 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 122 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 122 times.
✗ Branch 11 not taken.
 122 Eigen::Matrix4cd R = D * L.adjoint() * U1;
153
154 // We try conjugating the L and R circuits to see if we can reduce CX count.
155 122 std::optional<Circuit> best_circ;
156 122 std::optional<Complex> best_z0;
157 122 std::optional<unsigned> best_n_cx;
158
3/4
✓ Branch 1 taken 122 times.
✗ Branch 2 not taken.
✓ Branch 8 taken 732 times.
✓ Branch 9 taken 122 times.
0/1
? Decision couldn't be analyzed.
 854 for (const Eigen::Matrix4cd &u_conj : get_conj_unitaries()) {
159 // 4. Decompose R into a 2-CX circuit followed by a diagonal.
160
1/2
✓ Branch 1 taken 732 times.
✗ Branch 2 not taken.
 732 auto u_conj_adj = u_conj.adjoint();
161
3/6
✓ Branch 1 taken 732 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 732 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 732 times.
✗ Branch 8 not taken.
 732 auto [R_circ, w0] = decompose_2cx_DV(u_conj * R);
162
1/2
✓ Branch 3 taken 732 times.
✗ Branch 4 not taken.
 732 Eigen::Vector4cd w = {w0, std::conj(w0), std::conj(w0), w0};
163
1/2
✓ Branch 1 taken 732 times.
✗ Branch 2 not taken.
 732 auto wD = w.asDiagonal();
164
165 // 5. Construct the circuit.
166 // The diagonal from R's decomposition commutes forward through the controls
167 // on qubits 1 and 2, so can be absorbed into L.
168
1/2
✓ Branch 2 taken 732 times.
✗ Branch 3 not taken.
 732 Circuit circ(3);
169
5/10
✓ Branch 1 taken 732 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 732 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 732 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 732 times.
✗ Branch 12 not taken.
✓ Branch 16 taken 732 times.
✗ Branch 17 not taken.
 3660 unit_map_t qm{{Qubit(0), Qubit(1)}, {Qubit(1), Qubit(2)}};
170
1/2
✓ Branch 1 taken 732 times.
✗ Branch 2 not taken.
 732 circ.append_with_map(R_circ, qm);
171
5/10
✓ Branch 1 taken 732 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 732 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 732 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 732 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 732 times.
✗ Branch 14 not taken.
 732 circ.append(two_qubit_diag_adjoint_plex(u_conj * D * u_conj_adj));
172 732 Complex z0;
173
1/2
✓ Branch 1 taken 732 times.
✗ Branch 2 not taken.
 732 Circuit L_circ;
174
2/2
✓ Branch 0 taken 366 times.
✓ Branch 1 taken 366 times.
2/2
✓ Decision 'true' taken 366 times.
✓ Decision 'false' taken 366 times.
 732 if (extract_final_diagonal) {
175
5/10
✓ Branch 1 taken 366 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 366 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 366 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 366 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 366 times.
✗ Branch 15 not taken.
 366 std::tie(L_circ, z0) = decompose_2cx_DV(L * u_conj_adj * wD);
176 } else {
177
5/10
✓ Branch 1 taken 366 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 366 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 366 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 366 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 366 times.
✗ Branch 14 not taken.
 366 L_circ = two_qubit_canonical(L * u_conj_adj * wD, OpType::CX);
178 366 z0 = 1.;
179 }
180
1/2
✓ Branch 1 taken 732 times.
✗ Branch 2 not taken.
 732 circ.append_with_map(L_circ, qm);
181
182
1/2
✓ Branch 1 taken 732 times.
✗ Branch 2 not taken.
 732 unsigned n_cx = circ.count_gates(OpType::CX);
183
6/6
✓ Branch 1 taken 610 times.
✓ Branch 2 taken 122 times.
✓ Branch 6 taken 21 times.
✓ Branch 7 taken 589 times.
✓ Branch 8 taken 143 times.
✓ Branch 9 taken 589 times.
2/2
✓ Decision 'true' taken 143 times.
✓ Decision 'false' taken 589 times.
 732 if (!best_circ || best_n_cx > n_cx) {
184
1/2
✓ Branch 1 taken 143 times.
✗ Branch 2 not taken.
 143 best_circ = circ;
185 143 best_z0 = z0;
186 143 best_n_cx = n_cx;
187 }
188 732 }
189 TKET_ASSERT(best_circ);
190
1/2
✓ Branch 3 taken 122 times.
✗ Branch 4 not taken.
 244 return {*best_circ, *best_z0};
191 122 }
192
193 // Return a 3-qubit circuit which implements the unitary
194 // [ U0 ]
195 // [ U1 ]
196 // where U0 and U1 are 4x4 unitaries.
197 86 static Circuit two_qubit_plex_tk(
198 const Eigen::Matrix4cd &U0, const Eigen::Matrix4cd &U1) {
199 // 1. Decompose U0 U1* as L T L* where L and T are unitary and T is diagonal.
200
3/6
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
 86 Eigen::Matrix4cd U = U0 * U1.adjoint();
201
1/2
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
 86 Eigen::ComplexSchur<Eigen::Matrix4cd> schur(U);
202
1/2
✓ Branch 2 taken 86 times.
✗ Branch 3 not taken.
 86 Eigen::Matrix4cd L = schur.matrixU();
203
1/2
✓ Branch 2 taken 86 times.
✗ Branch 3 not taken.
 86 Eigen::Matrix4cd T = schur.matrixT();
204 // By construction T is unitary and upper-triangular, hence diagonal.
205 TKET_ASSERT(T.isDiagonal());
206 // 2. Let D = sqrt(T)
207
2/4
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
 86 Eigen::Matrix4cd D = Eigen::Matrix4cd::Zero();
208
2/2
✓ Branch 0 taken 344 times.
✓ Branch 1 taken 86 times.
2/2
✓ Decision 'true' taken 344 times.
✓ Decision 'false' taken 86 times.
 430 for (unsigned i = 0; i < 4; i++) {
209
2/4
✓ Branch 1 taken 344 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 344 times.
✗ Branch 6 not taken.
 344 D(i, i) = sqrt(T(i, i));
210 }
211 // 3. Compute R such that U0 = L D R and U1 = L D* R.
212
4/8
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 86 times.
✗ Branch 11 not taken.
 86 Eigen::Matrix4cd R = D * L.adjoint() * U1;
213 // 4. Construct the circuit.
214
1/2
✓ Branch 2 taken 86 times.
✗ Branch 3 not taken.
 86 Circuit circ(3);
215
5/10
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 86 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 86 times.
✗ Branch 12 not taken.
✓ Branch 16 taken 86 times.
✗ Branch 17 not taken.
 430 unit_map_t qm{{Qubit(0), Qubit(1)}, {Qubit(1), Qubit(2)}};
216
2/4
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
 86 circ.append_with_map(two_qubit_canonical(R), qm);
217
2/4
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
 86 circ.append(two_qubit_diag_adjoint_plex_tk(D));
218
2/4
✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
 86 circ.append_with_map(two_qubit_canonical(L), qm);
219 172 return circ;
220 86 }
221
222 // Return a 3-qubit circuit implementing the unitary
223 // [ C -S ]
224 // [ DS DC ]
225 // using TK2 and 1-qubit operations, where C and S are 4x4 real diagonal
226 // matrices, C^2 + S^2 = I, and D = diag(1,-1,1,-1).
227 //
228 // Note that
229 // [ C -S ] = U [ C -S ]
230 // [ DS DC ] [ S C ]
231 // where U represents a CZ on qubits 0 and 2. We convert the CZ operations to CX
232 // by adding Hadamards and simplifying H Ry(t) H to Ry(-t).
233 61 static Circuit two_qubit_modified_cossin_circ(
234 const Eigen::Matrix4d &C, const Eigen::Matrix4d &S) {
235 // Compute angles a_i such that C_ii = cos(pi/2 a_i) and S_ii = sin(pi/2 a_i).
236 static const double f = 2 / PI;
237 61 double a0 = f * atan2(S(0, 0), C(0, 0));
238 61 double a1 = f * atan2(S(1, 1), C(1, 1));
239 61 double a2 = f * atan2(S(2, 2), C(2, 2));
240 61 double a3 = f * atan2(S(3, 3), C(3, 3));
241 61 double t0 = (a0 + a1 + a2 + a3) / 4;
242 61 double t1 = (a0 + a1 - a2 - a3) / 4;
243 61 double t2 = (a0 - a1 - a2 + a3) / 4;
244 61 double t3 = (a0 - a1 + a2 - a3) / 4;
245
1/2
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
 61 Circuit circ(3);
246
3/6
✓ Branch 3 taken 61 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 61 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 61 times.
✗ Branch 10 not taken.
 61 circ.add_op<unsigned>(OpType::Ry, t0, {0});
247
2/4
✓ Branch 3 taken 61 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 61 times.
✗ Branch 7 not taken.
 61 circ.add_op<unsigned>(OpType::H, {0});
248
2/4
✓ Branch 3 taken 61 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 61 times.
✗ Branch 7 not taken.
 61 circ.add_op<unsigned>(OpType::CX, {1, 0});
249
3/6
✓ Branch 3 taken 61 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 61 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 61 times.
✗ Branch 10 not taken.
 61 circ.add_op<unsigned>(OpType::Ry, -t1, {0});
250
2/4
✓ Branch 3 taken 61 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 61 times.
✗ Branch 7 not taken.
 61 circ.add_op<unsigned>(OpType::CX, {2, 0});
251
3/6
✓ Branch 3 taken 61 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 61 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 61 times.
✗ Branch 10 not taken.
 61 circ.add_op<unsigned>(OpType::Ry, -t2, {0});
252
2/4
✓ Branch 3 taken 61 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 61 times.
✗ Branch 7 not taken.
 61 circ.add_op<unsigned>(OpType::CX, {1, 0});
253
2/4
✓ Branch 3 taken 61 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 61 times.
✗ Branch 7 not taken.
 61 circ.add_op<unsigned>(OpType::H, {0});
254
3/6
✓ Branch 3 taken 61 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 61 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 61 times.
✗ Branch 10 not taken.
 61 circ.add_op<unsigned>(OpType::Ry, t3, {0});
255 61 return circ;
256 }
257
258 // Return a 3-qubit circuit implementing the unitary
259 // [ C -S ]
260 // [ DS DC ]
261 // using TK2 and 1-qubit operations, where C and S are 4x4 real diagonal
262 // matrices, C^2 + S^2 = I, and D = diag(1,-1,1,-1).
263 //
264 // Note that
265 // [ C -S ] = U [ C -S ]
266 // [ DS DC ] [ S C ]
267 // where U represents a CZ on qubits 0 and 2. We convert the CZ operations to CX
268 // by adding Hadamards and simplifying H Ry(t) H to Ry(-t).
269 43 static Circuit two_qubit_modified_cossin_circ_tk(
270 const Eigen::Matrix4d &C, const Eigen::Matrix4d &S) {
271 // Compute angles a_i such that C_ii = cos(pi/2 a_i) and S_ii = sin(pi/2 a_i).
272 static const double f = 2 / PI;
273 43 double a0 = f * atan2(S(0, 0), C(0, 0));
274 43 double a1 = f * atan2(S(1, 1), C(1, 1));
275 43 double a2 = f * atan2(S(2, 2), C(2, 2));
276 43 double a3 = f * atan2(S(3, 3), C(3, 3));
277 43 double t0 = (a0 + a1 + a2 + a3) / 4;
278 43 double t1 = (a0 + a1 - a2 - a3) / 4;
279 43 double t2 = (a0 - a1 - a2 + a3) / 4;
280 43 double t3 = (a0 - a1 + a2 - a3) / 4;
281
1/2
✓ Branch 2 taken 43 times.
✗ Branch 3 not taken.
 43 Circuit circ(3);
282
3/6
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 43 times.
✗ Branch 10 not taken.
 43 circ.add_op<unsigned>(OpType::Ry, t0, {0});
283
2/4
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
 43 circ.add_op<unsigned>(OpType::H, {0});
284
5/10
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 86 times.
✓ Branch 15 taken 43 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 172 circ.append_with_map(
285
4/8
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43 times.
✗ Branch 11 not taken.
 129 CircPool::CX_using_TK2(), {{Qubit(0), Qubit(1)}, {Qubit(1), Qubit(0)}});
286
3/6
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 43 times.
✗ Branch 10 not taken.
 43 circ.add_op<unsigned>(OpType::Ry, -t1, {0});
287
5/10
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 86 times.
✓ Branch 15 taken 43 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 172 circ.append_with_map(
288
4/8
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43 times.
✗ Branch 11 not taken.
 129 CircPool::CX_using_TK2(), {{Qubit(0), Qubit(2)}, {Qubit(1), Qubit(0)}});
289
3/6
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 43 times.
✗ Branch 10 not taken.
 43 circ.add_op<unsigned>(OpType::Ry, -t2, {0});
290
5/10
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 86 times.
✓ Branch 15 taken 43 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 172 circ.append_with_map(
291
4/8
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43 times.
✗ Branch 11 not taken.
 129 CircPool::CX_using_TK2(), {{Qubit(0), Qubit(1)}, {Qubit(1), Qubit(0)}});
292
2/4
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
 43 circ.add_op<unsigned>(OpType::H, {0});
293
3/6
✓ Branch 3 taken 43 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 43 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 43 times.
✗ Branch 10 not taken.
 43 circ.add_op<unsigned>(OpType::Ry, t3, {0});
294 43 return circ;
295 }
296
297 // Given matrices U and V, if UV* is a scalar multiple of the identity, return
298 // the scalar.
299 751 static std::optional<Complex> id_coeff(
300 const Eigen::Matrix4cd &U, const Eigen::Matrix4cd &V) {
301
3/6
✓ Branch 1 taken 751 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 751 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 751 times.
✗ Branch 8 not taken.
 751 Eigen::Matrix4cd W = U * V.adjoint();
302
1/2
✓ Branch 1 taken 751 times.
✗ Branch 2 not taken.
 751 Complex w = W(0, 0);
303
5/8
✓ Branch 2 taken 751 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 751 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 751 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 417 times.
✓ Branch 11 taken 334 times.
2/2
✓ Decision 'true' taken 417 times.
✓ Decision 'false' taken 334 times.
 751 if (W.isApprox(w * Eigen::Matrix4cd::Identity())) {
304 417 return w;
305 }
306
3/4
✓ Branch 2 taken 334 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 13 times.
✓ Branch 5 taken 321 times.
2/2
✓ Decision 'true' taken 13 times.
✓ Decision 'false' taken 321 times.
 334 if (W.isZero()) {
307 // Eigen's isApprox compares multiplicatively, so doesn't catch this case.
308
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 13 return 0.;
309 }
310 321 return std::nullopt;
311 }
312
313 // If the given 8x8 unitary represents a circuit with no entanglement between
314 // qubit 0 and the other qubits, return a decomposition into a circuit on qubit
315 // 0 and a circuit on qubits 1 and 2. (The qubits on the second circuit are
316 // indexed with 0 and 1.)
317 330 static std::optional<std::pair<Circuit, Circuit>> separate_0_12(
318 const Eigen::MatrixXcd &U) {
319 // We want to check whether the unitary is of the form
320 // [ w_00 V w_01 V ]
321 // [ w_10 V w_11 V ]
322 // where W is a 2x2 unitary and V is a 4x4 unitary.
323 // W.l.o.g. we will assume w_00 (or w_10) is real and positive, compute the
324 // w_ij assuming the above form, and then check that the form is correct.
325
2/4
✓ Branch 1 taken 330 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 330 times.
✗ Branch 5 not taken.
 330 Eigen::Matrix4cd U00 = U.topLeftCorner(4, 4);
326
2/4
✓ Branch 1 taken 330 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 330 times.
✗ Branch 5 not taken.
 330 Eigen::Matrix4cd U01 = U.topRightCorner(4, 4);
327
2/4
✓ Branch 1 taken 330 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 330 times.
✗ Branch 5 not taken.
 330 Eigen::Matrix4cd U10 = U.bottomLeftCorner(4, 4);
328
2/4
✓ Branch 1 taken 330 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 330 times.
✗ Branch 5 not taken.
 330 Eigen::Matrix4cd U11 = U.bottomRightCorner(4, 4);
329 // If U is of the desired form, then U_ij U_kl* = w_ij W_kl* I for all
330 // i, j, k, l.
331
1/2
✓ Branch 1 taken 330 times.
✗ Branch 2 not taken.
 330 std::optional<Complex> w0000 = id_coeff(U00, U00); // |w_00|^2
332
2/2
✓ Branch 1 taken 171 times.
✓ Branch 2 taken 159 times.
2/2
✓ Decision 'true' taken 159 times.
✓ Decision 'false' taken 171 times.
 330 if (!w0000) return std::nullopt;
333
1/2
✓ Branch 1 taken 159 times.
✗ Branch 2 not taken.
 159 std::optional<Complex> w0101 = id_coeff(U01, U01); // |w_01|^2
334
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 159 times.
1/2
✓ Decision 'true' taken 159 times.
✗ Decision 'false' not taken.
 159 if (!w0101) return std::nullopt;
335 159 double x0000 = w0000->real(), y0000 = w0000->imag();
336 159 double x0101 = w0101->real(), y0101 = w0101->imag();
337
5/10
✓ Branch 1 taken 159 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 159 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 159 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 159 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 159 times.
 159 if (std::abs(y0000) > EPS || std::abs(y0101) > EPS || x0000 < -EPS ||
338 x0101 < -EPS) {
339 ✗ return std::nullopt;
340 }
341 // If negative because of rounding errors, clamp to 0:
342
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
 159 if (x0000 < 0.) x0000 = 0.;
343
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
 159 if (x0101 < 0.) x0101 = 0.;
344 // By unitarity of U, w0000 + w0101 = 1, and so x0000 + x0101 = 1.
345 // Choose the larger one to work with.
346 159 Complex w00, w01, w10, w11;
347
1/2
✓ Branch 1 taken 159 times.
✗ Branch 2 not taken.
 159 Eigen::Matrix4cd V;
348
2/2
✓ Branch 0 taken 107 times.
✓ Branch 1 taken 52 times.
 159 if (x0000 >= x0101) {
349 107 w00 = sqrt(x0000);
350
2/4
✓ Branch 1 taken 107 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 107 times.
✗ Branch 5 not taken.
 107 V = U00 / w00;
351
1/2
✓ Branch 1 taken 107 times.
✗ Branch 2 not taken.
 107 std::optional<Complex> w0001 = id_coeff(U00, U01);
352
2/2
✓ Branch 1 taken 63 times.
✓ Branch 2 taken 44 times.
 149 if (!w0001) return std::nullopt;
353
1/2
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
 44 std::optional<Complex> w0010 = id_coeff(U00, U10);
354
2/2
✓ Branch 1 taken 23 times.
✓ Branch 2 taken 21 times.
 44 if (!w0010) return std::nullopt;
355
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 21 std::optional<Complex> w0011 = id_coeff(U00, U11);
356
2/2
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 2 times.
 21 if (!w0011) return std::nullopt;
357 2 w01 = std::conj(*w0001) / w00;
358 2 w10 = std::conj(*w0010) / w00;
359 2 w11 = std::conj(*w0011) / w00;
360 } else {
361 52 w01 = sqrt(x0101);
362
2/4
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 52 times.
✗ Branch 5 not taken.
 52 V = U01 / w01;
363
1/2
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
 52 std::optional<Complex> w0100 = id_coeff(U01, U00);
364
2/2
✓ Branch 1 taken 21 times.
✓ Branch 2 taken 31 times.
 76 if (!w0100) return std::nullopt;
365
1/2
✓ Branch 1 taken 31 times.
✗ Branch 2 not taken.
 31 std::optional<Complex> w0110 = id_coeff(U01, U10);
366
2/2
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 7 times.
 31 if (!w0110) return std::nullopt;
367
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 std::optional<Complex> w0111 = id_coeff(U01, U11);
368
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 7 times.
 7 if (!w0111) return std::nullopt;
369 7 w00 = std::conj(*w0100) / w01;
370 7 w10 = std::conj(*w0110) / w01;
371 7 w11 = std::conj(*w0111) / w01;
372 }
373
374
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 Eigen::Matrix2cd W;
375
4/8
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 9 times.
✗ Branch 11 not taken.
 9 W << w00, w01, w10, w11;
376
3/6
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 9 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
 9 if (U.isApprox(Eigen::KroneckerProduct(W, V))) {
377
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 std::vector<double> angs = tk1_angles_from_unitary(W);
378
1/2
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
 9 Circuit c_1q(1);
379
8/16
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 9 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 9 times.
✗ Branch 16 not taken.
✓ Branch 19 taken 9 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 9 times.
✗ Branch 23 not taken.
✓ Branch 26 taken 27 times.
✓ Branch 27 taken 9 times.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
 36 c_1q.add_op<unsigned>(OpType::TK1, {angs[0], angs[1], angs[2]}, {0});
380
2/4
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 9 times.
✗ Branch 6 not taken.
 9 c_1q.add_phase(angs[3]);
381
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 Circuit c_2q = two_qubit_canonical(V);
382
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 18 return std::pair<Circuit, Circuit>(c_1q, c_2q);
383 9 }
384
385 ✗ return std::nullopt;
386 }
387
388 // Special cases worth handling. This is not necessary for correctness, but
389 // allows us to obtain circuits that are more amenable to later optimization.
390 113 static std::optional<Circuit> special_3q_synth(const Eigen::MatrixXcd &U) {
391 1 static const Eigen::PermutationMatrix<8> P1 = []() {
392
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Eigen::VectorXi V1(8);
393
8/16
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
 1 V1 << 0, 1, 4, 5, 2, 3, 6, 7;
394
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 return Eigen::PermutationMatrix<8>(V1);
395
4/8
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 112 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 114 }();
396 1 static const Eigen::PermutationMatrix<8> P2 = []() {
397
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Eigen::VectorXi V2(8);
398
8/16
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
 1 V2 << 0, 4, 2, 6, 1, 5, 3, 7;
399
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 return Eigen::PermutationMatrix<8>(V2);
400
4/8
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 112 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 114 }();
401
402 // Try separating qubit 0 from qubits 1 and 2:
403
1/2
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
 113 std::optional<std::pair<Circuit, Circuit>> c0 = separate_0_12(U);
404
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 110 times.
 113 if (c0) {
405
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 Circuit c_1q = c0->first;
406
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 Circuit c_2q = c0->second;
407
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 Circuit c(3);
408
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 c.append(c_1q);
409
8/16
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 3 times.
✗ Branch 12 not taken.
✓ Branch 16 taken 3 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✓ Branch 23 taken 6 times.
✓ Branch 24 taken 3 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
 9 c.append_with_map(c_2q, {{Qubit(0), Qubit(1)}, {Qubit(1), Qubit(2)}});
410
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 return c;
411 3 }
412
413 // Try separating qubit 1 from qubits 0 and 2:
414
4/8
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 110 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 110 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 110 times.
✗ Branch 11 not taken.
 110 std::optional<std::pair<Circuit, Circuit>> c1 = separate_0_12(P1 * U * P1);
415
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 107 times.
 110 if (c1) {
416
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 Circuit c_1q = c1->first;
417
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 Circuit c_2q = c1->second;
418
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 Circuit c(3);
419
6/12
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 3 times.
✓ Branch 17 taken 3 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
 6 c.append_with_map(c_1q, {{Qubit(0), Qubit(1)}});
420
6/12
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 3 times.
✓ Branch 17 taken 3 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
 6 c.append_with_map(c_2q, {{Qubit(1), Qubit(2)}});
421
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 return c;
422 3 }
423
424 // Try separating qubit 2 from qubits 0 and 1:
425
4/8
✓ Branch 1 taken 107 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 107 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 107 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 107 times.
✗ Branch 11 not taken.
 107 std::optional<std::pair<Circuit, Circuit>> c2 = separate_0_12(P2 * U * P2);
426
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 104 times.
 107 if (c2) {
427
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 Circuit c_1q = c2->first;
428
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 Circuit c_2q = c2->second;
429
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 Circuit c(3);
430
6/12
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 3 times.
✓ Branch 17 taken 3 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
 6 c.append_with_map(c_1q, {{Qubit(0), Qubit(2)}});
431
8/16
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 3 times.
✗ Branch 12 not taken.
✓ Branch 16 taken 3 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✓ Branch 23 taken 6 times.
✓ Branch 24 taken 3 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
 9 c.append_with_map(c_2q, {{Qubit(0), Qubit(1)}, {Qubit(1), Qubit(0)}});
432
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 return c;
433 3 }
434
435 104 return std::nullopt;
436 113 }
437
438 66 Circuit three_qubit_synthesis(const Eigen::MatrixXcd &U) {
439
3/6
✓ Branch 1 taken 66 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 66 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 66 times.
 66 if (U.rows() != 8 || U.cols() != 8) {
440 ✗ throw std::invalid_argument("Wrong-size matrix for three-qubit synthesis");
441 }
442
443
1/2
✓ Branch 1 taken 66 times.
✗ Branch 2 not taken.
 66 std::optional<Circuit> c_special = special_3q_synth(U);
444
3/4
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 61 times.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
 66 if (c_special) return *c_special;
445
446
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 auto [l0, l1, r0, r1, c, s] = CS_decomp(U);
447
3/6
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 61 times.
✗ Branch 8 not taken.
 61 auto [R_circ, z0] = two_qubit_plex(r0, r1, true);
448 61 Complex z1 = std::conj(z0);
449
1/2
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
 61 Circuit circ(3);
450
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 circ.append(R_circ);
451
4/8
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 61 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 61 times.
✗ Branch 11 not taken.
 61 circ.append(two_qubit_modified_cossin_circ(c, s));
452 // We chopped off the last CZ (on qubits 0 and 2) from the circuit
453 // implementing the CS decomposition. Account for this by changing the signs
454 // of columns 1 and 3 of l1.
455 //
456 // We also carried a diagonal from the earlier subcircuit, which commutes
457 // through the controls in the middle circuit and merges with l0 and l1.
458 //
459 // Together these imply the following adjustments:
460
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 61 l0.col(0) *= z0;
461
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 61 l0.col(1) *= z1;
462
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 61 l0.col(2) *= z1;
463
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 61 l0.col(3) *= z0;
464
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 61 l1.col(0) *= z0;
465
2/4
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 61 times.
✗ Branch 6 not taken.
 61 l1.col(1) *= -z1;
466
2/4
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
 61 l1.col(2) *= z1;
467
2/4
✓ Branch 2 taken 61 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 61 times.
✗ Branch 6 not taken.
 61 l1.col(3) *= -z0;
468
4/8
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 61 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 61 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 61 times.
✗ Branch 11 not taken.
 61 circ.append(two_qubit_plex(l0, l1, false).first);
469
1/2
✓ Branch 1 taken 61 times.
✗ Branch 2 not taken.
 61 return circ;
470 66 }
471
472 47 Circuit three_qubit_tk_synthesis(const Eigen::MatrixXcd &U) {
473
3/6
✓ Branch 1 taken 47 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 47 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 47 times.
 47 if (U.rows() != 8 || U.cols() != 8) {
474 ✗ throw std::invalid_argument("Wrong-size matrix for three-qubit synthesis");
475 }
476
477
1/2
✓ Branch 1 taken 47 times.
✗ Branch 2 not taken.
 47 std::optional<Circuit> c_special = special_3q_synth(U);
478
3/4
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 43 times.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 47 if (c_special) return *c_special;
479
480
1/2
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
 43 auto [l0, l1, r0, r1, c, s] = CS_decomp(U);
481
1/2
✓ Branch 2 taken 43 times.
✗ Branch 3 not taken.
 43 Circuit circ(3);
482
4/8
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43 times.
✗ Branch 11 not taken.
 43 circ.append(two_qubit_plex_tk(r0, r1));
483
4/8
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43 times.
✗ Branch 11 not taken.
 43 circ.append(two_qubit_modified_cossin_circ_tk(c, s));
484 // We chopped off the last CZ (on qubits 0 and 2) from the circuit
485 // implementing the CS decomposition. Account for this by changing the signs
486 // of columns 1 and 3 of l1.
487
2/4
✓ Branch 2 taken 43 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 43 times.
✗ Branch 6 not taken.
 43 l1.col(1) *= -1;
488
2/4
✓ Branch 2 taken 43 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 43 times.
✗ Branch 6 not taken.
 43 l1.col(3) *= -1;
489
4/8
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43 times.
✗ Branch 11 not taken.
 43 circ.append(two_qubit_plex_tk(l0, l1));
490
1/2
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
 43 return circ;
491 47 }
492
493 98 Eigen::MatrixXcd get_3q_unitary(const Circuit &c) {
494
2/4
✓ Branch 1 taken 98 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 98 times.
 98 if (c.n_qubits() != 3) {
495 ✗ throw CircuitInvalidity("Circuit in get_3q_unitary must have 3 qubits");
496 }
497
498 // Construct map from qubits to indices {0,1,2}.
499
1/2
✓ Branch 1 taken 98 times.
✗ Branch 2 not taken.
 98 qubit_vector_t all_qbs = c.all_qubits();
500 98 std::map<Qubit, unsigned> idx;
501
2/2
✓ Branch 0 taken 294 times.
✓ Branch 1 taken 98 times.
 392 for (unsigned i = 0; i < 3; i++) {
502
1/2
✓ Branch 2 taken 294 times.
✗ Branch 3 not taken.
 294 idx[all_qbs[i]] = i;
503 }
504
505 // Step through commands, building unitary as we go.
506
2/4
✓ Branch 1 taken 98 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 98 times.
✗ Branch 5 not taken.
 98 Eigen::MatrixXcd U = Eigen::MatrixXcd::Identity(8, 8);
507
6/10
✓ Branch 1 taken 98 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 98 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1189 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1189 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1189 times.
✓ Branch 14 taken 98 times.
 1287 for (const Command &cmd : c) {
508
1/2
✓ Branch 1 taken 1189 times.
✗ Branch 2 not taken.
 1189 qubit_vector_t qbs = cmd.get_qubits();
509 1189 Op_ptr op = cmd.get_op_ptr();
510
1/2
✓ Branch 2 taken 1189 times.
✗ Branch 3 not taken.
 1189 Gate_ptr gate = as_gate_ptr(op);
511
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1189 times.
 1189 if (!gate) {
512 ✗ throw CircuitInvalidity("Circuit in get_3q_unitary not unitary");
513 }
514
2/4
✓ Branch 1 taken 1189 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1189 times.
✗ Branch 5 not taken.
 1189 Eigen::MatrixXcd M = Eigen::MatrixXcd::Zero(8, 8);
515
2/3
✓ Branch 1 taken 722 times.
✓ Branch 2 taken 467 times.
✗ Branch 3 not taken.
 1189 switch (qbs.size()) {
516 722 case 1: {
517
1/2
✓ Branch 2 taken 722 times.
✗ Branch 3 not taken.
 722 std::vector<Expr> angles = gate->get_tk1_angles();
518
2/4
✓ Branch 1 taken 722 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 722 times.
✗ Branch 5 not taken.
 722 Eigen::Matrix2cd u = get_matrix_from_tk1_angles(angles);
519 // Construct the 8x8 matrix representing u.
520
4/6
✓ Branch 2 taken 722 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 272 times.
✓ Branch 5 taken 234 times.
✓ Branch 6 taken 216 times.
✗ Branch 7 not taken.
 722 unsigned i = idx[qbs[0]];
521 switch (i) {
522 272 case 0:
523
5/10
✓ Branch 1 taken 272 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 272 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 272 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 272 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 272 times.
✗ Branch 14 not taken.
 272 M(0, 0) = M(1, 1) = M(2, 2) = M(3, 3) = u(0, 0);
524
5/10
✓ Branch 1 taken 272 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 272 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 272 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 272 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 272 times.
✗ Branch 14 not taken.
 272 M(0, 4) = M(1, 5) = M(2, 6) = M(3, 7) = u(0, 1);
525
5/10
✓ Branch 1 taken 272 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 272 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 272 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 272 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 272 times.
✗ Branch 14 not taken.
 272 M(4, 0) = M(5, 1) = M(6, 2) = M(7, 3) = u(1, 0);
526
5/10
✓ Branch 1 taken 272 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 272 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 272 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 272 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 272 times.
✗ Branch 14 not taken.
 272 M(4, 4) = M(5, 5) = M(6, 6) = M(7, 7) = u(1, 1);
527 272 break;
528 234 case 1:
529
5/10
✓ Branch 1 taken 234 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 234 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 234 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 234 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 234 times.
✗ Branch 14 not taken.
 234 M(0, 0) = M(1, 1) = M(4, 4) = M(5, 5) = u(0, 0);
530
5/10
✓ Branch 1 taken 234 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 234 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 234 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 234 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 234 times.
✗ Branch 14 not taken.
 234 M(0, 2) = M(1, 3) = M(4, 6) = M(5, 7) = u(0, 1);
531
5/10
✓ Branch 1 taken 234 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 234 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 234 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 234 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 234 times.
✗ Branch 14 not taken.
 234 M(2, 0) = M(3, 1) = M(6, 4) = M(7, 5) = u(1, 0);
532
5/10
✓ Branch 1 taken 234 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 234 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 234 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 234 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 234 times.
✗ Branch 14 not taken.
 234 M(2, 2) = M(3, 3) = M(6, 6) = M(7, 7) = u(1, 1);
533 234 break;
534 216 case 2:
535
5/10
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 216 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 216 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 216 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 216 times.
✗ Branch 14 not taken.
 216 M(0, 0) = M(2, 2) = M(4, 4) = M(6, 6) = u(0, 0);
536
5/10
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 216 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 216 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 216 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 216 times.
✗ Branch 14 not taken.
 216 M(0, 1) = M(2, 3) = M(4, 5) = M(6, 7) = u(0, 1);
537
5/10
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 216 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 216 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 216 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 216 times.
✗ Branch 14 not taken.
 216 M(1, 0) = M(3, 2) = M(5, 4) = M(7, 6) = u(1, 0);
538
5/10
✓ Branch 1 taken 216 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 216 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 216 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 216 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 216 times.
✗ Branch 14 not taken.
 216 M(1, 1) = M(3, 3) = M(5, 5) = M(7, 7) = u(1, 1);
539 216 break;
540 ✗ default:
541 TKET_ASSERT(!"Invalid index");
542 }
543 722 break;
544 722 }
545 467 case 2: {
546
2/4
✓ Branch 2 taken 467 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 467 times.
✗ Branch 6 not taken.
 467 Eigen::Matrix4cd m = gate->get_unitary();
547 // Note reversal of indices here:
548
1/2
✓ Branch 2 taken 467 times.
✗ Branch 3 not taken.
 467 unsigned i = 2 - idx[qbs[0]];
549
1/2
✓ Branch 2 taken 467 times.
✗ Branch 3 not taken.
 467 unsigned j = 2 - idx[qbs[1]];
550 // Construct the 8x8 matrix representing the gate.
551 // Let k be the untouched index, so that {i,j,k} = {0,1,2}:
552 467 unsigned k = 3 - i - j;
553 467 unsigned t = 1 << k;
554
2/2
✓ Branch 0 taken 1868 times.
✓ Branch 1 taken 467 times.
 2335 for (unsigned s0 = 0; s0 < 4; s0++) {
555 1868 unsigned s0_ = ((s0 >> 1) << i) + ((s0 & 1) << j);
556
2/2
✓ Branch 0 taken 7472 times.
✓ Branch 1 taken 1868 times.
 9340 for (unsigned s1 = 0; s1 < 4; s1++) {
557 7472 unsigned s1_ = ((s1 >> 1) << i) + ((s1 & 1) << j);
558
3/6
✓ Branch 1 taken 7472 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7472 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7472 times.
✗ Branch 8 not taken.
 7472 M(s0_, s1_) = M(s0_ + t, s1_ + t) = m(s0, s1);
559 }
560 }
561 467 break;
562 }
563 ✗ default:
564 ✗ throw CircuitInvalidity(
565 ✗ "Circuit in get_3q_unitary contains gates with more than 2 qubits");
566 }
567
2/4
✓ Branch 1 taken 1189 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1189 times.
✗ Branch 5 not taken.
 1189 U = M * U;
568 1287 }
569
570
5/10
✓ Branch 1 taken 98 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 98 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 98 times.
✗ Branch 8 not taken.
✓ Branch 13 taken 98 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 98 times.
✗ Branch 17 not taken.
 294 return std::exp(i_ * PI * eval_expr(c.get_phase()).value()) * U;
571 98 }
572
573 } // namespace tket
574