|
| 1 | +import math |
| 2 | + |
| 3 | +from sympy import Quaternion, pi, cos, sin, N |
| 4 | + |
| 5 | + |
| 6 | +# ---- Gate → quaternion (SU(2) part) ---- |
| 7 | + |
| 8 | + |
| 9 | +def r_gate(theta, phi): |
| 10 | + return Quaternion( |
| 11 | + cos(theta / 2), sin(theta / 2) * cos(phi), sin(theta / 2) * sin(phi), 0 |
| 12 | + ) |
| 13 | + |
| 14 | + |
| 15 | +def u2_gate(phi, lam): |
| 16 | + return Quaternion.from_euler([phi, pi / 2, lam], "ZYZ") |
| 17 | + |
| 18 | + |
| 19 | +# ---- Euler extraction matching the C++ pass exactly ---- |
| 20 | + |
| 21 | + |
| 22 | +def normalize_angle(a): |
| 23 | + """Normalize angle to [-pi, pi], matching the pass's normalizeAngle.""" |
| 24 | + two_pi = 2 * math.pi |
| 25 | + return a - math.floor((a + math.pi) / two_pi) * two_pi |
| 26 | + |
| 27 | + |
| 28 | +def angles_from_quaternion(w, x, y, z): |
| 29 | + """ZYZ Euler angles from quaternion, matching anglesFromQuaternion in the pass. |
| 30 | +
|
| 31 | + Does NOT normalize the quaternion sign — uses the raw Hamilton product |
| 32 | + result, just like the C++ code. Returns U gate parameters (theta, phi, lambda). |
| 33 | + """ |
| 34 | + eps = 1e-12 |
| 35 | + |
| 36 | + # beta = acos(clamp(2*(w^2 + z^2) - 1, -1, 1)) |
| 37 | + arg = 2 * (w * w + z * z) - 1 |
| 38 | + arg = max(-1.0, min(1.0, arg)) |
| 39 | + beta = math.acos(arg) |
| 40 | + |
| 41 | + abs_beta = abs(beta) |
| 42 | + abs_beta_minus_pi = abs(beta - math.pi) |
| 43 | + |
| 44 | + safe1 = abs_beta >= eps # not near 0 |
| 45 | + safe2 = abs_beta_minus_pi >= eps # not near pi |
| 46 | + safe = safe1 and safe2 |
| 47 | + |
| 48 | + theta_plus = math.atan2(z, w) |
| 49 | + theta_minus = math.atan2(-x, y) |
| 50 | + |
| 51 | + if safe: |
| 52 | + alpha = theta_plus + theta_minus |
| 53 | + gamma = theta_plus - theta_minus |
| 54 | + elif not safe1: |
| 55 | + # beta near 0 |
| 56 | + alpha = 2 * theta_plus |
| 57 | + gamma = 0.0 |
| 58 | + else: |
| 59 | + # beta near pi |
| 60 | + alpha = 2 * theta_minus |
| 61 | + gamma = 0.0 |
| 62 | + |
| 63 | + alpha = normalize_angle(alpha) |
| 64 | + gamma = normalize_angle(gamma) |
| 65 | + |
| 66 | + # U gate convention: theta=beta, phi=alpha, lambda=gamma |
| 67 | + return beta, alpha, gamma |
| 68 | + |
| 69 | + |
| 70 | +# ---- Global phase per gate type ---- |
| 71 | + |
| 72 | + |
| 73 | +def global_phase(gate_type, *angles): |
| 74 | + """Returns the global phase contribution of a gate. |
| 75 | +
|
| 76 | + U = e^{i*phase} * SU(2), this returns 'phase'. |
| 77 | + """ |
| 78 | + if gate_type in ("RX", "RY", "RZ", "R"): |
| 79 | + return 0 |
| 80 | + elif gate_type == "P": |
| 81 | + return angles[0] / 2 |
| 82 | + elif gate_type == "U": |
| 83 | + # U(theta, phi, lambda): phase = (phi + lambda) / 2 |
| 84 | + theta, phi, lam = angles |
| 85 | + return (phi + lam) / 2 |
| 86 | + elif gate_type == "U2": |
| 87 | + # U2(phi, lambda): phase = (phi + lambda) / 2 |
| 88 | + phi, lam = angles |
| 89 | + return (phi + lam) / 2 |
| 90 | + |
| 91 | + |
| 92 | +def output_phase(phi, lam): |
| 93 | + """Intrinsic phase of the synthesized U(theta, phi, lambda).""" |
| 94 | + return (phi + lam) / 2 |
| 95 | + |
| 96 | + |
| 97 | +def gphase_correction(input_phase, phi, lam): |
| 98 | + """GPhaseOp correction = total_input_phase - output_UOp_phase.""" |
| 99 | + return input_phase - output_phase(phi, lam) |
| 100 | + |
| 101 | + |
| 102 | +# ---- Helper to compute merge + gphase for a chain ---- |
| 103 | + |
| 104 | + |
| 105 | +def compute_merge(chain): |
| 106 | + """ |
| 107 | + chain: list of (gate_type, quaternion, *angles) |
| 108 | + Returns (theta, phi, lambda, gphase) all as floats. |
| 109 | +
|
| 110 | + Uses our own Euler extraction that matches the C++ pass exactly: |
| 111 | + no quaternion sign normalization, same atan2/acos/clamp logic, |
| 112 | + same gimbal-lock handling, same angle normalization. |
| 113 | + """ |
| 114 | + _, q0, *a0 = chain[0] |
| 115 | + q = q0 |
| 116 | + total_input_phase = global_phase(chain[0][0], *a0) |
| 117 | + |
| 118 | + for entry in chain[1:]: |
| 119 | + gt, qi, *ai = entry |
| 120 | + q = qi.mul(q) # Hamilton product in circuit order |
| 121 | + total_input_phase += global_phase(gt, *ai) |
| 122 | + |
| 123 | + # Extract Euler angles matching the pass (no sign normalization) |
| 124 | + w, x, y, z = float(N(q.a)), float(N(q.b)), float(N(q.c)), float(N(q.d)) |
| 125 | + theta, phi, lam = angles_from_quaternion(w, x, y, z) |
| 126 | + |
| 127 | + corr = gphase_correction(total_input_phase, phi, lam) |
| 128 | + |
| 129 | + return theta, phi, lam, float(N(corr)) |
| 130 | + |
| 131 | + |
| 132 | +# ---- Build gates ---- |
| 133 | + |
| 134 | +rx = Quaternion.from_euler([1, 0, 0], "xyz") |
| 135 | +ry = Quaternion.from_euler([0, 1, 0], "xyz") |
| 136 | +rz = Quaternion.from_euler([0, 0, 1], "xyz") |
| 137 | +mx = Quaternion.from_euler([-1, 0, 0], "xyz") |
| 138 | +my = Quaternion.from_euler([0, -1, 0], "xyz") |
| 139 | +mz = Quaternion.from_euler([0, 0, -1], "xyz") |
| 140 | +px = Quaternion.from_euler([pi, 0, 0], "xyz") |
| 141 | +py = Quaternion.from_euler([0, pi, 0], "xyz") |
| 142 | +pz = Quaternion.from_euler([0, 0, pi], "xyz") |
| 143 | +smallx = Quaternion.from_euler([0.001, 0, 0], "xyz") |
| 144 | +smally = Quaternion.from_euler([0, 0.001, 0], "xyz") |
| 145 | + |
| 146 | +# P gate has same SU(2) quaternion as RZ |
| 147 | +p1 = Quaternion.from_euler([0, 0, 1], "xyz") # P(1) same rotation as RZ(1) |
| 148 | + |
| 149 | +u1 = Quaternion.from_euler([2, 1, 3], "ZYZ") # U(1,2,3) |
| 150 | +u2 = Quaternion.from_euler([5, 4, 6], "ZYZ") # U(4,5,6) |
| 151 | + |
| 152 | +u2_12 = u2_gate(1, 2) |
| 153 | +u2_34 = u2_gate(3, 4) |
| 154 | + |
| 155 | +r12 = r_gate(1, 2) |
| 156 | +r34 = r_gate(3, 4) |
| 157 | +r11 = r_gate(1, 1) |
| 158 | + |
| 159 | +cases = [ |
| 160 | + ("RX+RY", [("RX", rx), ("RY", ry)]), |
| 161 | + ("RX+RZ", [("RX", rx), ("RZ", rz)]), |
| 162 | + ("RY+RX", [("RY", ry), ("RX", rx)]), |
| 163 | + ("RY+RZ", [("RY", ry), ("RZ", rz)]), |
| 164 | + ("RZ+RX", [("RZ", rz), ("RX", rx)]), |
| 165 | + ("RZ+RY", [("RZ", rz), ("RY", ry)]), |
| 166 | + ("RX+RX", [("RX", rx), ("RX", rx)]), |
| 167 | + ("RY+RY", [("RY", ry), ("RY", ry)]), |
| 168 | + ("RZ+RZ", [("RZ", rz), ("RZ", rz)]), |
| 169 | + ("U+U", [("U", u1, 1.0, 2.0, 3.0), ("U", u2, 4.0, 5.0, 6.0)]), |
| 170 | + ("P+RX", [("P", p1, 1.0), ("RX", rx)]), |
| 171 | + ("U2+U2", [("U2", u2_12, 1.0, 2.0), ("U2", u2_34, 3.0, 4.0)]), |
| 172 | + ("R+R", [("R", r12, 1.0, 2.0), ("R", r34, 3.0, 4.0)]), |
| 173 | + ("R+R same", [("R", r11, 1.0, 1.0), ("R", r11, 1.0, 1.0)]), |
| 174 | + ("small RX+RY", [("RX", smallx), ("RY", smally)]), |
| 175 | + ("RX(pi)+RY(pi)", [("RX", px), ("RY", py)]), |
| 176 | + ("RZ+RY+RX pi", [("RZ", pz), ("RY", py), ("RX", px)]), |
| 177 | + ("RY+RZ+RZ-+RY-", [("RY", ry), ("RZ", rz), ("RZ", mz), ("RY", my)]), |
| 178 | +] |
| 179 | + |
| 180 | +if __name__ == "__main__": |
| 181 | + for name, chain in cases: |
| 182 | + theta, phi, lam, gphase = compute_merge(chain) |
| 183 | + print(f"{name}: U({theta}, {phi}, {lam}) gphase={gphase}") |
0 commit comments