Skip to content

Circuits

qsp_proc.circuits

Public circuit builders and utilities for QSP constructions.

build_controlled_qsp_circuit(result, oracle, control_qubits=1, num_system_qubits=None)

Build a controlled QSP circuit following Appendix E.1.

Qubit layout is: [controls(0..c-1), ancilla(c), system(c+1..c+n)].

All signal-processing unitaries (E_0, V_j, V_j^\dagger) are controlled on the extra control register. The signal oracle is applied only on [ancilla, system] (i.e., no additional controls), enabling oracle sharing across controlled branches.

Parameters:

Name Type Description Default
result DecompositionResult

Recursive-carving decomposition containing E_0, projectors, and optional carve_sides metadata.

 required
oracle QuantumCircuit | ndarray

Signal oracle U as either a Qiskit circuit (on system qubits) or a unitary matrix.

 required
control_qubits int

Number of additional control qubits.

 1
num_system_qubits int | None

Number of system qubits. If None, infer from oracle.

 None

Returns:

Type Description
QuantumCircuit

A Qiskit QuantumCircuit implementing the controlled QSP sequence.

Raises:

Type Description
ValueError

If input dimensions or qubit counts are invalid.
TypeError

If oracle type is unsupported.
Source code in src\qsp_proc\circuits\controlled.py 
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def build_controlled_qsp_circuit(
    result: DecompositionResult,
    oracle: QuantumCircuit | np.ndarray,
    control_qubits: int = 1,
    num_system_qubits: int | None = None,
) -> QuantumCircuit:
    r"""Build a controlled QSP circuit following Appendix E.1.

    Qubit layout is:
    ``[controls(0..c-1), ancilla(c), system(c+1..c+n)]``.

    All signal-processing unitaries (``E_0``, ``V_j``, ``V_j^\dagger``) are
    controlled on the extra control register. The signal oracle is applied only
    on ``[ancilla, system]`` (i.e., no additional controls), enabling oracle
    sharing across controlled branches.

    Args:
        result: Recursive-carving decomposition containing ``E_0``, projectors,
            and optional ``carve_sides`` metadata.
        oracle: Signal oracle ``U`` as either a Qiskit circuit (on system
            qubits) or a unitary matrix.
        control_qubits: Number of additional control qubits.
        num_system_qubits: Number of system qubits. If ``None``, infer from
            ``oracle``.

    Returns:
        A Qiskit ``QuantumCircuit`` implementing the controlled QSP sequence.

    Raises:
        ValueError: If input dimensions or qubit counts are invalid.
        TypeError: If ``oracle`` type is unsupported.
    """
    if control_qubits < 0:
        raise ValueError("control_qubits must be non-negative.")

    u_gate, inferred_system_qubits = _coerce_oracle_to_gate(oracle)
    system_qubit_count = (
        inferred_system_qubits if num_system_qubits is None else int(num_system_qubits)
    )
    if system_qubit_count != inferred_system_qubits:
        raise ValueError(
            f"num_system_qubits={num_system_qubits} does not match oracle size "
            f"({inferred_system_qubits})."
        )
    if system_qubit_count <= 0:
        raise ValueError("num_system_qubits must be positive.")

    # Global layout: [controls..., ancilla, system...]
    total_qubits = control_qubits + 1 + system_qubit_count
    qc = QuantumCircuit(total_qubits, name="controlled_qsp")
    control_wires = list(range(control_qubits))
    ancilla_wire = control_qubits
    system_wires = list(range(control_qubits + 1, total_qubits))
    controlled_signal_wires = control_wires + [ancilla_wire]
    oracle_wires = [ancilla_wire] + system_wires

    # Oracle blocks (no extra controls beyond ancilla).
    cu_gate = u_gate.control(1)
    cu_dagger_gate = u_gate.inverse().control(1)

    # Controlled E0.
    c_e0 = _controlled_signal_processing_gate(
        np.asarray(result.e0, dtype=np.complex128).T,
        num_controls=control_qubits,
        label="cE0",
    )
    qc.append(c_e0, controlled_signal_wires)

    projectors = [np.asarray(projector, dtype=np.complex128) for projector in result.projectors]
    carve_sides = getattr(result, "carve_sides", None)
    side_sequence = _normalized_carve_sides(carve_sides, expected_length=len(projectors))
    use_carve_sides = side_sequence is not None

    # Recursive carving stores projectors in extraction order; application to
    # realize the original polynomial requires reversed order.
    iter_projectors = list(reversed(projectors))
    iter_sides = list(reversed(side_sequence)) if use_carve_sides else None

    # Controlled signal-processing blocks with shared oracle calls.
    for step, projector in enumerate(iter_projectors):
        v = projector_to_su2(projector)[:, [1, 0]]
        c_v = _controlled_signal_processing_gate(v, num_controls=control_qubits, label="cV")
        c_v_dagger = _controlled_signal_processing_gate(
            v.conj().T,
            num_controls=control_qubits,
            label="cVdg",
        )

        qc.append(c_v, controlled_signal_wires)

        if use_carve_sides:
            # Match base QSP semantics: top-side carve uses U, bottom uses U^\dagger.
            use_dagger = not iter_sides[step]
        else:
            # Fallback alternation when side metadata is unavailable.
            use_dagger = (step % 2) == 1

        qc.append(cu_dagger_gate if use_dagger else cu_gate, oracle_wires)

        qc.append(c_v_dagger, controlled_signal_wires)

    return qc

build_lcu_circuit(results, oracle, num_system_qubits=None)

Build the general LCU QSP circuit (Theorem 1).

For L = len(results), this constructs an LCU register of l = ceil(log2(L)) ancillas, one QSP ancilla, and a system register. The circuit prepares a uniform superposition on the LCU ancillas, applies branch-local single-qubit basis changes conditioned on |k> with interleaved/shared oracle layers across branches, pads any unused LCU basis states with explicit zero-contribution branches, then uncomputes with Hadamards. Post-selection on |0^l> yields the subnormalized block (1 / 2^l) * sum_k P_k.

Parameters:

Name Type Description Default
results list[DecompositionResult]

Decomposition outputs, one per polynomial/QSP sequence in the LCU sum.

 required
oracle QuantumCircuit | ndarray

Signal oracle as a Qiskit circuit or unitary matrix.

 required
num_system_qubits int | None

Number of system qubits. If None, infer from the oracle where possible.

 None

Returns:

Type Description
QuantumCircuit

A Qiskit circuit implementing the LCU-composed QSP construction.

Raises:

Type Description
ValueError

If inputs are empty, dimensions are invalid, or qubit counts are inconsistent.
TypeError

