Repository: PennyLaneAI/pennylane
Pull request #5125 (2024-01-30 15:24:05)

Controlled operations rework Part 1
**Context:**
All controlled operations should inherit from the general Controlled class, and the decomposition of controlled operations is not consistent for custom and non-custom controlled operations. This is a continuation of https://github.com/PennyLaneAI/pennylane/pull/5069

This is the first PR out of two for this rework. The second PR will focus on making sure that all custom controlled operations inherit from Controlled for more consistent inheritance structure.

**Description of the Change:**
- Make `MultiControlledX` inherit from ControlledOp.
- `qml.ctrl` called on operators with custom controlled versions will return instances of the custom class.
- Special handling of `PauliX` based controlled operations (`PauliX`, `CNOT`, `Toffoli`, `MultiControlledX`)
  - Calling `qml.ctrl` on one of these operators will always resolve to the best option in `CNOT`, `Toffoli`, or `MultiControlledX` depending on the number of control wires and control values.
- `qml.ctrl` will flatten nested controlled operators to a single multi-controlled operation.
- Controlled operators with a custom controlled version decomposes like how their controlled counterpart decomposes, as opposed to decomposing into their controlled version.
  - Special handling of `PauliX` based controlled operations: e.g., `Controlled(CNOT([0, 1]), [2, 3])` will have the same decomposition behaviour as a `MultiControlledX([2, 3, 0, 1])`

**Benefits:**
Cleaner code and more consistent behaviour

**Possible Drawbacks:**
Change of decomposition behaviour may cause issues.
~For `MultiControlledX`, the `wires` attribute now refers to all wires, as in `control_wires + target_wire + work_wires`, to access only the `control_wires + target_wires`, use the `active_wires` attribute.~

**Related GitHub Issues:**
https://github.com/PennyLaneAI/pennylane/pull/5069
https://github.com/PennyLaneAI/pennylane/issues/1447

**Related Shortcut Stories**
[sc-55949]
[sc-55131]
[sc-55358]

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/5125?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
Attention: `3 lines` in your changes are missing coverage. Please review.
> Comparison is base [(`cfe2991`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/cfe29911a30c4515f9500e640867305fdea46433?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.69% compared to head [(`964e008`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5125?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.67%.
> Report is 1 commits behind head on master.

| [Files](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5125?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\_ops.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5125?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9vcF9tYXRoL2NvbnRyb2xsZWRfb3BzLnB5) | 96.10% | [3 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5125?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    #5125      +/-   ##
==========================================
- Coverage   99.69%   99.67%   -0.02%     
==========================================
  Files         393      393              
  Lines       35864    35629     -235     
==========================================
- Hits        35753    35514     -239     
- Misses        111      115       +4     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5125?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/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 16ea70792..9a33628a9 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -165,9 +165,17 @@
   and the codecov check itself would never execute.
   [(#5101)](https://github.com/PennyLaneAI/pennylane/pull/5101)
 
+* `qml.ctrl` called on operators with custom controlled versions will return instances
+  of the custom class, and it will also flatten nested controlled operators to a single
+  multi-controlled operation. For `PauliX`, `CNOT`, `Toffoli`, and `MultiControlledX`,
+  calling `qml.ctrl` will always resolve to the best option in `CNOT`, `Toffoli`, or
+  `MultiControlledX` depending on the number of control wires and control values.
+  [(#5125)](https://github.com/PennyLaneAI/pennylane/pull/5125/)
+
 * `qml.Identity()` can be initialized without wires. Measuring it is currently not possible though.
   [(#5106)](https://github.com/PennyLaneAI/pennylane/pull/5106)
 
+
 <h4>Community contributions 🥳</h4>
 
 * The transform `split_non_commuting` now accepts measurements of type `probs`, `sample` and `counts` which accept both wires and observables.
@@ -203,7 +211,12 @@
   (with potentially negative eigenvalues) has been implemented.
   [(#5048)](https://github.com/PennyLaneAI/pennylane/pull/5048)
 
-* The decomposition of an operator created with calling `qml.ctrl` on a parametric operator (specifically `RX`, `RY`, `RZ`, `Rot`, `PhaseShift`) with a single control wire will now be the full decomposition instead of a single controlled gate. For example:
+* Controlled operators with a custom controlled version decomposes like how their
+  controlled counterpart decomposes, as opposed to decomposing into their controlled version.   
+  [(#5069)](https://github.com/PennyLaneAI/pennylane/pull/5069)
+  [(#5125)](https://github.com/PennyLaneAI/pennylane/pull/5125/)
+  
+  For example:
   ```
   >>> qml.ctrl(qml.RX(0.123, wires=1), control=0).decomposition()
   [
@@ -215,7 +228,6 @@
     RZ(-1.5707963267948966, wires=[1])
   ]
   ```
-  [(#5069)](https://github.com/PennyLaneAI/pennylane/pull/5069)
 
 * `QuantumScript.is_sampled` and `QuantumScript.all_sampled` have been removed. Users should now
   validate these properties manually.
diff --git a/pennylane/devices/qubit/apply_operation.py b/pennylane/devices/qubit/apply_operation.py
index cd6e4949e..8f7c9a7a1 100644
--- a/pennylane/devices/qubit/apply_operation.py
+++ b/pennylane/devices/qubit/apply_operation.py
@@ -274,9 +274,7 @@ def apply_multicontrolledx(
         return _apply_operation_default(op, state, is_state_batched, debugger)
     ctrl_wires = [w + is_state_batched for w in op.control_wires]
     # apply x on all control wires with control value 0
-    roll_axes = [
-        w for val, w in zip(op.hyperparameters["control_values"], ctrl_wires) if val == "0"
-    ]
+    roll_axes = [w for val, w in zip(op.control_values, ctrl_wires) if val is False]
     for ax in roll_axes:
         state = math.roll(state, 1, ax)
 
diff --git a/pennylane/drawer/tape_mpl.py b/pennylane/drawer/tape_mpl.py
index ce8200a4e..97464edad 100644
--- a/pennylane/drawer/tape_mpl.py
+++ b/pennylane/drawer/tape_mpl.py
@@ -107,8 +107,7 @@ def _(op: ops.Toffoli, drawer, layer, _):
 
 @_add_operation_to_drawer.register
 def _(op: ops.MultiControlledX, drawer, layer, _):
-    control_values = [(i == "1") for i in op.hyperparameters["control_values"]]
-    drawer.CNOT(layer, op.wires, control_values=control_values)
+    drawer.CNOT(layer, op.active_wires, control_values=op.control_values)
 
 
 @_add_operation_to_drawer.register
diff --git a/pennylane/ops/functions/bind_new_parameters.py b/pennylane/ops/functions/bind_new_parameters.py
index 96b93b29c..dd8bff6f5 100644
--- a/pennylane/ops/functions/bind_new_parameters.py
+++ b/pennylane/ops/functions/bind_new_parameters.py
@@ -113,8 +113,9 @@ def bind_new_parameters_composite_op(op: CompositeOp, params: Sequence[TensorLik
 
 @bind_new_parameters.register(qml.CY)
 @bind_new_parameters.register(qml.CZ)
+@bind_new_parameters.register(qml.MultiControlledX)
 def bind_new_parameters_copy(
-    op: Union[qml.CY, qml.CZ], params: Sequence[TensorLike]
+    op: Union[qml.CY, qml.CZ, qml.MultiControlledX], params: Sequence[TensorLike]
 ):  # pylint:disable=unused-argument
     return copy.copy(op)
 
diff --git a/pennylane/ops/op_math/__init__.py b/pennylane/ops/op_math/__init__.py
index e15ab1cf6..20089d30d 100644
--- a/pennylane/ops/op_math/__init__.py
+++ b/pennylane/ops/op_math/__init__.py
@@ -109,6 +109,7 @@ from .controlled_ops import (
     CRZ,
     CY,
     CZ,
+    MultiControlledX,
 )
 from .decompositions import one_qubit_decomposition, two_qubit_decomposition, sk_decomposition
 from .evolution import Evolution
@@ -124,6 +125,7 @@ controlled_qubit_ops = {
     "ControlledQubitUnitary",
     "CY",
     "CZ",
+    "MultiControlledX",
     "CRX",
     "CRY",
     "CRZ",
diff --git a/pennylane/ops/op_math/controlled.py b/pennylane/ops/op_math/controlled.py
index 37bd06f96..b0ba87b20 100644
--- a/pennylane/ops/op_math/controlled.py
+++ b/pennylane/ops/op_math/controlled.py
@@ -15,6 +15,7 @@
 This submodule defines the symbolic operation that indicates the control of an operator.
 """
 import warnings
+import functools
 from copy import copy
 from functools import wraps
 from inspect import signature
@@ -127,36 +128,41 @@ def ctrl(op, control, control_values=None, work_wires=None):
         ops_loader = available_eps[active_jit]["ops"].load()
         return ops_loader.ctrl(op, control, control_values=control_values, work_wires=work_wires)
 
-    custom_ops = {
-        (qml.PauliZ, 1): qml.CZ,
-        (qml.PauliY, 1): qml.CY,
-        (qml.PauliX, 1): qml.CNOT,
-        (qml.PauliX, 2): qml.Toffoli,
-        (qml.RX, 1): qml.CRX,
-        (qml.RY, 1): qml.CRY,
-        (qml.RZ, 1): qml.CRZ,
-        (qml.Rot, 1): qml.CRot,
-        (qml.PhaseShift, 1): qml.ControlledPhaseShift,
-    }
-    control_values = [control_values] if isinstance(control_values, (int, bool)) else control_values
     control = qml.wires.Wires(control)
-    custom_key = (type(op), len(control))
+    if isinstance(control_values, (int, bool)):
+        control_values = [control_values]
+    elif control_values is None:
+        control_values = [True] * len(control)
 
-    if custom_key in custom_ops and (control_values is None or all(control_values)):
-        qml.QueuingManager.remove(op)
-        return custom_ops[custom_key](*op.data, control + op.wires)
-    if isinstance(op, qml.PauliX):
+    ctrl_op = _try_wrap_in_custom_ctrl_op(
+        op, control, control_values=control_values, work_wires=work_wires
+    )
+    if ctrl_op is not None:
+        return ctrl_op
+
+    pauli_x_based_ctrl_ops = _get_pauli_x_based_ops()
+
+    # Special handling for PauliX-based controlled operations
+    if isinstance(op, pauli_x_based_ctrl_ops):
         qml.QueuingManager.remove(op)
-        control_string = (
-            None if control_values is None else "".join([str(int(v)) for v in control_values])
-        )
-        return qml.MultiControlledX(
-            wires=control + op.wires, control_values=control_string, work_wires=work_wires
+        return _handle_pauli_x_based_controlled_ops(op, control, control_values, work_wires)
+
+    # Flatten nested controlled operations to a multi-controlled operation for better
+    # decomposition algorithms. This includes special cases like CRX, CRot, etc.
+    if isinstance(op, Controlled):
+        work_wires = work_wires or []
+        return ctrl(
+            op.base,
+            control=control + op.control_wires,
+            control_values=control_values + op.control_values,
+            work_wires=work_wires + op.work_wires,
         )
+
     if isinstance(op, Operator):
         return Controlled(
             op, control_wires=control, control_values=control_values, work_wires=work_wires
         )
+
     if not callable(op):
         raise ValueError(
             f"The object {op} of type {type(op)} is not an Operator or callable. "
@@ -190,6 +196,115 @@ def ctrl(op, control, control_values=None, work_wires=None):
     return wrapper
 
 
+@functools.lru_cache()
+def _get_special_ops():
+    """Gets a list of special operations with custom controlled versions.
+
+    This is placed inside a function to avoid circular imports.
+
+    """
+
+    ops_with_custom_ctrl_ops = {
+        (qml.PauliZ, 1): qml.CZ,
+        (qml.PauliZ, 2): qml.CCZ,
+        (qml.PauliY, 1): qml.CY,
+        (qml.CZ, 1): qml.CCZ,
+        (qml.SWAP, 1): qml.CSWAP,
+        (qml.Hadamard, 1): qml.CH,
+        (qml.RX, 1): qml.CRX,
+        (qml.RY, 1): qml.CRY,
+        (qml.RZ, 1): qml.CRZ,
+        (qml.Rot, 1): qml.CRot,
+        (qml.PhaseShift, 1): qml.ControlledPhaseShift,
+    }
+    return ops_with_custom_ctrl_ops
+
+
+@functools.lru_cache()
+def _get_pauli_x_based_ops():
+    """Gets a list of pauli-x based operations
+
+    This is placed inside a function to avoid circular imports.
+
+    """
+    return qml.PauliX, qml.CNOT, qml.Toffoli, qml.MultiControlledX
+
+
+def _try_wrap_in_custom_ctrl_op(op, control, control_values=None, work_wires=None):
+    """Wraps a controlled operation in custom ControlledOp, returns None if not applicable."""
+
+    ops_with_custom_ctrl_ops = _get_special_ops()
+    custom_key = (type(op), len(control))
+
+    if custom_key in ops_with_custom_ctrl_ops and all(control_values):
+        qml.QueuingManager.remove(op)
+        return ops_with_custom_ctrl_ops[custom_key](*op.data, control + op.wires)
+
+    if isinstance(op, qml.QubitUnitary):
+        return qml.ControlledQubitUnitary(
+            op, control_wires=control, control_values=control_values, work_wires=work_wires
+        )
+
+    # A controlled ControlledPhaseShift should not be compressed to a multi controlled
+    # PhaseShift because the decomposition of PhaseShift contains a GlobalPhase that we
+    # do not have a controlled version of.
+    # TODO: remove this special case when we support ControlledGlobalPhase (sc-44933)
+    if isinstance(op, qml.ControlledPhaseShift):
+        return Controlled(
+            op, control_wires=control, control_values=control_values, work_wires=work_wires
+        )
+    # Similarly, compress the bottom levels of a multi-controlled PhaseShift to a
+    # ControlledPhaseShift if possible to avoid dealing with a controlled GlobalPhase
+    # during decomposition. This should also be removed in the future.
+    if isinstance(op, qml.PhaseShift) and control_values[-1]:
+        op = qml.ControlledPhaseShift(*op.data, wires=control[-1:] + op.wires)
+        return Controlled(
+            op,
+            control_wires=control[:-1],
+            control_values=control_values[:-1],
+            work_wires=work_wires,
+        )
+
+    return None
+
+
+def _handle_pauli_x_based_controlled_ops(op, control, control_values, work_wires):
+    """Handles PauliX-based controlled operations."""
+
+    op_map = {
+        (qml.PauliX, 1): qml.CNOT,
+        (qml.PauliX, 2): qml.Toffoli,
+        (qml.CNOT, 1): qml.Toffoli,
+    }
+
+    custom_key = (type(op), len(control))
+    if custom_key in op_map and all(control_values):
+        qml.QueuingManager.remove(op)
+        return op_map[custom_key](wires=control + op.wires)
+
+    if isinstance(op, qml.PauliX):
+        return qml.MultiControlledX(
+            wires=control + op.wires, control_values=control_values, work_wires=work_wires
+        )
+
+    # TODO: remove special handling of CNOT and Toffoli when they inherit from Controlled
+    if isinstance(op, qml.CNOT):
+        return qml.MultiControlledX(
+            wires=control + op.wires, control_values=control_values + [1], work_wires=work_wires
+        )
+    if isinstance(op, qml.Toffoli):
+        return qml.MultiControlledX(
+            wires=control + op.wires, control_values=control_values + [1, 1], work_wires=work_wires
+        )
+
+    work_wires = work_wires or []
+    return qml.MultiControlledX(
+        wires=control + op.wires,
+        control_values=control_values + op.control_values,
+        work_wires=work_wires + op.work_wires,
+    )
+
+
 # pylint: disable=too-many-arguments, too-many-public-methods
 class Controlled(SymbolicOp):
     """Symbolic operator denoting a controlled operator.
@@ -531,7 +646,7 @@ class Controlled(SymbolicOp):
             return True
         if len(self.control_wires) == 1 and hasattr(self.base, "_controlled"):
             return True
-        if isinstance(self.base, qml.PauliX):
+        if isinstance(self.base, _get_pauli_x_based_ops()):
             return True
         if _is_single_qubit_special_unitary(self.base):
             return True
@@ -631,21 +746,57 @@ def _is_single_qubit_special_unitary(op):
     return qmlmath.allclose(det, 1)
 
 
-# pylint: disable=protected-access
-def _decompose_no_control_values(op: "operation.Operator") -> List["operation.Operator"]:
-    """Provides a decomposition without considering control values. Returns None if
-    no decomposition.
-    """
+def _decompose_pauli_x_based_no_control_values(op: Controlled):
+    """Decomposes a PauliX-based operation"""
+
+    if isinstance(op.base, qml.PauliX) and len(op.control_wires) == 1:
+        return [qml.CNOT(wires=op.active_wires)]
+
+    if isinstance(op.base, qml.PauliX) and len(op.control_wires) == 2:
+        return qml.Toffoli.compute_decomposition(wires=op.active_wires)
+
+    if isinstance(op.base, qml.CNOT) and len(op.control_wires) == 1:
+        return qml.Toffoli.compute_decomposition(wires=op.active_wires)
+
+    return qml.MultiControlledX.compute_decomposition(
+        wires=op.active_wires,
+        work_wires=op.work_wires,
+    )
+
+
+def _decompose_custom_ops(op: Controlled) -> List["operation.Operator"]:
+    """Custom handling for decomposing a controlled operation"""
+
+    pauli_x_based_ctrl_ops = _get_pauli_x_based_ops()
+    ops_with_custom_ctrl_ops = _get_special_ops()
+
+    custom_key = (type(op.base), len(op.control_wires))
+    if custom_key in ops_with_custom_ctrl_ops:
+        custom_op_cls = ops_with_custom_ctrl_ops[custom_key]
+        return custom_op_cls.compute_decomposition(*op.data, op.active_wires)
+    if type(op.base) in pauli_x_based_ctrl_ops:
+        # has some special case handling of its own for further decomposition
+        return _decompose_pauli_x_based_no_control_values(op)
+
+    # TODO: will be removed in the second part of the controlled rework
     if len(op.control_wires) == 1 and hasattr(op.base, "_controlled"):
-        result = op.base._controlled(op.control_wires[0])
+        result = op.base._controlled(op.control_wires[0])  # pylint: disable=protected-access
         # disallow decomposing to itself
         # pylint: disable=unidiomatic-typecheck
         if type(result) != type(op):
             return [result]
         qml.QueuingManager.remove(result)
-    if isinstance(op.base, qml.PauliX):
-        # has some special case handling of its own for further decomposition
-        return [qml.MultiControlledX(wires=op.active_wires, work_wires=op.work_wires)]
+
+    return None
+
+
+def _decompose_no_control_values(op: Controlled) -> List["operation.Operator"]:
+    """Decompose without considering control values. Returns None if no decomposition."""
+
+    decomp = _decompose_custom_ops(op)
+    if decomp is not None:
+        return decomp
+
     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)
@@ -663,7 +814,7 @@ def _decompose_no_control_values(op: "operation.Operator") -> List["operation.Op
             UserWarning,
         )
 
-    return [Controlled(newop, op.control_wires, work_wires=op.work_wires) for newop in base_decomp]
+    return [ctrl(newop, op.control_wires, work_wires=op.work_wires) for newop in base_decomp]
 
 
 class ControlledOp(Controlled, operation.Operation):
diff --git a/pennylane/ops/op_math/controlled_decompositions.py b/pennylane/ops/op_math/controlled_decompositions.py
index bbf14575a..b66a71386 100644
--- a/pennylane/ops/op_math/controlled_decompositions.py
+++ b/pennylane/ops/op_math/controlled_decompositions.py
@@ -199,15 +199,15 @@ def ctrl_decomp_zyz(target_operation: Operator, control_wires: Wires):
         decomp.extend(
             [
                 qml.RY(theta / 2, wires=target_wire),
-                qml.MultiControlledX(wires=control_wires + target_wire),
+                qml.ctrl(qml.PauliX(wires=target_wire), control=control_wires),
                 qml.RY(-theta / 2, wires=target_wire),
             ]
         )
     else:
-        decomp.append(qml.MultiControlledX(wires=control_wires + target_wire))
+        decomp.append(qml.ctrl(qml.PauliX(wires=target_wire), control=control_wires))
     if not qml.math.isclose(-(phi + omega) / 2, 0.0, atol=1e-6, rtol=0):
         decomp.append(qml.RZ(-(phi + omega) / 2, wires=target_wire))
-    decomp.append(qml.MultiControlledX(wires=control_wires + target_wire))
+    decomp.append(qml.ctrl(qml.PauliX(wires=target_wire), control=control_wires))
     if not qml.math.isclose((omega - phi) / 2, 0.0, atol=1e-8, rtol=0):
         decomp.append(qml.RZ((omega - phi) / 2, wires=target_wire))
 
@@ -256,9 +256,9 @@ def _ctrl_decomp_bisect_od(
 
     def component():
         return [
-            qml.MultiControlledX(wires=control_k1 + target_wire, work_wires=control_k2),
+            qml.ctrl(qml.PauliX(wires=target_wire), control=control_k1, work_wires=control_k2),
             qml.QubitUnitary(a, target_wire),
-            qml.MultiControlledX(wires=control_k2 + target_wire, work_wires=control_k1),
+            qml.ctrl(qml.PauliX(wires=target_wire), control=control_k2, work_wires=control_k1),
             qml.adjoint(qml.QubitUnitary(a, target_wire)),
         ]
 
@@ -350,9 +350,9 @@ def _ctrl_decomp_bisect_general(
 
     component = [
         qml.QubitUnitary(c2t, target_wire),
-        qml.MultiControlledX(wires=control_k2 + target_wire, work_wires=control_k1),
+        qml.ctrl(qml.PauliX(wires=target_wire), control=control_k2, work_wires=control_k1),
         qml.adjoint(qml.QubitUnitary(c1, target_wire)),
-        qml.MultiControlledX(wires=control_k1 + target_wire, work_wires=control_k2),
+        qml.ctrl(qml.PauliX(wires=target_wire), control=control_k1, work_wires=control_k2),
     ]
 
     od_decomp = _ctrl_decomp_bisect_od(d, target_wire, control_wires)
@@ -439,3 +439,76 @@ def ctrl_decomp_bisect(
         return _ctrl_decomp_bisect_md(target_matrix, target_wire, control_wires)
     # General algorithm - 20n+O(1) CNOTs
     return _ctrl_decomp_bisect_general(target_matrix, target_wire, control_wires)
+
+
+def decompose_mcx(control_wires, target_wire, work_wires):
+    """Decomposes the multi-controlled PauliX gate"""
+
+    num_work_wires_needed = len(control_wires) - 2
+
+    if len(work_wires) >= num_work_wires_needed:
+        return _decompose_mcx_with_many_workers(control_wires, target_wire, work_wires)
+
+    return _decompose_mcx_with_one_worker(control_wires, target_wire, work_wires[0])
+
+
+def _decompose_mcx_with_many_workers(control_wires, target_wire, work_wires):
+    """Decomposes the multi-controlled PauliX gate using the approach in Lemma 7.2 of
+    https://arxiv.org/abs/quant-ph/9503016, which requires a suitably large register of
+    work wires"""
+    num_work_wires_needed = len(control_wires) - 2
+    work_wires = work_wires[:num_work_wires_needed]
+
+    work_wires_reversed = list(reversed(work_wires))
+    control_wires_reversed = list(reversed(control_wires))
+
+    gates = []
+
+    for i in range(len(work_wires)):
+        ctrl1 = control_wires_reversed[i]
+        ctrl2 = work_wires_reversed[i]
+        t = target_wire if i == 0 else work_wires_reversed[i - 1]
+        gates.append(qml.Toffoli(wires=[ctrl1, ctrl2, t]))
+
+    gates.append(qml.Toffoli(wires=[*control_wires[:2], work_wires[0]]))
+
+    for i in reversed(range(len(work_wires))):
+        ctrl1 = control_wires_reversed[i]
+        ctrl2 = work_wires_reversed[i]
+        t = target_wire if i == 0 else work_wires_reversed[i - 1]
+        gates.append(qml.Toffoli(wires=[ctrl1, ctrl2, t]))
+
+    for i in range(len(work_wires) - 1):
+        ctrl1 = control_wires_reversed[i + 1]
+        ctrl2 = work_wires_reversed[i + 1]
+        t = work_wires_reversed[i]
+        gates.append(qml.Toffoli(wires=[ctrl1, ctrl2, t]))
+
+    gates.append(qml.Toffoli(wires=[*control_wires[:2], work_wires[0]]))
+
+    for i in reversed(range(len(work_wires) - 1)):
+        ctrl1 = control_wires_reversed[i + 1]
+        ctrl2 = work_wires_reversed[i + 1]
+        t = work_wires_reversed[i]
+        gates.append(qml.Toffoli(wires=[ctrl1, ctrl2, t]))
+
+    return gates
+
+
+def _decompose_mcx_with_one_worker(control_wires, target_wire, work_wire):
+    """Decomposes the multi-controlled PauliX gate using the approach in Lemma 7.3 of
+    https://arxiv.org/abs/quant-ph/9503016, which requires a single work wire"""
+    tot_wires = len(control_wires) + 2
+    partition = int(np.ceil(tot_wires / 2))
+
+    first_part = control_wires[:partition]
+    second_part = control_wires[partition:]
+
+    gates = [
+        qml.ctrl(qml.PauliX(work_wire), control=first_part, work_wires=second_part + target_wire),
+        qml.ctrl(qml.PauliX(target_wire), control=second_part + work_wire, work_wires=first_part),
+        qml.ctrl(qml.PauliX(work_wire), control=first_part, work_wires=second_part + target_wire),
+        qml.ctrl(qml.PauliX(target_wire), control=second_part + work_wire, work_wires=first_part),
+    ]
+
+    return gates
diff --git a/pennylane/ops/op_math/controlled_ops.py b/pennylane/ops/op_math/controlled_ops.py
index 31ab1c607..5601351bb 100644
--- a/pennylane/ops/op_math/controlled_ops.py
+++ b/pennylane/ops/op_math/controlled_ops.py
@@ -20,12 +20,17 @@ from typing import Iterable
 from functools import lru_cache
 
 import numpy as np
+from scipy.linalg import block_diag
 
 import pennylane as qml
-from pennylane.operation import AnyWires
+from pennylane.operation import AnyWires, Wires
 from pennylane.ops.qubit.matrix_ops import QubitUnitary
 from pennylane.ops.qubit.parametric_ops_single_qubit import stack_last
 from .controlled import ControlledOp
+from .controlled_decompositions import decompose_mcx
+
+
+INV_SQRT2 = 1 / qml.math.sqrt(2)
 
 
 # pylint: disable=too-few-public-methods
@@ -141,8 +146,8 @@ class ControlledQubitUnitary(ControlledOp):
         self._name = "ControlledQubitUnitary"
 
     def _controlled(self, wire):
-        ctrl_wires = self.control_wires + wire
-        values = None if self.control_values is None else self.control_values + [True]
+        ctrl_wires = wire + self.control_wires
+        values = None if self.control_values is None else [True] + self.control_values
         return ControlledQubitUnitary(
             self.base,
             control_wires=ctrl_wires,
@@ -339,6 +344,242 @@ class CZ(ControlledOp):
     def _controlled(self, wire):
         return qml.CCZ(wires=wire + self.wires)
 
+    @staticmethod
+    def compute_decomposition(wires):  # pylint: disable=arguments-differ
+        return [qml.ControlledPhaseShift(np.pi, wires=wires)]
+
+    def decomposition(self):
+        return self.compute_decomposition(self.wires)
+
+
+def _check_and_convert_control_values(control_values, control_wires):
+
+    if isinstance(control_values, str):
+        # Make sure all values are either 0 or 1
+        if not set(control_values).issubset({"1", "0"}):
+            raise ValueError("String of control values can contain only '0' or '1'.")
+
+        control_values = [int(x) for x in control_values]
+
+    if control_values is None:
+        return [1] * len(control_wires)
+
+    if len(control_values) != len(control_wires):
+        raise ValueError("Length of control values must equal number of control wires.")
+
+    return control_values
+
+
+class MultiControlledX(ControlledOp):
+    r"""MultiControlledX(control_wires, wires, control_values)
+    Apply a Pauli X gate controlled on an arbitrary computational basis state.
+
+    **Details:**
+
+    * Number of wires: Any (the operation can act on any number of wires)
+    * Number of parameters: 0
+    * Gradient recipe: None
+
+    Args:
+        control_wires (Union[Wires, Sequence[int], or int]): Deprecated way to indicate the control wires.
+            Now users should use "wires" to indicate both the control wires and the target wire.
+        wires (Union[Wires, Sequence[int], or int]): control wire(s) followed by a single target wire where
+            the operation acts on
+        control_values (Union[bool, list[bool], int, list[int]]): The value(s) the control wire(s)
+                should take. Integers other than 0 or 1 will be treated as ``int(bool(x))``.
+        work_wires (Union[Wires, Sequence[int], or int]): optional work wires used to decompose
+            the operation into a series of Toffoli gates
+
+
+    .. note::
+
+        If ``MultiControlledX`` is not supported on the targeted device, PennyLane will decompose
+        the operation into :class:`~.Toffoli` and/or :class:`~.CNOT` gates. When controlling on
+        three or more wires, the Toffoli-based decompositions described in Lemmas 7.2 and 7.3 of
+        `Barenco et al. <https://arxiv.org/abs/quant-ph/9503016>`__ will be used. These methods
+        require at least one work wire.
+
+        The number of work wires provided determines the decomposition method used and the resulting
+        number of Toffoli gates required. When ``MultiControlledX`` is controlling on :math:`n`
+        wires:
+
+        #. If at least :math:`n - 2` work wires are provided, the decomposition in Lemma 7.2 will be
+           applied using the first :math:`n - 2` work wires.
+        #. If fewer than :math:`n - 2` work wires are provided, a combination of Lemmas 7.3 and 7.2
+           will be applied using only the first work wire.
+
+        These methods present a tradeoff between qubit number and depth. The method in point 1
+        requires fewer Toffoli gates but a greater number of qubits.
+
+        Note that the state of the work wires before and after the decomposition takes place is
+        unchanged.
+
+    """
+
+    is_self_inverse = True
+    """bool: Whether or not the operator is self-inverse."""
+
+    num_wires = AnyWires
+    """int: Number of wires the operation acts on."""
+
+    num_params = 0
+    """int: Number of trainable parameters that the operator depends on."""
+
+    ndim_params = ()
+    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
+
+    name = "MultiControlledX"
+
+    def _flatten(self):
+        return (), (self.active_wires, tuple(self.control_values), self.work_wires)
+
+    @classmethod
+    def _unflatten(cls, _, metadata):
+        return cls(wires=metadata[0], control_values=metadata[1], work_wires=metadata[2])
+
+    # pylint: disable=too-many-arguments
+    def __init__(self, control_wires=None, wires=None, control_values=None, work_wires=None):
+        if wires is None:
+            raise ValueError("Must specify the wires where the operation acts on")
+        wires = wires if isinstance(wires, Wires) else Wires(wires)
+        if control_wires is not None:
+            warnings.warn(
+                "The control_wires keyword will be removed soon. Use wires = (control_wires, "
+                "target_wire) instead. See the documentation for more information.",
+                UserWarning,
+            )
+            if len(wires) != 1:
+                raise ValueError("MultiControlledX accepts a single target wire.")
+        else:
+            if len(wires) < 2:
+                raise ValueError(
+                    f"MultiControlledX: wrong number of wires. {len(wires)} wire(s) given. Need at least 2."
+                )
+            control_wires = wires[:-1]
+            wires = wires[-1:]
+
+        control_values = _check_and_convert_control_values(control_values, control_wires)
+
+        super().__init__(
+            qml.PauliX(wires=wires),
+            control_wires=control_wires,
+            control_values=control_values,
+            work_wires=work_wires,
+        )
+
+    def __repr__(self):
+        return f"MultiControlledX(wires={self.active_wires.tolist()}, control_values={self.control_values})"
+
+    @property
+    def wires(self):
+        return self.active_wires
+
+    # pylint: disable=unused-argument, arguments-differ
+    @staticmethod
+    def compute_matrix(control_wires, control_values=None, **kwargs):
+        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
+
+        The canonical matrix is the textbook matrix representation that does not consider wires.
+        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
+
+        .. seealso:: :meth:`~.MultiControlledX.matrix`
+
+        Args:
+            control_wires (Any or Iterable[Any]): wires to place controls on
+            control_values (Union[bool, list[bool], int, list[int]]): The value(s) the control wire(s)
+                should take. Integers other than 0 or 1 will be treated as ``int(bool(x))``.
+
+        Returns:
+            tensor_like: matrix representation
+
+        **Example**
+
+        >>> print(qml.MultiControlledX.compute_matrix([0], 1))
+        [[1. 0. 0. 0.]
+         [0. 1. 0. 0.]
+         [0. 0. 0. 1.]
+         [0. 0. 1. 0.]]
+        >>> print(qml.MultiControlledX.compute_matrix([1], 0))
+        [[0. 1. 0. 0.]
+         [1. 0. 0. 0.]
+         [0. 0. 1. 0.]
+         [0. 0. 0. 1.]]
+
+        """
+
+        control_values = _check_and_convert_control_values(control_values, control_wires)
+        padding_left = sum(2**i * int(val) for i, val in enumerate(reversed(control_values))) * 2
+        padding_right = 2 ** (len(control_wires) + 1) - 2 - padding_left
+        return block_diag(np.eye(padding_left), qml.PauliX.compute_matrix(), np.eye(padding_right))
+
+    def matrix(self, wire_order=None):
+        canonical_matrix = self.compute_matrix(self.control_wires, self.control_values)
+        wire_order = wire_order or self.wires
+        return qml.math.expand_matrix(
+            canonical_matrix, wires=self.active_wires, wire_order=wire_order
+        )
+
+    # pylint: disable=unused-argument, arguments-differ
+    @staticmethod
+    def compute_decomposition(wires=None, work_wires=None, control_values=None, **kwargs):
+        r"""Representation of the operator as a product of other operators (static method).
+
+        .. math:: O = O_1 O_2 \dots O_n.
+
+        .. seealso:: :meth:`~.MultiControlledX.decomposition`.
+
+        Args:
+            wires (Iterable[Any] or Wires): wires that the operation acts on
+            work_wires (Wires): optional work wires used to decompose
+                the operation into a series of Toffoli gates.
+            control_values (Union[bool, list[bool], int, list[int]]): The value(s) the control wire(s)
+                should take. Integers other than 0 or 1 will be treated as ``int(bool(x))``.
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+
+        **Example:**
+
+        >>> print(qml.MultiControlledX.compute_decomposition(
+        ...     wires=[0,1,2,3], control_values=[1,1,1], work_wires=qml.wires.Wires("aux")))
+        [Toffoli(wires=[2, 'aux', 3]),
+        Toffoli(wires=[0, 1, 'aux']),
+        Toffoli(wires=[2, 'aux', 3]),
+        Toffoli(wires=[0, 1, 'aux'])]
+
+        """
+
+        if len(wires) < 2:
+            raise ValueError(f"Wrong number of wires. {len(wires)} given. Need at least 2.")
+
+        target_wire = wires[-1]
+        control_wires = wires[:-1]
+
+        if control_values is None:
+            control_values = [True] * len(control_wires)
+
+        work_wires = work_wires or []
+        if len(control_wires) > 2 and len(work_wires) == 0:
+            raise ValueError(
+                "At least one work wire is required to decompose operation: MultiControlledX"
+            )
+
+        flips1 = [qml.PauliX(wires=w) for w, val in zip(control_wires, control_values) if not val]
+
+        if len(control_wires) == 1:
+            decomp = [qml.CNOT(wires=wires)]
+        elif len(control_wires) == 2:
+            decomp = qml.Toffoli.compute_decomposition(wires=wires)
+        else:
+            decomp = decompose_mcx(control_wires, target_wire, work_wires)
+
+        flips2 = [qml.PauliX(wires=w) for w, val in zip(control_wires, control_values) if not val]
+
+        return flips1 + decomp + flips2
+
+    def decomposition(self):
+        return self.compute_decomposition(self.active_wires, self.work_wires, self.control_values)
+
 
 class CRX(ControlledOp):
     r"""The controlled-RX operator
@@ -423,6 +664,7 @@ class CRX(ControlledOp):
                 [0.0+0.0j, 0.0+0.0j, 0.9689+0.0j, 0.0-0.2474j],
                 [0.0+0.0j, 0.0+0.0j, 0.0-0.2474j, 0.9689+0.0j]])
         """
+
         interface = qml.math.get_interface(theta)
 
         c = qml.math.cos(theta / 2)
diff --git a/pennylane/ops/qubit/__init__.py b/pennylane/ops/qubit/__init__.py
index d9484a05c..147f41009 100644
--- a/pennylane/ops/qubit/__init__.py
+++ b/pennylane/ops/qubit/__init__.py
@@ -88,7 +88,6 @@ __ops__ = {
     "QubitUnitary",
     "BlockEncode",
     "SpecialUnitary",
-    "MultiControlledX",
     "IntegerComparator",
     "DiagonalQubitUnitary",
     "SingleExcitation",
diff --git a/pennylane/ops/qubit/arithmetic_ops.py b/pennylane/ops/qubit/arithmetic_ops.py
index 27236b782..fa84ab73c 100644
--- a/pennylane/ops/qubit/arithmetic_ops.py
+++ b/pennylane/ops/qubit/arithmetic_ops.py
@@ -24,7 +24,6 @@ import pennylane as qml
 from pennylane.operation import AnyWires, Operation
 from pennylane.wires import Wires
 from pennylane.ops import Identity
-from pennylane.ops.qubit.non_parametric_ops import MultiControlledX
 
 
 class QubitCarry(Operation):
@@ -465,7 +464,7 @@ class IntegerComparator(Operation):
             control_values_list = [format(n, binary) for n in values]
             mat = np.eye(2 ** (len(control_wires) + 1))
             for control_values in control_values_list:
-                mat = mat @ MultiControlledX.compute_matrix(
+                mat = mat @ qml.MultiControlledX.compute_matrix(
                     control_wires, control_values=control_values
                 )
 
@@ -523,7 +522,7 @@ class IntegerComparator(Operation):
             gates = []
             for control_values in control_values_list:
                 gates.append(
-                    MultiControlledX(
+                    qml.MultiControlledX(
                         wires=control_wires + wires,
                         control_values=control_values,
                         work_wires=work_wires,
diff --git a/pennylane/ops/qubit/non_parametric_ops.py b/pennylane/ops/qubit/non_parametric_ops.py
index 6233444ee..4cc695f2d 100644
--- a/pennylane/ops/qubit/non_parametric_ops.py
+++ b/pennylane/ops/qubit/non_parametric_ops.py
@@ -17,17 +17,15 @@ not depend on any parameters.
 """
 # pylint:disable=abstract-method,arguments-differ,protected-access,invalid-overridden-method, no-member
 import cmath
-import warnings
 from copy import copy
 from functools import lru_cache
 
 import numpy as np
 
 from scipy import sparse
-from scipy.linalg import block_diag
 
 import pennylane as qml
-from pennylane.operation import AnyWires, Observable, Operation
+from pennylane.operation import Observable, Operation
 from pennylane.utils import pauli_eigs
 from pennylane.wires import Wires
 
@@ -158,12 +156,11 @@ class Hadamard(Observable, Operation):
         PhaseShift(1.5707963267948966, wires=[0])]
 
         """
-        decomp_ops = [
+        return [
             qml.PhaseShift(np.pi / 2, wires=wires),
             qml.RX(np.pi / 2, wires=wires),
             qml.PhaseShift(np.pi / 2, wires=wires),
         ]
-        return decomp_ops
 
     def _controlled(self, wire):
         return CH(wires=Wires(wire) + self.wires)
@@ -311,12 +308,11 @@ class PauliX(Observable, Operation):
         PhaseShift(1.5707963267948966, wires=[0])]
 
         """
-        decomp_ops = [
+        return [
             qml.PhaseShift(np.pi / 2, wires=wires),
             qml.RX(np.pi, wires=wires),
             qml.PhaseShift(np.pi / 2, wires=wires),
         ]
-        return decomp_ops
 
     def adjoint(self):
         return PauliX(wires=self.wires)
@@ -469,12 +465,11 @@ class PauliY(Observable, Operation):
         PhaseShift(1.5707963267948966, wires=[0])]
 
         """
-        decomp_ops = [
+        return [
             qml.PhaseShift(np.pi / 2, wires=wires),
             qml.RY(np.pi, wires=wires),
             qml.PhaseShift(np.pi / 2, wires=wires),
         ]
-        return decomp_ops
 
     def adjoint(self):
         return PauliY(wires=self.wires)
@@ -959,13 +954,12 @@ class SX(Operation):
         PhaseShift(1.5707963267948966, wires=[0])]
 
         """
-        decomp_ops = [
+        return [
             qml.RZ(np.pi / 2, wires=wires),
             qml.RY(np.pi / 2, wires=wires),
             qml.RZ(-np.pi, wires=wires),
             qml.PhaseShift(np.pi / 2, wires=wires),
         ]
-        return decomp_ops
 
     def pow(self, z):
         z_mod4 = z % 4
@@ -1224,12 +1218,11 @@ class SWAP(Operation):
         [CNOT(wires=[0, 1]), CNOT(wires=[1, 0]), CNOT(wires=[0, 1])]
 
         """
-        decomp_ops = [
+        return [
             qml.CNOT(wires=[wires[0], wires[1]]),
             qml.CNOT(wires=[wires[1], wires[0]]),
             qml.CNOT(wires=[wires[0], wires[1]]),
         ]
-        return decomp_ops
 
     def pow(self, z):
         return super().pow(z % 2)
@@ -1471,7 +1464,7 @@ class ISWAP(Operation):
         Hadamard(wires=[1])]
 
         """
-        decomp_ops = [
+        return [
             S(wires=wires[0]),
             S(wires=wires[1]),
             Hadamard(wires=wires[0]),
@@ -1479,7 +1472,6 @@ class ISWAP(Operation):
             CNOT(wires=[wires[1], wires[0]]),
             Hadamard(wires=wires[1]),
         ]
-        return decomp_ops
 
     def pow(self, z):
         z_mod2 = z % 2
@@ -1601,7 +1593,7 @@ class SISWAP(Operation):
         SX(wires=[1])]
 
         """
-        decomp_ops = [
+        return [
             SX(wires=wires[0]),
             qml.RZ(np.pi / 2, wires=wires[0]),
             CNOT(wires=[wires[0], wires[1]]),
@@ -1615,7 +1607,6 @@ class SISWAP(Operation):
             SX(wires=wires[0]),
             SX(wires=wires[1]),
         ]
-        return decomp_ops
 
     def pow(self, z):
         z_mod4 = z % 4
@@ -1865,7 +1856,7 @@ class CCZ(Operation):
             T(wires=wires[0]),
             qml.adjoint(T(wires=wires[1])),
             CNOT(wires=[wires[0], wires[1]]),
-            Hadamard(wires=[2]),
+            Hadamard(wires=wires[2]),
         ]
 
     def adjoint(self):
@@ -2024,329 +2015,3 @@ class Toffoli(Operation):
     @property
     def is_hermitian(self):
         return True
-
-
-class MultiControlledX(Operation):
-    r"""MultiControlledX(control_wires, wires, control_values)
-    Apply a Pauli X gate controlled on an arbitrary computational basis state.
-
-    **Details:**
-
-    * Number of wires: Any (the operation can act on any number of wires)
-    * Number of parameters: 0
-    * Gradient recipe: None
-
-    Args:
-        control_wires (Union[Wires, Sequence[int], or int]): Deprecated way to indicate the control wires.
-            Now users should use "wires" to indicate both the control wires and the target wire.
-        wires (Union[Wires, Sequence[int], or int]): control wire(s) followed by a single target wire where
-            the operation acts on
-        control_values (str): a string of bits representing the state of the control
-            wires to control on (default is the all 1s state)
-        work_wires (Union[Wires, Sequence[int], or int]): optional work wires used to decompose
-            the operation into a series of Toffoli gates
-
-
-    .. note::
-
-        If ``MultiControlledX`` is not supported on the targeted device, PennyLane will decompose
-        the operation into :class:`~.Toffoli` and/or :class:`~.CNOT` gates. When controlling on
-        three or more wires, the Toffoli-based decompositions described in Lemmas 7.2 and 7.3 of
-        `Barenco et al. <https://arxiv.org/abs/quant-ph/9503016>`__ will be used. These methods
-        require at least one work wire.
-
-        The number of work wires provided determines the decomposition method used and the resulting
-        number of Toffoli gates required. When ``MultiControlledX`` is controlling on :math:`n`
-        wires:
-
-        #. If at least :math:`n - 2` work wires are provided, the decomposition in Lemma 7.2 will be
-           applied using the first :math:`n - 2` work wires.
-        #. If fewer than :math:`n - 2` work wires are provided, a combination of Lemmas 7.3 and 7.2
-           will be applied using only the first work wire.
-
-        These methods present a tradeoff between qubit number and depth. The method in point 1
-        requires fewer Toffoli gates but a greater number of qubits.
-
-        Note that the state of the work wires before and after the decomposition takes place is
-        unchanged.
-
-    """
-
-    is_self_inverse = True
-    num_wires = AnyWires
-    num_params = 0
-    """int: Number of trainable parameters that the operator depends on."""
-
-    grad_method = None
-
-    def _flatten(self):
-        hyperparameters = (
-            ("wires", self.wires),
-            ("control_values", self.hyperparameters["control_values"]),
-            ("work_wires", self.hyperparameters["work_wires"]),
-        )
-        return tuple(), hyperparameters
-
-    @classmethod
-    def _unflatten(cls, _, metadata):
-        return cls(**dict(metadata))
-
-    # pylint: disable=too-many-arguments
-    def __init__(self, control_wires=None, wires=None, control_values=None, work_wires=None):
-        if wires is None:
-            raise ValueError("Must specify the wires where the operation acts on")
-        if control_wires is None:
-            if len(wires) > 1:
-                control_wires = Wires(wires[:-1])
-                wires = Wires(wires[-1])
-            else:
-                raise ValueError(
-                    "MultiControlledX: wrong number of wires. "
-                    f"{len(wires)} wire(s) given. Need at least 2."
-                )
-        else:
-            wires = Wires(wires)
-            control_wires = Wires(control_wires)
-
-            warnings.warn(
-                "The control_wires keyword will be removed soon. "
-                "Use wires = (control_wires, target_wire) instead. "
-                "See the documentation for more information.",
-                category=UserWarning,
-            )
-
-            if len(wires) != 1:
-                raise ValueError("MultiControlledX accepts a single target wire.")
-
-        work_wires = Wires([]) if work_wires is None else Wires(work_wires)
-        total_wires = control_wires + wires
-
-        if Wires.shared_wires([total_wires, work_wires]):
-            raise ValueError("The work wires must be different from the control and target wires")
-
-        if not control_values:
-            control_values = "1" * len(control_wires)
-
-        self.hyperparameters["control_wires"] = control_wires
-        self.hyperparameters["work_wires"] = work_wires
-        self.hyperparameters["control_values"] = control_values
-        self.total_wires = total_wires
-
-        super().__init__(wires=self.total_wires)
-
-    def __repr__(self):
-        return f'MultiControlledX(wires={list(self.total_wires._labels)}, control_values="{self.hyperparameters["control_values"]}")'
-
-    def label(self, decimals=None, base_label=None, cache=None):
-        return base_label or "X"
-
-    # pylint: disable=unused-argument
-    @staticmethod
-    def compute_matrix(
-        control_wires, control_values=None, **kwargs
-    ):  # pylint: disable=arguments-differ
-        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
-
-        The canonical matrix is the textbook matrix representation that does not consider wires.
-        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
-
-        .. seealso:: :meth:`~.MultiControlledX.matrix`
-
-        Args:
-            control_wires (Any or Iterable[Any]): wires to place controls on
-            control_values (str): string of bits determining the controls
-
-        Returns:
-           tensor_like: matrix representation
-
-        **Example**
-
-        >>> print(qml.MultiControlledX.compute_matrix([0], '1'))
-        [[1. 0. 0. 0.]
-         [0. 1. 0. 0.]
-         [0. 0. 0. 1.]
-         [0. 0. 1. 0.]]
-        >>> print(qml.MultiControlledX.compute_matrix([1], '0'))
-        [[0. 1. 0. 0.]
-         [1. 0. 0. 0.]
-         [0. 0. 1. 0.]
-         [0. 0. 0. 1.]]
-
-        """
-        if control_values is None:
-            control_values = "1" * len(control_wires)
-
-        if isinstance(control_values, str):
-            if len(control_values) != len(control_wires):
-                raise ValueError("Length of control bit string must equal number of control wires.")
-
-            # Make sure all values are either 0 or 1
-            if not set(control_values).issubset({"1", "0"}):
-                raise ValueError("String of control values can contain only '0' or '1'.")
-
-            control_int = int(control_values, 2)
-        else:
-            raise ValueError("Control values must be passed as a string.")
-
-        padding_left = control_int * 2
-        padding_right = 2 ** (len(control_wires) + 1) - 2 - padding_left
-        cx = block_diag(np.eye(padding_left), PauliX.compute_matrix(), np.eye(padding_right))
-        return cx
-
-    @property
-    def control_wires(self):
-        return self.wires[:~0]
-
-    def adjoint(self):
-        return MultiControlledX(
-            wires=self.wires,
-            control_values=self.hyperparameters["control_values"],
-        )
-
-    def pow(self, z):
-        return super().pow(z % 2)
-
-    @staticmethod
-    def compute_decomposition(wires=None, work_wires=None, control_values=None, **kwargs):
-        r"""Representation of the operator as a product of other operators (static method).
-
-        .. math:: O = O_1 O_2 \dots O_n.
-
-
-        .. seealso:: :meth:`~.MultiControlledX.decomposition`.
-
-        Args:
-            wires (Iterable[Any] or Wires): wires that the operation acts on
-            work_wires (Wires): optional work wires used to decompose
-                the operation into a series of Toffoli gates.
-            control_values (str): a string of bits representing the state of the control
-                wires to control on (default is the all 1s state)
-
-        Returns:
-            list[Operator]: decomposition into lower level operations
-
-        **Example:**
-
-        >>> print(qml.MultiControlledX.compute_decomposition(wires=[0,1,2,3],control_values="111", work_wires=qml.wires.Wires("aux")))
-        [Toffoli(wires=[2, 'aux', 3]),
-        Toffoli(wires=[0, 1, 'aux']),
-        Toffoli(wires=[2, 'aux', 3]),
-        Toffoli(wires=[0, 1, 'aux'])]
-
-        """
-
-        target_wire = wires[~0]
-        control_wires = wires[:~0]
-
-        if control_values is None:
-            control_values = "1" * len(control_wires)
-
-        if len(control_wires) > 2 and len(work_wires) == 0:
-            raise ValueError(
-                "At least one work wire is required to decompose operation: MultiControlledX"
-            )
-
-        flips1 = [
-            qml.PauliX(control_wires[i]) for i, val in enumerate(control_values) if val == "0"
-        ]
-
-        if len(control_wires) == 1:
-            decomp = [qml.CNOT(wires=[control_wires[0], target_wire])]
-        elif len(control_wires) == 2:
-            decomp = [qml.Toffoli(wires=[*control_wires, target_wire])]
-        else:
-            num_work_wires_needed = len(control_wires) - 2
-
-            if len(work_wires) >= num_work_wires_needed:
-                decomp = MultiControlledX._decomposition_with_many_workers(
-                    control_wires, target_wire, work_wires
-                )
-            else:
-                work_wire = work_wires[0]
-                decomp = MultiControlledX._decomposition_with_one_worker(
-                    control_wires, target_wire, work_wire
-                )
-
-        flips2 = [
-            qml.PauliX(control_wires[i]) for i, val in enumerate(control_values) if val == "0"
-        ]
-
-        return flips1 + decomp + flips2
-
-    @staticmethod
-    def _decomposition_with_many_workers(control_wires, target_wire, work_wires):
-        """Decomposes the multi-controlled PauliX gate using the approach in Lemma 7.2 of
-        https://arxiv.org/abs/quant-ph/9503016, which requires a suitably large register of
-        work wires"""
-        num_work_wires_needed = len(control_wires) - 2
-        work_wires = work_wires[:num_work_wires_needed]
-
-        work_wires_reversed = list(reversed(work_wires))
-        control_wires_reversed = list(reversed(control_wires))
-
-        gates = []
-
-        for i in range(len(work_wires)):
-            ctrl1 = control_wires_reversed[i]
-            ctrl2 = work_wires_reversed[i]
-            t = target_wire if i == 0 else work_wires_reversed[i - 1]
-            gates.append(qml.Toffoli(wires=[ctrl1, ctrl2, t]))
-
-        gates.append(qml.Toffoli(wires=[*control_wires[:2], work_wires[0]]))
-
-        for i in reversed(range(len(work_wires))):
-            ctrl1 = control_wires_reversed[i]
-            ctrl2 = work_wires_reversed[i]
-            t = target_wire if i == 0 else work_wires_reversed[i - 1]
-            gates.append(qml.Toffoli(wires=[ctrl1, ctrl2, t]))
-
-        for i in range(len(work_wires) - 1):
-            ctrl1 = control_wires_reversed[i + 1]
-            ctrl2 = work_wires_reversed[i + 1]
-            t = work_wires_reversed[i]
-            gates.append(qml.Toffoli(wires=[ctrl1, ctrl2, t]))
-
-        gates.append(qml.Toffoli(wires=[*control_wires[:2], work_wires[0]]))
-
-        for i in reversed(range(len(work_wires) - 1)):
-            ctrl1 = control_wires_reversed[i + 1]
-            ctrl2 = work_wires_reversed[i + 1]
-            t = work_wires_reversed[i]
-            gates.append(qml.Toffoli(wires=[ctrl1, ctrl2, t]))
-
-        return gates
-
-    @staticmethod
-    def _decomposition_with_one_worker(control_wires, target_wire, work_wire):
-        """Decomposes the multi-controlled PauliX gate using the approach in Lemma 7.3 of
-        https://arxiv.org/abs/quant-ph/9503016, which requires a single work wire"""
-        tot_wires = len(control_wires) + 2
-        partition = int(np.ceil(tot_wires / 2))
-
-        first_part = control_wires[:partition]
-        second_part = control_wires[partition:]
-
-        gates = [
-            MultiControlledX(
-                wires=first_part + work_wire,
-                work_wires=second_part + target_wire,
-            ),
-            MultiControlledX(
-                wires=second_part + work_wire + target_wire,
-                work_wires=first_part,
-            ),
-            MultiControlledX(
-                wires=first_part + work_wire,
-                work_wires=second_part + target_wire,
-            ),
-            MultiControlledX(
-                wires=second_part + work_wire + target_wire,
-                work_wires=first_part,
-            ),
-        ]
-
-        return gates
-
-    @property
-    def is_hermitian(self):
-        return True
diff --git a/pennylane/ops/qubit/parametric_ops_single_qubit.py b/pennylane/ops/qubit/parametric_ops_single_qubit.py
index 0bde74d6d..5ca6ff93f 100644
--- a/pennylane/ops/qubit/parametric_ops_single_qubit.py
+++ b/pennylane/ops/qubit/parametric_ops_single_qubit.py
@@ -655,12 +655,11 @@ class Rot(Operation):
         [RZ(1.2, wires=[0]), RY(2.3, wires=[0]), RZ(3.4, wires=[0])]
 
         """
-        decomp_ops = [
+        return [
             RZ(phi, wires=wires),
             RY(theta, wires=wires),
             RZ(omega, wires=wires),
         ]
-        return decomp_ops
 
     def adjoint(self):
         phi, theta, omega = self.parameters
@@ -927,12 +926,11 @@ class U2(Operation):
 
         """
         pi_half = qml.math.ones_like(delta) * (np.pi / 2)
-        decomp_ops = [
+        return [
             Rot(delta, pi_half, -delta, wires=wires),
             PhaseShift(delta, wires=wires),
             PhaseShift(phi, wires=wires),
         ]
-        return decomp_ops
 
     def adjoint(self):
         phi, delta = self.parameters
@@ -1084,12 +1082,11 @@ class U3(Operation):
         PhaseShift(2.34, wires=[0])]
 
         """
-        decomp_ops = [
+        return [
             Rot(delta, theta, -delta, wires=wires),
             PhaseShift(delta, wires=wires),
             PhaseShift(phi, wires=wires),
         ]
-        return decomp_ops
 
     def adjoint(self):
         theta, phi, delta = self.parameters
diff --git a/pennylane/transforms/qmc.py b/pennylane/transforms/qmc.py
index 995f89475..1f3405daf 100644
--- a/pennylane/transforms/qmc.py
+++ b/pennylane/transforms/qmc.py
@@ -122,7 +122,7 @@ def apply_controlled_Q(
         wires = Wires(wires)
         target_wire = Wires(target_wire)
         control_wire = Wires(control_wire)
-        work_wires = Wires(work_wires)
+        work_wires = Wires(work_wires) if work_wires is not None else Wires([])
 
         if not wires.contains_wires(target_wire):
             raise ValueError("The target wire must be contained within wires")


Test patch:
diff --git a/tests/drawer/test_tape_mpl.py b/tests/drawer/test_tape_mpl.py
index 2baab80b4..8af4d1697 100644
--- a/tests/drawer/test_tape_mpl.py
+++ b/tests/drawer/test_tape_mpl.py
@@ -23,6 +23,7 @@ import pennylane as qml
 
 from pennylane.drawer import tape_mpl
 from pennylane.tape import QuantumScript
+from pennylane.ops.op_math import Controlled
 
 mpl = pytest.importorskip("matplotlib")
 plt = pytest.importorskip("matplotlib.pyplot")
@@ -544,9 +545,9 @@ class TestControlledGates:
         when they are provided as a list of bools."""
 
         with qml.queuing.AnnotatedQueue() as q_tape:
-            qml.ctrl(
-                qml.ctrl(qml.PauliX(wires=4), control=[2, 3], control_values=[1, 0]),
-                control=[0, 1],
+            Controlled(
+                qml.ctrl(qml.PauliY(wires=4), control=[2, 3], control_values=[1, 0]),
+                control_wires=[0, 1],
                 control_values=[1, 0],
             )
 
diff --git a/tests/ops/op_math/test_controlled.py b/tests/ops/op_math/test_controlled.py
index 54d3a3ddc..3f804d153 100644
--- a/tests/ops/op_math/test_controlled.py
+++ b/tests/ops/op_math/test_controlled.py
@@ -16,18 +16,30 @@
 from copy import copy
 from functools import partial
 
-import numpy as onp
+import numpy as np
 import pytest
-from gate_data import CNOT, CSWAP, CZ, CRot3, CRotx, CRoty, CRotz, Toffoli
+from gate_data import (
+    CNOT,
+    CSWAP,
+    CY,
+    CZ,
+    CCZ,
+    CH,
+    CRot3,
+    CRotx,
+    CRoty,
+    CRotz,
+    Toffoli,
+    ControlledPhaseShift,
+)
 from scipy import sparse
 
 import pennylane as qml
-from pennylane import numpy as np
+from pennylane import numpy as pnp
 from pennylane.operation import DecompositionUndefinedError, Operation, Operator
 from pennylane.ops.op_math.controlled import (
     Controlled,
     ControlledOp,
-    _decompose_no_control_values,
     ctrl,
 )
 from pennylane.tape import QuantumScript
@@ -38,6 +50,7 @@ from pennylane.wires import Wires
 # pylint: disable=protected-access
 # pylint: disable=pointless-statement
 # pylint: disable=expression-not-assigned
+# pylint: disable=too-many-arguments
 
 
 def equal_list(lhs, rhs):
@@ -48,18 +61,6 @@ def equal_list(lhs, rhs):
     return len(lhs) == len(rhs) and all(qml.equal(l, r) for l, r in zip(lhs, rhs))
 
 
-base_num_control_mats = [
-    (qml.PauliX("a"), 1, CNOT),
-    (qml.PauliZ("a"), 1, CZ),
-    (qml.SWAP(("a", "b")), 1, CSWAP),
-    (qml.PauliX("a"), 2, Toffoli),
-    (qml.RX(1.234, "b"), 1, CRotx(1.234)),
-    (qml.RY(-0.432, "a"), 1, CRoty(-0.432)),
-    (qml.RZ(6.78, "a"), 1, CRotz(6.78)),
-    (qml.Rot(1.234, -0.432, 9.0, "a"), 1, CRot3(1.234, -0.432, 9.0)),
-]
-
-
 class TempOperator(Operator):
     num_wires = 1
 
@@ -68,6 +69,16 @@ class TempOperation(Operation):
     num_wires = 1
 
 
+class OpWithDecomposition(Operation):
+    @staticmethod
+    def compute_decomposition(*params, wires=None, **_):
+        return [
+            qml.Hadamard(wires=wires[0]),
+            qml.S(wires=wires[1]),
+            qml.RX(params[0], wires=wires[0]),
+        ]
+
+
 class TestControlledInheritance:
     """Test the inheritance structure modified through dynamic __new__ method."""
 
@@ -110,7 +121,7 @@ class TestControlledInheritance:
         assert type(op) is ControlledOp  # pylint: disable=unidiomatic-typecheck
 
 
-class TestInitialization:
+class TestControlledInit:
     """Test the initialization process and standard properties."""
 
     temp_op = TempOperator("a")
@@ -182,25 +193,25 @@ class TestInitialization:
             Controlled(self.temp_op, control_wires="b", work_wires="b")
 
 
-class TestProperties:
+class TestControlledProperties:
     """Test the properties of the ``Controlled`` symbolic operator."""
 
     def test_data(self):
         """Test that the base data can be get and set through Controlled class."""
 
-        x = np.array(1.234)
+        x = pnp.array(1.234)
 
         base = qml.RX(x, wires="a")
         op = Controlled(base, (0, 1))
 
         assert op.data == (x,)
 
-        x_new = (np.array(2.3454),)
+        x_new = (pnp.array(2.3454),)
         op.data = x_new
         assert op.data == (x_new,)
         assert base.data == (x_new,)
 
-        x_new2 = (np.array(3.456),)
+        x_new2 = (pnp.array(3.456),)
         base.data = x_new2
         assert op.data == (x_new2,)
         assert op.parameters == [x_new2]
@@ -346,7 +357,7 @@ class TestProperties:
         assert op.work_wires == Wires(("extra"))
 
 
-class TestMiscMethods:
+class TestControlledMiscMethods:
     """Test miscellaneous minor Controlled methods."""
 
     def test_repr(self):
@@ -416,7 +427,7 @@ class TestMiscMethods:
 
     def test_label_matrix_param(self):
         """Test that the label method simply returns the label of the base and updates the cache."""
-        U = np.eye(2)
+        U = pnp.eye(2)
         base = qml.QubitUnitary(U, wires=0)
         op = Controlled(base, ["a", "b"])
 
@@ -430,10 +441,10 @@ class TestMiscMethods:
         op = Controlled(base, (2, 3))
 
         mat = op.matrix()
-        mat_eigvals = np.sort(qml.math.linalg.eigvals(mat))
+        mat_eigvals = pnp.sort(qml.math.linalg.eigvals(mat))
 
         eigs = op.eigvals()
-        sort_eigs = np.sort(eigs)
+        sort_eigs = pnp.sort(eigs)
 
         assert qml.math.allclose(mat_eigvals, sort_eigs)
 
@@ -524,7 +535,7 @@ class TestMiscMethods:
         assert op7.hash != op1.hash
 
 
-class TestOperationProperties:
+class TestControlledOperationProperties:
     """Test ControlledOp specific properties."""
 
     # pylint:disable=no-member
@@ -590,7 +601,7 @@ class TestOperationProperties:
             op.parameter_frequencies
 
 
-class TestSimplify:
+class TestControlledSimplify:
     """Test qml.sum simplify method and depth property."""
 
     def test_depth_property(self):
@@ -634,7 +645,7 @@ class TestSimplify:
         assert simplified_op.arithmetic_depth == final_op.arithmetic_depth
 
 
-class TestQueuing:
+class TestControlledQueuing:
     """Test that Controlled operators queue and update base metadata."""
 
     def test_queuing(self):
@@ -659,27 +670,31 @@ class TestQueuing:
 
 base_num_control_mats = [
     (qml.PauliX("a"), 1, CNOT),
+    (qml.PauliX("a"), 2, Toffoli),
+    (qml.CNOT(["a", "b"]), 1, Toffoli),
+    (qml.PauliY("a"), 1, CY),
     (qml.PauliZ("a"), 1, CZ),
+    (qml.PauliZ("a"), 2, CCZ),
     (qml.SWAP(("a", "b")), 1, CSWAP),
-    (qml.PauliX("a"), 2, Toffoli),
+    (qml.Hadamard("a"), 1, CH),
     (qml.RX(1.234, "b"), 1, CRotx(1.234)),
     (qml.RY(-0.432, "a"), 1, CRoty(-0.432)),
     (qml.RZ(6.78, "a"), 1, CRotz(6.78)),
     (qml.Rot(1.234, -0.432, 9.0, "a"), 1, CRot3(1.234, -0.432, 9.0)),
+    (qml.PhaseShift(1.234, wires="a"), 1, ControlledPhaseShift(1.234)),
 ]
 
 
 class TestMatrix:
     """Tests of Controlled.matrix and Controlled.sparse_matrix"""
 
-    def test_correct_matrix_dimenions_with_batching(self):
+    def test_correct_matrix_dimensions_with_batching(self):
         """Test batching returns a matrix of the correct dimensions"""
-        x = np.array([1.0, 2.0, 3.0])
+
+        x = pnp.array([1.0, 2.0, 3.0])
         base = qml.RX(x, 0)
         op = Controlled(base, 1)
-
         matrix = op.matrix()
-
         assert matrix.shape == (3, 4, 4)
 
     @pytest.mark.parametrize("base, num_control, mat", base_num_control_mats)
@@ -781,63 +796,228 @@ class TestMatrix:
         assert isinstance(lil_mat, sparse.lil_matrix)
 
 
-class TestHelperMethod:
-    """Unittests for the _decompose_no_control_values helper function."""
+special_non_par_op_decomps = [
+    (qml.PauliY, [], [0], [1], qml.CY, [qml.CRY(np.pi, wires=[1, 0]), qml.S(1)]),
+    (qml.PauliZ, [], [1], [0], qml.CZ, [qml.ControlledPhaseShift(np.pi, wires=[0, 1])]),
+    (
+        qml.Hadamard,
+        [],
+        [1],
+        [0],
+        qml.CH,
+        [qml.RY(-np.pi / 4, wires=1), qml.CZ(wires=[0, 1]), qml.RY(np.pi / 4, wires=1)],
+    ),
+    (
+        qml.PauliZ,
+        [],
+        [0],
+        [2, 1],
+        qml.CCZ,
+        [
+            qml.CNOT(wires=[1, 0]),
+            qml.adjoint(qml.T(wires=0)),
+            qml.CNOT(wires=[2, 0]),
+            qml.T(wires=0),
+            qml.CNOT(wires=[1, 0]),
+            qml.adjoint(qml.T(wires=0)),
+            qml.CNOT(wires=[2, 0]),
+            qml.T(wires=0),
+            qml.T(wires=1),
+            qml.CNOT(wires=[2, 1]),
+            qml.Hadamard(wires=0),
+            qml.T(wires=2),
+            qml.adjoint(qml.T(wires=1)),
+            qml.CNOT(wires=[2, 1]),
+            qml.Hadamard(wires=0),
+        ],
+    ),
+    (
+        qml.CZ,
+        [],
+        [1, 2],
+        [0],
+        qml.CCZ,
+        [
+            qml.CNOT(wires=[1, 2]),
+            qml.adjoint(qml.T(wires=2)),
+            qml.CNOT(wires=[0, 2]),
+            qml.T(wires=2),
+            qml.CNOT(wires=[1, 2]),
+            qml.adjoint(qml.T(wires=2)),
+            qml.CNOT(wires=[0, 2]),
+            qml.T(wires=2),
+            qml.T(wires=1),
+            qml.CNOT(wires=[0, 1]),
+            qml.Hadamard(wires=2),
+            qml.T(wires=0),
+            qml.adjoint(qml.T(wires=1)),
+            qml.CNOT(wires=[0, 1]),
+            qml.Hadamard(wires=[2]),
+        ],
+    ),
+    (
+        qml.SWAP,
+        [],
+        [1, 2],
+        [0],
+        qml.CSWAP,
+        [qml.Toffoli(wires=[0, 2, 1]), qml.Toffoli(wires=[0, 1, 2]), qml.Toffoli(wires=[0, 2, 1])],
+    ),
+]
 
-    def test_crx(self):
-        """Test case with single control wire and defined _controlled"""
-        base = qml.RX(1.0, wires=0)
-        op = Controlled(base, 1)
-        decomp = _decompose_no_control_values(op)
-        assert len(decomp) == 1
-        assert qml.equal(decomp[0], qml.CRX(1.0, wires=(1, 0)))
-
-    def test_toffoli(self):
-        """Test case when PauliX with two controls."""
-        op = Controlled(qml.PauliX("c"), ("a", 2))
-        decomp = _decompose_no_control_values(op)
-        assert len(decomp) == 1
-        assert equal_list(decomp, qml.MultiControlledX(wires=("a", 2, "c")))
-
-    def test_multicontrolledx(self):
-        """Test case when PauliX has many controls."""
-        op = Controlled(qml.PauliX(4), (0, 1, 2, 3))
-        decomp = _decompose_no_control_values(op)
-        assert len(decomp) == 1
-        assert qml.equal(decomp[0], qml.MultiControlledX(wires=(0, 1, 2, 3, 4)))
-
-    def test_decomposes_target(self):
-        """Test that we decompose the target if we don't have a special case."""
-        target = qml.IsingXX(1.0, wires=(0, 1))
-        op = Controlled(target, (3, 4))
-
-        decomp = _decompose_no_control_values(op)
-        assert len(decomp) == 3
-
-        target_decomp = target.expand().circuit
-        for op1, target in zip(decomp, target_decomp):
-            assert isinstance(op1, Controlled)
-            assert op1.control_wires == (3, 4)
-
-            assert qml.equal(op1.base, target)
-
-    def test_None_default(self):
-        """Test that helper returns None if no special decomposition."""
-        op = Controlled(TempOperator(0), (1, 2))
-        assert _decompose_no_control_values(op) is None
+special_par_op_decomps = [
+    (
+        qml.RX,
+        [0.123],
+        [1],
+        [0],
+        qml.CRX,
+        [
+            qml.RZ(np.pi / 2, wires=1),
+            qml.RY(0.123 / 2, wires=1),
+            qml.CNOT(wires=[0, 1]),
+            qml.RY(-0.123 / 2, wires=1),
+            qml.CNOT(wires=[0, 1]),
+            qml.RZ(-np.pi / 2, wires=1),
+        ],
+    ),
+    (
+        qml.RY,
+        [0.123],
+        [1],
+        [0],
+        qml.CRY,
+        [
+            qml.RY(0.123 / 2, 1),
+            qml.CNOT(wires=(0, 1)),
+            qml.RY(-0.123 / 2, 1),
+            qml.CNOT(wires=(0, 1)),
+        ],
+    ),
+    (
+        qml.RZ,
+        [0.123],
+        [0],
+        [1],
+        qml.CRZ,
+        [
+            qml.PhaseShift(0.123 / 2, wires=0),
+            qml.CNOT(wires=[1, 0]),
+            qml.PhaseShift(-0.123 / 2, wires=0),
+            qml.CNOT(wires=[1, 0]),
+        ],
+    ),
+    (
+        qml.Rot,
+        [0.1, 0.2, 0.3],
+        [1],
+        [0],
+        qml.CRot,
+        [
+            qml.RZ((0.1 - 0.3) / 2, wires=1),
+            qml.CNOT(wires=[0, 1]),
+            qml.RZ(-(0.1 + 0.3) / 2, wires=1),
+            qml.RY(-0.2 / 2, wires=1),
+            qml.CNOT(wires=[0, 1]),
+            qml.RY(0.2 / 2, wires=1),
+            qml.RZ(0.3, wires=1),
+        ],
+    ),
+    (
+        qml.PhaseShift,
+        [0.123],
+        [1],
+        [0],
+        qml.ControlledPhaseShift,
+        [
+            qml.PhaseShift(0.123 / 2, wires=0),
+            qml.CNOT(wires=[0, 1]),
+            qml.PhaseShift(-0.123 / 2, wires=1),
+            qml.CNOT(wires=[0, 1]),
+            qml.PhaseShift(0.123 / 2, wires=1),
+        ],
+    ),
+]
+
+custom_ctrl_op_decomps = special_non_par_op_decomps + special_par_op_decomps
+
+pauli_x_based_op_decomps = [
+    (qml.PauliX, [0], [1], [qml.CNOT([1, 0])]),
+    (
+        qml.PauliX,
+        [2],
+        [0, 1],
+        qml.Toffoli.compute_decomposition(wires=[0, 1, 2]),
+    ),
+    (
+        qml.CNOT,
+        [1, 2],
+        [0],
+        qml.Toffoli.compute_decomposition(wires=[0, 1, 2]),
+    ),
+    (
+        qml.PauliX,
+        [3],
+        [0, 1, 2],
+        qml.MultiControlledX.compute_decomposition(wires=[0, 1, 2, 3], work_wires=Wires("aux")),
+    ),
+    (
+        qml.CNOT,
+        [2, 3],
+        [0, 1],
+        qml.MultiControlledX.compute_decomposition(wires=[0, 1, 2, 3], work_wires=Wires("aux")),
+    ),
+    (
+        qml.Toffoli,
+        [1, 2, 3],
+        [0],
+        qml.MultiControlledX.compute_decomposition(wires=[0, 1, 2, 3], work_wires=Wires("aux")),
+    ),
+]
+
+
+class TestDecomposition:
+    """Test decomposition of Controlled."""
+
+    @pytest.mark.parametrize(
+        "target, decomp",
+        [
+            (
+                OpWithDecomposition(0.123, wires=[0, 1]),
+                [
+                    qml.CH(wires=[2, 0]),
+                    Controlled(qml.S(wires=1), control_wires=2),
+                    qml.CRX(0.123, wires=[2, 0]),
+                ],
+            ),
+            (
+                qml.IsingXX(0.123, wires=[0, 1]),
+                [
+                    qml.Toffoli(wires=[2, 0, 1]),
+                    qml.CRX(0.123, wires=[2, 0]),
+                    qml.Toffoli(wires=[2, 0, 1]),
+                ],
+            ),
+        ],
+    )
+    def test_decomposition(self, target, decomp):
+        """Test that we decompose a normal controlled operation"""
+        op = Controlled(target, 2)
+        assert op.decomposition() == decomp
 
     def test_non_differentiable_one_qubit_special_unitary(self):
-        """Assert that a non differentiable on qubit special unitary uses the bisect decomposition."""
+        """Assert that a non-differentiable on qubit special unitary uses the bisect decomposition."""
+
         op = qml.ctrl(qml.RZ(1.2, wires=0), (1, 2, 3, 4))
         decomp = op.decomposition()
 
-        assert qml.equal(decomp[0], qml.MultiControlledX(wires=(1, 2, 0), work_wires=(3, 4)))
+        assert qml.equal(decomp[0], qml.Toffoli(wires=(1, 2, 0)))
         assert isinstance(decomp[1], qml.QubitUnitary)
-        assert qml.equal(decomp[2], qml.MultiControlledX(wires=(3, 4, 0), work_wires=(1, 2)))
+        assert qml.equal(decomp[2], qml.Toffoli(wires=(3, 4, 0)))
         assert isinstance(decomp[3].base, qml.QubitUnitary)
-        assert qml.equal(decomp[4], qml.MultiControlledX(wires=(1, 2, 0), work_wires=(3, 4)))
+        assert qml.equal(decomp[4], qml.Toffoli(wires=(1, 2, 0)))
         assert isinstance(decomp[5], qml.QubitUnitary)
-        assert qml.equal(decomp[6], qml.MultiControlledX(wires=(3, 4, 0), work_wires=(1, 2)))
+        assert qml.equal(decomp[6], qml.Toffoli(wires=(3, 4, 0)))
         assert isinstance(decomp[7].base, qml.QubitUnitary)
 
         decomp_mat = qml.matrix(op.decomposition, wire_order=op.wires)()
@@ -858,6 +1038,93 @@ class TestHelperMethod:
         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,
+    )
+    def test_decomposition_custom_ops(
+        self,
+        base_cls,
+        params,
+        base_wires,
+        ctrl_wires,
+        custom_ctrl_cls,
+        expected,
+        tol,
+    ):
+        """Tests decompositions of custom operations"""
+
+        active_wires = ctrl_wires + base_wires
+        base_op = base_cls(*params, wires=base_wires)
+        ctrl_op = Controlled(base_op, control_wires=ctrl_wires)
+        custom_ctrl_op = custom_ctrl_cls(*params, active_wires)
+
+        assert ctrl_op.decomposition() == expected
+        assert ctrl_op.expand().circuit == expected
+        assert custom_ctrl_op.decomposition() == expected
+        assert custom_ctrl_cls.compute_decomposition(*params, active_wires) == expected
+
+        mat = qml.matrix(ctrl_op.decomposition, wire_order=active_wires)()
+        assert np.allclose(mat, custom_ctrl_op.matrix(), atol=tol, rtol=0)
+
+    @pytest.mark.parametrize(
+        "base_cls, params, base_wires, ctrl_wires, custom_ctrl_cls, expected",
+        special_par_op_decomps,
+    )
+    def test_decomposition_custom_par_ops_broadcasted(
+        self,
+        base_cls,
+        params,
+        base_wires,
+        ctrl_wires,
+        custom_ctrl_cls,
+        expected,
+        tol,
+    ):
+        """Tests broadcasted decompositions of custom controlled ops"""
+        broad_casted_params = [np.array([p, p]) for p in params]
+        self.test_decomposition_custom_ops(
+            base_cls,
+            broad_casted_params,
+            base_wires,
+            ctrl_wires,
+            custom_ctrl_cls,
+            expected,
+            tol,
+        )
+
+    @pytest.mark.parametrize(
+        "base_cls, base_wires, ctrl_wires, expected",
+        pauli_x_based_op_decomps,
+    )
+    def test_decomposition_pauli_x(self, base_cls, base_wires, ctrl_wires, expected):
+        """Tests decompositions where the base is PauliX"""
+
+        base_op = base_cls(wires=base_wires)
+        ctrl_op = Controlled(base_op, control_wires=ctrl_wires, work_wires=Wires("aux"))
+
+        assert ctrl_op.decomposition() == expected
+        assert ctrl_op.expand().circuit == expected
+
+    def test_decomposition_nested(self):
+        """Tests decompositions of nested controlled operations"""
+
+        ctrl_op = Controlled(Controlled(qml.RZ(0.123, wires=0), control_wires=1), control_wires=2)
+        expected = [
+            qml.ControlledPhaseShift(0.123 / 2, wires=[2, 0]),
+            qml.Toffoli(wires=[2, 1, 0]),
+            qml.ControlledPhaseShift(-0.123 / 2, wires=[2, 0]),
+            qml.Toffoli(wires=[2, 1, 0]),
+        ]
+        assert ctrl_op.decomposition() == expected
+        assert ctrl_op.expand().circuit == expected
+
+    def test_decomposition_undefined(self):
+        """Tests error raised when decomposition is undefined"""
+        op = Controlled(TempOperator(0), (1, 2))
+        with pytest.raises(DecompositionUndefinedError):
+            op.decomposition()
+
     def test_global_phase_decomp_raises_warning(self):
         """Test that ctrl(GlobalPhase).decomposition() raises a warning."""
         op = qml.ctrl(qml.GlobalPhase(1.23), control=[0])
@@ -866,12 +1133,7 @@ class TestHelperMethod:
         ):
             assert op.decomposition() == []
 
-
-@pytest.mark.parametrize("test_expand", (False, True))
-class TestDecomposition:
-    """Test controlled's decomposition method."""
-
-    def test_control_values_no_special_decomp(self, test_expand):
+    def test_control_on_zero(self):
         """Test decomposition applies PauliX gates to flip any control-on-zero wires."""
 
         control_wires = (0, 1, 2)
@@ -880,57 +1142,45 @@ class TestDecomposition:
         base = TempOperator("a")
         op = Controlled(base, control_wires, control_values)
 
-        decomp = op.expand().circuit if test_expand else op.decomposition()
-
-        assert qml.equal(decomp[0], qml.PauliX(1))
-        assert qml.equal(decomp[1], qml.PauliX(2))
+        decomp1 = op.decomposition()
+        decomp2 = op.expand().circuit
 
-        assert isinstance(decomp[2], Controlled)
-        assert decomp[2].control_values == [True, True, True]
+        for decomp in [decomp1, decomp2]:
+            assert qml.equal(decomp[0], qml.PauliX(1))
+            assert qml.equal(decomp[1], qml.PauliX(2))
 
-        assert qml.equal(decomp[3], qml.PauliX(1))
-        assert qml.equal(decomp[4], qml.PauliX(2))
+            assert isinstance(decomp[2], Controlled)
+            assert decomp[2].control_values == [True, True, True]
 
-    def test_control_values_special_decomp(self, test_expand):
-        """Test decomposition when needs control_values flips and special decomp exists."""
+            assert qml.equal(decomp[3], qml.PauliX(1))
+            assert qml.equal(decomp[4], qml.PauliX(2))
 
-        base = qml.PauliX(2)
-        op = Controlled(base, (0, 1), (True, False))
-
-        decomp = op.expand().circuit if test_expand else op.decomposition()
-        expected = [qml.PauliX(1), qml.MultiControlledX(wires=(0, 1, 2)), qml.PauliX(1)]
-        assert equal_list(decomp, expected)
+    @pytest.mark.parametrize(
+        "base_cls, params, base_wires, ctrl_wires, _, expected",
+        custom_ctrl_op_decomps,
+    )
+    def test_control_on_zero_custom_ops(
+        self, base_cls, params, base_wires, ctrl_wires, _, expected
+    ):
+        """Tests that custom ops are not converted when wires are control-on-zero."""
 
-    def test_no_control_values_special_decomp(self, test_expand):
-        """Test a case with no control values but a special decomposition."""
-        base = qml.RX(1.0, 2)
-        op = Controlled(base, 1)
-        decomp = op.expand().circuit if test_expand else op.decomposition()
-        assert len(decomp) == 1
-        assert qml.equal(decomp[0], qml.CRX(1.0, (1, 2)))
-
-    def test_no_control_values_target_decomposition(self, test_expand):
-        """Tests a case with no control values and no special decomposition but
-        the ability to decompose the target."""
-        base = qml.IsingXX(1.23, wires=(0, 1))
-        op = Controlled(base, "a")
+        base_op = base_cls(*params, wires=base_wires)
+        op = Controlled(base_op, control_wires=ctrl_wires, control_values=[False] * len(ctrl_wires))
 
-        decomp = op.expand().circuit if test_expand else op.decomposition()
-        base_decomp = base.decomposition()
-        for cop, base_op in zip(decomp, base_decomp):
-            assert isinstance(cop, Controlled)
-            assert qml.equal(cop.base, base_op)
+        decomp = op.decomposition()
 
-    def test_no_control_values_no_special_decomp(self, test_expand):
-        """Test if all control_values are true and no special decomposition exists,
-        the method raises a DecompositionUndefinedError."""
+        i = 0
+        for ctrl_wire in ctrl_wires:
+            assert decomp[i] == qml.PauliX(wires=ctrl_wire)
+            i += 1
 
-        base = TempOperator("a")
-        op = Controlled(base, (0, 1, 2))
+        for exp in expected:
+            assert decomp[i] == exp
+            i += 1
 
-        with pytest.raises(DecompositionUndefinedError):
-            # pylint: disable=unused-variable
-            decomp = op.expand().circuit if test_expand else op.decomposition()
+        for ctrl_wire in ctrl_wires:
+            assert decomp[i] == qml.PauliX(wires=ctrl_wire)
+            i += 1
 
 
 class TestArithmetic:
@@ -991,7 +1241,7 @@ class TestDifferentiation:
         """Test differentiation using autograd"""
 
         dev = qml.device("default.qubit", wires=2)
-        init_state = np.array([1.0, -1.0], requires_grad=False) / np.sqrt(2)
+        init_state = pnp.array([1.0, -1.0], requires_grad=False) / pnp.sqrt(2)
 
         @qml.qnode(dev, diff_method=diff_method)
         def circuit(b):
@@ -999,11 +1249,11 @@ class TestDifferentiation:
             Controlled(qml.RY(b, wires=1), control_wires=0)
             return qml.expval(qml.PauliX(0))
 
-        b = np.array(0.123, requires_grad=True)
+        b = pnp.array(0.123, requires_grad=True)
         res = qml.grad(circuit)(b)
-        expected = np.sin(b / 2) / 2
+        expected = pnp.sin(b / 2) / 2
 
-        assert np.allclose(res, expected)
+        assert pnp.allclose(res, expected)
 
     @pytest.mark.torch
     def test_torch(self, diff_method):
@@ -1013,7 +1263,7 @@ class TestDifferentiation:
         dev = qml.device("default.qubit", wires=2)
         init_state = torch.tensor(
             [1.0, -1.0], requires_grad=False, dtype=torch.complex128
-        ) / np.sqrt(2)
+        ) / pnp.sqrt(2)
 
         @qml.qnode(dev, diff_method=diff_method)
         def circuit(b):
@@ -1026,9 +1276,9 @@ class TestDifferentiation:
         loss.backward()  # pylint:disable=no-member
 
         res = b.grad.detach()
-        expected = np.sin(b.detach() / 2) / 2
+        expected = pnp.sin(b.detach() / 2) / 2
 
-        assert np.allclose(res, expected)
+        assert pnp.allclose(res, expected)
 
     @pytest.mark.jax
     @pytest.mark.parametrize("jax_interface", ["auto", "jax", "jax-python"])
@@ -1045,16 +1295,16 @@ class TestDifferentiation:
 
         @qml.qnode(dev, diff_method=diff_method, interface=jax_interface)
         def circuit(b):
-            init_state = onp.array([1.0, -1.0]) / np.sqrt(2)
+            init_state = np.array([1.0, -1.0]) / pnp.sqrt(2)
             qml.StatePrep(init_state, wires=0)
             Controlled(qml.RY(b, wires=1), control_wires=0)
             return qml.expval(qml.PauliX(0))
 
         b = jnp.array(0.123)
         res = jax.grad(circuit)(b)
-        expected = np.sin(b / 2) / 2
+        expected = pnp.sin(b / 2) / 2
 
-        assert np.allclose(res, expected)
+        assert pnp.allclose(res, expected)
 
     @pytest.mark.tf
     def test_tf(self, diff_method):
@@ -1062,7 +1312,7 @@ class TestDifferentiation:
         import tensorflow as tf
 
         dev = qml.device("default.qubit", wires=2)
-        init_state = tf.constant([1.0, -1.0], dtype=tf.complex128) / np.sqrt(2)
+        init_state = tf.constant([1.0, -1.0], dtype=tf.complex128) / pnp.sqrt(2)
 
         @qml.qnode(dev, diff_method=diff_method)
         def circuit(b):
@@ -1076,9 +1326,9 @@ class TestDifferentiation:
             loss = circuit(b)
 
         res = tape.gradient(loss, b)
-        expected = np.sin(b / 2) / 2
+        expected = pnp.sin(b / 2) / 2
 
-        assert np.allclose(res, expected)
+        assert pnp.allclose(res, expected)
 
 
 class TestControlledSupportsBroadcasting:
@@ -1145,7 +1395,7 @@ class TestControlledSupportsBroadcasting:
         """Test that a Controlled operation whose base is a single-scalar-parameter operations
         on a single wire marked as supporting parameter broadcasting actually do support broadcasting.
         """
-        par = np.array([0.25, 2.1, -0.42])
+        par = pnp.array([0.25, 2.1, -0.42])
         wires = ["wire0"]
 
         cls = getattr(qml, name)
@@ -1162,7 +1412,7 @@ class TestControlledSupportsBroadcasting:
         """Test that a Controlled operation whose base is a single-scalar-parameter operations
         on multiple wires marked as supporting parameter broadcasting actually do support broadcasting.
         """
-        par = np.array([0.25, 2.1, -0.42])
+        par = pnp.array([0.25, 2.1, -0.42])
         cls = getattr(qml, name)
 
         # Provide up to 6 wires and take as many as the class requires
@@ -1181,7 +1431,7 @@ class TestControlledSupportsBroadcasting:
         """Test that a Controlled operation whose base is a two-scalar-parameter operations
         on a single wire marked as supporting parameter broadcasting actually do support broadcasting.
         """
-        par = (np.array([0.25, 2.1, -0.42]), np.array([-6.2, 0.12, 0.421]))
+        par = (pnp.array([0.25, 2.1, -0.42]), pnp.array([-6.2, 0.12, 0.421]))
         wires = ["wire0"]
 
         cls = getattr(qml, name)
@@ -1200,9 +1450,9 @@ class TestControlledSupportsBroadcasting:
         on a single wire marked as supporting parameter broadcasting actually do support broadcasting.
         """
         par = (
-            np.array([0.25, 2.1, -0.42]),
-            np.array([-6.2, 0.12, 0.421]),
-            np.array([0.2, 1.1, -5.2]),
+            pnp.array([0.25, 2.1, -0.42]),
+            pnp.array([-6.2, 0.12, 0.421]),
+            pnp.array([0.2, 1.1, -5.2]),
         )
         wires = ["wire0"]
 
@@ -1222,9 +1472,9 @@ class TestControlledSupportsBroadcasting:
         on multiple wires marked as supporting parameter broadcasting actually do support broadcasting.
         """
         par = (
-            np.array([0.25, 2.1, -0.42]),
-            np.array([-6.2, 0.12, 0.421]),
-            np.array([0.2, 1.1, -5.2]),
+            pnp.array([0.25, 2.1, -0.42]),
+            pnp.array([-6.2, 0.12, 0.421]),
+            pnp.array([0.2, 1.1, -5.2]),
         )
         wires = ["wire0", 214]
 
@@ -1241,7 +1491,7 @@ class TestControlledSupportsBroadcasting:
     def test_controlled_diagonal_qubit_unitary(self):
         """Test that a Controlled operation whose base is a DiagonalQubitUnitary, which is marked
         as supporting parameter broadcasting, actually does support broadcasting."""
-        diag = np.array([[1j, 1, 1, -1j], [-1j, 1j, 1, -1], [1j, -1j, 1.0, -1]])
+        diag = pnp.array([[1j, 1, 1, -1j], [-1j, 1j, 1, -1], [1j, -1j, 1.0, -1]])
         wires = ["a", 5]
 
         base = qml.DiagonalQubitUnitary(diag, wires=wires)
@@ -1260,7 +1510,7 @@ class TestControlledSupportsBroadcasting:
     def test_controlled_pauli_rot(self, pauli_word, wires):
         """Test that a Controlled operation whose base is PauliRot, which is marked as supporting
         parameter broadcasting, actually does support broadcasting."""
-        par = np.array([0.25, 2.1, -0.42])
+        par = pnp.array([0.25, 2.1, -0.42])
 
         base = qml.PauliRot(par, pauli_word, wires=wires)
         op = Controlled(base, "wire1")
@@ -1276,7 +1526,7 @@ class TestControlledSupportsBroadcasting:
     def test_controlled_multi_rz(self, wires):
         """Test that a Controlled operation whose base is MultiRZ, which is marked as supporting
         parameter broadcasting, actually does support broadcasting."""
-        par = np.array([0.25, 2.1, -0.42])
+        par = pnp.array([0.25, 2.1, -0.42])
 
         base = qml.MultiRZ(par, wires=wires)
         op = Controlled(base, "wire1")
@@ -1288,13 +1538,13 @@ class TestControlledSupportsBroadcasting:
 
     @pytest.mark.parametrize(
         "state_, num_wires",
-        [([1.0, 0.0], 1), ([0.5, -0.5j, 0.5, -0.5], 2), (np.ones(8) / np.sqrt(8), 3)],
+        [([1.0, 0.0], 1), ([0.5, -0.5j, 0.5, -0.5], 2), (pnp.ones(8) / pnp.sqrt(8), 3)],
     )
     def test_controlled_qubit_state_vector(self, state_, num_wires):
         """Test that StatePrep, which is marked as supporting parameter broadcasting,
         actually does support broadcasting."""
 
-        state = np.array([state_])
+        state = pnp.array([state_])
         base = qml.StatePrep(state, wires=list(range(num_wires)))
         op = Controlled(base, "wire1")
 
@@ -1302,7 +1552,7 @@ class TestControlledSupportsBroadcasting:
         qml.StatePrep.compute_decomposition(state, list(range(num_wires)))
         op.decomposition()
 
-        state = np.array([state_] * 3)
+        state = pnp.array([state_] * 3)
         base = qml.StatePrep(state, wires=list(range(num_wires)))
         op = Controlled(base, "wire1")
         assert op.batch_size == 3
@@ -1311,20 +1561,20 @@ class TestControlledSupportsBroadcasting:
 
     @pytest.mark.parametrize(
         "state, num_wires",
-        [([1.0, 0.0], 1), ([0.5, -0.5j, 0.5, -0.5], 2), (np.ones(8) / np.sqrt(8), 3)],
+        [([1.0, 0.0], 1), ([0.5, -0.5j, 0.5, -0.5], 2), (pnp.ones(8) / pnp.sqrt(8), 3)],
     )
     def test_controlled_amplitude_embedding(self, state, num_wires):
         """Test that AmplitudeEmbedding, which is marked as supporting parameter broadcasting,
         actually does support broadcasting."""
 
-        features = np.array([state])
+        features = pnp.array([state])
         base = qml.AmplitudeEmbedding(features, wires=list(range(num_wires)))
         op = Controlled(base, "wire1")
         assert op.batch_size == 1
         qml.AmplitudeEmbedding.compute_decomposition(features, list(range(num_wires)))
         op.decomposition()
 
-        features = np.array([state] * 3)
+        features = pnp.array([state] * 3)
         base = qml.AmplitudeEmbedding(features, wires=list(range(num_wires)))
         op = Controlled(base, "wire1")
         assert op.batch_size == 3
@@ -1334,9 +1584,9 @@ class TestControlledSupportsBroadcasting:
     @pytest.mark.parametrize(
         "angles, num_wires",
         [
-            (np.array([[0.5], [2.1]]), 1),
-            (np.array([[0.5, -0.5], [0.2, 1.5]]), 2),
-            (np.ones((2, 5)), 5),
+            (pnp.array([[0.5], [2.1]]), 1),
+            (pnp.array([[0.5, -0.5], [0.2, 1.5]]), 2),
+            (pnp.ones((2, 5)), 5),
         ],
     )
     def test_controlled_angle_embedding(self, angles, num_wires):
@@ -1352,9 +1602,9 @@ class TestControlledSupportsBroadcasting:
     @pytest.mark.parametrize(
         "features, num_wires",
         [
-            (np.array([[0.5], [2.1]]), 1),
-            (np.array([[0.5, -0.5], [0.2, 1.5]]), 2),
-            (np.ones((2, 5)), 5),
+            (pnp.array([[0.5], [2.1]]), 1),
+            (pnp.array([[0.5, -0.5], [0.2, 1.5]]), 2),
+            (pnp.ones((2, 5)), 5),
         ],
     )
     def test_controlled_iqp_embedding(self, features, num_wires):
@@ -1375,9 +1625,9 @@ class TestControlledSupportsBroadcasting:
     @pytest.mark.parametrize(
         "features, weights, num_wires, batch_size",
         [
-            (np.array([[0.5], [2.1]]), np.array([[0.61], [0.3]]), 1, 2),
-            (np.array([[0.5, -0.5], [0.2, 1.5]]), np.ones((2, 4, 3)), 2, 2),
-            (np.array([0.5, -0.5, 0.2]), np.ones((3, 2, 6)), 3, 3),
+            (pnp.array([[0.5], [2.1]]), pnp.array([[0.61], [0.3]]), 1, 2),
+            (pnp.array([[0.5, -0.5], [0.2, 1.5]]), pnp.ones((2, 4, 3)), 2, 2),
+            (pnp.array([0.5, -0.5, 0.2]), pnp.ones((3, 2, 6)), 3, 3),
         ],
     )
     def test_controlled_qaoa_embedding(self, features, weights, num_wires, batch_size):
@@ -1393,213 +1643,268 @@ class TestControlledSupportsBroadcasting:
         op.decomposition()
 
 
-##### TESTS FOR THE ctrl TRANSFORM #####
+custom_ctrl_ops = [
+    (qml.PauliY(wires=0), [1], qml.CY(wires=[1, 0])),
+    (qml.PauliZ(wires=0), [1], qml.CZ(wires=[1, 0])),
+    (qml.RX(0.123, wires=0), [1], qml.CRX(0.123, wires=[1, 0])),
+    (qml.RY(0.123, wires=0), [1], qml.CRY(0.123, wires=[1, 0])),
+    (qml.RZ(0.123, wires=0), [1], qml.CRZ(0.123, wires=[1, 0])),
+    (
+        qml.Rot(0.123, 0.234, 0.456, wires=0),
+        [1],
+        qml.CRot(0.123, 0.234, 0.456, wires=[1, 0]),
+    ),
+    (qml.PhaseShift(0.123, wires=0), [1], qml.ControlledPhaseShift(0.123, wires=[1, 0])),
+    (
+        qml.QubitUnitary(np.array([[0, 1], [1, 0]]), wires=0),
+        [1, 2],
+        qml.ControlledQubitUnitary(np.array([[0, 1], [1, 0]]), wires=0, control_wires=[1, 2]),
+    ),
+]
 
 
-def test_invalid_input_error():
-    """Test that a ValueError is raised upon invalid inputs."""
-    err = r"The object 1 of type <class 'int'> is not an Operator or callable."
-    with pytest.raises(ValueError, match=err):
-        qml.ctrl(1, control=2)
+class TestCtrl:
+    """Tests for the ctrl transform."""
+
+    def test_invalid_input_error(self):
+        """Test that a ValueError is raised upon invalid inputs."""
+        with pytest.raises(ValueError, match=r"<class 'int'> is not an Operator or callable."):
+            qml.ctrl(1, control=2)
+
+    @pytest.mark.parametrize("op, ctrl_wires, expected_op", custom_ctrl_ops)
+    def test_custom_controlled_ops(self, op, ctrl_wires, expected_op):
+        """Tests custom controlled operations are handled correctly."""
+        assert qml.ctrl(op, control=ctrl_wires) == expected_op
+
+    @pytest.mark.parametrize("op, ctrl_wires, _", custom_ctrl_ops)
+    def test_custom_controlled_ops_ctrl_on_zero(self, op, ctrl_wires, _):
+        """Tests custom controlled ops with control on zero are handled correctly."""
+
+        if isinstance(op, qml.QubitUnitary):
+            pytest.skip("ControlledQubitUnitary can accept any control values.")
+
+        ctrl_values = [False] * len(ctrl_wires)
+
+        if isinstance(op, Controlled):
+            expected = Controlled(
+                op.base,
+                control_wires=ctrl_wires + op.control_wires,
+                control_values=ctrl_values + op.control_values,
+            )
+        else:
+            expected = Controlled(op, control_wires=ctrl_wires, control_values=ctrl_values)
+
+        assert qml.ctrl(op, control=ctrl_wires, control_values=ctrl_values) == expected
+
+    @pytest.mark.parametrize("op, ctrl_wires, _", custom_ctrl_ops)
+    def test_custom_controlled_ops_wrong_wires(self, op, ctrl_wires, _):
+        """Tests custom controlled ops with wrong number of wires are handled correctly."""
+
+        ctrl_wires = ctrl_wires + ["a", "b", "c"]
+
+        if isinstance(op, qml.PhaseShift):
+            # TODO: remove this special case once ControlledGlobalPhase is implemented.
+            pytest.skip("PhaseShift has its temporary custom logic.")
+        if isinstance(op, qml.QubitUnitary):
+            pytest.skip("ControlledQubitUnitary can accept any number of control wires.")
+        elif isinstance(op, Controlled):
+            expected = Controlled(
+                op.base,
+                control_wires=ctrl_wires + op.control_wires,
+            )
+        else:
+            expected = Controlled(op, control_wires=ctrl_wires)
+
+        assert qml.ctrl(op, control=ctrl_wires) == expected
+
+    def test_nested_controls(self):
+        """Tests that nested controls are flattened correctly."""
+
+        op = qml.ctrl(
+            Controlled(
+                Controlled(qml.S(wires=[0]), control_wires=[1]),
+                control_wires=[2],
+                control_values=[0],
+            ),
+            control=[3],
+        )
+        expected = Controlled(
+            qml.S(wires=[0]),
+            control_wires=[3, 2, 1],
+            control_values=[1, 0, 1],
+        )
+        assert op == expected
 
+    @pytest.mark.parametrize("op, ctrl_wires, ctrl_op", custom_ctrl_ops)
+    def test_nested_custom_controls(self, op, ctrl_wires, ctrl_op):
+        """Tests that nested controls of custom controlled ops are flattened correctly."""
 
-def test_ctrl_sanity_check():
-    """Test that control works on a very standard usecase."""
+        if isinstance(ctrl_op, qml.ControlledQubitUnitary):
+            pytest.skip("ControlledQubitUnitary has its own logic")
 
-    def make_ops():
-        qml.RX(0.123, wires=0)
-        qml.RY(0.456, wires=2)
-        qml.RX(0.789, wires=0)
-        qml.Rot(0.111, 0.222, 0.333, wires=2)
-        qml.PauliX(wires=2)
-        qml.PauliY(wires=4)
-        qml.PauliZ(wires=0)
+        if isinstance(ctrl_op, qml.ControlledPhaseShift):
+            # TODO: remove this special case once ControlledGlobalPhase is implemented
+            pytest.skip("ControlledPhaseShift has its own logic")
 
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        cmake_ops = ctrl(make_ops, control=1)
-        # Execute controlled version.
-        cmake_ops()
+        expected_base = op.base if isinstance(op, Controlled) else op
+        base_ctrl_wires = (
+            ctrl_wires + op.control_wires if isinstance(op, Controlled) else ctrl_wires
+        )
+        ctrl_values = [1] * len(ctrl_wires)
+        base_ctrl_values = (
+            ctrl_values + op.control_values if isinstance(op, Controlled) else ctrl_values
+        )
 
-    tape = QuantumScript.from_queue(q_tape)
-    expanded_tape = tape.expand()
+        op = qml.ctrl(
+            Controlled(
+                ctrl_op,
+                control_wires=["b"],
+                control_values=[0],
+            ),
+            control=["a"],
+        )
+        expected = Controlled(
+            expected_base,
+            control_wires=["a", "b"] + base_ctrl_wires,
+            control_values=[1, 0] + base_ctrl_values,
+        )
+        assert op == expected
+
+    def test_nested_ctrl_qubit_unitaries(self):
+        """Tests that nested controlled qubit unitaries are flattened correctly."""
+
+        op = qml.ctrl(
+            Controlled(
+                qml.ControlledQubitUnitary(
+                    np.array([[0, 1], [1, 0]]), control_wires=[1], wires=[0]
+                ),
+                control_wires=[2],
+                control_values=[0],
+            ),
+            control=[3],
+        )
+        expected = qml.ControlledQubitUnitary(
+            np.array([[0, 1], [1, 0]]), control_wires=[3, 2, 1], wires=[0], control_values=[1, 0, 1]
+        )
+        assert op == expected
 
-    expected = [
-        *qml.CRX(0.123, wires=[1, 0]).decomposition(),
-        *qml.CRY(0.456, wires=[1, 2]).decomposition(),
-        *qml.CRX(0.789, wires=[1, 0]).decomposition(),
-        *qml.CRot(0.111, 0.222, 0.333, wires=[1, 2]).decomposition(),
-        qml.CNOT(wires=[1, 2]),
-        *qml.CY(wires=[1, 4]).decomposition(),
-        *qml.CZ(wires=[1, 0]).decomposition(),
-    ]
-    assert len(tape.operations) == 7
-    for op1, op2 in zip(expanded_tape, expected):
-        assert qml.equal(op1, op2)
-
-
-def test_adjoint_of_ctrl():
-    """Test adjoint(ctrl(fn)) and ctrl(adjoint(fn))"""
-
-    def my_op(a, b, c):
-        qml.RX(a, wires=2)
-        qml.RY(b, wires=3)
-        qml.RZ(c, wires=0)
-
-    with qml.queuing.AnnotatedQueue() as q1:
-        cmy_op_dagger = qml.simplify(qml.adjoint(ctrl(my_op, 5)))
-        # Execute controlled and adjointed version of my_op.
-        cmy_op_dagger(0.789, 0.123, c=0.456)
-
-    tape1 = QuantumScript.from_queue(q1)
-    with qml.queuing.AnnotatedQueue() as q2:
-        cmy_op_dagger = qml.simplify(ctrl(qml.adjoint(my_op), 5))
-        # Execute adjointed and controlled version of my_op.
-        cmy_op_dagger(0.789, 0.123, c=0.456)
-
-    tape2 = QuantumScript.from_queue(q2)
-    expected = [
-        *qml.CRZ(4 * onp.pi - 0.456, wires=[5, 0]).decomposition(),
-        *qml.CRY(4 * onp.pi - 0.123, wires=[5, 3]).decomposition(),
-        *qml.CRX(4 * onp.pi - 0.789, wires=[5, 2]).decomposition(),
-    ]
-    for tape in [tape1.expand(depth=1), tape2.expand(depth=1)]:
-        for op1, op2 in zip(tape, expected):
-            assert qml.equal(op1, op2)
-
-
-def test_nested_ctrl():
-    """Test nested use of control"""
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        CCS = ctrl(ctrl(qml.S, 7), 3)
-        CCS(wires=0)
-    tape = QuantumScript.from_queue(q_tape)
-    assert len(tape.operations) == 1
-    op = tape.operations[0]
-    assert isinstance(op, Controlled)
-    new_tape = tape.expand(depth=2)
-    assert qml.equal(new_tape[0], Controlled(qml.ControlledPhaseShift(np.pi / 2, [7, 0]), [3]))
-
-
-def test_multi_ctrl():
-    """Test control with a list of wires."""
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        CCS = ctrl(qml.S, control=[3, 7])
-        CCS(wires=0)
-    tape = QuantumScript.from_queue(q_tape)
-    assert len(tape.operations) == 1
-    op = tape.operations[0]
-    assert isinstance(op, Controlled)
-    new_tape = tape.expand(depth=1)
-    assert qml.equal(new_tape[0], Controlled(qml.PhaseShift(np.pi / 2, 0), [3, 7]))
-
-
-def test_ctrl_with_qnode():
-    """Test ctrl works when in a qnode cotext."""
-    dev = qml.device("default.qubit", wires=3)
-
-    def my_ansatz(params):
-        qml.RY(params[0], wires=0)
-        qml.RY(params[1], wires=1)
-        qml.CNOT(wires=[0, 1])
-        qml.RX(params[2], wires=1)
-        qml.RX(params[3], wires=0)
-        qml.CNOT(wires=[1, 0])
-
-    def controlled_ansatz(params):
-        qml.CRY(params[0], wires=[2, 0])
-        qml.CRY(params[1], wires=[2, 1])
-        qml.Toffoli(wires=[2, 0, 1])
-        qml.CRX(params[2], wires=[2, 1])
-        qml.CRX(params[3], wires=[2, 0])
-        qml.Toffoli(wires=[2, 1, 0])
-
-    def circuit(ansatz, params):
-        qml.RX(np.pi / 4.0, wires=2)
-        ansatz(params)
-        return qml.state()
-
-    params = [0.123, 0.456, 0.789, 1.345]
-    circuit1 = qml.qnode(dev)(partial(circuit, ansatz=ctrl(my_ansatz, 2)))
-    circuit2 = qml.qnode(dev)(partial(circuit, ansatz=controlled_ansatz))
-    res1 = circuit1(params=params)
-    res2 = circuit2(params=params)
-    assert qml.math.allclose(res1, res2)
-
-
-def test_ctrl_within_ctrl():
-    """Test using ctrl on a method that uses ctrl."""
-
-    def ansatz(params):
-        qml.RX(params[0], wires=0)
-        ctrl(qml.PauliX, control=0)(wires=1)
-        qml.RX(params[1], wires=0)
-
-    controlled_ansatz = ctrl(ansatz, 2)
-
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        controlled_ansatz([0.123, 0.456])
-
-    tape = QuantumScript.from_queue(q_tape)
-    tape = tape.expand(2, stop_at=lambda op: not isinstance(op, Controlled))
-
-    expected = [
-        *qml.CRX(0.123, wires=[2, 0]).decomposition(),
-        qml.Toffoli(wires=[2, 0, 1]),
-        *qml.CRX(0.456, wires=[2, 0]).decomposition(),
-    ]
-    for op1, op2 in zip(tape, expected):
-        assert qml.equal(op1, op2)
-
-
-def test_diagonal_ctrl():
-    """Test ctrl on diagonal gates."""
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        ctrl(qml.DiagonalQubitUnitary, 1)(onp.array([-1.0, 1.0j]), wires=0)
-    tape = QuantumScript.from_queue(q_tape)
-    tape = tape.expand(3, stop_at=lambda op: not isinstance(op, Controlled))
-    assert qml.equal(
-        tape[0], qml.DiagonalQubitUnitary(onp.array([1.0, 1.0, -1.0, 1.0j]), wires=[1, 0])
+    @pytest.mark.parametrize(
+        "op, ctrl_wires, ctrl_values, expected_op",
+        [
+            (qml.PauliX(wires=[0]), [1], [1], qml.CNOT([1, 0])),
+            (
+                qml.PauliX(wires=[2]),
+                [0, 1],
+                [1, 1],
+                qml.Toffoli(wires=[0, 1, 2]),
+            ),
+            (
+                qml.CNOT(wires=[1, 2]),
+                [0],
+                [1],
+                qml.Toffoli(wires=[0, 1, 2]),
+            ),
+            (
+                qml.PauliX(wires=[0]),
+                [1],
+                [0],
+                qml.MultiControlledX(wires=[1, 0], control_values=[0], work_wires=["aux"]),
+            ),
+            (
+                qml.PauliX(wires=[2]),
+                [0, 1],
+                [1, 0],
+                qml.MultiControlledX(wires=[0, 1, 2], control_values=[1, 0], work_wires=["aux"]),
+            ),
+            (
+                qml.CNOT(wires=[1, 2]),
+                [0],
+                [0],
+                qml.MultiControlledX(wires=[0, 1, 2], control_values=[0, 1], work_wires=["aux"]),
+            ),
+            (
+                qml.PauliX(wires=[3]),
+                [0, 1, 2],
+                [1, 1, 1],
+                qml.MultiControlledX(wires=[0, 1, 2, 3], work_wires=Wires("aux")),
+            ),
+            (
+                qml.CNOT(wires=[2, 3]),
+                [0, 1],
+                [1, 1],
+                qml.MultiControlledX(wires=[0, 1, 2, 3], work_wires=Wires("aux")),
+            ),
+            (
+                qml.Toffoli(wires=[1, 2, 3]),
+                [0],
+                [1],
+                qml.MultiControlledX(wires=[0, 1, 2, 3], work_wires=Wires("aux")),
+            ),
+        ],
     )
+    def test_pauli_x_based_ctrl_ops(self, op, ctrl_wires, ctrl_values, expected_op):
+        """Tests that PauliX-based ops are handled correctly."""
+        op = qml.ctrl(op, control=ctrl_wires, control_values=ctrl_values, work_wires=["aux"])
+        assert op == expected_op
+
+    def test_nested_pauli_x_based_ctrl_ops(self):
+        """Tests that nested PauliX-based ops are handled correctly."""
+
+        op = qml.ctrl(
+            Controlled(
+                qml.CNOT(wires=[1, 0]),
+                control_wires=[2],
+                control_values=[0],
+            ),
+            control=[3],
+        )
+        expected = qml.MultiControlledX(wires=[3, 2, 1, 0], control_values=[1, 0, 1])
+        assert op == expected
 
+    def test_controlled_phase_shift_special_logic(self):
+        """Tests special logic for PhaseShift"""
 
-@pytest.mark.parametrize(
-    "M",
-    [
-        qml.PauliX.compute_matrix(),
-        qml.PauliY.compute_matrix(),
-        qml.PauliZ.compute_matrix(),
-        qml.Hadamard.compute_matrix(),
-        np.array(
-            [
-                [1 + 2j, -3 + 4j],
-                [3 + 4j, 1 - 2j],
-            ]
+        # TODO: remove this special case once ControlledGlobalPhase is implemented.
+        # Tests special logic that wraps the phase shift in a controlled version.
+        op1 = qml.ctrl(qml.PhaseShift(0.123, wires=2), control=[0, 1])
+        expected = Controlled(qml.ControlledPhaseShift(0.123, wires=[1, 2]), control_wires=[0])
+        assert op1 == expected
+
+        op2 = qml.ctrl(qml.ControlledPhaseShift(0.123, wires=[1, 2]), control=[0])
+        assert op2 == expected
+
+        # Tests that special logic is not triggered when the control value is False.
+        op = qml.ctrl(qml.PhaseShift(0.123, wires=2), control=[0, 1], control_values=[True, False])
+        expected = Controlled(
+            qml.PhaseShift(0.123, wires=[2]), control_wires=[0, 1], control_values=[True, False]
         )
-        * 30**-0.5,
-    ],
-)
-def test_qubit_unitary(M):
-    """Test ctrl on QubitUnitary and ControlledQubitUnitary"""
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        ctrl(qml.QubitUnitary, 1)(M, wires=0)
+        assert op == expected
+
+
+class _Rot(Operation):
+    """A rotation operation that is not an instance of Rot
+
+    Used to test the default behaviour of expanding tapes without custom handling
+    of custom controlled operators (bypass automatic simplification of controlled
+    Rot to CRot gates in decompositions).
 
-    tape = QuantumScript.from_queue(q_tape)
-    expected = qml.ControlledQubitUnitary(M, control_wires=1, wires=0)
-    assert equal_list(list(tape), expected)
+    """
 
-    # causes decomposition into more basic operators
-    tape = tape.expand(3, stop_at=lambda op: not isinstance(op, Controlled))
-    assert not equal_list(list(tape), expected)
+    @staticmethod
+    def compute_decomposition(*params, wires=None):
+        return qml.Rot.compute_decomposition(*params, wires=wires)
 
+    def decomposition(self):
+        return self.compute_decomposition(*self.parameters, wires=self.wires)
 
-@pytest.mark.parametrize(
-    "M",
+
+unitaries = (
     [
-        pytest.param(qml.PauliX.compute_matrix(), marks=pytest.mark.xfail),
-        pytest.param(qml.PauliY.compute_matrix(), marks=pytest.mark.xfail),
-        pytest.param(qml.PauliZ.compute_matrix(), marks=pytest.mark.xfail),
-        pytest.param(qml.Hadamard.compute_matrix(), marks=pytest.mark.xfail),
-        np.array(
+        qml.PauliX.compute_matrix(),
+        qml.PauliY.compute_matrix(),
+        qml.PauliZ.compute_matrix(),
+        qml.Hadamard.compute_matrix(),
+        pnp.array(
             [
                 [1 + 2j, -3 + 4j],
                 [3 + 4j, 1 - 2j],
@@ -1608,146 +1913,288 @@ def test_qubit_unitary(M):
         * 30**-0.5,
     ],
 )
-def test_controlledqubitunitary(M):
-    """Test ctrl on ControlledQubitUnitary."""
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        ctrl(qml.ControlledQubitUnitary, 1)(M, control_wires=2, wires=0)
 
-    tape = QuantumScript.from_queue(q_tape)
-    # will immediately decompose according to selected decomposition algorithm
-    tape = tape.expand(3, stop_at=lambda op: not isinstance(op, Controlled))
 
-    expected = qml.ControlledQubitUnitary(M, control_wires=[2, 1], wires=0).decomposition()
-    assert equal_list(list(tape), expected)
+class TestTapeExpansionWithControlled:
+    """Tests expansion of tapes containing Controlled operations"""
+
+    def test_ctrl_values_sanity_check(self):
+        """Test that control works with control values on a very standard usecase."""
+
+        def make_ops():
+            qml.RX(0.123, wires=0)
+            qml.RY(0.456, wires=2)
+            qml.RX(0.789, wires=0)
+            qml.Rot(0.111, 0.222, 0.333, wires=2)
+            qml.PauliX(wires=2)
+            qml.PauliY(wires=4)
+            qml.PauliZ(wires=0)
+
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            ctrl(make_ops, control=1, control_values=0)()
+
+        tape = QuantumScript.from_queue(q_tape)
+        expected = [
+            qml.PauliX(wires=1),
+            *qml.CRX(0.123, wires=[1, 0]).decomposition(),
+            *qml.CRY(0.456, wires=[1, 2]).decomposition(),
+            *qml.CRX(0.789, wires=[1, 0]).decomposition(),
+            *qml.CRot(0.111, 0.222, 0.333, wires=[1, 2]).decomposition(),
+            qml.CNOT(wires=[1, 2]),
+            *qml.CY(wires=[1, 4]).decomposition(),
+            *qml.CZ(wires=[1, 0]).decomposition(),
+            qml.PauliX(wires=1),
+        ]
+        assert len(tape) == 9
+        expanded = tape.expand(stop_at=lambda obj: not isinstance(obj, Controlled))
+        assert expanded.circuit == expected
 
+    @pytest.mark.parametrize(
+        "op",
+        [
+            qml.ctrl(qml.ctrl(_Rot, 7), 3),  # nested control
+            qml.ctrl(_Rot, [3, 7]),  # multi-wire control
+        ],
+    )
+    def test_nested_ctrl(self, op, tol):
+        """Tests that nested controlled ops are expanded correctly"""
 
-def test_no_control_defined():
-    """Test a custom operation with no control transform defined."""
-    # QFT has no control rule defined.
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        ctrl(qml.templates.QFT, 2)(wires=[0, 1])
-    tape = QuantumScript.from_queue(q_tape)
-    tape = tape.expand(depth=3, stop_at=lambda op: not isinstance(op, Controlled))
-    assert len(tape.operations) == 8
-    # Check that all operations are updated to their controlled version.
-    for op in tape.operations:
-        assert type(op) in {qml.ControlledPhaseShift, qml.Toffoli, qml.CRX, qml.CSWAP, qml.CH}
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            op(0.1, 0.2, 0.3, wires=0)
 
+        tape = QuantumScript.from_queue(q_tape)
+        assert tape.expand(depth=1).circuit == [
+            Controlled(qml.RZ(0.1, 0), control_wires=[3, 7]),
+            Controlled(qml.RY(0.2, 0), control_wires=[3, 7]),
+            Controlled(qml.RZ(0.3, 0), control_wires=[3, 7]),
+        ]
 
-def test_decomposition_defined():
-    """Test that a controlled gate that has no control transform defined,
-    and a decomposition transformed defined, still works correctly"""
+        # Tests that the decomposition of the nested controlled _Rot gate is ultimately
+        # equivalent to the decomposition of the controlled CRot
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            for op_ in qml.CRot.compute_decomposition(0.1, 0.2, 0.3, wires=[7, 0]):
+                qml.ctrl(op_, control=3)
+        tape_expected = QuantumScript.from_queue(q_tape)
 
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        ctrl(qml.CY, 0)(wires=[1, 2])
+        def stopping_condition(o):
+            return not isinstance(o, Controlled) or not o.has_decomposition
 
-    tape = QuantumScript.from_queue(q_tape)
-    tape = tape.expand()
+        actual = tape.expand(depth=10, stop_at=stopping_condition)
+        expected = tape_expected.expand(depth=10, stop_at=stopping_condition)
+        actual_mat = qml.matrix(actual, wire_order=[3, 7, 0])
+        expected_mat = qml.matrix(expected, wire_order=[3, 7, 0])
+        assert qml.math.allclose(actual_mat, expected_mat, atol=tol, rtol=0)
 
-    assert len(tape.operations) == 2
+    @pytest.mark.parametrize(
+        "op",
+        [
+            qml.ctrl(qml.ctrl(qml.S, 7), 3),  # nested control
+            qml.ctrl(qml.S, [3, 7]),  # multi-wire control
+        ],
+    )
+    def test_nested_ctrl_containing_phase_shift(self, op):
+        """Test that nested controlled ops are expanded correctly when phase shift is involved
 
-    assert tape.operations[0].name == "C(CRY)"
-    assert tape.operations[1].name == "C(S)"
+        The decomposition of S gate contains a PhaseShift. In the nested case, we do not want
+        to apply control to the expanded PhaseShift like how it is typically done for other
+        operations, because the decomposition of PhaseShift contains a GlobalPhase, the controlled
+        version of which we do not have handling for.
 
+        TODO: remove this special case once ControlledGlobalPhase is implemented.
 
-def test_ctrl_template():
-    """Test that a controlled template correctly expands
-    on a device that doesn't support it"""
+        """
 
-    weights = np.ones([3, 2])
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            op(wires=0)
 
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        ctrl(qml.templates.BasicEntanglerLayers, 0)(weights, wires=[1, 2])
+        tape = QuantumScript.from_queue(q_tape)
+        assert tape.expand(depth=1).circuit == [
+            Controlled(qml.ControlledPhaseShift(np.pi / 2, wires=[7, 0]), control_wires=[3])
+        ]
 
-    tape = QuantumScript.from_queue(q_tape)
-    tape = expand_tape(tape, depth=2)
-    assert len(tape) == 9
-    assert all(o.name in {"CRX", "Toffoli"} for o in tape.operations)
+        assert tape.expand(depth=2).circuit == [
+            qml.ControlledPhaseShift(np.pi / 4, wires=[3, 7]),
+            qml.Toffoli(wires=[3, 7, 0]),
+            qml.ControlledPhaseShift(-np.pi / 4, wires=[3, 0]),
+            qml.Toffoli(wires=[3, 7, 0]),
+            qml.ControlledPhaseShift(np.pi / 4, wires=[3, 0]),
+        ]
 
+    def test_adjoint_of_ctrl(self):
+        """Tests that adjoint(ctrl(fn)) and ctrl(adjoint(fn)) are equivalent"""
+
+        def my_op(a, b, c):
+            qml.RX(a, wires=2)
+            qml.RY(b, wires=3)
+            qml.RZ(c, wires=0)
+
+        with qml.queuing.AnnotatedQueue() as q1:
+            # Execute controlled and adjoint version of my_op.
+            cmy_op_dagger = qml.simplify(qml.adjoint(ctrl(my_op, 5)))
+            cmy_op_dagger(0.789, 0.123, c=0.456)
+        tape1 = QuantumScript.from_queue(q1)
+
+        with qml.queuing.AnnotatedQueue() as q2:
+            # Execute adjoint and controlled version of my_op.
+            cmy_op_dagger = qml.simplify(ctrl(qml.adjoint(my_op), 5))
+            cmy_op_dagger(0.789, 0.123, c=0.456)
+        tape2 = QuantumScript.from_queue(q2)
+
+        expected = [
+            *qml.CRZ(4 * np.pi - 0.456, wires=[5, 0]).decomposition(),
+            *qml.CRY(4 * np.pi - 0.123, wires=[5, 3]).decomposition(),
+            *qml.CRX(4 * np.pi - 0.789, wires=[5, 2]).decomposition(),
+        ]
+        assert tape1.expand(depth=1).circuit == expected
+        assert tape2.expand(depth=1).circuit == expected
+
+    def test_ctrl_with_qnode(self):
+        """Test ctrl works when in a qnode cotext."""
+        dev = qml.device("default.qubit", wires=3)
+
+        def my_ansatz(params):
+            qml.RY(params[0], wires=0)
+            qml.RY(params[1], wires=1)
+            qml.CNOT(wires=[0, 1])
+            qml.RX(params[2], wires=1)
+            qml.RX(params[3], wires=0)
+            qml.CNOT(wires=[1, 0])
+
+        def controlled_ansatz(params):
+            qml.CRY(params[0], wires=[2, 0])
+            qml.CRY(params[1], wires=[2, 1])
+            qml.Toffoli(wires=[2, 0, 1])
+            qml.CRX(params[2], wires=[2, 1])
+            qml.CRX(params[3], wires=[2, 0])
+            qml.Toffoli(wires=[2, 1, 0])
+
+        def circuit(ansatz, params):
+            qml.RX(pnp.pi / 4.0, wires=2)
+            ansatz(params)
+            return qml.state()
+
+        params = [0.123, 0.456, 0.789, 1.345]
+        circuit1 = qml.qnode(dev)(partial(circuit, ansatz=ctrl(my_ansatz, 2)))
+        circuit2 = qml.qnode(dev)(partial(circuit, ansatz=controlled_ansatz))
+        res1 = circuit1(params=params)
+        res2 = circuit2(params=params)
+        assert qml.math.allclose(res1, res2)
+
+    def test_ctrl_within_ctrl(self):
+        """Test using ctrl on a method that uses ctrl."""
+
+        def ansatz(params):
+            qml.RX(params[0], wires=0)
+            ctrl(qml.PauliX, control=0)(wires=1)
+            qml.RX(params[1], wires=0)
+
+        controlled_ansatz = ctrl(ansatz, 2)
+
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            controlled_ansatz([0.123, 0.456])
+
+        tape = QuantumScript.from_queue(q_tape)
+        assert tape.expand(1).circuit == [
+            *qml.CRX(0.123, wires=[2, 0]).decomposition(),
+            *qml.Toffoli(wires=[2, 0, 1]).decomposition(),
+            *qml.CRX(0.456, wires=[2, 0]).decomposition(),
+        ]
 
-def test_ctrl_template_and_operations():
-    """Test that a combination of controlled templates and operations correctly expands
-    on a device that doesn't support it"""
+    @pytest.mark.parametrize("ctrl_values", [[0, 0], [0, 1], [1, 0], [1, 1]])
+    def test_multi_ctrl_values(self, ctrl_values):
+        """Test control with a list of wires and control values."""
+
+        def expected_ops(ctrl_val):
+            exp_op = []
+            ctrl_wires = [3, 7]
+            for i, j in enumerate(ctrl_val):
+                if not bool(j):
+                    exp_op.append(qml.PauliX(ctrl_wires[i]))
+            exp_op.append(Controlled(qml.ControlledPhaseShift(pnp.pi / 2, [7, 0]), 3))
+            for i, j in enumerate(ctrl_val):
+                if not bool(j):
+                    exp_op.append(qml.PauliX(ctrl_wires[i]))
+
+            return exp_op
+
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            ctrl(qml.S, control=[3, 7], control_values=ctrl_values)(wires=0)
+        tape = QuantumScript.from_queue(q_tape)
+        assert len(tape.operations) == 1
+        op = tape.operations[0]
+        assert isinstance(op, Controlled)
+        new_tape = expand_tape(tape, 1)
+        assert equal_list(list(new_tape), expected_ops(ctrl_values))
 
-    weights = np.ones([3, 2])
+    def test_diagonal_ctrl(self):
+        """Test ctrl on diagonal gates."""
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            qml.ctrl(qml.DiagonalQubitUnitary, 1)(np.array([-1.0, 1.0j]), wires=0)
+        tape = QuantumScript.from_queue(q_tape)
+        tape = tape.expand(3, stop_at=lambda op: not isinstance(op, Controlled))
+        assert tape[0] == qml.DiagonalQubitUnitary(np.array([1.0, 1.0, -1.0, 1.0j]), wires=[1, 0])
 
-    def ansatz(weights, wires):
-        qml.PauliX(wires=wires[0])
-        qml.templates.BasicEntanglerLayers(weights, wires=wires)
+    @pytest.mark.parametrize("M", unitaries)
+    def test_qubit_unitary(self, M):
+        """Test ctrl on QubitUnitary"""
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            ctrl(qml.QubitUnitary, 1)(M, wires=0)
 
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        ctrl(ansatz, 0)(weights, wires=[1, 2])
+        tape = QuantumScript.from_queue(q_tape)
+        expected = qml.ControlledQubitUnitary(M, control_wires=1, wires=0)
+        assert equal_list(list(tape), expected)
 
-    tape = QuantumScript.from_queue(q_tape)
-    tape = tape.expand(depth=2, stop_at=lambda obj: not isinstance(obj, Controlled))
-    assert len(tape.operations) == 10
-    assert all(o.name in {"CNOT", "CRX", "Toffoli"} for o in tape.operations)
+        # causes decomposition into more basic operators
+        tape = tape.expand(3, stop_at=lambda op: not isinstance(op, Controlled))
+        assert not equal_list(list(tape), expected)
 
+    @pytest.mark.parametrize("M", unitaries)
+    def test_controlled_qubit_unitary(self, M):
+        """Test ctrl on ControlledQubitUnitary."""
 
-custom_controlled_ops = [  # operators with their own controlled class
-    (qml.PauliX, 1, qml.CNOT),
-    (qml.PauliY, 1, qml.CY),
-    (qml.PauliZ, 1, qml.CZ),
-    (qml.PauliX, 2, qml.Toffoli),
-]
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            ctrl(qml.ControlledQubitUnitary, 1)(M, control_wires=2, wires=0)
 
+        tape = QuantumScript.from_queue(q_tape)
+        # will immediately decompose according to selected decomposition algorithm
+        tape = tape.expand(1, stop_at=lambda op: not isinstance(op, Controlled))
 
-class TestCtrlCustomOperator:
-    @pytest.mark.parametrize("op_cls, num_ctrl_wires, custom_op_cls", custom_controlled_ops)
-    def test_ctrl_custom_operators(self, op_cls, num_ctrl_wires, custom_op_cls):
-        """Test that ctrl returns operators with their own controlled class."""
-        ctrl_wires = list(range(1, num_ctrl_wires + 1))
-        op = op_cls(wires=0)
-        ctrl_op = qml.ctrl(op, control=ctrl_wires)
-        custom_op = custom_op_cls(wires=ctrl_wires + [0])
-        assert qml.equal(ctrl_op, custom_op)
-        assert ctrl_op.name == custom_op.name
-
-    @pytest.mark.parametrize("op_cls, _, custom_op_cls", custom_controlled_ops)
-    def test_no_ctrl_custom_operators_excess_wires(self, op_cls, _, custom_op_cls):
-        """Test that ctrl returns a `Controlled` class when there is an excess of control wires."""
-        if op_cls is qml.PauliX:
-            pytest.skip("ctrl(PauliX) becomes MultiControlledX, not Controlled")
-
-        control_wires = list(range(1, 6))
-        op = op_cls(wires=0)
-        ctrl_op = qml.ctrl(op, control=control_wires)
-        expected = Controlled(op, control_wires=control_wires)
-        assert not isinstance(ctrl_op, custom_op_cls)
-        assert qml.equal(ctrl_op, expected)
-
-    @pytest.mark.parametrize("op_cls, num_ctrl_wires, custom_op_cls", custom_controlled_ops)
-    def test_no_ctrl_custom_operators_control_values(self, op_cls, num_ctrl_wires, custom_op_cls):
-        """Test that ctrl returns a `Controlled` class when the control value is not `True`."""
-        if op_cls is qml.PauliX:
-            pytest.skip("ctrl(PauliX) becomes MultiControlledX, not Controlled")
-
-        ctrl_wires = list(range(1, num_ctrl_wires + 1))
-        op = op_cls(wires=0)
-        ctrl_op = qml.ctrl(op, ctrl_wires, control_values=[False] * num_ctrl_wires)
-        expected = Controlled(op, ctrl_wires, control_values=[False] * num_ctrl_wires)
-        assert not isinstance(ctrl_op, custom_op_cls)
-        assert qml.equal(ctrl_op, expected)
+        expected = qml.ControlledQubitUnitary(M, control_wires=[1, 2], wires=0).decomposition()
+        assert tape.circuit == expected
 
     @pytest.mark.parametrize(
-        "control_wires,control_values,expected_values",
+        "op, params, depth, expected",
         [
-            ([1], (False), "0"),
-            ([1, 2], (0, 1), "01"),
-            ([1, 2, 3], (True, True, True), "111"),
-            ([1, 2, 3], (True, True, False), "110"),
-            ([1, 2, 3], None, None),
+            (qml.templates.QFT, [], 2, 14),
+            (qml.templates.BasicEntanglerLayers, [pnp.ones([3, 2])], 1, 9),
         ],
     )
-    def test_ctrl_PauliX_MultiControlledX(self, control_wires, control_values, expected_values):
-        """Tests that ctrl(PauliX) with 3+ control wires or Falsy control values make a MCX"""
-        with qml.queuing.AnnotatedQueue() as q:
-            op = qml.ctrl(qml.PauliX(0), control_wires, control_values=control_values)
+    def test_ctrl_templates(self, op, params, depth, expected):
+        """Test ctrl on two different templates."""
 
-        expected = qml.MultiControlledX(wires=control_wires + [0], control_values=expected_values)
-        assert len(q) == 1
-        assert qml.equal(op, expected)
-        assert qml.equal(q.queue[0], expected)
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            ctrl(op, 2)(*params, wires=[0, 1])
+        tape = QuantumScript.from_queue(q_tape)
+        expanded_tape = tape.expand(depth=depth)
+        assert len(expanded_tape.operations) == expected
+
+    def test_ctrl_template_and_operations(self):
+        """Test that a combination of controlled templates and operations correctly expands
+        on a device that doesn't support it"""
+
+        weights = pnp.ones([3, 2])
+
+        def ansatz(weights, wires):
+            qml.PauliX(wires=wires[0])
+            qml.templates.BasicEntanglerLayers(weights, wires=wires)
+
+        with qml.queuing.AnnotatedQueue() as q_tape:
+            ctrl(ansatz, 0)(weights, wires=[1, 2])
+
+        tape = QuantumScript.from_queue(q_tape)
+        tape = tape.expand(depth=1, stop_at=lambda obj: not isinstance(obj, Controlled))
+        assert len(tape.operations) == 10
+        assert all(o.name in {"CNOT", "CRX", "Toffoli"} for o in tape.operations)
 
 
 @pytest.mark.parametrize("diff_method", ["backprop", "parameter-shift", "finite-diff"])
@@ -1759,7 +2206,7 @@ class TestCtrlTransformDifferentiation:
         """Test differentiation using autograd"""
 
         dev = qml.device("default.qubit", wires=2)
-        init_state = np.array([1.0, -1.0], requires_grad=False) / np.sqrt(2)
+        init_state = pnp.array([1.0, -1.0], requires_grad=False) / pnp.sqrt(2)
 
         @qml.qnode(dev, diff_method=diff_method)
         def circuit(b):
@@ -1767,11 +2214,11 @@ class TestCtrlTransformDifferentiation:
             qml.ctrl(qml.RY, control=0)(b, wires=[1])
             return qml.expval(qml.PauliX(0))
 
-        b = np.array(0.123, requires_grad=True)
+        b = pnp.array(0.123, requires_grad=True)
         res = qml.grad(circuit)(b)
-        expected = np.sin(b / 2) / 2
+        expected = pnp.sin(b / 2) / 2
 
-        assert np.allclose(res, expected)
+        assert pnp.allclose(res, expected)
 
     @pytest.mark.torch
     def test_torch(self, diff_method):
@@ -1781,7 +2228,7 @@ class TestCtrlTransformDifferentiation:
         dev = qml.device("default.qubit", wires=2)
         init_state = torch.tensor(
             [1.0, -1.0], requires_grad=False, dtype=torch.complex128
-        ) / np.sqrt(2)
+        ) / pnp.sqrt(2)
 
         @qml.qnode(dev, diff_method=diff_method)
         def circuit(b):
@@ -1794,9 +2241,9 @@ class TestCtrlTransformDifferentiation:
         loss.backward()  # pylint:disable=no-member
 
         res = b.grad.detach()
-        expected = np.sin(b.detach() / 2) / 2
+        expected = pnp.sin(b.detach() / 2) / 2
 
-        assert np.allclose(res, expected)
+        assert pnp.allclose(res, expected)
 
     @pytest.mark.jax
     @pytest.mark.parametrize("jax_interface", ["auto", "jax", "jax-python"])
@@ -1813,16 +2260,16 @@ class TestCtrlTransformDifferentiation:
 
         @qml.qnode(dev, diff_method=diff_method, interface=jax_interface)
         def circuit(b):
-            init_state = onp.array([1.0, -1.0]) / onp.sqrt(2)
+            init_state = np.array([1.0, -1.0]) / np.sqrt(2)
             qml.StatePrep(init_state, wires=0)
             qml.ctrl(qml.RY, control=0)(b, wires=[1])
             return qml.expval(qml.PauliX(0))
 
         b = jnp.array(0.123)
         res = jax.grad(circuit)(b)
-        expected = np.sin(b / 2) / 2
+        expected = pnp.sin(b / 2) / 2
 
-        assert np.allclose(res, expected)
+        assert pnp.allclose(res, expected)
 
     @pytest.mark.tf
     def test_tf(self, diff_method):
@@ -1830,7 +2277,7 @@ class TestCtrlTransformDifferentiation:
         import tensorflow as tf
 
         dev = qml.device("default.qubit", wires=2)
-        init_state = tf.constant([1.0, -1.0], dtype=tf.complex128) / np.sqrt(2)
+        init_state = tf.constant([1.0, -1.0], dtype=tf.complex128) / pnp.sqrt(2)
 
         @qml.qnode(dev, diff_method=diff_method)
         def circuit(b):
@@ -1844,69 +2291,6 @@ class TestCtrlTransformDifferentiation:
             loss = circuit(b)
 
         res = tape.gradient(loss, b)
-        expected = np.sin(b / 2) / 2
-
-        assert np.allclose(res, expected)
-
-
-def test_ctrl_values_sanity_check():
-    """Test that control works with control values on a very standard usecase."""
-
-    def make_ops():
-        qml.RX(0.123, wires=0)
-        qml.RY(0.456, wires=2)
-        qml.RX(0.789, wires=0)
-        qml.Rot(0.111, 0.222, 0.333, wires=2)
-        qml.PauliX(wires=2)
-        qml.PauliY(wires=4)
-        qml.PauliZ(wires=0)
-
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        cmake_ops = ctrl(make_ops, control=1, control_values=0)
-        # Execute controlled version.
-        cmake_ops()
-
-    tape = QuantumScript.from_queue(q_tape)
-    expected = [
-        qml.PauliX(wires=1),
-        *qml.CRX(0.123, wires=[1, 0]).decomposition(),
-        *qml.CRY(0.456, wires=[1, 2]).decomposition(),
-        *qml.CRX(0.789, wires=[1, 0]).decomposition(),
-        *qml.CRot(0.111, 0.222, 0.333, wires=[1, 2]).decomposition(),
-        qml.CNOT(wires=[1, 2]),
-        *qml.CY(wires=[1, 4]).decomposition(),
-        *qml.CZ(wires=[1, 0]).decomposition(),
-        qml.PauliX(wires=1),
-    ]
-    assert len(tape) == 9
-    expanded = tape.expand(stop_at=lambda obj: not isinstance(obj, Controlled))
-    for op1, op2 in zip(expanded, expected):
-        assert qml.equal(op1, op2)
-
-
-@pytest.mark.parametrize("ctrl_values", [[0, 0], [0, 1], [1, 0], [1, 1]])
-def test_multi_ctrl_values(ctrl_values):
-    """Test control with a list of wires and control values."""
-
-    def expected_ops(ctrl_val):
-        exp_op = []
-        ctrl_wires = [3, 7]
-        for i, j in enumerate(ctrl_val):
-            if not bool(j):
-                exp_op.append(qml.PauliX(ctrl_wires[i]))
-        exp_op.append(Controlled(qml.PhaseShift(np.pi / 2, 0), [3, 7]))
-        for i, j in enumerate(ctrl_val):
-            if not bool(j):
-                exp_op.append(qml.PauliX(ctrl_wires[i]))
-
-        return exp_op
-
-    with qml.queuing.AnnotatedQueue() as q_tape:
-        CCS = ctrl(qml.S, control=[3, 7], control_values=ctrl_values)
-        CCS(wires=0)
-    tape = QuantumScript.from_queue(q_tape)
-    assert len(tape.operations) == 1
-    op = tape.operations[0]
-    assert isinstance(op, Controlled)
-    new_tape = expand_tape(tape, 1)
-    assert equal_list(list(new_tape), expected_ops(ctrl_values))
+        expected = pnp.sin(b / 2) / 2
+
+        assert pnp.allclose(res, expected)
diff --git a/tests/ops/op_math/test_controlled_decompositions.py b/tests/ops/op_math/test_controlled_decompositions.py
index 0c64ce991..d0b2560c1 100644
--- a/tests/ops/op_math/test_controlled_decompositions.py
+++ b/tests/ops/op_math/test_controlled_decompositions.py
@@ -14,8 +14,10 @@
 """
 Tests for the controlled decompositions.
 """
-import pytest
 
+import itertools
+
+import pytest
 import numpy as np
 import pennylane as qml
 from pennylane.ops import ctrl_decomp_zyz
@@ -28,6 +30,8 @@ from pennylane.ops.op_math.controlled_decompositions import (
     _convert_to_su2,
     _bisect_compute_a,
     _bisect_compute_b,
+    _decompose_mcx_with_one_worker,
+    _decompose_mcx_with_many_workers,
 )
 from pennylane.ops.op_math.controlled import (
     Controlled,
@@ -193,14 +197,14 @@ class TestControlledDecompositionZYZ:
         decomp = ctrl_decomp_zyz(op, [1])
         expected = [
             qml.RZ(np.pi, wires=0),
-            qml.MultiControlledX(wires=[1, 0]),
+            qml.CNOT(wires=[1, 0]),
             qml.RZ(-np.pi / 2, wires=0),
-            qml.MultiControlledX(wires=[1, 0]),
+            qml.CNOT(wires=[1, 0]),
             qml.RZ(-np.pi / 2, wires=0),
         ]
 
         assert len(decomp) == 5
-        assert all(qml.equal(o, e) for o, e in zip(decomp, expected))
+        assert decomp == expected
 
     @pytest.mark.parametrize("test_expand", [False, True])
     def test_zyz_decomp_no_control_values(self, test_expand):
@@ -372,9 +376,7 @@ class TestControlledBisectOD:
         assert qml.equal(mcx1, op_seq[0])
         assert qml.equal(mcx1, op_seq[4])
 
-        mcx2 = qml.MultiControlledX(
-            control_wires=Wires([4, 5]), wires=Wires(0), work_wires=Wires([1, 2, 3])
-        )
+        mcx2 = qml.Toffoli(wires=[4, 5, 0])
         assert qml.equal(mcx2, op_seq[2])
         assert qml.equal(mcx2, op_seq[6])
 
@@ -664,6 +666,69 @@ class TestControlledBisectGeneral:
         assert np.allclose(res, expected, atol=tol, rtol=tol)
 
 
+class TestMCXDecomposition:
+    @pytest.mark.parametrize("n_ctrl_wires", range(3, 6))
+    def test_decomposition_with_many_workers(self, n_ctrl_wires):
+        """Test that the decomposed MultiControlledX gate performs the same unitary as the
+        matrix-based version by checking if U^dagger U applies the identity to each basis
+        state. This test focuses on the case where there are many work wires."""
+        # pylint: disable=protected-access
+        control_wires = range(n_ctrl_wires)
+        target_wire = n_ctrl_wires
+        work_wires = range(n_ctrl_wires + 1, 2 * n_ctrl_wires + 1)
+
+        dev = qml.device("default.qubit", wires=2 * n_ctrl_wires + 1)
+
+        with qml.queuing.AnnotatedQueue() as q:
+            _decompose_mcx_with_many_workers(control_wires, target_wire, work_wires)
+        tape = qml.tape.QuantumScript.from_queue(q)
+        assert all(isinstance(op, qml.Toffoli) for op in tape.operations)
+
+        @qml.qnode(dev)
+        def f(bitstring):
+            qml.BasisState(bitstring, wires=range(n_ctrl_wires + 1))
+            qml.MultiControlledX(wires=list(control_wires) + [target_wire])
+            for op in tape.operations:
+                op.queue()
+            return qml.probs(wires=range(n_ctrl_wires + 1))
+
+        u = np.array(
+            [f(np.array(b)) for b in itertools.product(range(2), repeat=n_ctrl_wires + 1)]
+        ).T
+        assert np.allclose(u, np.eye(2 ** (n_ctrl_wires + 1)))
+
+    @pytest.mark.parametrize("n_ctrl_wires", range(3, 6))
+    def test_decomposition_with_one_worker(self, n_ctrl_wires):
+        """Test that the decomposed MultiControlledX gate performs the same unitary as the
+        matrix-based version by checking if U^dagger U applies the identity to each basis
+        state. This test focuses on the case where there is one work wire."""
+
+        # pylint: disable=protected-access
+        control_wires = Wires(range(n_ctrl_wires))
+        target_wire = n_ctrl_wires
+        work_wires = n_ctrl_wires + 1
+
+        dev = qml.device("default.qubit", wires=n_ctrl_wires + 2)
+
+        with qml.queuing.AnnotatedQueue() as q:
+            _decompose_mcx_with_one_worker(control_wires, target_wire, work_wires)
+        tape = qml.tape.QuantumScript.from_queue(q)
+        tape = tape.expand(depth=1)
+
+        @qml.qnode(dev)
+        def f(bitstring):
+            qml.BasisState(bitstring, wires=range(n_ctrl_wires + 1))
+            qml.MultiControlledX(wires=list(control_wires) + [target_wire])
+            for op in tape.operations:
+                op.queue()
+            return qml.probs(wires=range(n_ctrl_wires + 1))
+
+        u = np.array(
+            [f(np.array(b)) for b in itertools.product(range(2), repeat=n_ctrl_wires + 1)]
+        ).T
+        assert np.allclose(u, np.eye(2 ** (n_ctrl_wires + 1)))
+
+
 def test_ControlledQubitUnitary_has_decomposition_correct():
     """Test that ControlledQubitUnitary reports has_decomposition=False if it is False"""
     U = qml.Toffoli(wires=[0, 1, 2]).matrix()
diff --git a/tests/ops/op_math/test_controlled_ops.py b/tests/ops/op_math/test_controlled_ops.py
index 2e46fa496..ad325fcc2 100644
--- a/tests/ops/op_math/test_controlled_ops.py
+++ b/tests/ops/op_math/test_controlled_ops.py
@@ -15,8 +15,6 @@
 Unit tests for Operators inheriting from ControlledOp.
 """
 
-import functools
-
 import numpy as np
 import pytest
 from scipy.linalg import fractional_matrix_power
@@ -695,215 +693,6 @@ class TestComputations:
         )
 
 
-def _dot_broadcasted(a, b):
-    return np.einsum("...ij,...jk->...ik", a, b)
-
-
-def _multi_dot_broadcasted(matrices):
-    return functools.reduce(_dot_broadcasted, matrices)
-
-
-class TestDecompositions:
-    @pytest.mark.parametrize("phi", [0.432, np.array([0.1, 2.1])])
-    def test_CRX(self, phi):
-        """Test the decomposition for CRX."""
-
-        ops1 = qml.CRX.compute_decomposition(phi, wires=[0, 1])
-        ops2 = qml.CRX(phi, wires=(0, 1)).decomposition()
-
-        expected_ops = [
-            qml.RZ(np.pi / 2, 1),
-            qml.RY(phi / 2, 1),
-            qml.CNOT(wires=(0, 1)),
-            qml.RY(-phi / 2, 1),
-            qml.CNOT(wires=(0, 1)),
-            qml.RZ(-np.pi / 2, 1),
-        ]
-
-        assert ops1 == expected_ops
-        assert ops2 == expected_ops
-
-    @pytest.mark.parametrize("phi", [0.432, np.array([2.1, 0.2])])
-    def test_CRY(self, phi):
-        """Test the decomposition for CRY."""
-
-        ops1 = qml.CRY.compute_decomposition(phi, wires=[0, 1])
-        ops2 = qml.CRY(phi, wires=(0, 1)).decomposition()
-
-        expected_ops = [
-            qml.RY(phi / 2, 1),
-            qml.CNOT(wires=(0, 1)),
-            qml.RY(-phi / 2, 1),
-            qml.CNOT(wires=(0, 1)),
-        ]
-
-        assert ops1 == expected_ops
-        assert ops2 == expected_ops
-
-    @pytest.mark.parametrize("phi", [0.321, np.array([0.6, 2.1])])
-    def test_CRZ(self, phi):
-        """Test the decomposition for CRZ."""
-
-        ops1 = qml.CRZ.compute_decomposition(phi, wires=[1, 0])
-        ops2 = qml.CRZ(phi, wires=(1, 0)).decomposition()
-
-        expected_ops = [
-            qml.PhaseShift(phi / 2, wires=0),
-            qml.CNOT(wires=[1, 0]),
-            qml.PhaseShift(-phi / 2, wires=0),
-            qml.CNOT(wires=[1, 0]),
-        ]
-
-        assert ops1 == expected_ops
-        assert ops2 == expected_ops
-
-    @pytest.mark.parametrize("phi, theta, omega", [[0.5, 0.6, 0.7], [0.1, -0.4, 0.7], [-10, 5, -1]])
-    def test_CRot(self, tol, phi, theta, omega):
-        """Tests that the decomposition of the CRot gate is correct"""
-
-        op = qml.CRot(phi, theta, omega, wires=[0, 1])
-        res = op.decomposition()
-
-        mats = []
-        for i in reversed(res):
-            if len(i.wires) == 1:
-                mats.append(np.kron(np.eye(2), i.matrix()))
-            else:
-                mats.append(i.matrix())
-
-        decomposed_matrix = np.linalg.multi_dot(mats)
-        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
-
-    @pytest.mark.parametrize(
-        "phi, theta, omega",
-        [
-            [np.array([0.1, 0.2]), np.array([-0.4, 2.19]), np.array([0.7, -0.7])],
-            [np.array([0.1, 0.2, 0.9]), -0.4, np.array([0.7, 0.0, -0.7])],
-        ],
-    )
-    def test_CRot_broadcasted(self, tol, phi, theta, omega):
-        """Tests that the decomposition of the broadcasted CRot gate is correct"""
-
-        op = qml.CRot(phi, theta, omega, wires=[0, 1])
-        res = op.decomposition()
-
-        mats = []
-        for i in reversed(res):
-            mat = i.matrix()
-            if len(i.wires) == 1:
-                I = np.eye(2)[np.newaxis] if qml.math.ndim(mat) == 3 else np.eye(2)
-                mats.append(np.kron(I, mat))
-            else:
-                mats.append(mat)
-
-        decomposed_matrix = _multi_dot_broadcasted(mats)
-        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
-
-    @pytest.mark.parametrize("phi", [-0.1, 0.2, 0.5])
-    def test_controlled_phase_shift(self, phi):
-        """Tests that the ControlledPhaseShift calculates the correct decomposition"""
-
-        op = qml.ControlledPhaseShift(phi, wires=[0, 2])
-        decomp = op.decomposition()
-
-        mats = []
-        for i in reversed(decomp):
-            if i.wires.tolist() == [0]:
-                mats.append(np.kron(i.matrix(), np.eye(4)))
-            elif i.wires.tolist() == [1]:
-                mats.append(np.kron(np.eye(2), np.kron(i.matrix(), np.eye(2))))
-            elif i.wires.tolist() == [2]:
-                mats.append(np.kron(np.eye(4), i.matrix()))
-            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 1]:
-                mats.append(np.kron(i.matrix(), np.eye(2)))
-            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 2]:
-                mats.append(
-                    np.array(
-                        [
-                            [1, 0, 0, 0, 0, 0, 0, 0],
-                            [0, 1, 0, 0, 0, 0, 0, 0],
-                            [0, 0, 1, 0, 0, 0, 0, 0],
-                            [0, 0, 0, 1, 0, 0, 0, 0],
-                            [0, 0, 0, 0, 0, 1, 0, 0],
-                            [0, 0, 0, 0, 1, 0, 0, 0],
-                            [0, 0, 0, 0, 0, 0, 0, 1],
-                            [0, 0, 0, 0, 0, 0, 1, 0],
-                        ]
-                    )
-                )
-
-        decomposed_matrix = np.linalg.multi_dot(mats)
-        lam = np.exp(1j * phi)
-        exp = np.array(
-            [
-                [1, 0, 0, 0, 0, 0, 0, 0],
-                [0, 1, 0, 0, 0, 0, 0, 0],
-                [0, 0, 1, 0, 0, 0, 0, 0],
-                [0, 0, 0, 1, 0, 0, 0, 0],
-                [0, 0, 0, 0, 1, 0, 0, 0],
-                [0, 0, 0, 0, 0, lam, 0, 0],
-                [0, 0, 0, 0, 0, 0, 1, 0],
-                [0, 0, 0, 0, 0, 0, 0, lam],
-            ]
-        )
-
-        assert np.allclose(decomposed_matrix, exp)
-
-    def test_controlled_phase_shift_broadcasted(self):
-        """Tests that the ControlledPhaseShift calculates the correct decomposition"""
-
-        phi = np.array([-0.2, 4.2, 1.8])
-        op = qml.ControlledPhaseShift(phi, wires=[0, 2])
-        decomp = op.decomposition()
-
-        mats = []
-        for i in reversed(decomp):
-            mat = i.matrix()
-            eye = np.eye(2)[np.newaxis] if np.ndim(mat) == 3 else np.eye(2)
-            if i.wires.tolist() == [0]:
-                mats.append(np.kron(mat, np.kron(eye, eye)))
-            elif i.wires.tolist() == [1]:
-                mats.append(np.kron(eye, np.kron(mat, eye)))
-            elif i.wires.tolist() == [2]:
-                mats.append(np.kron(np.kron(eye, eye), mat))
-            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 1]:
-                mats.append(np.kron(mat, eye))
-            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 2]:
-                mats.append(
-                    np.array(
-                        [
-                            [1, 0, 0, 0, 0, 0, 0, 0],
-                            [0, 1, 0, 0, 0, 0, 0, 0],
-                            [0, 0, 1, 0, 0, 0, 0, 0],
-                            [0, 0, 0, 1, 0, 0, 0, 0],
-                            [0, 0, 0, 0, 0, 1, 0, 0],
-                            [0, 0, 0, 0, 1, 0, 0, 0],
-                            [0, 0, 0, 0, 0, 0, 0, 1],
-                            [0, 0, 0, 0, 0, 0, 1, 0],
-                        ]
-                    )
-                )
-
-        decomposed_matrix = _multi_dot_broadcasted(mats)
-        lam = np.exp(1j * phi)
-        exp = np.array([np.diag([1, 1, 1, 1, 1, el, 1, el]) for el in lam])
-
-        assert np.allclose(decomposed_matrix, exp)
-
-    @pytest.mark.parametrize("ops, expected_ops", NON_PARAMETRIC_OPS_DECOMPOSITIONS)
-    def test_non_parametric_op_decompositions(self, ops, expected_ops, tol):
-        """Tests that decompositions of non-parametrized operations are correct"""
-
-        op = ops(wires=[0, 1])
-        decomps = op.decomposition()
-        decomposed_matrix = qml.matrix(op.decomposition)()
-
-        for gate, expected in zip(decomps, expected_ops):
-            assert qml.equal(gate, expected, atol=tol, rtol=0)
-
-        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
-
-
 def test_simplify_crot():
     """Simplify CRot operations with different parameters."""
 
diff --git a/tests/ops/op_math/test_pow_op.py b/tests/ops/op_math/test_pow_op.py
index 975dc53f0..f47ecb7e6 100644
--- a/tests/ops/op_math/test_pow_op.py
+++ b/tests/ops/op_math/test_pow_op.py
@@ -465,10 +465,11 @@ class TestSimplify:
     def test_simplify_method_with_controlled_operation(self):
         """Test simplify method with controlled operation."""
         pow_op = Pow(ControlledOp(base=qml.Hadamard(0), control_wires=1, id=3), z=3)
-        final_op = ControlledOp(base=qml.Hadamard(0), control_wires=1, id=3)
+        final_op = qml.CH([1, 0], id=3)
         simplified_op = pow_op.simplify()
 
-        assert isinstance(simplified_op, ControlledOp)
+        # TODO: uncomment this check when all controlled operations inherit from ControlledOp
+        # assert isinstance(simplified_op, ControlledOp)
         assert final_op.data == simplified_op.data
         assert final_op.wires == simplified_op.wires
         assert final_op.arithmetic_depth == simplified_op.arithmetic_depth
diff --git a/tests/ops/qubit/test_non_parametric_ops.py b/tests/ops/qubit/test_non_parametric_ops.py
index 56edb2717..ae27f8175 100644
--- a/tests/ops/qubit/test_non_parametric_ops.py
+++ b/tests/ops/qubit/test_non_parametric_ops.py
@@ -574,13 +574,7 @@ class TestMultiControlledX:
                 None,
                 [0, 1, 2],
                 "011",
-                "Length of control bit string must equal number of control wires.",
-            ),
-            (
-                None,
-                [0, 1, 2],
-                [0, 1],
-                "Control values must be passed as a string.",
+                "Length of control values must equal number of control wires.",
             ),
             (
                 None,
@@ -590,8 +584,7 @@ class TestMultiControlledX:
             ),
             ([0], None, "", "Must specify the wires where the operation acts on"),
             ([0, 1], 2, "ab", "String of control values can contain only '0' or '1'."),
-            ([0, 1], 2, "011", "Length of control bit string must equal number of control wires."),
-            ([0, 1], 2, [0, 1], "Control values must be passed as a string."),
+            ([0, 1], 2, "011", "Length of control values must equal number of control wires."),
             ([0, 1], [2, 3], "10", "MultiControlledX accepts a single target wire."),
         ],
     )
@@ -670,6 +663,11 @@ class TestMultiControlledX:
 
         assert np.allclose(mpmct_state, pauli_x_state)
 
+    def test_decomposition_not_enough_wires(self):
+        """Test that the decomposition raises an error if the number of wires"""
+        with pytest.raises(ValueError, match="Wrong number of wires"):
+            qml.MultiControlledX.compute_decomposition((0,), control_values=[1])
+
     def test_decomposition_no_control_values(self):
         """Test decomposition has default control values of all ones."""
         decomp1 = qml.MultiControlledX.compute_decomposition((0, 1, 2))
@@ -741,71 +739,6 @@ class TestMultiControlledX:
 
         assert np.allclose(mpmct_state, pauli_x_state)
 
-    @pytest.mark.parametrize("n_ctrl_wires", range(3, 6))
-    def test_decomposition_with_many_workers(self, n_ctrl_wires):
-        """Test that the decomposed MultiControlledX gate performs the same unitary as the
-        matrix-based version by checking if U^dagger U applies the identity to each basis
-        state. This test focuses on the case where there are many work wires."""
-        # pylint: disable=protected-access
-        control_wires = range(n_ctrl_wires)
-        target_wire = n_ctrl_wires
-        work_wires = range(n_ctrl_wires + 1, 2 * n_ctrl_wires + 1)
-
-        dev = qml.device("default.qubit", wires=2 * n_ctrl_wires + 1)
-
-        with qml.queuing.AnnotatedQueue() as q:
-            qml.MultiControlledX._decomposition_with_many_workers(
-                control_wires, target_wire, work_wires
-            )
-        tape = qml.tape.QuantumScript.from_queue(q)
-        assert all(isinstance(op, qml.Toffoli) for op in tape.operations)
-
-        @qml.qnode(dev)
-        def f(bitstring):
-            qml.BasisState(bitstring, wires=range(n_ctrl_wires + 1))
-            qml.MultiControlledX(wires=list(control_wires) + [target_wire])
-            for op in tape.operations:
-                op.queue()
-            return qml.probs(wires=range(n_ctrl_wires + 1))
-
-        u = np.array(
-            [f(np.array(b)) for b in itertools.product(range(2), repeat=n_ctrl_wires + 1)]
-        ).T
-        assert np.allclose(u, np.eye(2 ** (n_ctrl_wires + 1)))
-
-    @pytest.mark.parametrize("n_ctrl_wires", range(3, 6))
-    def test_decomposition_with_one_worker(self, n_ctrl_wires):
-        """Test that the decomposed MultiControlledX gate performs the same unitary as the
-        matrix-based version by checking if U^dagger U applies the identity to each basis
-        state. This test focuses on the case where there is one work wire."""
-        # pylint: disable=protected-access
-        control_wires = Wires(range(n_ctrl_wires))
-        target_wire = n_ctrl_wires
-        work_wires = n_ctrl_wires + 1
-
-        dev = qml.device("default.qubit", wires=n_ctrl_wires + 2)
-
-        with qml.queuing.AnnotatedQueue() as q:
-            qml.MultiControlledX._decomposition_with_one_worker(
-                control_wires, target_wire, work_wires
-            )
-        tape = qml.tape.QuantumScript.from_queue(q)
-        tape = tape.expand(depth=1)
-        assert all(isinstance(op, (qml.Toffoli, qml.CNOT)) for op in tape.operations)
-
-        @qml.qnode(dev)
-        def f(bitstring):
-            qml.BasisState(bitstring, wires=range(n_ctrl_wires + 1))
-            qml.MultiControlledX(wires=list(control_wires) + [target_wire])
-            for op in tape.operations:
-                op.queue()
-            return qml.probs(wires=range(n_ctrl_wires + 1))
-
-        u = np.array(
-            [f(np.array(b)) for b in itertools.product(range(2), repeat=n_ctrl_wires + 1)]
-        ).T
-        assert np.allclose(u, np.eye(2 ** (n_ctrl_wires + 1)))
-
     def test_not_enough_workers(self):
         """Test that a ValueError is raised when more than 2 control wires are to be decomposed with
         no work wires supplied"""
@@ -821,7 +754,8 @@ class TestMultiControlledX:
         control_target_wires = range(4)
         work_wires = range(2)
         with pytest.raises(
-            ValueError, match="The work wires must be different from the control and target wires"
+            ValueError,
+            match="Work wires must be different the control_wires and base operation wires.",
         ):
             qml.MultiControlledX(wires=control_target_wires, work_wires=work_wires)
 
@@ -934,9 +868,9 @@ class TestMultiControlledX:
     def test_repr(self):
         """Test ``__repr__`` method that shows ``control_values``"""
         wires = [0, 1, 2]
-        control_values = "01"
+        control_values = [False, True]
         op_repr = qml.MultiControlledX(wires=wires, control_values=control_values).__repr__()
-        assert op_repr == f'MultiControlledX(wires={wires}, control_values="{control_values}")'
+        assert op_repr == f"MultiControlledX(wires={wires}, control_values={control_values})"
 
 
 period_two_ops = (
diff --git a/tests/transforms/test_optimization/test_undo_swaps.py b/tests/transforms/test_optimization/test_undo_swaps.py
index a5ff7ec33..7e9a1c878 100644
--- a/tests/transforms/test_optimization/test_undo_swaps.py
+++ b/tests/transforms/test_optimization/test_undo_swaps.py
@@ -32,7 +32,7 @@ class TestUndoSwaps:
         ops = [
             qml.adjoint(qml.S(0)),
             qml.PauliRot(1.2, "XY", wires=(0, 2)),
-            qml.ctrl(qml.PauliX(0), 2, control_values=[0, 0]),
+            qml.ctrl(qml.PauliX(0), [2, 3], control_values=[0, 0]),
             qml.SWAP((0, 1)),
         ]
 
@@ -42,7 +42,7 @@ class TestUndoSwaps:
         expected_ops = [
             qml.adjoint(qml.S(1)),
             qml.PauliRot(1.2, "XY", wires=(1, 2)),
-            qml.ctrl(qml.PauliX(1), 2, control_values=[0, 0]),
+            qml.ctrl(qml.PauliX(1), [2, 3], control_values=[0, 0]),
         ]
         assert len(batch) == 1
         assert batch[0].shots == tape.shots