Repository: PennyLaneAI/pennylane
Pull request #6042 (2024-07-26 14:26:59)

Add support for multi-controlled zyz
### Before submitting

Please complete the following checklist when submitting a PR:

- [x] All new features must include a unit test.
      If you've fixed a bug or added code that should be tested, add a test to the
      test directory!

- [x] All new functions and code must be clearly commented and documented.
      If you do make documentation changes, make sure that the docs build and
      render correctly by running `make docs`.

- [x] Ensure that the test suite passes, by running `make test`.

- [x] Add a new entry to the `doc/releases/changelog-dev.md` file, summarizing the
      change, and including a link back to the PR.

- [x] The PennyLane source code conforms to
      [PEP8 standards](https://www.python.org/dev/peps/pep-0008/).
      We check all of our code against [Pylint](https://www.pylint.org/).
      To lint modified files, simply `pip install pylint`, and then
      run `pylint pennylane/path/to/file.py`.

When all the above are checked, delete everything above the dashed
line and fill in the pull request template.

------------------------------------------------------------------------------------------------------------

**Context:**
This PR implements the decomposition of zyz for special unitaries with multiple control wires defined in Lemma 7.9 of https://arxiv.org/pdf/quant-ph/9503016.

**Benefits:**
Support decomposition of zyz for special unitaries.

**Possible Drawbacks:**

**Related GitHub Issues:**
[sc-68471]

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6042?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
Attention: Patch coverage is `97.91667%` with `1 line` in your changes missing coverage. Please review.
> Project coverage is 99.57%. Comparing base [(`b7675ca`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/b7675ca42f2e237d396788b287de256254a71958?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`27432b7`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/27432b7ddae51145ea93ae789afb67cd62acfd57?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> Report is 1 commits behind head on master.

| [Files](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6042?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Patch % | Lines |
|---|---|---|
| [pennylane/ops/op\_math/controlled\_decompositions.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6042?src=pr&el=tree&filepath=pennylane%2Fops%2Fop_math%2Fcontrolled_decompositions.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9vcF9tYXRoL2NvbnRyb2xsZWRfZGVjb21wb3NpdGlvbnMucHk=) | 97.77% | [1 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6042?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #6042      +/-   ##
==========================================
- Coverage   99.66%   99.57%   -0.09%     
==========================================
  Files         428      428              
  Lines       41182    40992     -190     
==========================================
- Hits        41042    40819     -223     
- Misses        140      173      +33     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/6042?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


Patch:
diff --git a/Makefile b/Makefile
index 5f7186e45..7f36c24c6 100644
--- a/Makefile
+++ b/Makefile
@@ -70,17 +70,17 @@ coverage:
 .PHONY:format
 format:
 ifdef check
-	isort --py 311 --profile black -l 100 -o autoray -p ./pennylane --skip __init__.py --filter-files ./pennylane ./tests --check
-	black -t py39 -t py310 -t py311 -l 100 ./pennylane ./tests --check
+	$(PYTHON) -m isort --py 311 --profile black -l 100 -o autoray -p ./pennylane --skip __init__.py --filter-files ./pennylane ./tests --check
+	$(PYTHON) -m black -t py39 -t py310 -t py311 -l 100 ./pennylane ./tests --check
 else
-	isort --py 311 --profile black -l 100 -o autoray -p ./pennylane --skip __init__.py --filter-files ./pennylane ./tests
-	black -t py39 -t py310 -t py311 -l 100 ./pennylane ./tests
+	$(PYTHON) -m isort --py 311 --profile black -l 100 -o autoray -p ./pennylane --skip __init__.py --filter-files ./pennylane ./tests
+	$(PYTHON) -m black -t py39 -t py310 -t py311 -l 100 ./pennylane ./tests
 endif
 
 .PHONY: lint
 lint:
-	pylint pennylane --rcfile .pylintrc
+	$(PYTHON) -m pylint pennylane --rcfile .pylintrc
 
 .PHONY: lint-test
 lint-test:
-	pylint tests pennylane/devices/tests --rcfile tests/.pylintrc
+	$(PYTHON) -m pylint tests pennylane/devices/tests --rcfile tests/.pylintrc
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 1e313abe1..2084b0599 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -216,6 +216,9 @@
 
 <h4>Other improvements</h4>
 
+* Added the decomposition of zyz for special unitaries with multiple control wires.
+  [(#6042)](https://github.com/PennyLaneAI/pennylane/pull/6042)
+
 * A new method `process_density_matrix` has been added to the `ProbabilityMP` and `DensityMatrixMP`
   classes, allowing for more efficient handling of quantum density matrices, particularly with batch
   processing support. This method simplifies the calculation of probabilities from quantum states
@@ -395,6 +398,7 @@ This release contains contributions from (in alphabetical order):
 
 Tarun Kumar Allamsetty,
 Guillermo Alonso,
+Ali Asadi,
 Utkarsh Azad,
 Tonmoy T. Bhattacharya,
 Gabriel Bottrill,
diff --git a/pennylane/ops/op_math/controlled.py b/pennylane/ops/op_math/controlled.py
index df95dcd83..4bc0690a9 100644
--- a/pennylane/ops/op_math/controlled.py
+++ b/pennylane/ops/op_math/controlled.py
@@ -864,7 +864,7 @@ def _decompose_custom_ops(op: Controlled) -> list["operation.Operator"]:
     return None
 
 
-def _decompose_no_control_values(op: Controlled) -> list["operation.Operator"]:
+def _decompose_no_control_values(op: Controlled) -> Optional[list["operation.Operator"]]:
     """Decompose without considering control values. Returns None if no decomposition."""
 
     decomp = _decompose_custom_ops(op)
@@ -874,7 +874,7 @@ def _decompose_no_control_values(op: Controlled) -> list["operation.Operator"]:
     if _is_single_qubit_special_unitary(op.base):
         if len(op.control_wires) >= 2 and qmlmath.get_interface(*op.data) == "numpy":
             return ctrl_decomp_bisect(op.base, op.control_wires)
-        return ctrl_decomp_zyz(op.base, op.control_wires)
+        return ctrl_decomp_zyz(op.base, op.control_wires, work_wires=op.work_wires)
 
     if not op.base.has_decomposition:
         return None
diff --git a/pennylane/ops/op_math/controlled_decompositions.py b/pennylane/ops/op_math/controlled_decompositions.py
index 17b8405c0..26de3761c 100644
--- a/pennylane/ops/op_math/controlled_decompositions.py
+++ b/pennylane/ops/op_math/controlled_decompositions.py
@@ -16,6 +16,7 @@ This submodule defines functions to decompose controlled operations
 """
 
 from copy import copy
+from typing import Optional
 
 import numpy as np
 import numpy.linalg as npl
@@ -134,12 +135,95 @@ def _bisect_compute_b(u: np.ndarray):
     return _param_su2(c, d, b, 0)
 
 
-def ctrl_decomp_zyz(target_operation: Operator, control_wires: Wires):
+def _multi_controlled_zyz(
+    rot_angles,
+    global_phase,
+    target_wire: Wires,
+    control_wires: Wires,
+    work_wires: Optional[Wires] = None,
+) -> list[Operator]:
+    # The decomposition of zyz for special unitaries with multiple control wires
+    # defined in Lemma 7.9 of https://arxiv.org/pdf/quant-ph/9503016
+
+    if not qml.math.allclose(0.0, global_phase, atol=1e-6, rtol=0):
+        raise ValueError(f"The global_phase should be zero, instead got: {global_phase}.")
+
+    # Unpack the rotation angles
+    phi, theta, omega = rot_angles
+
+    # We use the conditional statements to account when decomposition is ran within a queue
+    decomp = []
+
+    cop_wires = (control_wires[-1], target_wire[0])
+
+    # Add A operator
+    if not qml.math.allclose(0.0, phi, atol=1e-8, rtol=0):
+        decomp.append(qml.CRZ(phi, wires=cop_wires))
+    if not qml.math.allclose(0.0, theta / 2, atol=1e-8, rtol=0):
+        decomp.append(qml.CRY(theta / 2, wires=cop_wires))
+
+    decomp.append(qml.ctrl(qml.X(target_wire), control=control_wires[:-1], work_wires=work_wires))
+
+    # Add B operator
+    if not qml.math.allclose(0.0, theta / 2, atol=1e-8, rtol=0):
+        decomp.append(qml.CRY(-theta / 2, wires=cop_wires))
+    if not qml.math.allclose(0.0, -(phi + omega) / 2, atol=1e-6, rtol=0):
+        decomp.append(qml.CRZ(-(phi + omega) / 2, wires=cop_wires))
+
+    decomp.append(qml.ctrl(qml.X(target_wire), control=control_wires[:-1], work_wires=work_wires))
+
+    # Add C operator
+    if not qml.math.allclose(0.0, (omega - phi) / 2, atol=1e-8, rtol=0):
+        decomp.append(qml.CRZ((omega - phi) / 2, wires=cop_wires))
+
+    return decomp
+
+
+def _single_control_zyz(rot_angles, global_phase, target_wire, control_wires: Wires):
+    # The zyz decomposition of a general unitary with single control wire
+    # defined in Lemma 7.9 of https://arxiv.org/pdf/quant-ph/9503016
+
+    # Unpack the rotation angles
+    phi, theta, omega = rot_angles
+    # We use the conditional statements to account when decomposition is ran within a queue
+    decomp = []
+    # Add negative of global phase. Compare definition of qml.GlobalPhase and Ph(delta) from section 4.1 of Barenco et al.
+    if not qml.math.allclose(0.0, global_phase, atol=1e-8, rtol=0):
+        decomp.append(
+            qml.ctrl(qml.GlobalPhase(phi=-global_phase, wires=target_wire), control=control_wires)
+        )
+    # Add A operator
+    if not qml.math.allclose(0.0, phi, atol=1e-8, rtol=0):
+        decomp.append(qml.RZ(phi, wires=target_wire))
+    if not qml.math.allclose(0.0, theta / 2, atol=1e-8, rtol=0):
+        decomp.append(qml.RY(theta / 2, wires=target_wire))
+
+    decomp.append(qml.ctrl(qml.X(target_wire), control=control_wires))
+
+    # Add B operator
+    if not qml.math.allclose(0.0, theta / 2, atol=1e-8, rtol=0):
+        decomp.append(qml.RY(-theta / 2, wires=target_wire))
+    if not qml.math.allclose(0.0, -(phi + omega) / 2, atol=1e-6, rtol=0):
+        decomp.append(qml.RZ(-(phi + omega) / 2, wires=target_wire))
+
+    decomp.append(qml.ctrl(qml.X(target_wire), control=control_wires))
+
+    # Add C operator
+    if not qml.math.allclose(0.0, (omega - phi) / 2, atol=1e-8, rtol=0):
+        decomp.append(qml.RZ((omega - phi) / 2, wires=target_wire))
+
+    return decomp
+
+
+def ctrl_decomp_zyz(
+    target_operation: Operator, control_wires: Wires, work_wires: Optional[Wires] = None
+) -> list[Operator]:
     """Decompose the controlled version of a target single-qubit operation
 
-    This function decomposes a controlled single-qubit target operation with one
-    single control using the decomposition defined in Lemma 4.3 and Lemma 5.1 of
-    `Barenco et al. (1995) <https://arxiv.org/abs/quant-ph/9503016>`_.
+    This function decomposes both single and multiple controlled single-qubit
+    target operations using the decomposition defined in Lemma 4.3 and Lemma 5.1
+    for single `controlled_wires`, and Lemma 7.9 for multiple `controlled_wires`
+    from `Barenco et al. (1995) <https://arxiv.org/abs/quant-ph/9503016>`_.
 
     Args:
         target_operation (~.operation.Operator): the target operation to decompose
@@ -190,53 +274,24 @@ def ctrl_decomp_zyz(target_operation: Operator, control_wires: Wires):
             f"got {target_operation.__class__.__name__}."
         )
     control_wires = Wires(control_wires)
-    if len(control_wires) > 1:
-        raise ValueError(
-            f"The control_wires should be a single wire, instead got: {len(control_wires)} wires."
-        )
 
     target_wire = target_operation.wires
 
     if isinstance(target_operation, Operation):
         try:
-            phi, theta, omega = target_operation.single_qubit_rot_angles()
+            rot_angles = target_operation.single_qubit_rot_angles()
         except NotImplementedError:
-            phi, theta, omega = _get_single_qubit_rot_angles_via_matrix(
-                qml.matrix(target_operation)
-            )
+            rot_angles = _get_single_qubit_rot_angles_via_matrix(qml.matrix(target_operation))
     else:
-        phi, theta, omega = _get_single_qubit_rot_angles_via_matrix(qml.matrix(target_operation))
+        rot_angles = _get_single_qubit_rot_angles_via_matrix(qml.matrix(target_operation))
 
     _, global_phase = _convert_to_su2(qml.matrix(target_operation), return_global_phase=True)
 
-    # We use the conditional statements to account when decomposition is ran within a queue
-    decomp = []
-    # Add negative of global phase. Compare definition of qml.GlobalPhase and Ph(delta) from section 4.1 of Barenco et al.
-    if not qml.math.allclose(0.0, global_phase, atol=1e-8, rtol=0):
-        decomp.append(
-            qml.ctrl(qml.GlobalPhase(phi=-global_phase, wires=target_wire), control=control_wires)
-        )
-    # Add A operator
-    if not qml.math.allclose(0.0, phi, atol=1e-8, rtol=0):
-        decomp.append(qml.RZ(phi, wires=target_wire))
-    if not qml.math.allclose(0.0, theta / 2, atol=1e-8, rtol=0):
-        decomp.append(qml.RY(theta / 2, wires=target_wire))
-
-    decomp.append(qml.ctrl(qml.X(target_wire), control=control_wires))
-
-    # Add B operator
-    if not qml.math.allclose(0.0, theta / 2, atol=1e-8, rtol=0):
-        decomp.append(qml.RY(-theta / 2, wires=target_wire))
-    if not qml.math.allclose(0.0, -(phi + omega) / 2, atol=1e-6, rtol=0):
-        decomp.append(qml.RZ(-(phi + omega) / 2, wires=target_wire))
-
-    decomp.append(qml.ctrl(qml.PauliX(wires=target_wire), control=control_wires))
-
-    # Add C operator
-    if not qml.math.allclose(0.0, (omega - phi) / 2, atol=1e-8, rtol=0):
-        decomp.append(qml.RZ((omega - phi) / 2, wires=target_wire))
-
-    return decomp
+    return (
+        _multi_controlled_zyz(rot_angles, global_phase, target_wire, control_wires, work_wires)
+        if len(control_wires) > 1
+        else _single_control_zyz(rot_angles, global_phase, target_wire, control_wires)
+    )
 
 
 def _ctrl_decomp_bisect_od(


Test patch:
diff --git a/tests/ops/op_math/test_controlled.py b/tests/ops/op_math/test_controlled.py
index d73616a29..40e315449 100644
--- a/tests/ops/op_math/test_controlled.py
+++ b/tests/ops/op_math/test_controlled.py
@@ -1050,20 +1050,39 @@ class TestDecomposition:
         decomp_mat = qml.matrix(op.decomposition, wire_order=op.wires)()
         assert qml.math.allclose(op.matrix(), decomp_mat)
 
-    def test_differentiable_one_qubit_special_unitary(self):
-        """Assert that a differentiable qubit special unitary uses the zyz decomposition."""
+    def test_differentiable_one_qubit_special_unitary_single_ctrl(self):
+        """
+        Assert that a differentiable qubit special unitary uses the zyz decomposition with a single controlled wire.
+        """
 
-        op = qml.ctrl(qml.RZ(qml.numpy.array(1.2), 0), (1))
+        theta = 1.2
+        op = qml.ctrl(qml.RZ(qml.numpy.array(theta), 0), (1))
         decomp = op.decomposition()
 
-        qml.assert_equal(decomp[0], qml.PhaseShift(qml.numpy.array(1.2 / 2), 0))
+        qml.assert_equal(decomp[0], qml.PhaseShift(qml.numpy.array(theta / 2), 0))
         qml.assert_equal(decomp[1], qml.CNOT(wires=(1, 0)))
-        qml.assert_equal(decomp[2], qml.PhaseShift(qml.numpy.array(-1.2 / 2), 0))
+        qml.assert_equal(decomp[2], qml.PhaseShift(qml.numpy.array(-theta / 2), 0))
         qml.assert_equal(decomp[3], qml.CNOT(wires=(1, 0)))
 
         decomp_mat = qml.matrix(op.decomposition, wire_order=op.wires)()
         assert qml.math.allclose(op.matrix(), decomp_mat)
 
+    def test_differentiable_one_qubit_special_unitary_multiple_ctrl(self):
+        """Assert that a differentiable qubit special unitary uses the zyz decomposition with multiple controlled wires."""
+
+        theta = 1.2
+        op = qml.ctrl(qml.RZ(qml.numpy.array(theta), 0), (1, 2, 3, 4))
+        decomp = op.decomposition()
+
+        assert qml.equal(decomp[0], qml.CRZ(qml.numpy.array(theta), [4, 0]))
+        assert qml.equal(decomp[1], qml.MultiControlledX(wires=[1, 2, 3, 0]))
+        assert qml.equal(decomp[2], qml.CRZ(qml.numpy.array(-theta / 2), wires=[4, 0]))
+        assert qml.equal(decomp[3], qml.MultiControlledX(wires=[1, 2, 3, 0]))
+        assert qml.equal(decomp[4], qml.CRZ(qml.numpy.array(-theta / 2), wires=[4, 0]))
+
+        decomp_mat = qml.matrix(op.decomposition, wire_order=op.wires)()
+        assert qml.math.allclose(op.matrix(), decomp_mat)
+
     @pytest.mark.parametrize(
         "base_cls, params, base_wires, ctrl_wires, custom_ctrl_cls, expected",
         custom_ctrl_op_decomps,
@@ -1730,7 +1749,6 @@ class TestCtrl:
 
         if isinstance(op, qml.QubitUnitary):
             pytest.skip("ControlledQubitUnitary can accept any number of control wires.")
-            expected = None  # to pass pylint(possibly-used-before-assignment) error
         elif isinstance(op, Controlled):
             expected = Controlled(
                 op.base,
diff --git a/tests/ops/op_math/test_controlled_decompositions.py b/tests/ops/op_math/test_controlled_decompositions.py
index c8cd1e4e3..6c1adda9f 100644
--- a/tests/ops/op_math/test_controlled_decompositions.py
+++ b/tests/ops/op_math/test_controlled_decompositions.py
@@ -87,14 +87,6 @@ class TestControlledDecompositionZYZ:
         ):
             _ = ctrl_decomp_zyz(qml.CNOT([0, 1]), [2])
 
-    def test_invalid_num_controls(self):
-        """Tests that an error is raised when an invalid number of control wires is passed"""
-        with pytest.raises(
-            ValueError,
-            match="The control_wires should be a single wire, instead got: 2",
-        ):
-            _ = ctrl_decomp_zyz(qml.X([1]), [0, 1])
-
     su2_ops = [
         qml.RX(0.123, wires=0),
         qml.RY(0.123, wires=0),
@@ -102,11 +94,7 @@ class TestControlledDecompositionZYZ:
         qml.Rot(0.123, 0.456, 0.789, wires=0),
     ]
 
-    unitary_ops = [
-        qml.Hadamard(0),
-        qml.PauliZ(0),
-        qml.S(0),
-        qml.PhaseShift(1.5, wires=0),
+    general_unitary_ops = [
         qml.QubitUnitary(
             np.array(
                 [
@@ -117,13 +105,17 @@ class TestControlledDecompositionZYZ:
             wires=0,
         ),
         qml.DiagonalQubitUnitary(np.array([1, -1]), wires=0),
+        qml.Hadamard(0),
+        qml.PauliZ(0),
+        qml.S(0),
+        qml.PhaseShift(1.5, wires=0),
     ]
 
-    @pytest.mark.parametrize("op", su2_ops + unitary_ops)
+    @pytest.mark.parametrize("op", su2_ops + general_unitary_ops)
     @pytest.mark.parametrize("control_wires", ([1], [2], [3]))
-    def test_decomposition_circuit(self, op, control_wires, tol):
+    def test_decomposition_circuit_general_ops(self, op, control_wires, tol):
         """Tests that the controlled decomposition of a single-qubit operation
-        behaves as expected in a quantum circuit"""
+        behaves as expected in a quantum circuit for general_unitary_ops"""
         dev = qml.device("default.qubit", wires=4)
 
         @qml.qnode(dev)
@@ -143,7 +135,26 @@ class TestControlledDecompositionZYZ:
 
         assert np.allclose(res, expected, atol=tol, rtol=0)
 
-    @pytest.mark.parametrize("control_wires", ([1], [2], [3]))
+    @pytest.mark.parametrize("op", general_unitary_ops)
+    @pytest.mark.parametrize("control_wires", ([1, 2], [1, 2, 3]))
+    def test_decomposition_circuit_general_ops_error(self, op, control_wires):
+        """Tests that the controlled decomposition of a single-qubit operation
+        with multiple controlled wires raises a ValueError for general_unitary_ops"""
+        dev = qml.device("default.qubit", wires=4)
+
+        @qml.qnode(dev)
+        def decomp_circuit():
+            qml.broadcast(unitary=qml.Hadamard, pattern="single", wires=control_wires)
+            ctrl_decomp_zyz(op, Wires(control_wires))
+            return qml.probs()
+
+        with pytest.raises(
+            ValueError,
+            match="The global_phase should be zero",
+        ):
+            decomp_circuit()
+
+    @pytest.mark.parametrize("control_wires", ([1], [1, 2], [1, 2, 3]))
     def test_decomposition_circuit_gradient(self, control_wires, tol):
         """Tests that the controlled decomposition of a single-qubit operation
         behaves as expected in a quantum circuit"""
@@ -193,23 +204,23 @@ class TestControlledDecompositionZYZ:
         """Test that the operations in the decomposition are correct."""
         phi, theta, omega = 0.123, 0.456, 0.789
         op = qml.Rot(phi, theta, omega, wires=0)
-        control_wires = [1]
+        control_wires = [1, 2, 3]
         decomps = ctrl_decomp_zyz(op, Wires(control_wires))
 
         expected_ops = [
-            qml.RZ(0.123, wires=0),
-            qml.RY(0.456 / 2, wires=0),
-            qml.CNOT(wires=control_wires + [0]),
-            qml.RY(-0.456 / 2, wires=0),
-            qml.RZ(-(0.123 + 0.789) / 2, wires=0),
-            qml.CNOT(wires=control_wires + [0]),
-            qml.RZ((0.789 - 0.123) / 2, wires=0),
+            qml.CRZ(0.123, wires=[3, 0]),
+            qml.CRY(0.456 / 2, wires=[3, 0]),
+            qml.Toffoli(wires=control_wires[:-1] + [0]),
+            qml.CRY(-0.456 / 2, wires=[3, 0]),
+            qml.CRZ(-(0.123 + 0.789) / 2, wires=[3, 0]),
+            qml.Toffoli(wires=control_wires[:-1] + [0]),
+            qml.CRZ((0.789 - 0.123) / 2, wires=[3, 0]),
         ]
         for decomp_op, expected_op in zip(decomps, expected_ops):
             qml.assert_equal(decomp_op, expected_op)
         assert len(decomps) == 7
 
-    @pytest.mark.parametrize("op", su2_ops + unitary_ops)
+    @pytest.mark.parametrize("op", su2_ops + general_unitary_ops)
     @pytest.mark.parametrize("control_wires", ([1], [2], [3]))
     def test_decomp_queues_correctly(self, op, control_wires, tol):
         """Test that any incorrect operations aren't queued when using
@@ -896,6 +907,30 @@ class TestMultiControlledUnitary:
         res = _decompose_multicontrolled_unitary(op, Wires(control_wires))
         assert equal_list(res, expected)
 
+    @pytest.mark.parametrize(
+        "op, controlled_wires, work_wires",
+        [
+            (qml.RX(0.123, wires=1), [0, 2], [3, 4, 5]),
+            (qml.Rot(0.123, 0.456, 0.789, wires=0), [1, 2, 3], [4, 5]),
+        ],
+    )
+    def test_with_many_workers(self, op, controlled_wires, work_wires):
+        """Tests ctrl_decomp_zyz with multiple workers"""
+
+        dev = qml.device("default.qubit", wires=6)
+
+        @qml.qnode(dev)
+        def decomp_circuit(op):
+            ctrl_decomp_zyz(op, controlled_wires, work_wires=work_wires)
+            return qml.probs()
+
+        @qml.qnode(dev)
+        def expected_circuit(op):
+            qml.ctrl(op, controlled_wires, work_wires=work_wires)
+            return qml.probs()
+
+        assert np.allclose(decomp_circuit(op), expected_circuit(op))
+
     controlled_wires = tuple(list(range(2, 1 + n)) for n in range(3, 7))
 
     @pytest.mark.parametrize("op", gen_ops + su2_gen_ops)