If oracle type is unsupported.
Source code in src\qsp_proc\circuits\lcu.py 
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
def build_lcu_circuit(
    results: list[DecompositionResult],
    oracle: QuantumCircuit | np.ndarray,
    num_system_qubits: int | None = None,
) -> QuantumCircuit:
    """Build the general LCU QSP circuit (Theorem 1).

    For ``L = len(results)``, this constructs an LCU register of
    ``l = ceil(log2(L))`` ancillas, one QSP ancilla, and a system register.
    The circuit prepares a uniform superposition on the LCU ancillas, applies
    branch-local single-qubit basis changes conditioned on ``|k>`` with
    interleaved/shared oracle layers across branches, pads any unused LCU basis
    states with explicit zero-contribution branches, then uncomputes with
    Hadamards. Post-selection on ``|0^l>`` yields the subnormalized block
    ``(1 / 2^l) * sum_k P_k``.

    Args:
        results: Decomposition outputs, one per polynomial/QSP sequence in the
            LCU sum.
        oracle: Signal oracle as a Qiskit circuit or unitary matrix.
        num_system_qubits: Number of system qubits. If ``None``, infer from the
            oracle where possible.

    Returns:
        A Qiskit circuit implementing the LCU-composed QSP construction.

    Raises:
        ValueError: If inputs are empty, dimensions are invalid, or qubit counts
            are inconsistent.
        TypeError: If ``oracle`` type is unsupported.
    """
    from qsp_proc.circuits.qsp_circuit import projector_to_su2

    if len(results) == 0:
        raise ValueError("results must contain at least one decomposition result.")

    num_system_qubits = _infer_system_qubits(
        oracle=oracle,
        num_system_qubits=num_system_qubits,
    )
    if num_system_qubits < 1:
        raise ValueError("num_system_qubits must be at least 1.")

    # Build oracle and controlled-oracle gates (controlled by QSP ancilla only).
    u_gate, u_dagger_gate = _oracle_gate_pair(
        oracle=oracle,
        num_system_qubits=num_system_qubits,
    )
    u_dagger_gate.label = "Udg"
    cu_gate = u_gate.control(1, label="CU")
    cu_dagger_gate = u_dagger_gate.control(1, label="CUdg")

    lcu_ancilla_count = int(np.ceil(np.log2(len(results))))
    total_qubits = lcu_ancilla_count + 1 + int(num_system_qubits)
    qsp_ancilla_wire = lcu_ancilla_count
    lcu_wires = list(range(lcu_ancilla_count))
    oracle_wires = [qsp_ancilla_wire] + list(range(lcu_ancilla_count + 1, total_qubits))

    qc = QuantumCircuit(total_qubits, name="LCU-QSP")

    def _append_branch_controlled_single_qubit(
        unitary: np.ndarray,
        branch: int,
        label: str,
    ) -> None:
        """Append a single-qubit unitary on QSP ancilla controlled on LCU |branch>."""
        gate = UnitaryGate(np.asarray(unitary, dtype=np.complex128), label=label)
        if lcu_ancilla_count == 0:
            qc.append(gate, [qsp_ancilla_wire])
            return
        ctrl_state = format(branch, f"0{lcu_ancilla_count}b")
        qc.append(
            gate.control(lcu_ancilla_count, ctrl_state=ctrl_state),
            lcu_wires + [qsp_ancilla_wire],
        )

    # Prepare uniform superposition over LCU ancillas.
    for wire in lcu_wires:
        qc.h(wire)

    total_branches = 1 << lcu_ancilla_count
    branch_steps: list[list[tuple[np.ndarray, bool]]] = []

    # Branch-local setup and per-layer data extraction.
    for k, result in enumerate(results):
        e0 = np.asarray(result.e0, dtype=np.complex128)
        if e0.shape != (2, 2):
            raise ValueError(f"results[{k}].e0 must have shape (2, 2), got {e0.shape}.")
        _append_branch_controlled_single_qubit(e0.T, k, label="E0")

        projectors = [np.asarray(projector, dtype=np.complex128) for projector in result.projectors]
        side_sequence = _normalized_carve_sides(
            getattr(result, "carve_sides", None),
            expected_length=len(projectors),
            context=f"results[{k}]",
        )
        use_carve_sides = side_sequence is not None

        # Recursive carving stores extraction order; apply reversed for synthesis.
        iter_projectors = list(reversed(projectors))
        iter_sides = list(reversed(side_sequence)) if use_carve_sides else None

        steps_for_branch: list[tuple[np.ndarray, bool]] = []
        for step, projector in enumerate(iter_projectors):
            v = projector_to_su2(projector)[:, [1, 0]]
            if use_carve_sides:
                use_dagger = not iter_sides[step]
            else:
                use_dagger = (step % 2) == 1
            steps_for_branch.append((v, use_dagger))
        branch_steps.append(steps_for_branch)

    # Pad non-power-of-two LCU spaces with explicit zero-contribution branches.
    # A branch-local X on the QSP ancilla makes the ancilla-|0> postselected block
    # exactly zero, avoiding a phantom identity contribution from "do nothing".
    for k in range(len(results), total_branches):
        _append_branch_controlled_single_qubit(
            np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex128),
            k,
            label="E0pad",
        )

    # Interleave branch layers so oracle calls are shared across branches.
    max_steps = max((len(steps) for steps in branch_steps), default=0)
    side_aligned = True
    for global_step in range(max_steps):
        step_sides = {
            steps[global_step][1] for steps in branch_steps if global_step < len(steps)
        }
        if len(step_sides) > 1:
            side_aligned = False
            break

    if side_aligned:
        for global_step in range(max_steps):
            active_branches = [
                branch for branch, steps in enumerate(branch_steps) if global_step < len(steps)
            ]
            if not active_branches:
                continue
            use_dagger = branch_steps[active_branches[0]][global_step][1]
            for branch in active_branches:
                v, _ = branch_steps[branch][global_step]
                _append_branch_controlled_single_qubit(v, branch, label="V")
            qc.append(cu_dagger_gate if use_dagger else cu_gate, oracle_wires)
            for branch in active_branches:
                v, _ = branch_steps[branch][global_step]
                _append_branch_controlled_single_qubit(v.conj().T, branch, label="Vdg")
    else:
        # Safety fallback for custom, mismatched side schedules.
        # Without side alignment, a fully shared oracle layer can change branch
        # dynamics. Revert to per-branch sequencing to preserve behavior.
        for branch, steps in enumerate(branch_steps):
            for v, use_dagger in steps:
                _append_branch_controlled_single_qubit(v, branch, label="V")
                qc.append(cu_dagger_gate if use_dagger else cu_gate, oracle_wires)
                _append_branch_controlled_single_qubit(v.conj().T, branch, label="Vdg")

    # Uncompute LCU superposition.
    for wire in lcu_wires:
        qc.h(wire)

    return qc

build_real_extraction_circuit(result, result_conj, oracle, num_system_qubits=None)

Build the Appendix E.2 real-part extraction circuit (P + P*) / 2.

Qubit layout is fixed as: - qubit 0: LCU ancilla, - qubit 1: QSP ancilla, - qubits 2..: system register.

The circuit implements the two-branch LCU pattern (Figure 10): H on LCU ancilla, controlled-P on |0>, controlled-P* on |1>, final H, and measurement of the LCU ancilla.

Parameters:

Name Type Description Default
result DecompositionResult

Decomposition result for P.

 required
result_conj DecompositionResult

Decomposition result for the conjugate branch P*.

 required
oracle QuantumCircuit | ndarray

Signal oracle as a Qiskit circuit or unitary matrix.

 required
num_system_qubits int | None

Number of system qubits. If None, infer from the oracle where possible.

 None

Returns:

Type Description
QuantumCircuit

A Qiskit circuit implementing real-part extraction via post-selection
QuantumCircuit

on measured LCU ancilla |0>.

Raises:

Type Description
ValueError

