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 5376e4d..cbc0e15 100644 --- a/pygridsynth/gridsynth.py +++ b/pygridsynth/gridsynth.py @@ -5,6 +5,7 @@ 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 @@ -17,14 +18,14 @@ 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 @@ -82,31 +83,22 @@ def generate_complex_unitary(sol): ) -def generate_target_Rz(theta): - return mpmath.matrix( - [[mpmath.exp(-1.0j * theta / 2), 0], [0, mpmath.exp(1.0j * 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 + 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"u_approx={u_approx.to_matrix}") print(f"{e=}") @@ -150,7 +142,6 @@ def gridsynth( if measure_time: time_of_solve_TDGP += time.time() - start start = time.time() - for z in sol: if (z * z.conj).residue == 0: continue @@ -183,9 +174,11 @@ def gridsynth( k += 1 -def gridsynth_gates( +def gridsynth_circuit( theta, epsilon, + wires=[0], + decompose_phase_gate=True, loop_controller=None, verbose=False, measure_time=False, @@ -202,8 +195,33 @@ def gridsynth_gates( ) start = time.time() if measure_time else 0.0 - 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 + + +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/normal_form.py b/pygridsynth/normal_form.py index 403f251..634b5b3 100644 --- a/pygridsynth/normal_form.py +++ b/pygridsynth/normal_form.py @@ -1,5 +1,7 @@ from enum import Enum +from .quantum_gate import HGate, QuantumCircuit, SGate, SXGate, TGate, WGate + class Axis(Enum): I = 0 @@ -127,6 +129,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 ( @@ -182,15 +197,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): @@ -204,13 +232,21 @@ 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: @@ -230,21 +266,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 063bd50..307d5e0 100644 --- a/pygridsynth/odgp.py +++ b/pygridsynth/odgp.py @@ -1,25 +1,27 @@ +import itertools + from .mymath import SQRT2, ceil, floor, floorlog, pow_sqrt2 from .ring import LAMBDA, DRootTwo, ZRootTwo 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 @@ -27,46 +29,43 @@ def _solve_ODGP_internal(I, J): 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..c0d0432 --- /dev/null +++ b/pygridsynth/quantum_gate.py @@ -0,0 +1,331 @@ +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.0j * theta / 2), 0], [0, mpmath.exp(1.0j * theta / 2)]] + ) + + +def Rx(theta): + 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: + 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 __str__(self) -> str: + return f"QuantumGate({self.matrix}, {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 __str__(self) -> str: + return f"SingleQubitGate({self.matrix}, {self.target_qubit})" + + 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) + + def __str__(self) -> str: + return f"HGate({self.target_qubit})" + + def to_simple_str(self) -> str: + return "H" + + +class TGate(SingleQubitGate): + def __init__(self, target_qubit): + matrix = mpmath.matrix([[1, 0], [0, mpmath.exp(1.0j * mpmath.pi / 4)]]) + super().__init__(matrix, target_qubit) + + def __str__(self) -> str: + return f"TGate({self.target_qubit})" + + def to_simple_str(self) -> str: + return "T" + + +class SGate(SingleQubitGate): + def __init__(self, target_qubit): + matrix = mpmath.matrix([[1, 0], [0, 1.0j]]) + super().__init__(matrix, target_qubit) + + def __str__(self) -> str: + return f"SGate({self.target_qubit})" + + def to_simple_str(self) -> str: + return "S" + + +class WGate(SingleQubitGate): + def __init__(self): + matrix = mpmath.exp(1.0j * W_PHASE) + super().__init__(matrix, []) + + def __str__(self) -> str: + return "WGate()" + + def to_simple_str(self) -> str: + return "W" + + 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) + + def __str__(self) -> str: + return f"SXGate({self.target_qubit})" + + def to_simple_str(self) -> str: + return "X" + + +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 __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) + 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 __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) + 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 __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) + 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}) * " + " * ".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): + 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.0j * 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 cb39d6e..45bc26b 100644 --- a/pygridsynth/synthesis_of_cliffordT.py +++ b/pygridsynth/synthesis_of_cliffordT.py @@ -1,13 +1,28 @@ from .normal_form import NormalForm from .ring import OMEGA_POWER from .unitary import DOmegaUnitary +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 @@ -19,45 +34,53 @@ def _reduce_denomexp(unitary): 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 75df50f..a8185f5 100644 --- a/pygridsynth/tdgp.py +++ b/pygridsynth/tdgp.py @@ -1,3 +1,5 @@ +import itertools + from .myplot import plot_sol from .odgp import solve_scaled_ODGP, solve_scaled_ODGP_with_parity from .region import Interval @@ -16,42 +18,43 @@ def solve_TDGP( 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)] + 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 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: + if show_graph: plot_sol( [sol_transformed, sol], setA.ellipse, diff --git a/pygridsynth/to_upright.py b/pygridsynth/to_upright.py index 6663f31..500644c 100644 --- a/pygridsynth/to_upright.py +++ b/pygridsynth/to_upright.py @@ -108,7 +108,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 736bab8..dad7980 100644 --- a/pygridsynth/unitary.py +++ b/pygridsynth/unitary.py @@ -2,6 +2,7 @@ import mpmath +from .quantum_gate import HGate, SGate, SXGate, TGate, WGate from .ring import DOmega @@ -129,6 +130,24 @@ def reduce_denomexp(self): def identity(cls): return cls(DOmega.from_int(1), DOmega.from_int(0), 0) + @classmethod + def from_circuit(cls, circuit): + unitary = cls.identity() + 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): + unitary = unitary.mul_by_T_from_left() + elif isinstance(g, SGate): + unitary = unitary.mul_by_S_from_left() + elif isinstance(g, SXGate): + unitary = unitary.mul_by_X_from_left() + elif isinstance(g, WGate): + unitary = unitary.mul_by_W_from_left() + else: + raise ValueError + return unitary.reduce_denomexp() + @classmethod def from_gates(cls, gates): unitary = cls.identity() @@ -143,4 +162,6 @@ def from_gates(cls, gates): unitary = unitary.mul_by_X_from_left() elif g == "W": unitary = unitary.mul_by_W_from_left() + else: + raise ValueError return unitary.reduce_denomexp() 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