From d2f4b35f00a41d3879901f11e7c915b65d2ff6ea Mon Sep 17 00:00:00 2001 From: shun0923 Date: Thu, 7 Aug 2025 22:56:05 +0900 Subject: [PATCH 1/3] Update files --- pygridsynth/gridsynth.py | 47 +++-- pygridsynth/normal_form.py | 74 +++++-- pygridsynth/odgp.py | 33 ++-- pygridsynth/quantum_gate.py | 270 ++++++++++++++++++++++++++ pygridsynth/synthesis_of_cliffordT.py | 51 +++-- pygridsynth/tdgp.py | 63 +++--- pygridsynth/to_upright.py | 2 +- pygridsynth/unitary.py | 13 +- 8 files changed, 444 insertions(+), 109 deletions(-) create mode 100644 pygridsynth/quantum_gate.py diff --git a/pygridsynth/gridsynth.py b/pygridsynth/gridsynth.py index 8ea87d5..7def454 100644 --- a/pygridsynth/gridsynth.py +++ b/pygridsynth/gridsynth.py @@ -9,20 +9,21 @@ from .diophantine import NO_SOLUTION, diophantine_dyadic from .unitary import DOmegaUnitary from .synthesis_of_cliffordT import decompose_domega_unitary +from .quantum_gate import Rz class EpsilonRegion(ConvexSet): def __init__(self, theta, epsilon): self._theta = theta self._epsilon = epsilon - self._d = 1 - epsilon ** 2 / 2 + self._d = sqrt(1 - epsilon ** 2 / 4) self._z_x = mpmath.cos(-theta / 2) self._z_y = mpmath.sin(-theta / 2) D_1 = mpmath.matrix([[self._z_x, -self._z_y], [self._z_y, self._z_x]]) - D_2 = mpmath.matrix([[4 * (1 / epsilon) ** 4, 0], [0, (1 / epsilon) ** 2]]) + D_2 = mpmath.matrix([[64 * (1 / epsilon) ** 4, 0], [0, 4 * (1 / epsilon) ** 2]]) D_3 = mpmath.matrix([[self._z_x, self._z_y], [-self._z_y, self._z_x]]) p = mpmath.matrix([self._d * self._z_x, self._d * self._z_y]) - ellipse = Ellipse(D_1 * D_2 * D_3, p) + ellipse = Ellipse(D_1 @ D_2 @ D_3, p) super().__init__(ellipse) @property @@ -79,36 +80,28 @@ def generate_complex_unitary(sol): [t.to_complex, u.conj.to_complex]]) -def generate_target_Rz(theta): - return mpmath.matrix([[mpmath.exp(- 1.j * theta / 2), 0], - [0, mpmath.exp(1.j * theta / 2)]]) - - def error(theta, gates): - Rz = generate_target_Rz(theta) - U = DOmegaUnitary.from_gates(gates).to_complex_matrix - E = U - Rz - return sqrt(mpmath.fabs(E[0, 0] * E[1, 1] - E[0, 1] * E[1, 0])) + u_target = Rz(theta) + u_approx = DOmegaUnitary.from_gates(gates).to_complex_matrix + return 2 * sqrt(mpmath.re(1 - mpmath.re(mpmath.conj(u_target[0, 0]) * u_approx[0, 0]) ** 2)) def check(theta, gates): - t_count = gates.count("T") - h_count = gates.count("H") - U_decomp = DOmegaUnitary.from_gates(gates) - # Rz = generate_target_Rz(theta) - # U = U_decomp.to_complex_matrix + # t_count = gates.count("T") + # h_count = gates.count("H") + u_approx = DOmegaUnitary.from_gates(gates) e = error(theta, gates) print(f"{gates=}") - print(f"{t_count=}, {h_count=}") - # print(f"{Rz=}") - print(f"U_decomp={U_decomp.to_matrix}") - # print(f"{U=}") + # print(f"{t_count=}, {h_count=}") + print(f"u_approx={u_approx.to_matrix}") print(f"{e=}") def gridsynth(theta, epsilon, diophantine_timeout=200, factoring_timeout=50, verbose=False, measure_time=False, show_graph=False): + theta = mpmath.mpmathify(theta) + epsilon = mpmath.mpmathify(epsilon) epsilon_region = EpsilonRegion(theta, epsilon) unit_disk = UnitDisk() k = 0 @@ -132,7 +125,6 @@ def gridsynth(theta, epsilon, if measure_time: time_of_solve_TDGP += time.time() - start start = time.time() - for z in sol: if (z * z.conj).residue == 0: continue @@ -164,9 +156,12 @@ def gridsynth(theta, epsilon, k += 1 -def gridsynth_gates(theta, epsilon, +def gridsynth_gates(theta, epsilon, wires=[0], + decompose_phase_gate=True, diophantine_timeout=200, factoring_timeout=50, verbose=False, measure_time=False, show_graph=False): + theta = mpmath.mpmathify(theta) + epsilon = mpmath.mpmathify(epsilon) if measure_time: start_total = time.time() u_approx = gridsynth(theta=theta, epsilon=epsilon, @@ -175,8 +170,10 @@ def gridsynth_gates(theta, epsilon, verbose=verbose, measure_time=measure_time, show_graph=show_graph) if measure_time: start = time.time() - gates = decompose_domega_unitary(u_approx) + circuit = decompose_domega_unitary(u_approx, wires=wires, decompose_phase_gate=decompose_phase_gate) + if measure_time: print(f"time of decompose_domega_unitary: {(time.time() - start) * 1000} ms") print(f"total time: {(time.time() - start_total) * 1000} ms") - return gates + + return circuit diff --git a/pygridsynth/normal_form.py b/pygridsynth/normal_form.py index 1e8373a..c3dd9e6 100644 --- a/pygridsynth/normal_form.py +++ b/pygridsynth/normal_form.py @@ -1,4 +1,5 @@ from enum import Enum +from .quantum_gate import QuantumCircuit, HGate, TGate, SGate, SXGate, WGate class Axis(Enum): @@ -82,6 +83,19 @@ def from_str(cls, g): else: raise ValueError + @classmethod + def from_gate(cls, g): + if isinstance(g, HGate): + return CLIFFORD_H + elif isinstance(g, SGate): + return CLIFFORD_S + elif isinstance(g, SXGate): + return CLIFFORD_X + elif isinstance(g, WGate): + return CLIFFORD_W + else: + raise ValueError + def __eq__(self, other): if isinstance(other, self.__class__): return (self._a == other.a and self._b == other.b @@ -133,15 +147,28 @@ def decompose_tconj(self): axis, c1, d1 = self.__class__._tconj(self._a, self._b) return axis, self.__class__(0, self._b, c1 + self._c, d1 + self._d) - def to_gates(self): + def to_circuit(self, wires): axis, c = self.decompose_coset() - return ("" if axis == Axis.I else axis.name) + "X" * c.b + "S" * c.c + "W" * c.d + circuit = QuantumCircuit() + if axis == Axis.H: + circuit.append(HGate(target_qubit=wires[0])) + elif axis == Axis.SH: + circuit.append(SGate(target_qubit=wires[0])) + circuit.append(HGate(target_qubit=wires[0])) + for _ in range(c.b): + circuit.append(SXGate(target_qubit=wires[0])) + for _ in range(c.c): + circuit.append(SGate(target_qubit=wires[0])) + for _ in range(c.d): + circuit.append(WGate()) + return circuit class NormalForm(): - def __init__(self, syllables, c): + def __init__(self, syllables, c, phase=0): self._syllables = syllables self._c = c + self._phase = phase @property def syllables(self): @@ -155,13 +182,19 @@ def c(self): def c(self, c): self._c = c + @property + def phase(self): + return self._phase + + @phase.setter + def phase(self, phase): + self._phase = phase + def __repr__(self): - return f"NormalForm({repr(self._syllables)}, {repr(self._c)})" + return f"NormalForm({repr(self._syllables)}, {repr(self._c)}, {repr(self._phase)})" def _append_gate(self, g): - if g in ["H", "S", "X", "W"]: - self.c *= Clifford.from_str(g) - elif g == "T": + if isinstance(g, TGate): axis, new_c = self.c.decompose_tconj() if axis == Axis.I: if len(self._syllables) == 0: @@ -181,21 +214,30 @@ def _append_gate(self, g): elif axis == Axis.SH: self._syllables.append(Syllable.SHT) self.c = new_c + else: + self.c *= Clifford.from_gate(g) @classmethod - def from_gates(cls, gates): - normal_form = NormalForm([], CLIFFORD_I) - for g in gates: + def from_circuit(cls, circuit): + normal_form = NormalForm([], CLIFFORD_I, phase=circuit.phase) + for g in circuit: normal_form._append_gate(g) return normal_form - def to_gates(self): - gates = "" + def to_circuit(self, wires): + circuit = QuantumCircuit(phase=self.phase) for syllable in self._syllables: - if syllable != Syllable.I: - gates += syllable.name - gates += self._c.to_gates() - return "I" if gates == "" else gates + if syllable == Syllable.T: + circuit.append(TGate(target_qubit=wires[0])) + elif syllable == Syllable.HT: + circuit.append(HGate(target_qubit=wires[0])) + circuit.append(TGate(target_qubit=wires[0])) + elif syllable == Syllable.SHT: + circuit.append(SGate(target_qubit=wires[0])) + circuit.append(HGate(target_qubit=wires[0])) + circuit.append(TGate(target_qubit=wires[0])) + circuit += self._c.to_circuit(wires=wires) + return circuit CLIFFORD_I = Clifford(0, 0, 0, 0) diff --git a/pygridsynth/odgp.py b/pygridsynth/odgp.py index 2078332..ed0236c 100644 --- a/pygridsynth/odgp.py +++ b/pygridsynth/odgp.py @@ -1,66 +1,67 @@ +import itertools + from .mymath import SQRT2, floor, ceil, pow_sqrt2, floorlog from .ring import ZRootTwo, DRootTwo, LAMBDA def _solve_ODGP_internal(I, J): if I.width < 0 or J.width < 0: - return [] + return iter([]) elif I.width > 0 and J.width <= 0: sol = _solve_ODGP_internal(J, I) - return [beta.conj_sq2 for beta in sol] + return map(lambda beta: beta.conj_sq2, sol) else: (n, _) = (0, 0) if J.width <= 0 else floorlog(J.width, LAMBDA.to_real) if n == 0: - sol = [] a_min = ceil((I.l + J.l) / 2) a_max = floor((I.r + J.r) / 2) - for a in range(a_min, a_max + 1): + + def gen_sol(a): b_min = ceil(SQRT2() * (a - J.r) / 2) b_max = floor(SQRT2() * (a - J.l) / 2) - for b in range(b_min, b_max + 1): - sol.append(ZRootTwo(a, b)) - return sol + return map(lambda b: ZRootTwo(a, b), range(b_min, b_max + 1)) + + return itertools.chain.from_iterable(map(gen_sol, range(a_min, a_max + 1))) else: lambda_n = LAMBDA ** n lambda_inv_n = LAMBDA ** -n lambda_conj_sq2_n = LAMBDA.conj_sq2 ** n sol = _solve_ODGP_internal(I * lambda_n.to_real, J * lambda_conj_sq2_n.to_real) - sol = [beta * lambda_inv_n for beta in sol] - return sol + return map(lambda beta: beta * lambda_inv_n, sol) def solve_ODGP(I, J): if I.width < 0 or J.width < 0: - return [] + return iter([]) a = floor((I.l + J.l) / 2) b = floor(SQRT2() * (I.l - J.l) / 4) alpha = ZRootTwo(a, b) sol = _solve_ODGP_internal(I - alpha.to_real, J - alpha.conj_sq2.to_real) - sol = [beta + alpha for beta in sol] - sol = [beta for beta in sol if I.within(beta.to_real) and J.within(beta.conj_sq2.to_real)] + sol = map(lambda beta: beta + alpha, sol) + sol = filter(lambda beta: I.within(beta.to_real) and J.within(beta.conj_sq2.to_real), sol) return sol def solve_ODGP_with_parity(I, J, beta): p = beta.parity sol = solve_ODGP((I - p) * SQRT2() / 2, (J - p) * (- SQRT2()) / 2) - sol = [alpha * ZRootTwo(0, 1) + p for alpha in sol] + sol = map(lambda alpha: alpha * ZRootTwo(0, 1) + p, sol) return sol def solve_scaled_ODGP(I, J, k): scale = pow_sqrt2(k) sol = solve_ODGP(I * scale, -J * scale if k & 1 else J * scale) - return [DRootTwo(alpha, k) for alpha in sol] + return map(lambda alpha: DRootTwo(alpha, k), sol) def solve_scaled_ODGP_with_parity(I, J, k, beta): if k == 0: sol = solve_ODGP_with_parity(I, J, beta.renew_denomexp(0)) - return [DRootTwo.from_zroottwo(alpha) for alpha in sol] + return map(lambda alpha: DRootTwo.from_zroottwo(alpha), sol) else: p = beta.renew_denomexp(k).parity offset = DRootTwo.from_int(0) if p == 0 else DRootTwo.power_of_inv_sqrt2(k) sol = solve_scaled_ODGP(I - offset.to_real, J - offset.conj_sq2.to_real, k - 1) - return [alpha + offset for alpha in sol] + return map(lambda alpha: alpha + offset, sol) diff --git a/pygridsynth/quantum_gate.py b/pygridsynth/quantum_gate.py new file mode 100644 index 0000000..e5ec713 --- /dev/null +++ b/pygridsynth/quantum_gate.py @@ -0,0 +1,270 @@ +import mpmath + +CNOT01 = mpmath.matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]) +W_PHASE = mpmath.mp.pi / 4 + + +def Rz(theta): + return mpmath.matrix([[mpmath.exp(- 1.j * theta / 2), 0], + [0, mpmath.exp(1.j * theta / 2)]]) + + +def Rx(theta): + return mpmath.matrix([[mpmath.cos(- theta / 2), 1.j * mpmath.sin(- theta / 2)], + [1.j * mpmath.sin(- theta / 2), mpmath.cos(- theta / 2)]]) + + +class QuantumGate(): + def __init__(self, matrix, wires): + self._matrix = matrix + self._wires = wires + + @property + def matrix(self): + return self._matrix + + @matrix.setter + def matrix(self, matrix): + self._matrix = matrix + + @property + def wires(self): + return self._wires + + def to_whole_matrix(self, num_qubits): + whole_matrix = mpmath.matrix(2 ** num_qubits) + m = len(self._wires) + target_qubit_sorted = sorted(self._wires, reverse=True) + + for i in range(2 ** m): + for j in range(2 ** m): + k = 0 + while k < 2 ** num_qubits: + i2 = 0 + j2 = 0 + for l in range(m): + t = num_qubits - target_qubit_sorted[l] - 1 + if k & 1 << t: + k += 1 << t + i2 |= (i >> l & 1) << t + j2 |= (j >> l & 1) << t + if k >= 2 ** num_qubits: + break + whole_matrix[k | i2, k | j2] = self._matrix[i, j] + k += 1 + + return whole_matrix + + +class SingleQubitGate(QuantumGate): + def __init__(self, matrix, target_qubit): + self._matrix = matrix + self._wires = [target_qubit] + + @property + def target_qubit(self): + return self._wires[0] + + def to_whole_matrix(self, num_qubits): + whole_matrix = mpmath.matrix(2 ** num_qubits) + t = num_qubits - self.target_qubit - 1 + + for i in range(2): + for j in range(2): + k = 0 + while k < 2 ** num_qubits: + if k & 1 << t: + k += 1 << t + if k >= 2 ** num_qubits: + break + whole_matrix[k | (i << t), k | (j << t)] = self._matrix[i, j] + k += 1 + + return whole_matrix + + +class HGate(SingleQubitGate): + def __init__(self, target_qubit): + matrix = mpmath.sqrt(2) / 2 * mpmath.matrix([[1, 1], [1, -1]]) + super().__init__(matrix, target_qubit) + + +class TGate(SingleQubitGate): + def __init__(self, target_qubit): + matrix = mpmath.matrix([[1, 0], [0, mpmath.exp(1.j * mpmath.pi / 4)]]) + super().__init__(matrix, target_qubit) + + +class SGate(SingleQubitGate): + def __init__(self, target_qubit): + matrix = mpmath.matrix([[1, 0], [0, 1.j]]) + super().__init__(matrix, target_qubit) + + +class WGate(SingleQubitGate): + def __init__(self): + matrix = mpmath.exp(1.j * W_PHASE) + super().__init__(matrix, []) + + def to_whole_matrix(self, num_qubits): + return self._matrix + + +class SXGate(SingleQubitGate): + def __init__(self, target_qubit): + matrix = mpmath.matrix([[0, 1], [1, 0]]) + super().__init__(matrix, target_qubit) + + +class RzGate(SingleQubitGate): + def __init__(self, theta, target_qubit): + self._theta = theta + matrix = Rz(theta) + super().__init__(matrix, target_qubit) + + @property + def theta(self): + return self._theta + + def __mul__(self, other): + if isinstance(other, Rz) and self._target_qubit == other.target_qubit: + return RzGate(self._theta + other.theta) + else: + return NotImplemented + + +class RxGate(SingleQubitGate): + def __init__(self, theta, target_qubit): + self._theta = theta + matrix = Rx(theta) + super().__init__(matrix, target_qubit) + + @property + def theta(self): + return self._theta + + def __mul__(self, other): + if isinstance(other, Rx) and self._target_qubit == other.target_qubit: + return RxGate(self._theta + other.theta) + else: + return NotImplemented + + +class CxGate(QuantumGate): + def __init__(self, control_qubit, target_qubit): + matrix = CNOT01 + super().__init__(matrix, [control_qubit, target_qubit]) + + @property + def control_qubit(self): + return self.wires[0] + + @property + def target_qubit(self): + return self.wires[1] + + def to_whole_matrix(self, num_qubits): + whole_matrix = mpmath.matrix(2 ** num_qubits) + c = num_qubits - self.control_qubit - 1 + t = num_qubits - self.target_qubit - 1 + for i0 in range(2): + for j0 in range(2): + for i1 in range(2): + for j1 in range(2): + k = 0 + while k < 2 ** num_qubits: + if k & (1 << min(c, t)): + k += 1 << min(c, t) + if k & (1 << max(c, t)): + k += 1 << max(c, t) + if k >= 2 ** num_qubits: + break + i2 = k | (i0 << c) | (i1 << t) + j2 = k | (j0 << c) | (j1 << t) + whole_matrix[i2, j2] = self._matrix[i0 * 2 + i1, j0 * 2 + j1] + k += 1 + return whole_matrix + + +class QuantumCircuit(list): + def __init__(self, phase=0, args=[]): + self._phase = phase % (2 * mpmath.mp.pi) + super().__init__(args) + + @classmethod + def from_list(cls, gates): + return cls(args=gates) + + @property + def phase(self): + return self._phase + + @phase.setter + def phase(self, phase): + self._phase = phase + + def __str__(self): + return f"exp(1.j * {self._phase}) * \n" + "* \n".join(self) + + def __add__(self, other): + if isinstance(other, QuantumCircuit): + return QuantumCircuit(self.phase + other.phase, super().__add__(other)) + elif isinstance(other, list): + return QuantumCircuit(self.phase, super().__add__(other)) + else: + return NotImplemented + + def __radd__(self, other): + if isinstance(other, QuantumCircuit): + return QuantumCircuit(other.phase + self.phase, super().__radd__(other)) + elif isinstance(other, list): + return QuantumCircuit(self.phase, other.__add__(self)) + else: + return NotImplemented + + def __iadd__(self, other): + if isinstance(other, QuantumCircuit): + self.phase += other.phase + super().__iadd__(other) + return self + elif isinstance(other, list): + super().__iadd__(other) + return self + else: + return NotImplemented + + def decompose_phase_gate(self): + self._phase %= 2 * mpmath.mp.pi + + for _ in range(round(float(self._phase / W_PHASE))): + self.append(WGate()) + + self._phase = 0 + + def to_matrix(self, num_qubits): + U = mpmath.eye(2 ** num_qubits) * mpmath.exp(1.j * self._phase) + for g in self: + U2 = mpmath.zeros(2 ** num_qubits) + if isinstance(g, WGate): + U2 = U * g.matrix + elif isinstance(g, CxGate): + for i in range(2 ** num_qubits): + for j in range(2 ** num_qubits): + for b0 in range(2): + for b1 in range(2): + t = num_qubits - g.target_qubit - 1 + c = num_qubits - g.control_qubit - 1 + j2 = j ^ b0 << c ^ b1 << t + U2[i, j2] += U[i, j] * g.matrix[(j >> c & 1) << 1 | (j >> t & 1), (j2 >> c & 1) << 1 | (j2 >> t & 1)] + elif isinstance(g, SingleQubitGate): + t = num_qubits - g.target_qubit - 1 + for i in range(2 ** num_qubits): + for j in range(2 ** num_qubits): + U2[i, j] += U[i, j] * g.matrix[j >> t & 1, j >> t & 1] + j2 = j ^ 1 << t + U2[i, j2] += U[i, j] * g.matrix[j >> t & 1, j2 >> t & 1] + else: + raise ValueError + U = U2 + + return U diff --git a/pygridsynth/synthesis_of_cliffordT.py b/pygridsynth/synthesis_of_cliffordT.py index d2c1ab8..6f51dbf 100644 --- a/pygridsynth/synthesis_of_cliffordT.py +++ b/pygridsynth/synthesis_of_cliffordT.py @@ -1,13 +1,24 @@ from .ring import OMEGA_POWER from .unitary import DOmegaUnitary from .normal_form import NormalForm +from .quantum_gate import W_PHASE, QuantumCircuit, HGate, TGate, SGate, SXGate, WGate + BIT_SHIFT = [0, 0, 1, 0, 2, 0, 1, 3, 3, 3, 0, 2, 2, 1, 0, 0] BIT_COUNT = [0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4] -def _reduce_denomexp(unitary): - T_POWER_and_H = ["H", "TH", "SH", "TSH"] +def _reduce_denomexp(unitary, wires): + def t_power_and_h(m): + if m == 0: + return [HGate(target_qubit=wires[0])] + elif m == 1: + return [TGate(target_qubit=wires[0]), HGate(target_qubit=wires[0])] + elif m == 2: + return [SGate(target_qubit=wires[0]), HGate(target_qubit=wires[0])] + elif m == 3: + return [TGate(target_qubit=wires[0]), SGate(target_qubit=wires[0]), HGate(target_qubit=wires[0])] + residue_z = unitary.z.residue residue_w = unitary.w.residue residue_squared_z = (unitary.z.u * unitary.z.conj.u).residue @@ -17,41 +28,49 @@ def _reduce_denomexp(unitary): m += 4 if residue_squared_z == 0b0000: unitary = unitary.mul_by_H_and_T_power_from_left(0).renew_denomexp(unitary.k - 1) - return T_POWER_and_H[0], unitary + return t_power_and_h(0), unitary elif residue_squared_z == 0b1010: unitary = unitary.mul_by_H_and_T_power_from_left(-m).renew_denomexp(unitary.k - 1) - return T_POWER_and_H[m], unitary + return t_power_and_h(m), unitary elif residue_squared_z == 0b0001: if BIT_COUNT[residue_z] == BIT_COUNT[residue_w]: unitary = unitary.mul_by_H_and_T_power_from_left(-m).renew_denomexp(unitary.k - 1) - return T_POWER_and_H[m], unitary + return t_power_and_h(m), unitary else: unitary = unitary.mul_by_H_and_T_power_from_left(-m) - return T_POWER_and_H[m], unitary + return t_power_and_h(m), unitary -def decompose_domega_unitary(unitary): - gates = "" +def decompose_domega_unitary(unitary, wires, decompose_phase_gate=True): + circuit = QuantumCircuit() while unitary.k > 0: - g, unitary = _reduce_denomexp(unitary) - gates += g + g, unitary = _reduce_denomexp(unitary, wires=wires) + circuit += g if unitary.n & 1: - gates += "T" + circuit.append(TGate(target_qubit=wires[0])) unitary = unitary.mul_by_T_inv_from_left() if unitary.z == 0: - gates += "X" + circuit.append(SXGate(target_qubit=wires[0])) unitary = unitary.mul_by_X_from_left() + + m_W = 0 for m in range(8): if unitary.z == OMEGA_POWER[m]: m_W = m unitary = unitary.mul_by_W_power_from_left(-m_W) break + m_S = unitary.n >> 1 - gates += "S" * m_S + for _ in range(m_S): + circuit.append(SGate(target_qubit=wires[0])) unitary = unitary.mul_by_S_power_from_left(-m_S) - gates += "W" * m_W + if decompose_phase_gate: + for _ in range(m_W): + circuit.append(WGate()) + else: + circuit.phase = m_W * W_PHASE assert unitary == DOmegaUnitary.identity(), "decomposition failed..." - gates = NormalForm.from_gates(gates).to_gates() - return gates + circuit = NormalForm.from_circuit(circuit).to_circuit(wires=wires) + return circuit diff --git a/pygridsynth/tdgp.py b/pygridsynth/tdgp.py index 55366fd..f65cd7c 100644 --- a/pygridsynth/tdgp.py +++ b/pygridsynth/tdgp.py @@ -1,3 +1,5 @@ +import itertools + from .ring import DRootTwo, DOmega from .region import Interval from .odgp import solve_scaled_ODGP, solve_scaled_ODGP_with_parity @@ -6,40 +8,41 @@ def solve_TDGP(setA, setB, opG, ellipseA_upright, ellipseB_upright, bboxA, bboxB, k, verbose=False, show_graph=False): - sol_sufficient = [] + sol_sufficient = iter([]) sol_x = solve_scaled_ODGP(bboxA.I_x, bboxB.I_x, k + 1) + try: + alpha0 = next(sol_x) + except StopIteration: + return iter([]) + sol_y = solve_scaled_ODGP(bboxA.I_y.fatten(bboxA.I_y.width * 1e-4), bboxB.I_y.fatten(bboxB.I_y.width * 1e-4), k + 1) - if len(sol_x) <= 0 or len(sol_y) <= 0: - sol_sufficient = [] - else: - alpha0 = sol_x[0] - for beta in sol_y: - dx = DRootTwo.power_of_inv_sqrt2(k) - z0 = opG.inv * DOmega.from_droottwo_vector(alpha0, beta, k + 1) - v = opG.inv * DOmega.from_droottwo_vector(dx, DRootTwo.from_int(0), k) - t_A = setA.intersect(z0, v) - t_B = setB.intersect(z0.conj_sq2, v.conj_sq2) - if t_A is None or t_B is None: - continue - - parity = (beta - alpha0).mul_by_sqrt2_power_renewing_denomexp(k) - intA, intB = Interval(*t_A), Interval(*t_B) - dtA = 10 / max(10, (1 << k) * intB.width) - dtB = 10 / max(10, (1 << k) * intA.width) - intA, intB = intA.fatten(dtA), intB.fatten(dtB) - sol_t = solve_scaled_ODGP_with_parity(intA, intB, 1, parity) - sol_x = [alpha * dx + alpha0 for alpha in sol_t] - for alpha in sol_x: - sol_sufficient.append(DOmega.from_droottwo_vector(alpha, beta, k)) - sol_transformed = [opG.inv * z for z in sol_sufficient] - sol = [z for z in sol_transformed if setA.inside(z) and setB.inside(z.conj_sq2)] - - if verbose and len(sol_sufficient) > 0: - print(f"{k=}") - print(f"size of sol_sufficient: {len(sol_sufficient)}, size of sol: {len(sol)}") - if show_graph and len(sol_sufficient) > 0: + + def gen_sol_sufficient(beta): + dx = DRootTwo.power_of_inv_sqrt2(k) + z0 = opG.inv * DOmega.from_droottwo_vector(alpha0, beta, k + 1) + v = opG.inv * DOmega.from_droottwo_vector(dx, DRootTwo.from_int(0), k) + t_A = setA.intersect(z0, v) + t_B = setB.intersect(z0.conj_sq2, v.conj_sq2) + if t_A is None or t_B is None: + return iter([]) + + parity = (beta - alpha0).mul_by_sqrt2_power_renewing_denomexp(k) + intA, intB = Interval(*t_A), Interval(*t_B) + dtA = 10 / max(10, (1 << k) * intB.width) + dtB = 10 / max(10, (1 << k) * intA.width) + intA, intB = intA.fatten(dtA), intB.fatten(dtB) + sol_t = solve_scaled_ODGP_with_parity(intA, intB, 1, parity) + sol_x = map(lambda alpha: alpha * dx + alpha0, sol_t) + return map(lambda alpha: DOmega.from_droottwo_vector(alpha, beta, k), sol_x) + + sol_sufficient = itertools.chain.from_iterable(map(gen_sol_sufficient, sol_y)) + + sol_transformed = map(lambda z: opG.inv * z, sol_sufficient) + sol = filter(lambda z: setA.inside(z) and setB.inside(z.conj_sq2), sol_transformed) + + if show_graph: plot_sol([sol_transformed, sol], setA.ellipse, setB.ellipse, None, None, color_list=['limegreen', 'blue'], size_list=[5, 10]) diff --git a/pygridsynth/to_upright.py b/pygridsynth/to_upright.py index ae02577..b67b541 100644 --- a/pygridsynth/to_upright.py +++ b/pygridsynth/to_upright.py @@ -102,7 +102,7 @@ def _to_upright_ellipse_pair(ellipseA, ellipseB, verbose=False): return opG_l * opG_r -def to_upright_set_pair(setA, setB, show_graph, verbose=False): +def to_upright_set_pair(setA, setB, show_graph=False, verbose=False): opG = _to_upright_ellipse_pair(setA.ellipse, setB.ellipse, verbose=verbose) ellipse_pair = opG * EllipsePair(setA.ellipse, setB.ellipse) ellipseA_upright = ellipse_pair.A diff --git a/pygridsynth/unitary.py b/pygridsynth/unitary.py index c792667..0482f21 100644 --- a/pygridsynth/unitary.py +++ b/pygridsynth/unitary.py @@ -2,6 +2,7 @@ import mpmath from .ring import DOmega +from .quantum_gate import HGate, TGate, SGate, SXGate, WGate class DOmegaUnitary(): @@ -113,14 +114,16 @@ def identity(cls): def from_gates(cls, gates): unitary = cls.identity() for g in reversed(gates): - if g == "H": + if isinstance(g, HGate): unitary = unitary.renew_denomexp(unitary.k + 1).mul_by_H_from_left() - elif g == "T": + elif isinstance(g, TGate): unitary = unitary.mul_by_T_from_left() - elif g == "S": + elif isinstance(g, SGate): unitary = unitary.mul_by_S_from_left() - elif g == "X": + elif isinstance(g, SXGate): unitary = unitary.mul_by_X_from_left() - elif g == "W": + elif isinstance(g, WGate): unitary = unitary.mul_by_W_from_left() + else: + raise ValueError return unitary.reduce_denomexp() From d5008da01d3e58a2f4454b438b925d70d522f424 Mon Sep 17 00:00:00 2001 From: shun0923 Date: Fri, 15 Aug 2025 05:54:43 +0900 Subject: [PATCH 2/3] format with black and isort --- LICENSE | 2 +- pygridsynth/gridsynth.py | 2 +- pygridsynth/normal_form.py | 3 +- pygridsynth/quantum_gate.py | 101 ++++++++++++++++++++++++++---------- pygridsynth/unitary.py | 2 +- requirements.txt | 2 +- 6 files changed, 79 insertions(+), 33 deletions(-) diff --git a/LICENSE b/LICENSE index 40c8fd7..061e14a 100644 --- a/LICENSE +++ b/LICENSE @@ -18,4 +18,4 @@ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. \ No newline at end of file +SOFTWARE. diff --git a/pygridsynth/gridsynth.py b/pygridsynth/gridsynth.py index 7a5192c..8a0f6ff 100644 --- a/pygridsynth/gridsynth.py +++ b/pygridsynth/gridsynth.py @@ -5,10 +5,10 @@ from .diophantine import NO_SOLUTION, diophantine_dyadic from .loop_controller import LoopController from .mymath import solve_quadratic, sqrt +from .quantum_gate import Rz from .region import ConvexSet, Ellipse from .ring import DRootTwo from .synthesis_of_cliffordT import decompose_domega_unitary -from .quantum_gate import Rz from .tdgp import solve_TDGP from .to_upright import to_upright_set_pair from .unitary import DOmegaUnitary diff --git a/pygridsynth/normal_form.py b/pygridsynth/normal_form.py index 3cb37d6..634b5b3 100644 --- a/pygridsynth/normal_form.py +++ b/pygridsynth/normal_form.py @@ -1,5 +1,6 @@ from enum import Enum -from .quantum_gate import QuantumCircuit, HGate, TGate, SGate, SXGate, WGate + +from .quantum_gate import HGate, QuantumCircuit, SGate, SXGate, TGate, WGate class Axis(Enum): diff --git a/pygridsynth/quantum_gate.py b/pygridsynth/quantum_gate.py index e5ec713..fd2576a 100644 --- a/pygridsynth/quantum_gate.py +++ b/pygridsynth/quantum_gate.py @@ -5,16 +5,21 @@ def Rz(theta): - return mpmath.matrix([[mpmath.exp(- 1.j * theta / 2), 0], - [0, mpmath.exp(1.j * theta / 2)]]) + return mpmath.matrix( + [[mpmath.exp(-1.0j * theta / 2), 0], [0, mpmath.exp(1.0j * theta / 2)]] + ) def Rx(theta): - return mpmath.matrix([[mpmath.cos(- theta / 2), 1.j * mpmath.sin(- theta / 2)], - [1.j * mpmath.sin(- theta / 2), mpmath.cos(- theta / 2)]]) + return mpmath.matrix( + [ + [mpmath.cos(-theta / 2), 1.0j * mpmath.sin(-theta / 2)], + [1.0j * mpmath.sin(-theta / 2), mpmath.cos(-theta / 2)], + ] + ) -class QuantumGate(): +class QuantumGate: def __init__(self, matrix, wires): self._matrix = matrix self._wires = wires @@ -31,15 +36,18 @@ def matrix(self, matrix): def wires(self): return self._wires + def __str__(self) -> str: + return f"QuantumGate({self.matrix}, {self.wires})" + def to_whole_matrix(self, num_qubits): - whole_matrix = mpmath.matrix(2 ** num_qubits) + whole_matrix = mpmath.matrix(2**num_qubits) m = len(self._wires) target_qubit_sorted = sorted(self._wires, reverse=True) - for i in range(2 ** m): - for j in range(2 ** m): + for i in range(2**m): + for j in range(2**m): k = 0 - while k < 2 ** num_qubits: + while k < 2**num_qubits: i2 = 0 j2 = 0 for l in range(m): @@ -48,7 +56,7 @@ def to_whole_matrix(self, num_qubits): k += 1 << t i2 |= (i >> l & 1) << t j2 |= (j >> l & 1) << t - if k >= 2 ** num_qubits: + if k >= 2**num_qubits: break whole_matrix[k | i2, k | j2] = self._matrix[i, j] k += 1 @@ -65,17 +73,20 @@ def __init__(self, matrix, target_qubit): def target_qubit(self): return self._wires[0] + def __str__(self) -> str: + return f"QuantumGate({self.matrix}, {self.target_qubit})" + def to_whole_matrix(self, num_qubits): - whole_matrix = mpmath.matrix(2 ** num_qubits) + whole_matrix = mpmath.matrix(2**num_qubits) t = num_qubits - self.target_qubit - 1 for i in range(2): for j in range(2): k = 0 - while k < 2 ** num_qubits: + while k < 2**num_qubits: if k & 1 << t: k += 1 << t - if k >= 2 ** num_qubits: + if k >= 2**num_qubits: break whole_matrix[k | (i << t), k | (j << t)] = self._matrix[i, j] k += 1 @@ -88,24 +99,36 @@ def __init__(self, target_qubit): matrix = mpmath.sqrt(2) / 2 * mpmath.matrix([[1, 1], [1, -1]]) super().__init__(matrix, target_qubit) + def __str__(self) -> str: + return "H" + class TGate(SingleQubitGate): def __init__(self, target_qubit): - matrix = mpmath.matrix([[1, 0], [0, mpmath.exp(1.j * mpmath.pi / 4)]]) + matrix = mpmath.matrix([[1, 0], [0, mpmath.exp(1.0j * mpmath.pi / 4)]]) super().__init__(matrix, target_qubit) + def __str__(self) -> str: + return "T" + class SGate(SingleQubitGate): def __init__(self, target_qubit): - matrix = mpmath.matrix([[1, 0], [0, 1.j]]) + matrix = mpmath.matrix([[1, 0], [0, 1.0j]]) super().__init__(matrix, target_qubit) + def __str__(self) -> str: + return "S" + class WGate(SingleQubitGate): def __init__(self): - matrix = mpmath.exp(1.j * W_PHASE) + matrix = mpmath.exp(1.0j * W_PHASE) super().__init__(matrix, []) + def __str__(self) -> str: + return "W" + def to_whole_matrix(self, num_qubits): return self._matrix @@ -115,6 +138,9 @@ def __init__(self, target_qubit): matrix = mpmath.matrix([[0, 1], [1, 0]]) super().__init__(matrix, target_qubit) + def __str__(self) -> str: + return "X" + class RzGate(SingleQubitGate): def __init__(self, theta, target_qubit): @@ -126,6 +152,9 @@ def __init__(self, theta, target_qubit): def theta(self): return self._theta + def __str__(self) -> str: + return f"RzGate({self.theta}, {self.target_qubit})" + def __mul__(self, other): if isinstance(other, Rz) and self._target_qubit == other.target_qubit: return RzGate(self._theta + other.theta) @@ -143,6 +172,9 @@ def __init__(self, theta, target_qubit): def theta(self): return self._theta + def __str__(self) -> str: + return f"RxGate({self.theta}, {self.target_qubit})" + def __mul__(self, other): if isinstance(other, Rx) and self._target_qubit == other.target_qubit: return RxGate(self._theta + other.theta) @@ -163,8 +195,11 @@ def control_qubit(self): def target_qubit(self): return self.wires[1] + def __str__(self) -> str: + return f"CxGate({self.control_qubit}, {self.target_qubit})" + def to_whole_matrix(self, num_qubits): - whole_matrix = mpmath.matrix(2 ** num_qubits) + whole_matrix = mpmath.matrix(2**num_qubits) c = num_qubits - self.control_qubit - 1 t = num_qubits - self.target_qubit - 1 for i0 in range(2): @@ -172,16 +207,18 @@ def to_whole_matrix(self, num_qubits): for i1 in range(2): for j1 in range(2): k = 0 - while k < 2 ** num_qubits: + while k < 2**num_qubits: if k & (1 << min(c, t)): k += 1 << min(c, t) if k & (1 << max(c, t)): k += 1 << max(c, t) - if k >= 2 ** num_qubits: + if k >= 2**num_qubits: break i2 = k | (i0 << c) | (i1 << t) j2 = k | (j0 << c) | (j1 << t) - whole_matrix[i2, j2] = self._matrix[i0 * 2 + i1, j0 * 2 + j1] + whole_matrix[i2, j2] = self._matrix[ + i0 * 2 + i1, j0 * 2 + j1 + ] k += 1 return whole_matrix @@ -204,7 +241,9 @@ def phase(self, phase): self._phase = phase def __str__(self): - return f"exp(1.j * {self._phase}) * \n" + "* \n".join(self) + return f"exp(1.j * {self._phase}) * \n" + "* \n".join( + str(gate) for gate in self + ) def __add__(self, other): if isinstance(other, QuantumCircuit): @@ -242,24 +281,30 @@ def decompose_phase_gate(self): self._phase = 0 def to_matrix(self, num_qubits): - U = mpmath.eye(2 ** num_qubits) * mpmath.exp(1.j * self._phase) + U = mpmath.eye(2**num_qubits) * mpmath.exp(1.0j * self._phase) for g in self: - U2 = mpmath.zeros(2 ** num_qubits) + U2 = mpmath.zeros(2**num_qubits) if isinstance(g, WGate): U2 = U * g.matrix elif isinstance(g, CxGate): - for i in range(2 ** num_qubits): - for j in range(2 ** num_qubits): + for i in range(2**num_qubits): + for j in range(2**num_qubits): for b0 in range(2): for b1 in range(2): t = num_qubits - g.target_qubit - 1 c = num_qubits - g.control_qubit - 1 j2 = j ^ b0 << c ^ b1 << t - U2[i, j2] += U[i, j] * g.matrix[(j >> c & 1) << 1 | (j >> t & 1), (j2 >> c & 1) << 1 | (j2 >> t & 1)] + U2[i, j2] += ( + U[i, j] + * g.matrix[ + (j >> c & 1) << 1 | (j >> t & 1), + (j2 >> c & 1) << 1 | (j2 >> t & 1), + ] + ) elif isinstance(g, SingleQubitGate): t = num_qubits - g.target_qubit - 1 - for i in range(2 ** num_qubits): - for j in range(2 ** num_qubits): + for i in range(2**num_qubits): + for j in range(2**num_qubits): U2[i, j] += U[i, j] * g.matrix[j >> t & 1, j >> t & 1] j2 = j ^ 1 << t U2[i, j2] += U[i, j] * g.matrix[j >> t & 1, j2 >> t & 1] diff --git a/pygridsynth/unitary.py b/pygridsynth/unitary.py index 064f105..daa0b83 100644 --- a/pygridsynth/unitary.py +++ b/pygridsynth/unitary.py @@ -2,8 +2,8 @@ import mpmath +from .quantum_gate import HGate, SGate, SXGate, TGate, WGate from .ring import DOmega -from .quantum_gate import HGate, TGate, SGate, SXGate, WGate class DOmegaUnitary: diff --git a/requirements.txt b/requirements.txt index 64d3e72..dda7c27 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1 +1 @@ -mpmath \ No newline at end of file +mpmath From 8db2c3ed2208ec8772169465680d966201226a7d Mon Sep 17 00:00:00 2001 From: shun0923 Date: Fri, 15 Aug 2025 07:46:53 +0900 Subject: [PATCH 3/3] Add feature to convert circuit to string --- pygridsynth/gridsynth.py | 30 ++++++++++++++++++++++++++---- pygridsynth/quantum_gate.py | 24 ++++++++++++++++++++---- pygridsynth/unitary.py | 22 ++++++++++++++++++++-- 3 files changed, 66 insertions(+), 10 deletions(-) diff --git a/pygridsynth/gridsynth.py b/pygridsynth/gridsynth.py index 8a0f6ff..cbc0e15 100644 --- a/pygridsynth/gridsynth.py +++ b/pygridsynth/gridsynth.py @@ -92,12 +92,12 @@ def error(theta, gates): def check(theta, gates): - # t_count = gates.count("T") - # h_count = gates.count("H") + t_count = gates.count("T") + h_count = gates.count("H") u_approx = DOmegaUnitary.from_gates(gates) e = error(theta, gates) print(f"{gates=}") - # print(f"{t_count=}, {h_count=}") + print(f"{t_count=}, {h_count=}") print(f"u_approx={u_approx.to_matrix}") print(f"{e=}") @@ -174,7 +174,7 @@ def gridsynth( k += 1 -def gridsynth_gates( +def gridsynth_circuit( theta, epsilon, wires=[0], @@ -203,3 +203,25 @@ def gridsynth_gates( print(f"total time: {(time.time() - start_total) * 1000} ms") return circuit + + +def gridsynth_gates( + theta, + epsilon, + decompose_phase_gate=True, + loop_controller=None, + verbose=False, + measure_time=False, + show_graph=False, +): + circuit = gridsynth_circuit( + theta=theta, + epsilon=epsilon, + wires=[0], + decompose_phase_gate=decompose_phase_gate, + loop_controller=loop_controller, + verbose=verbose, + measure_time=measure_time, + show_graph=show_graph, + ) + return circuit.to_simple_str() diff --git a/pygridsynth/quantum_gate.py b/pygridsynth/quantum_gate.py index fd2576a..c0d0432 100644 --- a/pygridsynth/quantum_gate.py +++ b/pygridsynth/quantum_gate.py @@ -74,7 +74,7 @@ def target_qubit(self): return self._wires[0] def __str__(self) -> str: - return f"QuantumGate({self.matrix}, {self.target_qubit})" + return f"SingleQubitGate({self.matrix}, {self.target_qubit})" def to_whole_matrix(self, num_qubits): whole_matrix = mpmath.matrix(2**num_qubits) @@ -100,6 +100,9 @@ def __init__(self, target_qubit): super().__init__(matrix, target_qubit) def __str__(self) -> str: + return f"HGate({self.target_qubit})" + + def to_simple_str(self) -> str: return "H" @@ -109,6 +112,9 @@ def __init__(self, target_qubit): super().__init__(matrix, target_qubit) def __str__(self) -> str: + return f"TGate({self.target_qubit})" + + def to_simple_str(self) -> str: return "T" @@ -118,6 +124,9 @@ def __init__(self, target_qubit): super().__init__(matrix, target_qubit) def __str__(self) -> str: + return f"SGate({self.target_qubit})" + + def to_simple_str(self) -> str: return "S" @@ -127,6 +136,9 @@ def __init__(self): super().__init__(matrix, []) def __str__(self) -> str: + return "WGate()" + + def to_simple_str(self) -> str: return "W" def to_whole_matrix(self, num_qubits): @@ -139,6 +151,9 @@ def __init__(self, target_qubit): super().__init__(matrix, target_qubit) def __str__(self) -> str: + return f"SXGate({self.target_qubit})" + + def to_simple_str(self) -> str: return "X" @@ -241,9 +256,10 @@ def phase(self, phase): self._phase = phase def __str__(self): - return f"exp(1.j * {self._phase}) * \n" + "* \n".join( - str(gate) for gate in self - ) + return f"exp(1.j * {self._phase}) * " + " * ".join(str(gate) for gate in self) + + def to_simple_str(self): + return "".join(g.to_simple_str() for g in self) def __add__(self, other): if isinstance(other, QuantumCircuit): diff --git a/pygridsynth/unitary.py b/pygridsynth/unitary.py index daa0b83..dad7980 100644 --- a/pygridsynth/unitary.py +++ b/pygridsynth/unitary.py @@ -131,9 +131,9 @@ def identity(cls): return cls(DOmega.from_int(1), DOmega.from_int(0), 0) @classmethod - def from_gates(cls, gates): + def from_circuit(cls, circuit): unitary = cls.identity() - for g in reversed(gates): + for g in reversed(circuit): if isinstance(g, HGate): unitary = unitary.renew_denomexp(unitary.k + 1).mul_by_H_from_left() elif isinstance(g, TGate): @@ -147,3 +147,21 @@ def from_gates(cls, gates): else: raise ValueError return unitary.reduce_denomexp() + + @classmethod + def from_gates(cls, gates): + unitary = cls.identity() + for g in reversed(gates): + if g == "H": + unitary = unitary.renew_denomexp(unitary.k + 1).mul_by_H_from_left() + elif g == "T": + unitary = unitary.mul_by_T_from_left() + elif g == "S": + unitary = unitary.mul_by_S_from_left() + elif g == "X": + unitary = unitary.mul_by_X_from_left() + elif g == "W": + unitary = unitary.mul_by_W_from_left() + else: + raise ValueError + return unitary.reduce_denomexp()