If num_system_qubits cannot be inferred or is invalid.
Source code in src\qsp_proc\circuits\lcu.py 
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
def build_real_extraction_circuit(
    result: DecompositionResult,
    result_conj: DecompositionResult,
    oracle: QuantumCircuit | np.ndarray,
    num_system_qubits: int | None = None,
) -> QuantumCircuit:
    """Build the Appendix E.2 real-part extraction circuit ``(P + P*) / 2``.

    Qubit layout is fixed as:
    - qubit ``0``: LCU ancilla,
    - qubit ``1``: QSP ancilla,
    - qubits ``2..``: system register.

    The circuit implements the two-branch LCU pattern (Figure 10):
    ``H`` on LCU ancilla, controlled-``P`` on ``|0>``, controlled-``P*`` on
    ``|1>``, final ``H``, and measurement of the LCU ancilla.

    Args:
        result: Decomposition result for ``P``.
        result_conj: Decomposition result for the conjugate branch ``P*``.
        oracle: Signal oracle as a Qiskit circuit or unitary matrix.
        num_system_qubits: Number of system qubits. If ``None``, infer from the
            oracle where possible.

    Returns:
        A Qiskit circuit implementing real-part extraction via post-selection
        on measured LCU ancilla ``|0>``.

    Raises:
        ValueError: If ``num_system_qubits`` cannot be inferred or is invalid.
    """
    from qsp_proc.circuits.controlled import build_controlled_qsp_circuit

    num_system_qubits = _infer_system_qubits(
        oracle=oracle,
        num_system_qubits=num_system_qubits,
    )

    if num_system_qubits < 1:
        raise ValueError("num_system_qubits must be at least 1.")

    total_qubits = 2 + num_system_qubits
    full_wires = list(range(total_qubits))

    # Controlled QSP branches (control qubit is wire 0 in each subcircuit).
    controlled_p = build_controlled_qsp_circuit(
        result=result,
        oracle=oracle,
        control_qubits=1,
        num_system_qubits=num_system_qubits,
    )
    controlled_p_conj = build_controlled_qsp_circuit(
        result=result_conj,
        oracle=oracle,
        control_qubits=1,
        num_system_qubits=num_system_qubits,
    )

    circuit = QuantumCircuit(total_qubits, 1, name="real_extraction_lcu")

    # Prepare equal LCU superposition.
    circuit.h(0)

    # Apply controlled-P on |0> branch (flip control basis with X).
    circuit.x(0)
    circuit.compose(controlled_p, qubits=full_wires, inplace=True)
    circuit.x(0)

    # Apply controlled-P* on |1> branch.
    circuit.compose(controlled_p_conj, qubits=full_wires, inplace=True)

    # Interfere branches and measure LCU ancilla for post-selection.
    circuit.h(0)
    circuit.measure(0, 0)

    return circuit

build_qsp_circuit(result, oracle, num_system_qubits=None)

Build the basic QSP circuit from a decomposition result.

Constructs the Figure-8-style ancilla + system circuit using the decomposition data (e0, projectors, carve_sides) and a signal oracle U.

Parameters:

Name Type Description Default
result DecompositionResult

Decomposition output containing e0, projector sequence, and carving side metadata.

 required
oracle QuantumCircuit | ndarray

Signal operator U as either a :class:qiskit.QuantumCircuit or a unitary matrix.

 required
num_system_qubits int | None

Number of system qubits. Inferred from oracle when oracle is a :class:qiskit.QuantumCircuit.

 None

Returns:

Name Type Description
A QuantumCircuit

class:qiskit.QuantumCircuit with one ancilla qubit plus system qubits.
Source code in src\qsp_proc\circuits\qsp_circuit.py 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
def build_qsp_circuit(
    result: DecompositionResult,
    oracle: QuantumCircuit | np.ndarray,
    num_system_qubits: int | None = None,
) -> QuantumCircuit:
    """Build the basic QSP circuit from a decomposition result.

    Constructs the Figure-8-style ancilla + system circuit using the decomposition
    data ``(e0, projectors, carve_sides)`` and a signal oracle ``U``.

    Args:
        result: Decomposition output containing ``e0``, projector sequence, and
            carving side metadata.
        oracle: Signal operator ``U`` as either a :class:`qiskit.QuantumCircuit`
            or a unitary matrix.
        num_system_qubits: Number of system qubits. Inferred from ``oracle`` when
            ``oracle`` is a :class:`qiskit.QuantumCircuit`.

    Returns:
        A :class:`qiskit.QuantumCircuit` with one ancilla qubit plus system qubits.
    """
    oracle_circuit = _coerce_oracle_circuit(
        oracle=oracle,
        num_system_qubits=num_system_qubits,
    )
    inferred_qubits = oracle_circuit.num_qubits
    oracle_dagger_circuit = oracle_circuit.inverse()
    system_wires = list(range(1, inferred_qubits + 1))
    ancilla_wire = 0

    qc = QuantumCircuit(1 + inferred_qubits, name="qsp")

    e0 = np.asarray(result.e0, dtype=np.complex128)
    if e0.shape != (2, 2):
        raise ValueError(f"result.e0 must have shape (2, 2), got {e0.shape}.")
    # Decomposition matrices are represented in coefficient-space convention;
    # transpose maps them to the ancilla gate action convention.
    qc.unitary(e0.T, [ancilla_wire], label="E0")

    projectors = [np.asarray(projector, dtype=np.complex128) for projector in result.projectors]
    side_sequence = _normalized_carve_sides(
        getattr(result, "carve_sides", None),
        expected_length=len(projectors),
    )
    use_carve_sides = side_sequence is not None

    # Recursive carving stores projectors in extraction order; applying factors to
    # realize the original polynomial requires reverse order (matching reconstruction).
    iter_projectors = list(reversed(projectors))
    iter_sides = list(reversed(side_sequence)) if use_carve_sides else None

    cu_oracle = oracle_circuit.control(1)
    cu_oracle_dagger = oracle_dagger_circuit.control(1)

    for step, projector in enumerate(iter_projectors):
        # projector_to_su2 maps P to |0><0|; controlled-U triggers on |1>, so
        # swap columns to map the signal projector to ancilla |1>.
        vj = projector_to_su2(projector)[:, [1, 0]]
        qc.unitary(vj, [ancilla_wire], label=f"V{step + 1}")

        if use_carve_sides:
            from_top = bool(iter_sides[step])
            # Top-side carve corresponds to E_P(z)=Q+zP; bottom carve uses z^{-1}.
            use_dagger = not from_top
        else:
            # Fallback alternation when side metadata is unavailable.
            use_dagger = (step % 2) == 1

        qc.append(
            cu_oracle_dagger if use_dagger else cu_oracle,
            [ancilla_wire] + system_wires,
        )

        qc.unitary(vj.conj().T, [ancilla_wire], label=f"V{step + 1}dg")

    return qc

verify_qsp_circuit(circuit, target_poly, thetas=None, num_samples=64, tol=1e-06, system_basis_index=0)

Verify a QSP circuit against a target Laurent polynomial on sampled angles.

The circuit is evaluated pointwise in theta (binding its single parameter if present), and the achieved polynomial is compared to target_poly(e^{i theta}).

Parameters:

Name Type Description Default
circuit QuantumCircuit

QSP circuit to verify (optionally parameterized by theta).

 required
target_poly LaurentPoly

Target Laurent polynomial.

 required
thetas ndarray | None

Optional explicit theta grid (radians). If None, a uniform grid on [0, 2Ï€) is used.

 None
num_samples int

Number of uniform samples when thetas is None.

 64
tol float

Non-negative tolerance used for input validation and optional logging.

 1e-06
