From 8cde44bd26ad3cbf4a3c0d79854e7704856bb2e1 Mon Sep 17 00:00:00 2001 From: Ondra Lengal Date: Tue, 13 Oct 2026 10:37:25 +0200 Subject: [PATCH 1/5] fix: bound the FiveTuple rescaling by the depth, not by the input's max_k add_rescaling covered k = floor(|k_out - k_target| / 2) only below Synthesizer.max_k, which is computed once from the input. Past the input depth (incremental search, -d, restricted gate sets, custom vector pairs) the rescaled vectors were left unconstrained, so any circuit passed the equivalence check: t;x;t;x with -d 1 gave h;h;h;h. encode_equivalence now covers every gap the encoded depth can produce (each layer raises k by at most 2), and add_rescaling asserts k < d, so a gap outside the case split is UNSAT instead of unconstrained. Part of #65 (the MILP counterpart, F15, is not addressed here). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_0177AhE9SXm3iJ1KP8GTSk6o --- src/generator.py | 5 ++- src/synth.py | 12 +++++- tests/test_rescaling_bound.py | 74 +++++++++++++++++++++++++++++++++++ 3 files changed, 88 insertions(+), 3 deletions(-) create mode 100644 tests/test_rescaling_bound.py diff --git a/src/generator.py b/src/generator.py index d9c9acd..aa7a020 100644 --- a/src/generator.py +++ b/src/generator.py @@ -751,8 +751,9 @@ def multiply_by_identity_unscaled(result_vec, source_vec): return self.And(*eqs) if self.logic == "QF_LIA" or self.logic == "QF_NIA": - if (d // 2) == 0: - d = 2 + # the case split below covers k in [0, d); a larger k would leave r1 and r2 unconstrained, and any + # circuit would pass the equivalence check (#65) + self.add_assertion(self.LT(k, self.Int(d))) for i in range(d): pow2 = self.Int(2 ** i) r1_multiply_by_m = self.And(multiply_by_m_scaled(r1, v1, pow2), multiply_by_identity_unscaled(r2, v2)) diff --git a/src/synth.py b/src/synth.py index 081ca93..c49bdce 100644 --- a/src/synth.py +++ b/src/synth.py @@ -27,6 +27,9 @@ from logger import Logger import numpy as np +# the largest increase of k (FiveTuple) caused by one layer: sx, sxdg, csx and sqrtswap add 2, h and ch add 1 +MAX_K_INCREMENT = 2 + def rus_recovery_pairs(output_vectors, vector_pairs, measured_qubits): # (input, output) pairs for the RUS recovery operation: it maps the state after a failed # attempt (measured qubits projected to 1) back to the original input state @@ -54,6 +57,7 @@ def __init__(self) -> None: self.circuit_q = None # qubits the synthesized circuit acts on, the first circuit_q of the q self.d = None self.max_k = 0 + self.pair_k = [] self.encoding_method = None self.curr_depth = 0 self.v = 0 @@ -719,7 +723,11 @@ def encode_equivalence(self, vectors1, vectors2, depth, already_measured : bool if self.gen.mode == "milp" or self.gen.mode == "gurobi": self.gen.add_milp_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, self.max_k) else: - self.gen.add_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, self.max_k) + # after `depth` layers, k lies in [k_in, k_in + MAX_K_INCREMENT*depth]; cover every possible + # floor(|k - k_target| / 2) (self.max_k is fixed from the input and too small past its depth) + k_in, k_target = self.pair_k[pair_idx] + max_gap = max(k_in + MAX_K_INCREMENT * depth - k_target, k_target - k_in) + self.gen.add_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, max_gap // 2 + 1) for i in range(2**self.q): self.gen.add_assertion(self.gen.Equals(rescaled1[i], rescaled2[i])) @@ -1005,6 +1013,8 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool: # encode vector bounds inter_states = [] target_states = [] + # k of every (input, target) pair, used to bound the FiveTuple rescaling + self.pair_k = [(input_state.k, output_state.k) for input_state, output_state in vector_pairs] bounds = [1] * (self.d + 1) if self.complex_representation != Complex: for d in range(self.d + 1): diff --git a/tests/test_rescaling_bound.py b/tests/test_rescaling_bound.py new file mode 100644 index 0000000..34a09df --- /dev/null +++ b/tests/test_rescaling_bound.py @@ -0,0 +1,74 @@ +"""FiveTuple equivalence when the k gap between the output and the target exceeds the input's max_k (#65).""" + +import numpy as np +import pytest +from qiskit import QuantumCircuit +from qiskit.quantum_info import Operator + +from complex.fivetuple import FiveTuple +from complex.vector import Vector +from gates import GateSet +from synth import Synthesizer + +SMT = dict(solving="smt", solver="yices2", complex_representation="FiveTuple") + + +def write_qasm(path, gates): + with open(path, "w") as f: + f.write('OPENQASM 3.0;\ninclude "stdgates.inc";\nqubit[1] q;\n') + f.writelines(f"{gate} q[0];\n" for gate in gates) + return str(path) + + +def operator(gates): + # reference unitary computed by Qiskit, independently of MILQ's simulator + qc = QuantumCircuit(1) + for gate in gates: + getattr(qc, gate)(0) + return Operator(qc).data + + +def synthesized_gates(circuit): + return [gate.name for gate in circuit.gates if gate is not None and gate.name != "id"] + + +def test_shorter_depth_bound_still_needs_an_exact_circuit(workdir): + # t x t x = w*I; with d=1 the search goes past the input depth, where the old bound left the equivalence + # unconstrained and h h h h (= I) was accepted + qasm_file = write_qasm(workdir / "in.qasm", ["t", "x", "t", "x"]) + res, circuit, _ = Synthesizer().synthesis(qasm_file=qasm_file, vectors="all", mode="incremental", d=1, + output_qasm="out.qasm", **SMT) + assert res + np.testing.assert_allclose(operator(synthesized_gates(circuit)), operator(["t", "x", "t", "x"]), atol=1e-9) + + +def test_unreachable_target_is_not_synthesized(workdir): + # X is not a product of H gates, at any depth + qasm_file = write_qasm(workdir / "in.qasm", ["x"]) + res, circuit, _ = Synthesizer().synthesis(qasm_file=qasm_file, vectors="all", mode="incremental", + gate_set=GateSet(["h"]), max_depth=6, output_qasm="out.qasm", **SMT) + assert (res, circuit) == (False, None) + + +def test_restricted_gate_set_gives_an_exact_circuit(workdir): + # Z = h sx sx h; the old bound accepted sx sx (= X) + qasm_file = write_qasm(workdir / "in.qasm", ["z"]) + res, circuit, _ = Synthesizer().synthesis(qasm_file=qasm_file, vectors="all", mode="incremental", + gate_set=GateSet(["h", "sx"]), output_qasm="out.qasm", **SMT) + assert res + gates = synthesized_gates(circuit) + assert len(gates) == 4 + np.testing.assert_allclose(operator(gates), operator(["z"]), atol=1e-9) + + +@pytest.mark.parametrize("mode", ["basic", "incremental"]) +def test_custom_target_with_a_large_k(workdir, mode): + # |0> -> |1>, with the target written as 4/sqrt(2)^4: the old bound only covered the gap of the input depth + inp = Vector(q=2, element_representation=FiveTuple, k=0) + out = Vector(q=2, element_representation=FiveTuple, k=4) + inp[0] = FiveTuple(a=1, b=0, c=0, d=0) + out[1] = FiveTuple(a=4, b=0, c=0, d=0) + res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=[(inp, out)], q=1, d=1, + gate_set=GateSet(["x", "h"]), mode=mode, output_qasm="out.qasm", **SMT) + assert res + assert synthesized_gates(circuit) == ["x"] From 6e5e6cee7ffc622bc54f7c8760f8d6c116e4886e Mon Sep 17 00:00:00 2001 From: Ondra Lengal Date: Tue, 13 Oct 2026 11:59:05 +0200 Subject: [PATCH 2/5] fix: reject a non-integral FiveTuple k in vector pairs The rescaling bound is computed from the pairs' k, so a custom pair with k=4.0 raised a TypeError in range(). Integral floats are now accepted and anything else is rejected with a ValueError. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_0177AhE9SXm3iJ1KP8GTSk6o --- src/synth.py | 9 ++++++++- tests/test_rescaling_bound.py | 22 +++++++++++++++++----- 2 files changed, 25 insertions(+), 6 deletions(-) diff --git a/src/synth.py b/src/synth.py index c49bdce..c343ea0 100644 --- a/src/synth.py +++ b/src/synth.py @@ -30,6 +30,12 @@ # the largest increase of k (FiveTuple) caused by one layer: sx, sxdg, csx and sqrtswap add 2, h and ch add 1 MAX_K_INCREMENT = 2 +def integral_k(k): + # a FiveTuple vector is divided by sqrt(2)^k, k has to be an integer (4.0 is accepted) + if k is None or int(k) != k: + raise ValueError(f"k of a FiveTuple vector must be an integer, got {k!r}") + return int(k) + def rus_recovery_pairs(output_vectors, vector_pairs, measured_qubits): # (input, output) pairs for the RUS recovery operation: it maps the state after a failed # attempt (measured qubits projected to 1) back to the original input state @@ -1014,7 +1020,8 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool: inter_states = [] target_states = [] # k of every (input, target) pair, used to bound the FiveTuple rescaling - self.pair_k = [(input_state.k, output_state.k) for input_state, output_state in vector_pairs] + if self.complex_representation != Complex: + self.pair_k = [(integral_k(input_state.k), integral_k(output_state.k)) for input_state, output_state in vector_pairs] bounds = [1] * (self.d + 1) if self.complex_representation != Complex: for d in range(self.d + 1): diff --git a/tests/test_rescaling_bound.py b/tests/test_rescaling_bound.py index 34a09df..ee490ba 100644 --- a/tests/test_rescaling_bound.py +++ b/tests/test_rescaling_bound.py @@ -61,14 +61,26 @@ def test_restricted_gate_set_gives_an_exact_circuit(workdir): np.testing.assert_allclose(operator(gates), operator(["z"]), atol=1e-9) -@pytest.mark.parametrize("mode", ["basic", "incremental"]) -def test_custom_target_with_a_large_k(workdir, mode): - # |0> -> |1>, with the target written as 4/sqrt(2)^4: the old bound only covered the gap of the input depth +def zero_to_one(k_target): + # |0> -> |1>, with the target written as 4/sqrt(2)^k_target inp = Vector(q=2, element_representation=FiveTuple, k=0) - out = Vector(q=2, element_representation=FiveTuple, k=4) + out = Vector(q=2, element_representation=FiveTuple, k=k_target) inp[0] = FiveTuple(a=1, b=0, c=0, d=0) out[1] = FiveTuple(a=4, b=0, c=0, d=0) - res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=[(inp, out)], q=1, d=1, + return [(inp, out)] + + +@pytest.mark.parametrize("k_target", [4, 4.0]) +@pytest.mark.parametrize("mode", ["basic", "incremental"]) +def test_custom_target_with_a_large_k(workdir, mode, k_target): + # the old bound only covered the gap of the input depth, and the empty circuit was accepted + res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=zero_to_one(k_target), q=1, d=1, gate_set=GateSet(["x", "h"]), mode=mode, output_qasm="out.qasm", **SMT) assert res assert synthesized_gates(circuit) == ["x"] + + +def test_custom_target_with_a_fractional_k(workdir): + with pytest.raises(ValueError, match="k of a FiveTuple vector must be an integer, got 4.5"): + Synthesizer().synthesis(vectors="custom", vector_pairs=zero_to_one(4.5), q=1, d=1, gate_set=GateSet(["x", "h"]), + mode="incremental", output_qasm="out.qasm", **SMT) From eede8f612ad5e6f9b16c8af2fd5cc2a31d1f49d1 Mon Sep 17 00:00:00 2001 From: Jakub Havlik Date: Sat, 17 Oct 2026 19:32:19 +0200 Subject: [PATCH 3/5] fix: integral_k check before synthesis even begins -- crashed pysmt before --- src/synth.py | 5 ++++- tests/test_rescaling_bound.py | 9 +++++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/src/synth.py b/src/synth.py index c343ea0..9651208 100644 --- a/src/synth.py +++ b/src/synth.py @@ -915,6 +915,9 @@ def synthesis(self, qasm_file=None, matrix=None, vectors=None, vector_pairs=None self.gate_set.add_gate("id", weight=0, qubits=1) self.circuit_q = self.q // 2 if vectors == "jamiolkowski" else self.q self.v = len(vector_pairs) + if self.complex_representation != Complex: + for state in (state for pair in vector_pairs for state in pair): + state.k = integral_k(state.k) self.curr_depth = self.d max = 0 for pair_idx, (input_vector, output_vector) in enumerate(vector_pairs): @@ -1021,7 +1024,7 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool: target_states = [] # k of every (input, target) pair, used to bound the FiveTuple rescaling if self.complex_representation != Complex: - self.pair_k = [(integral_k(input_state.k), integral_k(output_state.k)) for input_state, output_state in vector_pairs] + self.pair_k = [(input_state.k, output_state.k) for input_state, output_state in vector_pairs] bounds = [1] * (self.d + 1) if self.complex_representation != Complex: for d in range(self.d + 1): diff --git a/tests/test_rescaling_bound.py b/tests/test_rescaling_bound.py index ee490ba..046d198 100644 --- a/tests/test_rescaling_bound.py +++ b/tests/test_rescaling_bound.py @@ -80,6 +80,15 @@ def test_custom_target_with_a_large_k(workdir, mode, k_target): assert synthesized_gates(circuit) == ["x"] +def test_float_k_reaches_pysmt_as_an_integer(workdir): + # yices2 tolerates a k of 4.0 in the formula; pySMT crashed on it, so the vectors must carry the int + res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=zero_to_one(4.0), q=1, d=1, + gate_set=GateSet(["x", "h"]), mode="basic", output_qasm="out.qasm", + solving="pysmt", solver="z3", complex_representation="FiveTuple") + assert res + assert synthesized_gates(circuit) == ["x"] + + def test_custom_target_with_a_fractional_k(workdir): with pytest.raises(ValueError, match="k of a FiveTuple vector must be an integer, got 4.5"): Synthesizer().synthesis(vectors="custom", vector_pairs=zero_to_one(4.5), q=1, d=1, gate_set=GateSet(["x", "h"]), From 79bfe6b47a72c187f96ea739267743a15a9b523e Mon Sep 17 00:00:00 2001 From: Jakub Havlik Date: Sun, 18 Oct 2026 00:44:42 +0200 Subject: [PATCH 4/5] fix: bound the MILP FiveTuple rescaling by the encoded depth (#65) add_milp_rescaling took its k bound and both big-Ms from Synthesizer.max_k, which is fixed from the input. Past the input depth, or with custom pairs whose k is not reduced, the correct circuit was infeasible: CBC found no sx;sx for X from {sx, z}, nor h;sx;sx;h for Z from {h, sx}. - The MILP rescaling uses the same per-depth gap bound as the SMT one, with a parity big-M of 2*max_gap + 1 and a value big-M from the vectors' bounds. - k of the rescaling has no ub: incremental search gets back the first depth's PuLP variable with its stale bound, the one-hot bounds k anyway. - The intermediate states' k bound includes the input's k. - The per-layer increase of k comes from a k_increments table in gates.py, used by encode_layer and by the bounds, instead of a constant 2. --- src/gates.py | 6 +++++ src/generator.py | 13 +++++----- src/synth.py | 44 +++++++++++++++++---------------- tests/test_rescaling_bound.py | 46 ++++++++++++++++++++++++++++++++++- 4 files changed, 81 insertions(+), 28 deletions(-) diff --git a/src/gates.py b/src/gates.py index 2e82e88..b2f203d 100644 --- a/src/gates.py +++ b/src/gates.py @@ -15,6 +15,12 @@ self_adjoints = ['h', 'x', 'y', 'z', 'cx', 'cz', 'cy', 'ch', 'swap', 'ccx', 'ccz'] +# FiveTuple: how much a gate raises k (the sqrt(2)^k denominator), the other gates keep it +k_increments = {'h': 1, 'ch': 1, 'sx': 2, 'sxdg': 2, 'csx': 2, 'sqrtswap': 2} + +def max_k_increment(gate_set) -> int: + return max((k_increments.get(gate, 0) for gate in gate_set), default=0) + # OpenQASM 3 definitions of supported gates that are not part of stdgates.inc nonstandard_gate_definitions = { 'sxdg': "gate sxdg a { inv @ sx a; }", diff --git a/src/generator.py b/src/generator.py index aa7a020..86a5e49 100644 --- a/src/generator.py +++ b/src/generator.py @@ -768,14 +768,15 @@ def multiply_by_identity_unscaled(result_vec, source_vec): else: raise ValueError("rescaling with wrong logic") - def add_milp_rescaling(self, r1, r2, v1, v2, pair_idx, d): - # first encode k as the result of the operation floor(abs(k1 - k2)/2) + def add_milp_rescaling(self, r1, r2, v1, v2, pair_idx, max_gap, value_bigM): + # first encode k as the result of the operation floor(abs(k1 - k2)/2), where abs(k1 - k2) <= max_gap + # value_bigM bounds |v - r| for every entry of v1, v2 and r1, r2 complex_representation = v1.element_representation if self.mode == "milp": - k = self.declare_integer(f"k{pair_idx}", lb=0, ub=d) + d = max_gap // 2 + k = self.declare_integer(f"k{pair_idx}", lb=0) q = self.declare_bool(f"q{pair_idx}") - - bigM = 2*d + 1 + bigM = 2*max_gap + 1 sleq = self.declare_bool(f"sleq{pair_idx}") self.add_assertion((v1.k - v2.k) - (2*k + q) <= bigM * sleq) self.add_assertion((2*k + q) - (v1.k - v2.k) <= bigM * sleq) @@ -789,7 +790,7 @@ def add_milp_rescaling(self, r1, r2, v1, v2, pair_idx, d): self.add_assertion(k == self.Sum([constants[i] * i for i in range(d+1)])) # now sleq, si, even encode all case splits needed for the rescaling # encode all combinations to assign to r1, r2 - bigM = (2**(d+1)) + 1 + bigM = value_bigM # for every si, I,M is multiplied by 2^i for i in range(len(v1)): for j, s in enumerate(constants): diff --git a/src/synth.py b/src/synth.py index 9651208..cb819ed 100644 --- a/src/synth.py +++ b/src/synth.py @@ -19,7 +19,7 @@ from pysmt.shortcuts import Portfolio, Symbol, Real, And, Equals, Plus, GT, LT, get_env, Int, Or, Not, Implies, GE, LE, Ite, Minus, Div, Times, write_smtlib, Solver from pysmt.typing import REAL, INT from pareto import Pareto -from gates import GateSet, supported_gates, check_supported, Circuit +from gates import GateSet, check_supported, Circuit, k_increments, max_k_increment from pulp import * from gurobipy import GRB from solvers import PortfolioSMTSolver, SMTSolver @@ -27,9 +27,6 @@ from logger import Logger import numpy as np -# the largest increase of k (FiveTuple) caused by one layer: sx, sxdg, csx and sqrtswap add 2, h and ch add 1 -MAX_K_INCREMENT = 2 - def integral_k(k): # a FiveTuple vector is divided by sqrt(2)^k, k has to be an integer (4.0 is accepted) if k is None or int(k) != k: @@ -64,6 +61,7 @@ def __init__(self) -> None: self.d = None self.max_k = 0 self.pair_k = [] + self.k_increment = 0 self.encoding_method = None self.curr_depth = 0 self.v = 0 @@ -195,7 +193,7 @@ def add_k_incr(sel, out_k, inp_k, increment): propagate_identities[pos].append(bool_var) propagate_identities[other].append(bool_var) if self.complex_representation == FiveTuple: - add_k_incr(bool_var, out.k, inp.k, 1) + add_k_incr(bool_var, out.k, inp.k, k_increments[gate]) elif gate == 's': for pos in range(size): one_flag = (pos >> q) & 1 @@ -261,7 +259,7 @@ def add_k_incr(sel, out_k, inp_k, increment): propagate_identities[pos].append(bool_var) propagate_identities[other].append(bool_var) if self.complex_representation == FiveTuple: - add_k_incr(bool_var, out.k, inp.k, 2) + add_k_incr(bool_var, out.k, inp.k, k_increments[gate]) elif gate == 'sxdg': modified_positions = [] for pos in range(size): @@ -274,7 +272,7 @@ def add_k_incr(sel, out_k, inp_k, increment): propagate_identities[pos].append(bool_var) propagate_identities[other].append(bool_var) if self.complex_representation == FiveTuple: - add_k_incr(bool_var, out.k, inp.k, 2) + add_k_incr(bool_var, out.k, inp.k, k_increments[gate]) elif gate == 'y': modified_positions = [] for pos in range(size): @@ -373,7 +371,7 @@ def add_k_incr(sel, out_k, inp_k, increment): elif cannot_propagate: saved_operations[pos].append((bool_var, inp[pos])) if self.complex_representation == FiveTuple: - add_k_incr(bool_var, out.k, inp.k, 1) + add_k_incr(bool_var, out.k, inp.k, k_increments[gate]) elif gate == 'csx': modified_positions = [] for pos in range(size): @@ -394,7 +392,7 @@ def add_k_incr(sel, out_k, inp_k, increment): elif cannot_propagate: saved_operations[pos].append((bool_var, inp[pos])) if self.complex_representation == FiveTuple: - add_k_incr(bool_var, out.k, inp.k, 2) + add_k_incr(bool_var, out.k, inp.k, k_increments[gate]) elif gate == 'cy': modified_positions = [] for pos in range(2**self.q): @@ -510,7 +508,7 @@ def add_k_incr(sel, out_k, inp_k, increment): elif cannot_propagate: saved_operations[pos].append((bool_var, inp[pos])) if self.complex_representation == FiveTuple: - add_k_incr(bool_var, out.k, inp.k, 2) + add_k_incr(bool_var, out.k, inp.k, k_increments[gate]) elif len(rest) == 4: # three qubit q1 = rest[1] q2 = rest[2] @@ -706,7 +704,7 @@ def encode_equivalence(self, vectors1, vectors2, depth, already_measured : bool for i in range(8): selectors.append(self.gen.declare_bool(f"global_phase_sel_{i}")) self.gen.ExactlyOne(*selectors) - vec = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Multiplied_by_global_phase_{vec2.name}", bound=2**self.max_k) + vec = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Multiplied_by_global_phase_{vec2.name}", bound=vec2.bound) vec_copy = vec2.copy() for i in range(8): # add selection shift @@ -724,15 +722,18 @@ def encode_equivalence(self, vectors1, vectors2, depth, already_measured : bool vec_target = vec # fivetuple rescaling - rescaled1 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled1_{vec1.name}", bound=2**self.max_k) - rescaled2 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled2_{vec2.name}", bound=2**self.max_k) + # after `depth` layers, k lies in [k_in, k_in + self.k_increment*depth], which bounds |k - k_target| + # (self.max_k is fixed from the input and too small past its depth) + k_in, k_target = self.pair_k[pair_idx] + max_gap = max(k_in + self.k_increment * depth - k_target, k_target - k_in) + # the rescaling multiplies by at most 2^(max_gap // 2), and an odd gap adds two components + value_bound = max(vec1.bound, vec_target.bound) + rescaled_bound = 2 * value_bound * 2**(max_gap // 2) + rescaled1 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled1_{vec1.name}", bound=rescaled_bound) + rescaled2 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled2_{vec2.name}", bound=rescaled_bound) if self.gen.mode == "milp" or self.gen.mode == "gurobi": - self.gen.add_milp_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, self.max_k) + self.gen.add_milp_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, max_gap, value_bound + rescaled_bound) else: - # after `depth` layers, k lies in [k_in, k_in + MAX_K_INCREMENT*depth]; cover every possible - # floor(|k - k_target| / 2) (self.max_k is fixed from the input and too small past its depth) - k_in, k_target = self.pair_k[pair_idx] - max_gap = max(k_in + MAX_K_INCREMENT * depth - k_target, k_target - k_in) self.gen.add_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, max_gap // 2 + 1) for i in range(2**self.q): self.gen.add_assertion(self.gen.Equals(rescaled1[i], rescaled2[i])) @@ -1022,9 +1023,10 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool: # encode vector bounds inter_states = [] target_states = [] - # k of every (input, target) pair, used to bound the FiveTuple rescaling + # k of every (input, target) pair and its largest increase per layer, used to bound the FiveTuple k if self.complex_representation != Complex: self.pair_k = [(input_state.k, output_state.k) for input_state, output_state in vector_pairs] + self.k_increment = max_k_increment(self.gate_set) bounds = [1] * (self.d + 1) if self.complex_representation != Complex: for d in range(self.d + 1): @@ -1049,7 +1051,7 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool: for d in range(self.d + 1): Inter = None if self.complex_representation == FiveTuple: - Inter = Vector(q=2**self.q, name=f"I_{pair_idx}_{d}", generator=self.gen, element_representation=self.complex_representation, k=input_state.k if d == 0 else 0, n = input_state.n, bound = bounds[d], k_bound = 2 * (d+1)) + Inter = Vector(q=2**self.q, name=f"I_{pair_idx}_{d}", generator=self.gen, element_representation=self.complex_representation, k=input_state.k if d == 0 else 0, n = input_state.n, bound = bounds[d], k_bound = input_state.k + self.k_increment * d) else: Inter = Vector(q=2**self.q, name=f"I_{pair_idx}_{d}", generator=self.gen, element_representation=self.complex_representation, bound = bounds[d], k_bound = 2 * (d+1)) @@ -1152,7 +1154,7 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool: for pair_idx in range(len(vector_pairs)): if self.curr_depth > self.d: if self.complex_representation == FiveTuple: - state = Vector(q=2**self.q, name=f"I_{pair_idx}_{self.curr_depth}", generator=self.gen, element_representation=self.complex_representation, k=0, n = vector_pairs[pair_idx][0].n, bound = self.bound, k_bound = 2 * (self.curr_depth+1)) + state = Vector(q=2**self.q, name=f"I_{pair_idx}_{self.curr_depth}", generator=self.gen, element_representation=self.complex_representation, k=0, n = vector_pairs[pair_idx][0].n, bound = self.bound, k_bound = vector_pairs[pair_idx][0].k + self.k_increment * self.curr_depth) else: state = Vector(q=2**self.q, name=f"I_{pair_idx}_{self.curr_depth}", generator=self.gen, element_representation=self.complex_representation, bound = self.bound, k_bound = 2 * (self.curr_depth+1)) inter_states[pair_idx].append(state) diff --git a/tests/test_rescaling_bound.py b/tests/test_rescaling_bound.py index 046d198..8e8d767 100644 --- a/tests/test_rescaling_bound.py +++ b/tests/test_rescaling_bound.py @@ -7,7 +7,7 @@ from complex.fivetuple import FiveTuple from complex.vector import Vector -from gates import GateSet +from gates import GateSet, max_k_increment from synth import Synthesizer SMT = dict(solving="smt", solver="yices2", complex_representation="FiveTuple") @@ -93,3 +93,47 @@ def test_custom_target_with_a_fractional_k(workdir): with pytest.raises(ValueError, match="k of a FiveTuple vector must be an integer, got 4.5"): Synthesizer().synthesis(vectors="custom", vector_pairs=zero_to_one(4.5), q=1, d=1, gate_set=GateSet(["x", "h"]), mode="incremental", output_qasm="out.qasm", **SMT) + + +MILP = dict(solving="milp", solver="cbc", complex_representation="FiveTuple") + + +@pytest.mark.parametrize("target, gate_set, expected", [ + # past the input's max_k, the stale k bound and big-M made the correct circuit infeasible (#65 F15) + (["x"], ["sx", "z"], ["sx", "sx"]), + (["z"], ["h", "sx"], ["h", "sx", "sx", "h"]), +]) +def test_milp_finds_circuits_past_the_input_depth(workdir, target, gate_set, expected): + qasm_file = write_qasm(workdir / "in.qasm", target) + res, circuit, _ = Synthesizer().synthesis(qasm_file=qasm_file, vectors="all", mode="incremental", + gate_set=GateSet(gate_set), max_depth=6, output_qasm="out.qasm", **MILP) + assert res + gates = synthesized_gates(circuit) + assert len(gates) == len(expected) + np.testing.assert_allclose(operator(gates), operator(target), atol=1e-9) + + +def test_milp_custom_target_with_a_large_k(workdir): + res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=zero_to_one(4), q=1, d=1, + gate_set=GateSet(["x", "h"]), mode="basic", output_qasm="out.qasm", **MILP) + assert res + assert synthesized_gates(circuit) == ["x"] + + +@pytest.mark.parametrize("mode", ["basic", "incremental"]) +def test_milp_input_with_a_large_k(workdir, mode): + # the intermediate states' k was capped at 2*(depth+1), below the input's own k of 6 + inp = Vector(q=2, element_representation=FiveTuple, k=6) + out = Vector(q=2, element_representation=FiveTuple, k=6) + inp[0] = FiveTuple(a=1, b=0, c=0, d=0) + out[1] = FiveTuple(a=1, b=0, c=0, d=0) + res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=[(inp, out)], q=1, d=1, + gate_set=GateSet(["x", "h"]), mode=mode, output_qasm="out.qasm", **MILP) + assert res + assert synthesized_gates(circuit) == ["x"] + + +@pytest.mark.parametrize("gates, expected", [(["x", "z", "t"], 0), (["h", "t"], 1), (["h", "sx"], 2), (["cx", "csx"], 2)]) +def test_k_increment_follows_the_gate_set(gates, expected): + # the k bounds of the encoding come from this, a gate missing from k_increments would make them too small + assert max_k_increment(GateSet(gates)) == expected From 594b5424f4561a403cebaa50a51a579eba84a282 Mon Sep 17 00:00:00 2001 From: Ondra Lengal Date: Thu, 22 Oct 2026 14:23:35 +0200 Subject: [PATCH 5/5] fix: declare the target's rescaled vector at every depth The FiveTuple rescaling of the target was named after the target alone, so in incremental MILP search, every depth got back the PuLP variables of the first depth, with that depth's bound. Since the empty circuit is tried first (#67), that is depth 0, whose gap is the smallest: CBC no longer found h;sx;sx;h for Z from {h, sx}. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_0177AhE9SXm3iJ1KP8GTSk6o --- src/synth.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/synth.py b/src/synth.py index cb819ed..c820ca4 100644 --- a/src/synth.py +++ b/src/synth.py @@ -730,7 +730,9 @@ def encode_equivalence(self, vectors1, vectors2, depth, already_measured : bool value_bound = max(vec1.bound, vec_target.bound) rescaled_bound = 2 * value_bound * 2**(max_gap // 2) rescaled1 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled1_{vec1.name}", bound=rescaled_bound) - rescaled2 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled2_{vec2.name}", bound=rescaled_bound) + # vec2 is the target, the same at every depth: without the depth in the name, MILP would get back the + # first depth's PuLP variables, with that depth's (smallest) bound + rescaled2 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled2_{vec2.name}_{depth}", bound=rescaled_bound) if self.gen.mode == "milp" or self.gen.mode == "gurobi": self.gen.add_milp_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, max_gap, value_bound + rescaled_bound) else: