GCC Code Coverage Report

Directory: ./
File: Gate/Rotation.cpp
Date: 2022-10-15 05:10:18
Warnings: 7 unchecked decisions!
Exec Total Coverage
Lines: 240 276 87.0%
Functions: 10 11 90.9%
Branches: 436 816 53.4%
Decisions: 90 140 64.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 "Rotation.hpp"
16
17 #include "OpType/OpDesc.hpp"
18 #include "OpType/OpType.hpp"
19 #include "Utils/Expression.hpp"
20 #include "Utils/Symbols.hpp"
21 #include "symengine/constants.h"
22
23 namespace tket {
24
25 5966 static Expr atan2_bypi(const Expr &a, const Expr &b) {
26
1/2
✓ Branch 1 taken 5966 times.
✗ Branch 2 not taken.
 5966 std::optional<double> va = eval_expr(a);
27
1/2
✓ Branch 1 taken 5966 times.
✗ Branch 2 not taken.
 5966 std::optional<double> vb = eval_expr(b);
28
6/6
✓ Branch 1 taken 5915 times.
✓ Branch 2 taken 51 times.
✓ Branch 4 taken 5914 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 5914 times.
✓ Branch 7 taken 52 times.
2/2
✓ Decision 'true' taken 5914 times.
✓ Decision 'false' taken 52 times.
 5966 if (va && vb) {
29
1/2
✓ Branch 1 taken 5914 times.
✗ Branch 2 not taken.
 5914 double vva = va.value();
30
1/2
✓ Branch 1 taken 5914 times.
✗ Branch 2 not taken.
 5914 double vvb = vb.value();
31
4/8
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 5868 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 46 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 5914 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
1/2
✓ Decision 'true' taken 5914 times.
✗ Decision 'false' not taken.
 5914 if (std::abs(vva) < EPS && std::abs(vvb) < EPS) return Expr(0.);
32
1/2
✓ Branch 1 taken 5914 times.
✗ Branch 2 not taken.
 5914 return atan2(vva, vvb) / PI;
33 } else {
34 // Convert symbolic zero to 0. This is a workaround for
35 // https://github.com/symengine/symengine/issues/1875 .
36
2/4
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 52 times.
✗ Branch 5 not taken.
 52 Expr a1 = a, b1 = b;
37
5/8
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 52 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 51 times.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
1/2
✓ Decision 'true' taken 52 times.
✗ Decision 'false' not taken.
 52 if (a1 == SymEngine::zero) a1 = 0.;
38
3/8
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 52 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 52 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
0/1
? Decision couldn't be analyzed.
 52 if (b1 == SymEngine::zero) b1 = 0.;
39
4/8
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 52 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 52 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 52 times.
✗ Branch 13 not taken.
 104 return SymEngine::div(SymEngine::atan2(a1, b1), SymEngine::pi);
40 52 }
41 }
42
43 741 static Expr acos_bypi(const Expr &a) {
44
1/2
✓ Branch 1 taken 741 times.
✗ Branch 2 not taken.
 741 std::optional<double> va = eval_expr(a);
45
2/2
✓ Branch 1 taken 731 times.
✓ Branch 2 taken 10 times.
2/2
✓ Decision 'true' taken 731 times.
✓ Decision 'false' taken 10 times.
 741 if (va) {
46
1/2
✓ Branch 1 taken 731 times.
✗ Branch 2 not taken.
 731 double vva = va.value();
47 // avoid undefined values due to rounding
48
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 731 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
1/2
✓ Decision 'true' taken 731 times.
✗ Decision 'false' not taken.
 731 if (vva >= 1.) return 0.;
49
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 731 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
1/2
✓ Decision 'true' taken 731 times.
✗ Decision 'false' not taken.
 731 if (vva <= -1.) return 1.;
50
1/2
✓ Branch 1 taken 731 times.
✗ Branch 2 not taken.
 731 return acos(vva) / PI;
51 } else {
52
4/8
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 10 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 10 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 10 times.
✗ Branch 12 not taken.
 20 return SymEngine::div(SymEngine::acos(a), SymEngine::pi);
53 }
54 }
55
56 // SymEngine::div does not always spot when the numerator is a scalar multiple
57 // of the denominator, for example in expressions like (-a + b) / (a - b) where
58 // a and b are symbolic. This function picks out the common cases.
59 1150 static Expr expr_div(const Expr &num, const Expr &den) {
60
4/6
✓ Branch 2 taken 1150 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1150 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 504 times.
✓ Branch 10 taken 646 times.
2/2
✓ Decision 'true' taken 646 times.
✓ Decision 'false' taken 504 times.
 1150 if (approx_0(SymEngine::expand(num - den))) return 1;
61
4/6
✓ Branch 2 taken 646 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 646 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 294 times.
✓ Branch 10 taken 352 times.
0/1
? Decision couldn't be analyzed.
 646 if (approx_0(SymEngine::expand(num + den))) return -1;
62
1/2
✓ Branch 4 taken 352 times.
✗ Branch 5 not taken.
 704 return SymEngine::div(num, den);
63 }
64
65 3600 static std::tuple<Expr, Expr, Expr> xyx_angles_from_coeffs(
66 const Expr &s, const Expr &i, const Expr &j, const Expr &k) {
67 // Handle exceptional cases first.
68
1/2
✓ Branch 1 taken 3600 times.
✗ Branch 2 not taken.
 3600 bool s_zero = approx_0(s);
69
3/6
✓ Branch 1 taken 3600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3600 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3600 times.
✗ Branch 8 not taken.
 3600 bool s_one = approx_0(s - 1);
70
1/2
✓ Branch 1 taken 3600 times.
✗ Branch 2 not taken.
 3600 bool i_zero = approx_0(i);
71
3/6
✓ Branch 1 taken 3600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3600 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3600 times.
✗ Branch 8 not taken.
 3600 bool i_one = approx_0(i - 1);
72
1/2
✓ Branch 1 taken 3600 times.
✗ Branch 2 not taken.
 3600 bool j_zero = approx_0(j);
73
3/6
✓ Branch 1 taken 3600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3600 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3600 times.
✗ Branch 8 not taken.
 3600 bool j_one = approx_0(j - 1);
74
1/2
✓ Branch 1 taken 3600 times.
✗ Branch 2 not taken.
 3600 bool k_zero = approx_0(k);
75
3/6
✓ Branch 1 taken 3600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3600 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3600 times.
✗ Branch 8 not taken.
 3600 bool k_one = approx_0(k - 1);
76
5/6
✓ Branch 0 taken 503 times.
✓ Branch 1 taken 3097 times.
✓ Branch 2 taken 358 times.
✓ Branch 3 taken 145 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 358 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 3600 times.
 3600 if (i_zero && j_zero && k_zero) {
77
0/2
✗ Decision 'true' not taken.
✗ Decision 'false' not taken.
 ✗ if (s_one)
78 ✗ return {0, 0, 0};
79 else
80 ✗ return {2, 0, 0};
81 }
82
5/6
✓ Branch 0 taken 1398 times.
✓ Branch 1 taken 2202 times.
✓ Branch 2 taken 268 times.
✓ Branch 3 taken 1130 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 268 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 3600 times.
 3600 if (s_zero && j_zero && k_zero) {
83
0/2
✗ Decision 'true' not taken.
✗ Decision 'false' not taken.
 ✗ if (i_one)
84 ✗ return {1, 0, 0};
85 else
86 ✗ return {3, 0, 0};
87 }
88
5/6
✓ Branch 0 taken 1398 times.
✓ Branch 1 taken 2202 times.
✓ Branch 2 taken 152 times.
✓ Branch 3 taken 1246 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 152 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 3600 times.
 3600 if (s_zero && i_zero && k_zero) {
89
0/2
✗ Decision 'true' not taken.
✗ Decision 'false' not taken.
 ✗ if (j_one)
90 ✗ return {0, 1, 0};
91 else
92 ✗ return {0, 3, 0};
93 }
94
6/6
✓ Branch 0 taken 1398 times.
✓ Branch 1 taken 2202 times.
✓ Branch 2 taken 152 times.
✓ Branch 3 taken 1246 times.
✓ Branch 4 taken 43 times.
✓ Branch 5 taken 109 times.
2/2
✓ Decision 'true' taken 43 times.
✓ Decision 'false' taken 3557 times.
 3600 if (s_zero && i_zero && j_zero) {
95
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 27 times.
2/2
✓ Decision 'true' taken 32 times.
✓ Decision 'false' taken 11 times.
 43 if (k_one)
96
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 32 return {3, 1, 0};
97 else
98
1/2
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
 54 return {1, 1, 0};
99 }
100
8/12
✓ Branch 0 taken 1355 times.
✓ Branch 1 taken 2202 times.
✓ Branch 2 taken 109 times.
✓ Branch 3 taken 1246 times.
✓ Branch 5 taken 109 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 109 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 109 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 109 times.
✗ Branch 15 not taken.
0/1
? Decision couldn't be analyzed.
 3666 if (s_zero && i_zero) return {-2 * atan2_bypi(k, j), 1, 0};
101
8/12
✓ Branch 0 taken 1246 times.
✓ Branch 1 taken 2202 times.
✓ Branch 2 taken 225 times.
✓ Branch 3 taken 1021 times.
✓ Branch 5 taken 225 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 225 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 225 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 225 times.
✗ Branch 15 not taken.
0/1
? Decision couldn't be analyzed.
 3673 if (s_zero && j_zero) return {0, 2 * atan2_bypi(k, i), 1};
102
8/12
✓ Branch 0 taken 1021 times.
✓ Branch 1 taken 2202 times.
✓ Branch 2 taken 1017 times.
✓ Branch 3 taken 4 times.
✓ Branch 5 taken 1017 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1017 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1017 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 1017 times.
✗ Branch 15 not taken.
2/2
✓ Decision 'true' taken 2521 times.
✓ Decision 'false' taken 1719 times.
 4240 if (s_zero && k_zero) return {0.5, 2 * atan2_bypi(j, i), 0.5};
103
8/12
✓ Branch 0 taken 351 times.
✓ Branch 1 taken 1855 times.
✓ Branch 2 taken 315 times.
✓ Branch 3 taken 36 times.
✓ Branch 5 taken 315 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 315 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 315 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 315 times.
✗ Branch 15 not taken.
2/2
✓ Decision 'true' taken 1891 times.
✓ Decision 'false' taken 630 times.
 2521 if (i_zero && j_zero) return {-0.5, 2 * atan2_bypi(k, s), 0.5};
104
3/12
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 1855 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 36 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
1/2
✓ Decision 'true' taken 1891 times.
✗ Decision 'false' not taken.
 1891 if (i_zero && k_zero) return {0, 2 * atan2_bypi(j, s), 0};
105
3/12
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1885 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
1/2
✓ Decision 'true' taken 1891 times.
✗ Decision 'false' not taken.
 1891 if (j_zero && k_zero) return {2 * atan2_bypi(i, s), 0, 0};
106
107 // This is a (partial) workaround for
108 // https://github.com/symengine/symengine/issues/1806
109 // (since it avoids the use of atan2 with proportional symbolics).
110 // Explanation: When the quaternion is of the form
111 // A + uA i + B j -/+ uB k
112 // where u is a pure number (with no free symbols) but A and B are symbolic,
113 // symengine wrongly simplifies atan2(uA, A) and atan2(uB, B), ignoring the
114 // symbols (whose sign ought to affect the result). So the general formula
115 // below does not work. We therefore handle these cases differently, noting
116 // that the quaternion factorizes as either
117 // (A + Bj) (1 + ui) or (1 + ui) (A + Bj)
118 // -- that is, an Rx followed by an Ry or vice versa.
119 // We will analyse the first case; the second is similar.
120 // Let alpha = atan(u). Note that -pi/2 < alpha < pi/2, so cos(alpha) > 0.
121 // The quaternion is then
122 // (A/cos(alpha) + B/cos(alpha) j) (cos(alpha) + sin(alpha) i)
123 // So we can take the angle of the Rx as 2*alpha, and the angle of the Ry as
124 // 2 * atan2(B, A).
125 // Finally, note that u must be well-defined because we have already dealt
126 // with all cases where s = 0.
127
7/12
✓ Branch 1 taken 1891 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1891 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1891 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1891 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1891 times.
✗ Branch 14 not taken.
✓ Branch 19 taken 456 times.
✓ Branch 20 taken 1435 times.
2/2
✓ Decision 'true' taken 456 times.
✓ Decision 'false' taken 1435 times.
 1891 if (approx_0(SymEngine::expand(i * j + s * k))) {
128
1/2
✓ Branch 1 taken 456 times.
✗ Branch 2 not taken.
 456 Expr u = expr_div(i, s);
129
3/6
✓ Branch 1 taken 456 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 456 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 456 times.
✗ Branch 9 not taken.
1/2
✓ Decision 'true' taken 456 times.
✗ Decision 'false' not taken.
 456 if (SymEngine::free_symbols(u).empty()) {
130
2/4
✓ Branch 2 taken 456 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 456 times.
✗ Branch 6 not taken.
 456 Expr a = SymEngine::atan(u);
131
3/6
✓ Branch 1 taken 456 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 456 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 456 times.
✗ Branch 8 not taken.
 912 Expr q = 2 * atan2_bypi(j, s);
132
5/10
✓ Branch 1 taken 456 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 456 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 456 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 456 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 456 times.
✗ Branch 15 not taken.
 912 Expr two_a_by_pi = SymEngine::div(2 * a, SymEngine::pi);
133
1/2
✓ Branch 1 taken 456 times.
✗ Branch 2 not taken.
 456 return std::tuple<Expr, Expr, Expr>(two_a_by_pi, q, 0);
134 456 }
135
8/14
✗ Branch 1 not taken.
✓ Branch 2 taken 456 times.
✓ Branch 4 taken 1435 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1435 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1435 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1435 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1435 times.
✗ Branch 17 not taken.
✓ Branch 22 taken 694 times.
✓ Branch 23 taken 741 times.
2/2
✓ Decision 'true' taken 694 times.
✓ Decision 'false' taken 1197 times.
 1891 } else if (approx_0(SymEngine::expand(i * j - s * k))) {
136
1/2
✓ Branch 1 taken 694 times.
✗ Branch 2 not taken.
 694 Expr u = expr_div(i, s);
137
3/6
✓ Branch 1 taken 694 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 694 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 694 times.
✗ Branch 9 not taken.
1/2
✓ Decision 'true' taken 694 times.
✗ Decision 'false' not taken.
 694 if (SymEngine::free_symbols(u).empty()) {
138
2/4
✓ Branch 2 taken 694 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 694 times.
✗ Branch 6 not taken.
 694 Expr a = SymEngine::atan(u);
139
3/6
✓ Branch 1 taken 694 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 694 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 694 times.
✗ Branch 8 not taken.
 1388 Expr q = 2 * atan2_bypi(j, s);
140
5/10
✓ Branch 1 taken 694 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 694 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 694 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 694 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 694 times.
✗ Branch 15 not taken.
 1388 Expr two_a_by_pi = SymEngine::div(2 * a, SymEngine::pi);
141
1/2
✓ Branch 1 taken 694 times.
✗ Branch 2 not taken.
 694 return std::tuple<Expr, Expr, Expr>(0, q, two_a_by_pi);
142 694 }
143
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 694 times.
 694 }
144
145 // Now the general case.
146
1/2
✓ Branch 1 taken 741 times.
✗ Branch 2 not taken.
 741 Expr a = atan2_bypi(i, s);
147
1/2
✓ Branch 1 taken 741 times.
✗ Branch 2 not taken.
 741 Expr b = atan2_bypi(k, j);
148
9/18
✓ Branch 1 taken 741 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 741 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 741 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 741 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 741 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 741 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 741 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 741 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 741 times.
✗ Branch 26 not taken.
 1482 Expr q = acos_bypi(SymEngine::expand(s * s + i * i - j * j - k * k));
149
3/6
✓ Branch 1 taken 741 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 741 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 741 times.
✗ Branch 8 not taken.
 1482 return std::tuple<Expr, Expr, Expr>(a - b, q, a + b);
150 741 }
151
152 188902 Rotation::Rotation(OpType optype, Expr a)
153
4/8
✓ Branch 2 taken 188902 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 188902 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 188902 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 188902 times.
✗ Branch 12 not taken.
 188902 : i_(0), j_(0), k_(0), optype_(optype), a_(a) {
154
3/4
✓ Branch 1 taken 188902 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 81264 times.
✓ Branch 4 taken 107638 times.
2/2
✓ Decision 'true' taken 81264 times.
✓ Decision 'false' taken 107638 times.
 188902 if (equiv_0(a, 4)) {
155 81264 rep_ = Rep::id;
156
1/2
✓ Branch 1 taken 81264 times.
✗ Branch 2 not taken.
 81264 s_ = 1;
157
3/6
✓ Branch 1 taken 81264 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 81264 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 81264 times.
✗ Branch 9 not taken.
 81264 i_ = j_ = k_ = 0;
158
5/8
✓ Branch 1 taken 107638 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 107638 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 107638 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 634 times.
✓ Branch 12 taken 107004 times.
2/2
✓ Decision 'true' taken 634 times.
✓ Decision 'false' taken 107004 times.
 107638 } else if (equiv_0(a - 2, 4)) {
159 634 rep_ = Rep::minus_id;
160
1/2
✓ Branch 1 taken 634 times.
✗ Branch 2 not taken.
 634 s_ = -1;
161
3/6
✓ Branch 1 taken 634 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 634 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 634 times.
✗ Branch 9 not taken.
 634 i_ = j_ = k_ = 0;
162 } else {
163 107004 rep_ = Rep::orth_rot;
164
1/2
✓ Branch 1 taken 107004 times.
✗ Branch 2 not taken.
 107004 s_ = cos_halfpi_times(a);
165
1/2
✓ Branch 1 taken 107004 times.
✗ Branch 2 not taken.
 107004 Expr t = sin_halfpi_times(a);
166
3/4
✓ Branch 0 taken 13945 times.
✓ Branch 1 taken 23220 times.
✓ Branch 2 taken 69839 times.
✗ Branch 3 not taken.
 107004 switch (optype) {
167
1/1
✓ Decision 'true' taken 13945 times.
 13945 case OpType::Rx:
168
1/2
✓ Branch 1 taken 13945 times.
✗ Branch 2 not taken.
 13945 i_ = t;
169 13945 break;
170
1/1
✓ Decision 'true' taken 23220 times.
 23220 case OpType::Ry:
171
1/2
✓ Branch 1 taken 23220 times.
✗ Branch 2 not taken.
 23220 j_ = t;
172 23220 break;
173
1/1
✓ Decision 'true' taken 69839 times.
 69839 case OpType::Rz:
174
1/2
✓ Branch 1 taken 69839 times.
✗ Branch 2 not taken.
 69839 k_ = t;
175 69839 break;
176
0/1
✗ Decision 'true' not taken.
 ✗ default:
177 ✗ throw std::logic_error(
178 "Quaternions can only be constructed "
179 ✗ "from Rx, Ry or Rz rotations");
180 }
181 107004 }
182 188902 }
183
184 64404 std::optional<Expr> Rotation::angle(OpType optype) const {
185
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 64404 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 64404 times.
 64404 if (rep_ == Rep::id) {
186 ✗ return Expr(0);
187
2/2
✓ Branch 0 taken 697 times.
✓ Branch 1 taken 63707 times.
2/2
✓ Decision 'true' taken 697 times.
✓ Decision 'false' taken 63707 times.
 64404 } else if (rep_ == Rep::minus_id) {
188 697 return Expr(2);
189
3/4
✓ Branch 0 taken 63707 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 48823 times.
✓ Branch 3 taken 14884 times.
2/2
✓ Decision 'true' taken 48823 times.
✓ Decision 'false' taken 14884 times.
 63707 } else if (rep_ == Rep::orth_rot && optype_ == optype) {
190 48823 return a_;
191 } else {
192 14884 return std::nullopt;
193 }
194 }
195
196 37351 std::tuple<Expr, Expr, Expr> Rotation::to_pqp(OpType p, OpType q) const {
197
2/2
✓ Branch 0 taken 9791 times.
✓ Branch 1 taken 27560 times.
2/2
✓ Decision 'true' taken 9791 times.
✓ Decision 'false' taken 27560 times.
 37351 if (rep_ == Rep::id) {
198
2/4
✓ Branch 2 taken 9791 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 9791 times.
✗ Branch 6 not taken.
 9791 return {Expr(0), Expr(0), Expr(0)};
199
2/2
✓ Branch 0 taken 137 times.
✓ Branch 1 taken 27423 times.
2/2
✓ Decision 'true' taken 137 times.
✓ Decision 'false' taken 27423 times.
 27560 } else if (rep_ == Rep::minus_id) {
200
2/4
✓ Branch 2 taken 137 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 137 times.
✗ Branch 6 not taken.
 137 return {Expr(2), Expr(0), Expr(0)};
201
2/2
✓ Branch 0 taken 24181 times.
✓ Branch 1 taken 3242 times.
2/2
✓ Decision 'true' taken 24181 times.
✓ Decision 'false' taken 3242 times.
 27423 } else if (rep_ == Rep::orth_rot) {
202
2/2
✓ Branch 0 taken 370 times.
✓ Branch 1 taken 23811 times.
2/2
✓ Decision 'true' taken 740 times.
✓ Decision 'false' taken 23441 times.
 24181 if (optype_ == p) {
203
2/4
✓ Branch 2 taken 370 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 370 times.
✗ Branch 6 not taken.
 740 return {a_, Expr(0), Expr(0)};
204
2/2
✓ Branch 0 taken 23453 times.
✓ Branch 1 taken 358 times.
0/1
? Decision couldn't be analyzed.
 23811 } else if (optype_ == q) {
205
2/4
✓ Branch 2 taken 23453 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 23453 times.
✗ Branch 6 not taken.
 46906 return {Expr(0), a_, Expr(0)};
206 }
207 }
208
3/4
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 3550 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 50 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 3600 times.
 3600 if (p == OpType::Rx && q == OpType::Ry) {
209 ✗ return xyx_angles_from_coeffs(s_, i_, j_, k_);
210
3/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3599 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 3600 times.
 3600 } else if (p == OpType::Ry && q == OpType::Rx) {
211 ✗ return xyx_angles_from_coeffs(s_, j_, i_, -k_);
212
3/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3599 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
2/2
✓ Decision 'true' taken 1 times.
✓ Decision 'false' taken 3599 times.
 3600 } else if (p == OpType::Ry && q == OpType::Rz) {
213 1 return xyx_angles_from_coeffs(s_, j_, k_, i_);
214
4/4
✓ Branch 0 taken 3549 times.
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 3364 times.
✓ Branch 3 taken 185 times.
0/1
? Decision couldn't be analyzed.
 3599 } else if (p == OpType::Rz && q == OpType::Ry) {
215
1/2
✓ Branch 2 taken 3364 times.
✗ Branch 3 not taken.
 6728 return xyx_angles_from_coeffs(s_, k_, j_, -i_);
216
3/4
✓ Branch 0 taken 185 times.
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 185 times.
✗ Branch 3 not taken.
2/2
✓ Decision 'true' taken 185 times.
✓ Decision 'false' taken 50 times.
 235 } else if (p == OpType::Rz && q == OpType::Rx) {
217 185 return xyx_angles_from_coeffs(s_, k_, i_, j_);
218
2/4
✓ Branch 0 taken 50 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
0/1
? Decision couldn't be analyzed.
 50 } else if (p == OpType::Rx && q == OpType::Rz) {
219
1/2
✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
 100 return xyx_angles_from_coeffs(s_, i_, k_, -j_);
220 } else {
221 ✗ throw std::logic_error("Axes must be a pair of X, Y, Z.");
222 }
223 }
224
225 // Table of compositions
226 static const std::map<
227 std::tuple<OpType, int, OpType, int>, std::pair<OpType, int>>
228 product = {
229 {{OpType::Rx, +1, OpType::Ry, +1}, {OpType::Rz, -1}},
230 {{OpType::Rx, +1, OpType::Rz, +1}, {OpType::Ry, +1}},
231 {{OpType::Rx, +1, OpType::Ry, -1}, {OpType::Rz, +1}},
232 {{OpType::Rx, +1, OpType::Rz, -1}, {OpType::Ry, -1}},
233 {{OpType::Rx, -1, OpType::Ry, +1}, {OpType::Rz, +1}},
234 {{OpType::Rx, -1, OpType::Rz, +1}, {OpType::Ry, -1}},
235 {{OpType::Rx, -1, OpType::Ry, -1}, {OpType::Rz, -1}},
236 {{OpType::Rx, -1, OpType::Rz, -1}, {OpType::Ry, +1}},
237 {{OpType::Ry, +1, OpType::Rz, +1}, {OpType::Rx, -1}},
238 {{OpType::Ry, +1, OpType::Rx, +1}, {OpType::Rz, +1}},
239 {{OpType::Ry, +1, OpType::Rz, -1}, {OpType::Rx, +1}},
240 {{OpType::Ry, +1, OpType::Rx, -1}, {OpType::Rz, -1}},
241 {{OpType::Ry, -1, OpType::Rz, +1}, {OpType::Rx, +1}},
242 {{OpType::Ry, -1, OpType::Rx, +1}, {OpType::Rz, -1}},
243 {{OpType::Ry, -1, OpType::Rz, -1}, {OpType::Rx, -1}},
244 {{OpType::Ry, -1, OpType::Rx, -1}, {OpType::Rz, +1}},
245 {{OpType::Rz, +1, OpType::Rx, +1}, {OpType::Ry, -1}},
246 {{OpType::Rz, +1, OpType::Ry, +1}, {OpType::Rx, +1}},
247 {{OpType::Rz, +1, OpType::Rx, -1}, {OpType::Ry, +1}},
248 {{OpType::Rz, +1, OpType::Ry, -1}, {OpType::Rx, -1}},
249 {{OpType::Rz, -1, OpType::Rx, +1}, {OpType::Ry, +1}},
250 {{OpType::Rz, -1, OpType::Ry, +1}, {OpType::Rx, -1}},
251 {{OpType::Rz, -1, OpType::Rx, -1}, {OpType::Ry, -1}},
252 {{OpType::Rz, -1, OpType::Ry, -1}, {OpType::Rx, +1}},
253 };
254
255 55248 void Rotation::apply(const Rotation &other) {
256
2/2
✓ Branch 0 taken 583 times.
✓ Branch 1 taken 54665 times.
2/2
✓ Decision 'true' taken 54665 times.
✓ Decision 'false' taken 28551 times.
 83216 if (other.rep_ == Rep::id) return;
257
258
2/2
✓ Branch 0 taken 27625 times.
✓ Branch 1 taken 27040 times.
2/2
✓ Decision 'true' taken 27625 times.
✓ Decision 'false' taken 27040 times.
 54665 if (rep_ == Rep::id) {
259 27625 rep_ = other.rep_;
260
1/2
✓ Branch 1 taken 27625 times.
✗ Branch 2 not taken.
 27625 s_ = other.s_;
261
1/2
✓ Branch 1 taken 27625 times.
✗ Branch 2 not taken.
 27625 i_ = other.i_;
262
1/2
✓ Branch 1 taken 27625 times.
✗ Branch 2 not taken.
 27625 j_ = other.j_;
263
1/2
✓ Branch 1 taken 27625 times.
✗ Branch 2 not taken.
 27625 k_ = other.k_;
264 27625 optype_ = other.optype_;
265
1/2
✓ Branch 1 taken 27625 times.
✗ Branch 2 not taken.
 27625 a_ = other.a_;
266 27625 return;
267 }
268
269
2/2
✓ Branch 0 taken 343 times.
✓ Branch 1 taken 26697 times.
2/2
✓ Decision 'true' taken 343 times.
✓ Decision 'false' taken 26697 times.
 27040 if (rep_ == Rep::minus_id) {
270
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 338 times.
2/2
✓ Decision 'true' taken 5 times.
✓ Decision 'false' taken 338 times.
 343 if (other.rep_ == Rep::minus_id) {
271 5 rep_ = Rep::id;
272
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 s_ = 1;
273
3/6
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 5 times.
✗ Branch 9 not taken.
 5 i_ = j_ = k_ = 0;
274
1/2
✓ Branch 0 taken 338 times.
✗ Branch 1 not taken.
1/2
✓ Decision 'true' taken 338 times.
✗ Decision 'false' not taken.
 338 } else if (other.rep_ == Rep::orth_rot) {
275 338 rep_ = Rep::orth_rot;
276 338 optype_ = other.optype_;
277
2/4
✓ Branch 1 taken 338 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 338 times.
✗ Branch 5 not taken.
 338 a_ = other.a_ + 2;
278
1/2
✓ Branch 1 taken 338 times.
✗ Branch 2 not taken.
 338 s_ = -other.s_;
279
1/2
✓ Branch 1 taken 338 times.
✗ Branch 2 not taken.
 338 i_ = -other.i_;
280
2/4
✓ Branch 1 taken 338 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 338 times.
✗ Branch 6 not taken.
 338 j_ = -other.j_, k_ = -other.k_;
281 } else {
282 ✗ rep_ = Rep::quat;
283 ✗ s_ = -other.s_;
284 ✗ i_ = -other.i_;
285 ✗ j_ = -other.j_, k_ = -other.k_;
286 }
287 343 return;
288 }
289
290
4/4
✓ Branch 0 taken 19205 times.
✓ Branch 1 taken 7492 times.
✓ Branch 2 taken 19112 times.
✓ Branch 3 taken 93 times.
2/2
✓ Decision 'true' taken 19112 times.
✓ Decision 'false' taken 7585 times.
 26697 if (rep_ == Rep::orth_rot && other.rep_ == Rep::orth_rot) {
291
2/2
✓ Branch 0 taken 14096 times.
✓ Branch 1 taken 5016 times.
2/2
✓ Decision 'true' taken 14096 times.
✓ Decision 'false' taken 5016 times.
 19112 if (optype_ == other.optype_) {
292
1/2
✓ Branch 1 taken 14096 times.
✗ Branch 2 not taken.
 14096 a_ += other.a_;
293
3/4
✓ Branch 1 taken 14096 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3513 times.
✓ Branch 4 taken 10583 times.
2/2
✓ Decision 'true' taken 3513 times.
✓ Decision 'false' taken 10583 times.
 14096 if (equiv_0(a_, 4)) {
294 3513 rep_ = Rep::id;
295
3/4
✓ Branch 1 taken 10583 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 854 times.
✓ Branch 4 taken 9729 times.
2/2
✓ Decision 'true' taken 854 times.
✓ Decision 'false' taken 9729 times.
 10583 } else if (equiv_0(a_, 2)) {
296 854 rep_ = Rep::minus_id;
297 }
298 5016 } else if (
299
8/10
✓ Branch 1 taken 5016 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4811 times.
✓ Branch 4 taken 205 times.
✓ Branch 6 taken 4811 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 756 times.
✓ Branch 9 taken 4055 times.
✓ Branch 10 taken 199 times.
✓ Branch 11 taken 4817 times.
 5977 (equiv_val(a_, 1., 4) || equiv_val(a_, -1., 4)) &&
300
6/8
✓ Branch 1 taken 961 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 949 times.
✓ Branch 4 taken 12 times.
✓ Branch 6 taken 949 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 187 times.
✓ Branch 9 taken 762 times.
 961 (equiv_val(other.a_, 1., 4) || equiv_val(other.a_, -1., 4))) {
301 // We are in a subgroup of order 8
302
3/4
✓ Branch 1 taken 199 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
✓ Branch 4 taken 181 times.
 199 int m0 = equiv_val(a_, 1., 4) ? 1 : -1;
303
3/4
✓ Branch 1 taken 199 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 187 times.
 199 int m1 = equiv_val(other.a_, 1., 4) ? 1 : -1;
304
2/4
✓ Branch 2 taken 199 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 199 times.
✗ Branch 7 not taken.
 199 std::tie(optype_, a_) = product.at({optype_, m0, other.optype_, m1});
305 } else
306 4817 rep_ = Rep::quat;
307 19112 } else
308 7585 rep_ = Rep::quat;
309
310
7/14
✓ Branch 1 taken 26697 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26697 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26697 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 26697 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 26697 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 26697 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 26697 times.
✗ Branch 20 not taken.
 53394 Expr s1 = other.s_ * s_ - other.i_ * i_ - other.j_ * j_ - other.k_ * k_;
311
7/14
✓ Branch 1 taken 26697 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26697 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26697 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 26697 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 26697 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 26697 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 26697 times.
✗ Branch 20 not taken.
 53394 Expr i1 = other.s_ * i_ + other.i_ * s_ + other.j_ * k_ - other.k_ * j_;
312
7/14
✓ Branch 1 taken 26697 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26697 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26697 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 26697 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 26697 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 26697 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 26697 times.
✗ Branch 20 not taken.
 53394 Expr j1 = other.s_ * j_ - other.i_ * k_ + other.j_ * s_ + other.k_ * i_;
313
7/14
✓ Branch 1 taken 26697 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26697 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26697 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 26697 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 26697 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 26697 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 26697 times.
✗ Branch 20 not taken.
 53394 Expr k1 = other.s_ * k_ + other.i_ * j_ - other.j_ * i_ + other.k_ * s_;
314
1/2
✓ Branch 1 taken 26697 times.
✗ Branch 2 not taken.
 26697 s_ = SymEngine::expand(s1);
315
1/2
✓ Branch 1 taken 26697 times.
✗ Branch 2 not taken.
 26697 i_ = SymEngine::expand(i1);
316
1/2
✓ Branch 1 taken 26697 times.
✗ Branch 2 not taken.
 26697 j_ = SymEngine::expand(j1);
317
1/2
✓ Branch 1 taken 26697 times.
✗ Branch 2 not taken.
 26697 k_ = SymEngine::expand(k1);
318
319
2/2
✓ Branch 0 taken 12402 times.
✓ Branch 1 taken 14295 times.
2/2
✓ Decision 'true' taken 12402 times.
✓ Decision 'false' taken 14295 times.
 26697 if (rep_ == Rep::quat) {
320 // See if we can simplify the representation.
321
1/2
✓ Branch 1 taken 12402 times.
✗ Branch 2 not taken.
 12402 bool i_zero = approx_0(i_);
322
1/2
✓ Branch 1 taken 12402 times.
✗ Branch 2 not taken.
 12402 bool j_zero = approx_0(j_);
323
1/2
✓ Branch 1 taken 12402 times.
✗ Branch 2 not taken.
 12402 bool k_zero = approx_0(k_);
324
5/6
✓ Branch 0 taken 2347 times.
✓ Branch 1 taken 10055 times.
✓ Branch 2 taken 491 times.
✓ Branch 3 taken 1856 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 491 times.
1/2
✗ Decision 'true' not taken.
✓ Decision 'false' taken 12402 times.
 12402 if (i_zero && j_zero && k_zero) {
325
0/2
✗ Decision 'true' not taken.
✗ Decision 'false' not taken.
 ✗ if (approx_0(s_ - 1)) {
326 ✗ rep_ = Rep::id;
327 ✗ s_ = 1;
328 ✗ i_ = j_ = k_ = 0;
329 } else {
330 ✗ rep_ = Rep::minus_id;
331 ✗ s_ = -1;
332 ✗ i_ = j_ = k_ = 0;
333 }
334
4/4
✓ Branch 0 taken 2282 times.
✓ Branch 1 taken 10120 times.
✓ Branch 2 taken 804 times.
✓ Branch 3 taken 1478 times.
2/2
✓ Decision 'true' taken 804 times.
✓ Decision 'false' taken 11598 times.
 12402 } else if (j_zero && k_zero) {
335 804 rep_ = Rep::orth_rot;
336 804 optype_ = OpType::Rx;
337
3/6
✓ Branch 1 taken 804 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 804 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 804 times.
✗ Branch 8 not taken.
 804 a_ = 2 * atan2_bypi(i_, s_);
338
2/4
✓ Branch 1 taken 804 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 804 times.
✗ Branch 6 not taken.
 804 j_ = k_ = 0;
339
4/4
✓ Branch 0 taken 1462 times.
✓ Branch 1 taken 10136 times.
✓ Branch 2 taken 373 times.
✓ Branch 3 taken 1089 times.
2/2
✓ Decision 'true' taken 373 times.
✓ Decision 'false' taken 11225 times.
 11598 } else if (k_zero && i_zero) {
340 373 rep_ = Rep::orth_rot;
341 373 optype_ = OpType::Ry;
342
3/6
✓ Branch 1 taken 373 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 373 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 373 times.
✗ Branch 8 not taken.
 373 a_ = 2 * atan2_bypi(j_, s_);
343
2/4
✓ Branch 1 taken 373 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 373 times.
✗ Branch 6 not taken.
 373 k_ = i_ = 0;
344
4/4
✓ Branch 0 taken 1974 times.
✓ Branch 1 taken 9251 times.
✓ Branch 2 taken 491 times.
✓ Branch 3 taken 1483 times.
2/2
✓ Decision 'true' taken 491 times.
✓ Decision 'false' taken 10734 times.
 11225 } else if (i_zero && j_zero) {
345 491 rep_ = Rep::orth_rot;
346 491 optype_ = OpType::Rz;
347
3/6
✓ Branch 1 taken 491 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 491 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 491 times.
✗ Branch 8 not taken.
 491 a_ = 2 * atan2_bypi(k_, s_);
348
2/4
✓ Branch 1 taken 491 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 491 times.
✗ Branch 6 not taken.
 491 i_ = j_ = 0;
349 }
350 }
351 26697 }
352
353 ✗ std::ostream &operator<<(std::ostream &os, const Rotation &q) {
354
0/2
✗ Decision 'true' not taken.
✗ Decision 'false' not taken.
 ✗ if (q.rep_ == Rotation::Rep::id) {
355 ✗ return os << "I";
356
0/2
✗ Decision 'true' not taken.
✗ Decision 'false' not taken.
 ✗ } else if (q.rep_ == Rotation::Rep::minus_id) {
357 ✗ return os << "-I";
358
0/2
✗ Decision 'true' not taken.
✗ Decision 'false' not taken.
 ✗ } else if (q.rep_ == Rotation::Rep::orth_rot) {
359 ✗ return os << OpDesc(q.optype_).name() << "(" << q.a_ << ")";
360 } else {
361 ✗ return os << q.s_ << " + " << q.i_ << " i + " << q.j_ << " j + " << q.k_
362 ✗ << " k";
363 }
364 }
365
366 15751 std::vector<double> tk1_angles_from_unitary(const Eigen::Matrix2cd &U) {
367 // clang-format off
368 //
369 // Assume U = e^{i pi p} TK1(a,b,c)
370 // = e^{i pi p} Rz(a) Rx(b) Rz(c)
371 // | e^{-i pi (a+c)/2} cos(pi b/2) -i e^{-i pi (a-c)/2} sin(pi b/2) |
372 // = e^{i pi p} | |
373 // | -i e^{ i pi (a-c)/2} sin(pi b/2) e^{ i pi (a+c)/2} cos(pi b/2) |
374 //
375 // clang-format on
376
377 double a, b, c, p; // to be found satisfying the above equation
378
379
3/6
✓ Branch 1 taken 15751 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15751 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 15751 times.
✗ Branch 8 not taken.
 15751 std::complex s = 0.5 * (U(0, 0) + U(1, 1));
380
4/8
✓ Branch 1 taken 15751 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15751 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 15751 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 15751 times.
✗ Branch 12 not taken.
 15751 std::complex x = 0.5 * i_ * (U(1, 0) + U(0, 1));
381
3/6
✓ Branch 1 taken 15751 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15751 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 15751 times.
✗ Branch 8 not taken.
 15751 std::complex y = 0.5 * (U(1, 0) - U(0, 1));
382
4/8
✓ Branch 1 taken 15751 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15751 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 15751 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 15751 times.
✗ Branch 12 not taken.
 15751 std::complex z = 0.5 * i_ * (U(0, 0) - U(1, 1));
383
384 // s = e^{i pi p} cos(pi b/2) cos(pi (a+c)/2)
385 // x = e^{i pi p} sin(pi b/2) cos(pi (a-c)/2)
386 // y = e^{i pi p} sin(pi b/2) sin(pi (a-c)/2)
387 // z = e^{i pi p} cos(pi b/2) sin(pi (a+c)/2)
388
389 // s, x, y and z all have phase e^{i pi p}. Extract it from the one with
390 // largest absolute value, to minimize numerical instability. Note that
391 // there are two possible values (p and p+1), but the choice is w.l.o.g.
392 // because (p,b) --> (p+1,b+2) does not change the value of the unitary.
393
394 15751 std::complex w = s;
395
2/2
✓ Branch 2 taken 5854 times.
✓ Branch 3 taken 9897 times.
1/2
✓ Decision 'true' taken 15751 times.
✗ Decision 'false' not taken.
 15751 if (std::abs(x) > std::abs(w)) w = x;
396
2/2
✓ Branch 2 taken 3779 times.
✓ Branch 3 taken 11972 times.
1/2
✓ Decision 'true' taken 15751 times.
✗ Decision 'false' not taken.
 15751 if (std::abs(y) > std::abs(w)) w = y;
397
2/2
✓ Branch 2 taken 2737 times.
✓ Branch 3 taken 13014 times.
1/2
✓ Decision 'true' taken 15751 times.
✗ Decision 'false' not taken.
 15751 if (std::abs(z) > std::abs(w)) w = z;
398 15751 std::complex eip = w / std::abs(w);
399
400 // eip = e^{i pi p}
401 15751 p = std::arg(eip) / PI;
402
403 // Now we've got the phase, factor it out.
404
405 15751 std::complex emip = std::conj(eip);
406
1/2
✓ Branch 1 taken 15751 times.
✗ Branch 2 not taken.
 15751 double s0 = std::real(emip * s);
407
1/2
✓ Branch 1 taken 15751 times.
✗ Branch 2 not taken.
 15751 double x0 = std::real(emip * x);
408
1/2
✓ Branch 1 taken 15751 times.
✗ Branch 2 not taken.
 15751 double y0 = std::real(emip * y);
409
1/2
✓ Branch 1 taken 15751 times.
✗ Branch 2 not taken.
 15751 double z0 = std::real(emip * z);
410
411 // s0 = cos(pi b/2) cos(pi (a+c)/2)
412 // x0 = sin(pi b/2) cos(pi (a-c)/2)
413 // y0 = sin(pi b/2) sin(pi (a-c)/2)
414 // z0 = cos(pi b/2) sin(pi (a+c)/2)
415
416 15751 std::complex u = {s0, z0}, v = {x0, y0};
417
418 // u = cos(pi b/2) e^{i pi (a+c)/2}
419 // v = sin(pi b/2) e^{i pi (a-c)/2}
420 // |u|^2 + |v|^2 = 1
421
422 // Note that (a,b) --> (a+2, b+2) does not change the value of the unitary, so
423 // we are free to choose either solution.
424 //
425 // We treat the two special cases u=0 and v=0 separately, then the general
426 // case.
427
428
2/2
✓ Branch 1 taken 484 times.
✓ Branch 2 taken 15267 times.
2/2
✓ Decision 'true' taken 484 times.
✓ Decision 'false' taken 15267 times.
 15751 if (std::abs(u) < EPS) { // special case, b = 1 or 3
429 // v = +/- e^{i pi (a-c)/2}
430
431 // We may as well choose c = 0.
432 484 c = 0;
433
434 // Assume b=1 and compute a.
435 // (b'=3, a'=a+2) is the other possibility but the unitary is the same.
436 484 b = 1;
437 484 a = 2 * std::arg(v) / PI;
438
2/2
✓ Branch 1 taken 1573 times.
✓ Branch 2 taken 13694 times.
2/2
✓ Decision 'true' taken 1573 times.
✓ Decision 'false' taken 13694 times.
 15267 } else if (std::abs(v) < EPS) { // special case, b = 0 or 2
439 // u = e^{i pi (a+c)/2}
440
441 // We may as well choose c = 0.
442 1573 c = 0;
443
444 // Assume b=0 and compute a.
445 // (b'=2, a'=a+2) is the other possibility but the unitary is the same.
446 1573 b = 0;
447 1573 a = 2 * std::arg(u) / PI;
448 } else { // general case
449 // s0^2 + z0^2 - x0^2 - y0^2 = cos(pi b)
450 13694 double t = s0 * s0 + z0 * z0 - x0 * x0 - y0 * y0;
451 // Rounding errors may mean t is outside the domain of acos. Fix this.
452
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13694 times.
1/2
✓ Decision 'true' taken 13694 times.
✗ Decision 'false' not taken.
 13694 if (t > +1.) t = +1.;
453
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13694 times.
1/2
✓ Decision 'true' taken 13694 times.
✗ Decision 'false' not taken.
 13694 if (t < -1.) t = -1.;
454 13694 b = std::acos(t) / PI;
455
456 // w.l.o.g. b is in the range (-1,+1).
457 13694 double ac0 = std::arg(u), ac1 = std::arg(v);
458 13694 a = (ac0 + ac1) / PI;
459 13694 c = (ac0 - ac1) / PI;
460 }
461
462
1/2
✓ Branch 2 taken 15751 times.
✗ Branch 3 not taken.
 15751 return {a, b, c, p};
463 }
464
465 49017 Eigen::Matrix2cd get_matrix_from_tk1_angles(std::vector<Expr> params) {
466
2/4
✓ Branch 2 taken 49017 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 49017 times.
✗ Branch 6 not taken.
 49017 double alpha = eval_expr(params[0]).value();
467
2/4
✓ Branch 2 taken 49017 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 49017 times.
✗ Branch 6 not taken.
 49017 double beta = eval_expr(params[1]).value();
468
2/4
✓ Branch 2 taken 49017 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 49017 times.
✗ Branch 6 not taken.
 49017 double gamma = eval_expr(params[2]).value();
469
2/4
✓ Branch 2 taken 49017 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 49017 times.
✗ Branch 6 not taken.
 49017 double t = eval_expr(params[3]).value();
470
1/2
✓ Branch 1 taken 49017 times.
✗ Branch 2 not taken.
 49017 Eigen::Matrix2cd m;
471 49017 alpha *= PI;
472 49017 beta *= PI;
473 49017 gamma *= PI;
474 49017 t *= PI;
475 49017 double c = cos(0.5 * beta);
476 49017 double s = sin(0.5 * beta);
477
1/2
✓ Branch 5 taken 49017 times.
✗ Branch 6 not taken.
 49017 m << exp(-0.5 * i_ * (alpha + gamma)) * c,
478
2/4
✓ Branch 5 taken 49017 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 49017 times.
✗ Branch 10 not taken.
 49017 -i_ * exp(0.5 * i_ * (gamma - alpha)) * s,
479
2/4
✓ Branch 5 taken 49017 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 49017 times.
✗ Branch 10 not taken.
 49017 -i_ * exp(0.5 * i_ * (alpha - gamma)) * s,
480
1/2
✓ Branch 5 taken 49017 times.
✗ Branch 6 not taken.
 49017 exp(0.5 * i_ * (alpha + gamma)) * c;
481
2/4
✓ Branch 3 taken 49017 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 49017 times.
✗ Branch 7 not taken.
 98034 return exp(i_ * t) * m;
482 }
483
484 } // namespace tket
485