system_basis_index int

Computational-basis index within the postselected ancilla-|0> system block used for scalar extraction.

 0

Returns:

Type Description
float

Maximum absolute error over sampled points.
Source code in src\qsp_proc\circuits\qsp_circuit.py 
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
def verify_qsp_circuit(
    circuit: QuantumCircuit,
    target_poly: LaurentPoly,
    thetas: np.ndarray | None = None,
    num_samples: int = 64,
    tol: float = 1e-6,
    system_basis_index: int = 0,
) -> float:
    """Verify a QSP circuit against a target Laurent polynomial on sampled angles.

    The circuit is evaluated pointwise in ``theta`` (binding its single parameter if
    present), and the achieved polynomial is compared to ``target_poly(e^{i theta})``.

    Args:
        circuit: QSP circuit to verify (optionally parameterized by ``theta``).
        target_poly: Target Laurent polynomial.
        thetas: Optional explicit theta grid (radians). If ``None``, a uniform grid
            on ``[0, 2Ï€)`` is used.
        num_samples: Number of uniform samples when ``thetas`` is ``None``.
        tol: Non-negative tolerance used for input validation and optional logging.
        system_basis_index: Computational-basis index within the postselected
            ancilla-``|0>`` system block used for scalar extraction.

    Returns:
        Maximum absolute error over sampled points.
    """
    if num_samples <= 0:
        raise ValueError("num_samples must be positive.")
    if tol < 0.0:
        raise ValueError("tol must be non-negative.")

    if thetas is None:
        theta_array = np.linspace(0.0, 2.0 * np.pi, num_samples, endpoint=False, dtype=np.float64)
    else:
        theta_array = np.asarray(thetas, dtype=np.float64).reshape(-1)
        if theta_array.size == 0:
            raise ValueError("thetas must contain at least one sample.")

    achieved = extract_polynomial_from_circuit(
        circuit=circuit,
        thetas=theta_array,
        system_basis_index=system_basis_index,
    )

    z = np.exp(1j * theta_array)
    target_vals = np.asarray(target_poly(z), dtype=np.complex128).reshape(-1)
    if target_vals.size == 1 and achieved.size > 1:
        target_vals = np.full_like(achieved, target_vals.item(), dtype=np.complex128)
    if target_vals.shape != achieved.shape:
        raise ValueError(
            "target_poly evaluation shape mismatch: "
            f"got {target_vals.shape}, expected {achieved.shape}."
        )

    max_error = float(np.max(np.abs(achieved - target_vals)))
    return max_error

projector_to_su2(projector)

Build the ancilla-basis unitary V associated with a rank-1 projector.

Given P = |p><p| on C^2, this returns the 2x2 unitary V = [p, p_perp] (columns), where p_perp is any normalized vector orthogonal to p. Therefore V |0> = |p> and P = V |0><0| V^\\dagger.

Parameters:

Name Type Description Default
projector ndarray

Rank-1 projector with shape (2, 2).

 required

Returns:

Type Description
ndarray

A 2x2 unitary matrix V as np.ndarray.

Raises:

Type Description
ValueError

If projector is not a valid rank-1 2x2 projector.
Source code in src\qsp_proc\circuits\qsp_circuit.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def projector_to_su2(projector: np.ndarray) -> np.ndarray:
    r"""Build the ancilla-basis unitary ``V`` associated with a rank-1 projector.

    Given ``P = |p><p|`` on ``C^2``, this returns the ``2x2`` unitary
    ``V = [p, p_perp]`` (columns), where ``p_perp`` is any normalized vector
    orthogonal to ``p``. Therefore ``V |0> = |p>`` and
    ``P = V |0><0| V^\\dagger``.

    Args:
        projector: Rank-1 projector with shape ``(2, 2)``.

    Returns:
        A ``2x2`` unitary matrix ``V`` as ``np.ndarray``.

    Raises:
        ValueError: If ``projector`` is not a valid rank-1 ``2x2`` projector.
    """
    proj = np.asarray(projector, dtype=np.complex128)

    if proj.shape != (2, 2):
        raise ValueError("projector must have shape (2, 2).")

    tol = 1e-10
    if not np.allclose(proj, proj.conj().T, atol=tol):
        raise ValueError("projector must be Hermitian.")
    if not np.allclose(proj @ proj, proj, atol=tol):
        raise ValueError("projector must be idempotent (P @ P = P).")

    evals, evecs = np.linalg.eigh(proj)
    idx = int(np.argmax(evals.real))
    if not np.isclose(evals[idx].real, 1.0, atol=1e-8):
        raise ValueError("projector must be rank-1 with eigenvalue spectrum {0, 1}.")

    p = evecs[:, idx]
    p = p / np.linalg.norm(p)

    # Canonical orthogonal complement in C^2.
    p_perp = np.array([-np.conj(p[1]), np.conj(p[0])], dtype=np.complex128)
    p_perp = p_perp / np.linalg.norm(p_perp)

    v = np.column_stack((p, p_perp))
    if not np.allclose(v.conj().T @ v, np.eye(2, dtype=np.complex128), atol=1e-10):
        raise ValueError("failed to construct a unitary basis from projector.")

    return v

qsp_proc.circuits.qsp_circuit

Basic QSP circuit construction (Figure 8, Skelton 2025).

This module provides interfaces for converting decomposition outputs (projectors and constant residue) into an executable ancilla-system QSP circuit.

projector_to_su2(projector)

Build the ancilla-basis unitary V associated with a rank-1 projector.

Given P = |p><p| on C^2, this returns the 2x2 unitary V = [p, p_perp] (columns), where p_perp is any normalized vector orthogonal to p. Therefore V |0> = |p> and P = V |0><0| V^\\dagger.

Parameters:

Name Type Description Default
projector ndarray

Rank-1 projector with shape (2, 2).

 required

Returns:

Type Description
ndarray

A 2x2 unitary matrix V as np.ndarray.

Raises:

Type Description
ValueError

If projector is not a valid rank-1 2x2 projector.
Source code in src\qsp_proc\circuits\qsp_circuit.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def projector_to_su2(projector: np.ndarray) -> np.ndarray:
    r"""Build the ancilla-basis unitary ``V`` associated with a rank-1 projector.

    Given ``P = |p><p|`` on ``C^2``, this returns the ``2x2`` unitary
    ``V = [p, p_perp]`` (columns), where ``p_perp`` is any normalized vector
    orthogonal to ``p``. Therefore ``V |0> = |p>`` and
    ``P = V |0><0| V^\\dagger``.

    Args:
        projector: Rank-1 projector with shape ``(2, 2)``.

    Returns:
        A ``2x2`` unitary matrix ``V`` as ``np.ndarray``.

    Raises:
        ValueError: If ``projector`` is not a valid rank-1 ``2x2`` projector.
    """
    proj = np.asarray(projector, dtype=np.complex128)

    if proj.shape != (2, 2):
        raise ValueError("projector must have shape (2, 2).")

    tol = 1e-10
    if not np.allclose(proj, proj.conj().T, atol=tol):
        raise ValueError("projector must be Hermitian.")
    if not np.allclose(proj @ proj, proj, atol=tol):
        raise ValueError("projector must be idempotent (P @ P = P).")

    evals, evecs = np.linalg.eigh(proj)
    idx = int(np.argmax(evals.real))
    if not np.isclose(evals[idx].real, 1.0, atol=1e-8):
        raise ValueError("projector must be rank-1 with eigenvalue spectrum {0, 1}.")

    p = evecs[:, idx]
    p = p / np.linalg.norm(p)

    # Canonical orthogonal complement in C^2.
    p_perp = np.array([-np.conj(p[1]), np.conj(p[0])], dtype=np.complex128)
    p_perp = p_perp / np.linalg.norm(p_perp)

    v = np.column_stack((p, p_perp))
    if not np.allclose(v.conj().T @ v, np.eye(2, dtype=np.complex128), atol=1e-10):
        raise ValueError("failed to construct a unitary basis from projector.")

    return v

build_qsp_circuit(result, oracle, num_system_qubits=None)

Build the basic QSP circuit from a decomposition result.

Constructs the Figure-8-style ancilla + system circuit using the decomposition data (e0, projectors, carve_sides) and a signal oracle U.

Parameters:

Name Type Description Default
result DecompositionResult

Decomposition output containing e0, projector sequence, and carving side metadata.

 required
oracle QuantumCircuit | ndarray

Signal operator U as either a :class:qiskit.QuantumCircuit or a unitary matrix.

 required
num_system_qubits int | None

Number of system qubits. Inferred from oracle when oracle is a :class:qiskit.QuantumCircuit.

 None

Returns:

Name Type Description
A QuantumCircuit

class:qiskit.QuantumCircuit with one ancilla qubit plus system qubits.
Source code in src\qsp_proc\circuits\qsp_circuit.py 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
def build_qsp_circuit(
    result: DecompositionResult,
    oracle: QuantumCircuit | np.ndarray,
    num_system_qubits: int | None = None,
) -> QuantumCircuit:
    """Build the basic QSP circuit from a decomposition result.

    Constructs the Figure-8-style ancilla + system circuit using the decomposition
    data ``(e0, projectors, carve_sides)`` and a signal oracle ``U``.

    Args:
        result: Decomposition output containing ``e0``, projector sequence, and
            carving side metadata.
        oracle: Signal operator ``U`` as either a :class:`qiskit.QuantumCircuit`
            or a unitary matrix.
        num_system_qubits: Number of system qubits. Inferred from ``oracle`` when
            ``oracle`` is a :class:`qiskit.QuantumCircuit`.

    Returns:
        A :class:`qiskit.QuantumCircuit` with one ancilla qubit plus system qubits.
    """
    oracle_circuit = _coerce_oracle_circuit(
        oracle=oracle,
        num_system_qubits=num_system_qubits,
    )
    inferred_qubits = oracle_circuit.num_qubits
    oracle_dagger_circuit = oracle_circuit.inverse()
    system_wires = list(range(1, inferred_qubits + 1))
    ancilla_wire = 0

    qc = QuantumCircuit(1 + inferred_qubits, name="qsp")

    e0 = np.asarray(result.e0, dtype=np.complex128)
    if e0.shape != (2, 2):
        raise ValueError(f"result.e0 must have shape (2, 2), got {e0.shape}.")
    # Decomposition matrices are represented in coefficient-space convention;
    # transpose maps them to the ancilla gate action convention.
    qc.unitary(e0.T, [ancilla_wire], label="E0")

    projectors = [np.asarray(projector, dtype=np.complex128) for projector in result.projectors]
    side_sequence = _normalized_carve_sides(
        getattr(result, "carve_sides", None),
        expected_length=len(projectors),
    )
    use_carve_sides = side_sequence is not None

    # Recursive carving stores projectors in extraction order; applying factors to
    # realize the original polynomial requires reverse order (matching reconstruction).
    iter_projectors = list(reversed(projectors))
    iter_sides = list(reversed(side_sequence)) if use_carve_sides else None

    cu_oracle = oracle_circuit.control(1)
    cu_oracle_dagger = oracle_dagger_circuit.control(1)

    for step, projector in enumerate(iter_projectors):
        # projector_to_su2 maps P to |0><0|; controlled-U triggers on |1>, so
        # swap columns to map the signal projector to ancilla |1>.
        vj = projector_to_su2(projector)[:, [1, 0]]
        qc.unitary(vj, [ancilla_wire], label=f"V{step + 1}")

        if use_carve_sides:
            from_top = bool(iter_sides[step])
            # Top-side carve corresponds to E_P(z)=Q+zP; bottom carve uses z^{-1}.
            use_dagger = not from_top
        else:
            # Fallback alternation when side metadata is unavailable.
            use_dagger = (step % 2) == 1

        qc.append(
            cu_oracle_dagger if use_dagger else cu_oracle,
            [ancilla_wire] + system_wires,
        )

        qc.unitary(vj.conj().T, [ancilla_wire], label=f"V{step + 1}dg")

    return qc

extract_polynomial_from_circuit(circuit, thetas, system_basis_index=0)

Extract achieved QSP polynomial values from a (possibly parameterized) circuit.

For each theta, this function binds the circuit's single free parameter (if present), computes the full unitary via :class:qiskit.quantum_info.Operator, extracts the ancilla |0> -> |0> block, and returns one diagonal system-basis entry from that block. This avoids trace-averaging across different oracle eigenvalues, which is incorrect for non-symmetric Laurent polynomials.

Parameters:

Name Type Description Default
circuit QuantumCircuit

QSP circuit, optionally parameterized by a single angle parameter.

 required
thetas ndarray

Sample points in radians.

 required
system_basis_index int

Computational-basis index within the postselected system block to read. Must satisfy 0 <= system_basis_index < 2^n for n system qubits.

 0

Returns:

Type Description
ndarray

Complex array of achieved polynomial values at each sample.
Source code in src\qsp_proc\circuits\qsp_circuit.py 
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
def extract_polynomial_from_circuit(
    circuit: QuantumCircuit,
    thetas: np.ndarray,
    system_basis_index: int = 0,
) -> np.ndarray:
    """Extract achieved QSP polynomial values from a (possibly parameterized) circuit.

    For each ``theta``, this function binds the circuit's single free parameter
    (if present), computes the full unitary via :class:`qiskit.quantum_info.Operator`,
    extracts the ancilla ``|0> -> |0>`` block, and returns one diagonal system-basis
    entry from that block. This avoids trace-averaging across different oracle
    eigenvalues, which is incorrect for non-symmetric Laurent polynomials.

    Args:
        circuit: QSP circuit, optionally parameterized by a single angle parameter.
        thetas: Sample points in radians.
        system_basis_index: Computational-basis index within the postselected system
            block to read. Must satisfy ``0 <= system_basis_index < 2^n`` for
            ``n`` system qubits.

    Returns:
        Complex array of achieved polynomial values at each sample.
    """
    theta_input = np.asarray(thetas, dtype=np.float64)
    if theta_input.size == 0:
        raise ValueError("thetas must contain at least one sample.")
    theta_array = theta_input.reshape(-1)

    params = tuple(circuit.parameters)
    if len(params) > 1:
        raise ValueError(
            "extract_polynomial_from_circuit supports circuits with at most one free "
            f"parameter, got {len(params)}."
        )

    if circuit.num_qubits <= 0:
        raise ValueError("circuit must contain at least one qubit (the ancilla).")

    num_system_qubits = circuit.num_qubits - 1
    dim_system = 1 << num_system_qubits
    if system_basis_index < 0 or system_basis_index >= dim_system:
        raise ValueError(
            "system_basis_index out of range: "
            f"got {system_basis_index}, valid range is [0, {dim_system - 1}]."
        )

    values = np.empty(theta_array.shape[0], dtype=np.complex128)

    for idx, theta in enumerate(theta_array):
        if len(params) == 1:
            bound_circuit = circuit.assign_parameters({params[0]: float(theta)}, inplace=False)
        else:
            bound_circuit = circuit

        unitary = Operator(bound_circuit).data
        dim = unitary.shape[0]

        # Ancilla is conventionally qubit 0 in this module. In Qiskit's basis
        # indexing, qubit 0 corresponds to the least-significant bit.
        basis_indices = np.arange(dim)
        ancilla_zero = basis_indices[(basis_indices & 1) == 0]
        block_00 = unitary[np.ix_(ancilla_zero, ancilla_zero)]

        values[idx] = block_00[system_basis_index, system_basis_index]

    return values.reshape(theta_input.shape)

verify_qsp_circuit(circuit, target_poly, thetas=None, num_samples=64, tol=1e-06, system_basis_index=0)

Verify a QSP circuit against a target Laurent polynomial on sampled angles.

The circuit is evaluated pointwise in theta (binding its single parameter if present), and the achieved polynomial is compared to target_poly(e^{i theta}).

Parameters:

Name Type Description Default
circuit QuantumCircuit

QSP circuit to verify (optionally parameterized by theta).

 required
target_poly LaurentPoly

Target Laurent polynomial.

 required
thetas ndarray | None

Optional explicit theta grid (radians). If None, a uniform grid on [0, 2Ï€) is used.

 None
num_samples int

Number of uniform samples when thetas is None.

 64
tol float

Non-negative tolerance used for input validation and optional logging.

 1e-06
system_basis_index int

Computational-basis index within the postselected ancilla-|0> system block used for scalar extraction.

 0

Returns:

Type Description
float

Maximum absolute error over sampled points.
Source code in src\qsp_proc\circuits\qsp_circuit.py 
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
def verify_qsp_circuit(
    circuit: QuantumCircuit,
    target_poly: LaurentPoly,
    thetas: np.ndarray | None = None,
    num_samples: int = 64,
    tol: float = 1e-6,
    system_basis_index: int = 0,
) -> float:
    """Verify a QSP circuit against a target Laurent polynomial on sampled angles.

    The circuit is evaluated pointwise in ``theta`` (binding its single parameter if
    present), and the achieved polynomial is compared to ``target_poly(e^{i theta})``.

    Args:
        circuit: QSP circuit to verify (optionally parameterized by ``theta``).
        target_poly: Target Laurent polynomial.
        thetas: Optional explicit theta grid (radians). If ``None``, a uniform grid
            on ``[0, 2Ï€)`` is used.
        num_samples: Number of uniform samples when ``thetas`` is ``None``.
        tol: Non-negative tolerance used for input validation and optional logging.
        system_basis_index: Computational-basis index within the postselected
            ancilla-``|0>`` system block used for scalar extraction.

    Returns:
        Maximum absolute error over sampled points.
    """
    if num_samples <= 0:
        raise ValueError("num_samples must be positive.")
    if tol < 0.0:
        raise ValueError("tol must be non-negative.")

    if thetas is None:
        theta_array = np.linspace(0.0, 2.0 * np.pi, num_samples, endpoint=False, dtype=np.float64)
    else:
        theta_array = np.asarray(thetas, dtype=np.float64).reshape(-1)
        if theta_array.size == 0:
            raise ValueError("thetas must contain at least one sample.")

    achieved = extract_polynomial_from_circuit(
        circuit=circuit,
        thetas=theta_array,
        system_basis_index=system_basis_index,
    )

    z = np.exp(1j * theta_array)
    target_vals = np.asarray(target_poly(z), dtype=np.complex128).reshape(-1)
    if target_vals.size == 1 and achieved.size > 1:
        target_vals = np.full_like(achieved, target_vals.item(), dtype=np.complex128)
    if target_vals.shape != achieved.shape:
        raise ValueError(
            "target_poly evaluation shape mismatch: "
            f"got {target_vals.shape}, expected {achieved.shape}."
        )

    max_error = float(np.max(np.abs(achieved - target_vals)))
    return max_error

qsp_proc.circuits.controlled

Controlled QSP circuit construction utilities (Skelton 2025, Appendix E.1).

This module defines interfaces for building controlled-QSP circuits where the QSP signal-processing operators (V_j and E_0) are controlled by additional control qubits, while the signal oracle is left unchanged with respect to those extra controls (it remains controlled only by the QSP ancilla, as in the base construction). This pattern supports parallel execution of multiple QSP sequences under shared oracle calls.

build_controlled_qsp_circuit(result, oracle, control_qubits=1, num_system_qubits=None)

Build a controlled QSP circuit following Appendix E.1.

Qubit layout is: [controls(0..c-1), ancilla(c), system(c+1..c+n)].

All signal-processing unitaries (E_0, V_j, V_j^\dagger) are controlled on the extra control register. The signal oracle is applied only on [ancilla, system] (i.e., no additional controls), enabling oracle sharing across controlled branches.

Parameters:

Name Type Description Default
result DecompositionResult

Recursive-carving decomposition containing E_0, projectors, and optional carve_sides metadata.

 required
oracle QuantumCircuit | ndarray

Signal oracle U as either a Qiskit circuit (on system qubits) or a unitary matrix.

 required
control_qubits int

Number of additional control qubits.

 1
num_system_qubits int | None

Number of system qubits. If None, infer from oracle.

 None

Returns:

Type Description
QuantumCircuit

A Qiskit QuantumCircuit implementing the controlled QSP sequence.

Raises:

Type Description
ValueError

If input dimensions or qubit counts are invalid.
TypeError

If oracle type is unsupported.
Source code in src\qsp_proc\circuits\controlled.py 
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def build_controlled_qsp_circuit(
    result: DecompositionResult,
    oracle: QuantumCircuit | np.ndarray,
    control_qubits: int = 1,
    num_system_qubits: int | None = None,
) -> QuantumCircuit:
    r"""Build a controlled QSP circuit following Appendix E.1.

    Qubit layout is:
    ``[controls(0..c-1), ancilla(c), system(c+1..c+n)]``.

    All signal-processing unitaries (``E_0``, ``V_j``, ``V_j^\dagger``) are
    controlled on the extra control register. The signal oracle is applied only
    on ``[ancilla, system]`` (i.e., no additional controls), enabling oracle
    sharing across controlled branches.

    Args:
        result: Recursive-carving decomposition containing ``E_0``, projectors,
            and optional ``carve_sides`` metadata.
        oracle: Signal oracle ``U`` as either a Qiskit circuit (on system
            qubits) or a unitary matrix.
        control_qubits: Number of additional control qubits.
        num_system_qubits: Number of system qubits. If ``None``, infer from
            ``oracle``.

    Returns:
        A Qiskit ``QuantumCircuit`` implementing the controlled QSP sequence.

    Raises:
        ValueError: If input dimensions or qubit counts are invalid.
        TypeError: If ``oracle`` type is unsupported.
    """
    if control_qubits < 0:
        raise ValueError("control_qubits must be non-negative.")

    u_gate, inferred_system_qubits = _coerce_oracle_to_gate(oracle)
    system_qubit_count = (
        inferred_system_qubits if num_system_qubits is None else int(num_system_qubits)
    )
    if system_qubit_count != inferred_system_qubits:
        raise ValueError(
            f"num_system_qubits={num_system_qubits} does not match oracle size "
            f"({inferred_system_qubits})."
        )
    if system_qubit_count <= 0:
        raise ValueError("num_system_qubits must be positive.")

    # Global layout: [controls..., ancilla, system...]
    total_qubits = control_qubits + 1 + system_qubit_count
    qc = QuantumCircuit(total_qubits, name="controlled_qsp")
    control_wires = list(range(control_qubits))
    ancilla_wire = control_qubits
    system_wires = list(range(control_qubits + 1, total_qubits))
    controlled_signal_wires = control_wires + [ancilla_wire]
    oracle_wires = [ancilla_wire] + system_wires

    # Oracle blocks (no extra controls beyond ancilla).
    cu_gate = u_gate.control(1)
    cu_dagger_gate = u_gate.inverse().control(1)

    # Controlled E0.
    c_e0 = _controlled_signal_processing_gate(
        np.asarray(result.e0, dtype=np.complex128).T,
        num_controls=control_qubits,
        label="cE0",
    )
    qc.append(c_e0, controlled_signal_wires)

    projectors = [np.asarray(projector, dtype=np.complex128) for projector in result.projectors]
    carve_sides = getattr(result, "carve_sides", None)
    side_sequence = _normalized_carve_sides(carve_sides, expected_length=len(projectors))
    use_carve_sides = side_sequence is not None

    # Recursive carving stores projectors in extraction order; application to
    # realize the original polynomial requires reversed order.
    iter_projectors = list(reversed(projectors))
    iter_sides = list(reversed(side_sequence)) if use_carve_sides else None

    # Controlled signal-processing blocks with shared oracle calls.
    for step, projector in enumerate(iter_projectors):
        v = projector_to_su2(projector)[:, [1, 0]]
        c_v = _controlled_signal_processing_gate(v, num_controls=control_qubits, label="cV")
        c_v_dagger = _controlled_signal_processing_gate(
            v.conj().T,
            num_controls=control_qubits,
            label="cVdg",
        )

        qc.append(c_v, controlled_signal_wires)

        if use_carve_sides:
            # Match base QSP semantics: top-side carve uses U, bottom uses U^\dagger.
            use_dagger = not iter_sides[step]
        else:
            # Fallback alternation when side metadata is unavailable.
            use_dagger = (step % 2) == 1

        qc.append(cu_dagger_gate if use_dagger else cu_gate, oracle_wires)

        qc.append(c_v_dagger, controlled_signal_wires)

    return qc

qsp_proc.circuits.lcu

LCU circuit-construction interfaces for QSP composition (Skelton 2025, Theorem 1).

This module declares circuit builders for linear-combination-of-unitaries (LCU) composition over multiple QSP sequences and a specialized real-part extraction construction from Appendix E.2 (Figure 10).

build_lcu_circuit(results, oracle, num_system_qubits=None)

Build the general LCU QSP circuit (Theorem 1).

For L = len(results), this constructs an LCU register of l = ceil(log2(L)) ancillas, one QSP ancilla, and a system register. The circuit prepares a uniform superposition on the LCU ancillas, applies branch-local single-qubit basis changes conditioned on |k> with interleaved/shared oracle layers across branches, pads any unused LCU basis states with explicit zero-contribution branches, then uncomputes with Hadamards. Post-selection on |0^l> yields the subnormalized block (1 / 2^l) * sum_k P_k.

Parameters:

Name Type Description Default
results list[DecompositionResult]

Decomposition outputs, one per polynomial/QSP sequence in the LCU sum.

 required
oracle QuantumCircuit | ndarray

Signal oracle as a Qiskit circuit or unitary matrix.

 required
num_system_qubits int | None

Number of system qubits. If None, infer from the oracle where possible.

 None

Returns:

Type Description
QuantumCircuit

A Qiskit circuit implementing the LCU-composed QSP construction.

Raises:

Type Description
ValueError

If inputs are empty, dimensions are invalid, or qubit counts are inconsistent.
TypeError

If oracle type is unsupported.
Source code in src\qsp_proc\circuits\lcu.py 
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
def build_lcu_circuit(
    results: list[DecompositionResult],
    oracle: QuantumCircuit | np.ndarray,
    num_system_qubits: int | None = None,
) -> QuantumCircuit:
    """Build the general LCU QSP circuit (Theorem 1).

    For ``L = len(results)``, this constructs an LCU register of
    ``l = ceil(log2(L))`` ancillas, one QSP ancilla, and a system register.
    The circuit prepares a uniform superposition on the LCU ancillas, applies
    branch-local single-qubit basis changes conditioned on ``|k>`` with
    interleaved/shared oracle layers across branches, pads any unused LCU basis
    states with explicit zero-contribution branches, then uncomputes with
    Hadamards. Post-selection on ``|0^l>`` yields the subnormalized block
    ``(1 / 2^l) * sum_k P_k``.

    Args:
        results: Decomposition outputs, one per polynomial/QSP sequence in the
            LCU sum.
        oracle: Signal oracle as a Qiskit circuit or unitary matrix.
        num_system_qubits: Number of system qubits. If ``None``, infer from the
            oracle where possible.

    Returns:
        A Qiskit circuit implementing the LCU-composed QSP construction.

    Raises:
        ValueError: If inputs are empty, dimensions are invalid, or qubit counts
            are inconsistent.
        TypeError: If ``oracle`` type is unsupported.
    """
    from qsp_proc.circuits.qsp_circuit import projector_to_su2

    if len(results) == 0:
        raise ValueError("results must contain at least one decomposition result.")

    num_system_qubits = _infer_system_qubits(
        oracle=oracle,
        num_system_qubits=num_system_qubits,
    )
    if num_system_qubits < 1:
        raise ValueError("num_system_qubits must be at least 1.")

    # Build oracle and controlled-oracle gates (controlled by QSP ancilla only).
    u_gate, u_dagger_gate = _oracle_gate_pair(
        oracle=oracle,
        num_system_qubits=num_system_qubits,
    )
    u_dagger_gate.label = "Udg"
    cu_gate = u_gate.control(1, label="CU")
    cu_dagger_gate = u_dagger_gate.control(1, label="CUdg")

    lcu_ancilla_count = int(np.ceil(np.log2(len(results))))
    total_qubits = lcu_ancilla_count + 1 + int(num_system_qubits)
    qsp_ancilla_wire = lcu_ancilla_count
    lcu_wires = list(range(lcu_ancilla_count))
    oracle_wires = [qsp_ancilla_wire] + list(range(lcu_ancilla_count + 1, total_qubits))

    qc = QuantumCircuit(total_qubits, name="LCU-QSP")

    def _append_branch_controlled_single_qubit(
        unitary: np.ndarray,
        branch: int,
        label: str,
    ) -> None:
        """Append a single-qubit unitary on QSP ancilla controlled on LCU |branch>."""
        gate = UnitaryGate(np.asarray(unitary, dtype=np.complex128), label=label)
        if lcu_ancilla_count == 0:
            qc.append(gate, [qsp_ancilla_wire])
            return
        ctrl_state = format(branch, f"0{lcu_ancilla_count}b")
        qc.append(
            gate.control(lcu_ancilla_count, ctrl_state=ctrl_state),
            lcu_wires + [qsp_ancilla_wire],
        )

    # Prepare uniform superposition over LCU ancillas.
    for wire in lcu_wires:
        qc.h(wire)

    total_branches = 1 << lcu_ancilla_count
    branch_steps: list[list[tuple[np.ndarray, bool]]] = []

    # Branch-local setup and per-layer data extraction.
    for k, result in enumerate(results):
        e0 = np.asarray(result.e0, dtype=np.complex128)
        if e0.shape != (2, 2):
            raise ValueError(f"results[{k}].e0 must have shape (2, 2), got {e0.shape}.")
        _append_branch_controlled_single_qubit(e0.T, k, label="E0")

        projectors = [np.asarray(projector, dtype=np.complex128) for projector in result.projectors]
        side_sequence = _normalized_carve_sides(
            getattr(result, "carve_sides", None),
            expected_length=len(projectors),
            context=f"results[{k}]",
        )
        use_carve_sides = side_sequence is not None

        # Recursive carving stores extraction order; apply reversed for synthesis.
        iter_projectors = list(reversed(projectors))
        iter_sides = list(reversed(side_sequence)) if use_carve_sides else None

        steps_for_branch: list[tuple[np.ndarray, bool]] = []
        for step, projector in enumerate(iter_projectors):
            v = projector_to_su2(projector)[:, [1, 0]]
            if use_carve_sides:
                use_dagger = not iter_sides[step]
            else:
                use_dagger = (step % 2) == 1
            steps_for_branch.append((v, use_dagger))
        branch_steps.append(steps_for_branch)

    # Pad non-power-of-two LCU spaces with explicit zero-contribution branches.
    # A branch-local X on the QSP ancilla makes the ancilla-|0> postselected block
    # exactly zero, avoiding a phantom identity contribution from "do nothing".
    for k in range(len(results), total_branches):
        _append_branch_controlled_single_qubit(
            np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex128),
            k,
            label="E0pad",
        )

    # Interleave branch layers so oracle calls are shared across branches.
    max_steps = max((len(steps) for steps in branch_steps), default=0)
    side_aligned = True
    for global_step in range(max_steps):
        step_sides = {
            steps[global_step][1] for steps in branch_steps if global_step < len(steps)
        }
        if len(step_sides) > 1:
            side_aligned = False
            break

    if side_aligned:
        for global_step in range(max_steps):
            active_branches = [
                branch for branch, steps in enumerate(branch_steps) if global_step < len(steps)
            ]
            if not active_branches:
                continue
            use_dagger = branch_steps[active_branches[0]][global_step][1]
            for branch in active_branches:
                v, _ = branch_steps[branch][global_step]
                _append_branch_controlled_single_qubit(v, branch, label="V")
            qc.append(cu_dagger_gate if use_dagger else cu_gate, oracle_wires)
            for branch in active_branches:
                v, _ = branch_steps[branch][global_step]
                _append_branch_controlled_single_qubit(v.conj().T, branch, label="Vdg")
    else:
        # Safety fallback for custom, mismatched side schedules.
        # Without side alignment, a fully shared oracle layer can change branch
        # dynamics. Revert to per-branch sequencing to preserve behavior.
        for branch, steps in enumerate(branch_steps):
            for v, use_dagger in steps:
                _append_branch_controlled_single_qubit(v, branch, label="V")
                qc.append(cu_dagger_gate if use_dagger else cu_gate, oracle_wires)
                _append_branch_controlled_single_qubit(v.conj().T, branch, label="Vdg")

    # Uncompute LCU superposition.
    for wire in lcu_wires:
        qc.h(wire)

    return qc

build_real_extraction_circuit(result, result_conj, oracle, num_system_qubits=None)

Build the Appendix E.2 real-part extraction circuit (P + P*) / 2.

Qubit layout is fixed as: - qubit 0: LCU ancilla, - qubit 1: QSP ancilla, - qubits 2..: system register.

The circuit implements the two-branch LCU pattern (Figure 10): H on LCU ancilla, controlled-P on |0>, controlled-P* on |1>, final H, and measurement of the LCU ancilla.

Parameters:

Name Type Description Default
result DecompositionResult

Decomposition result for P.

 required
result_conj DecompositionResult

Decomposition result for the conjugate branch P*.

 required
oracle QuantumCircuit | ndarray

Signal oracle as a Qiskit circuit or unitary matrix.

 required
num_system_qubits int | None

Number of system qubits. If None, infer from the oracle where possible.

 None

Returns:

Type Description
QuantumCircuit

A Qiskit circuit implementing real-part extraction via post-selection
QuantumCircuit

on measured LCU ancilla |0>.

Raises:

Type Description
ValueError

If num_system_qubits cannot be inferred or is invalid.
Source code in src\qsp_proc\circuits\lcu.py 
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
def build_real_extraction_circuit(
    result: DecompositionResult,
    result_conj: DecompositionResult,
    oracle: QuantumCircuit | np.ndarray,
    num_system_qubits: int | None = None,
) -> QuantumCircuit:
    """Build the Appendix E.2 real-part extraction circuit ``(P + P*) / 2``.

    Qubit layout is fixed as:
    - qubit ``0``: LCU ancilla,
    - qubit ``1``: QSP ancilla,
    - qubits ``2..``: system register.

    The circuit implements the two-branch LCU pattern (Figure 10):
    ``H`` on LCU ancilla, controlled-``P`` on ``|0>``, controlled-``P*`` on
    ``|1>``, final ``H``, and measurement of the LCU ancilla.

    Args:
        result: Decomposition result for ``P``.
        result_conj: Decomposition result for the conjugate branch ``P*``.
        oracle: Signal oracle as a Qiskit circuit or unitary matrix.
        num_system_qubits: Number of system qubits. If ``None``, infer from the
            oracle where possible.

    Returns:
        A Qiskit circuit implementing real-part extraction via post-selection
        on measured LCU ancilla ``|0>``.

    Raises:
        ValueError: If ``num_system_qubits`` cannot be inferred or is invalid.
    """
    from qsp_proc.circuits.controlled import build_controlled_qsp_circuit

    num_system_qubits = _infer_system_qubits(
        oracle=oracle,
        num_system_qubits=num_system_qubits,
    )

    if num_system_qubits < 1:
        raise ValueError("num_system_qubits must be at least 1.")

    total_qubits = 2 + num_system_qubits
    full_wires = list(range(total_qubits))

    # Controlled QSP branches (control qubit is wire 0 in each subcircuit).
    controlled_p = build_controlled_qsp_circuit(
        result=result,
        oracle=oracle,
        control_qubits=1,
        num_system_qubits=num_system_qubits,
    )
    controlled_p_conj = build_controlled_qsp_circuit(
        result=result_conj,
        oracle=oracle,
        control_qubits=1,
        num_system_qubits=num_system_qubits,
    )

    circuit = QuantumCircuit(total_qubits, 1, name="real_extraction_lcu")

    # Prepare equal LCU superposition.
    circuit.h(0)

    # Apply controlled-P on |0> branch (flip control basis with X).
    circuit.x(0)
    circuit.compose(controlled_p, qubits=full_wires, inplace=True)
    circuit.x(0)

    # Apply controlled-P* on |1> branch.
    circuit.compose(controlled_p_conj, qubits=full_wires, inplace=True)

    # Interfere branches and measure LCU ancilla for post-selection.
    circuit.h(0)
    circuit.measure(0, 0)

    return circuit