Repository: PennyLaneAI/pennylane
Pull request #7591 (2025-06-03 16:36:54)

[Decomposition] decomposition of PCPhase as a quantum function
**Context:**

**Description of the Change:**

Translates the decomposition of `PCPhase` into a quantum function.

**Benefits:**

**Possible Drawbacks:**

**Related GitHub Issues:**

[sc-89611]

mudit2812: Could you remove `[WIP]` from the PR title?

Patch:
diff --git a/.github/stable/all_interfaces.txt b/.github/stable/all_interfaces.txt
index 764a02376..71bf05eed 100644
--- a/.github/stable/all_interfaces.txt
+++ b/.github/stable/all_interfaces.txt
@@ -25,14 +25,14 @@ execnet==2.1.1
 filelock==3.18.0
 flaky==3.8.1
 flatbuffers==25.2.10
-fonttools==4.58.1
+fonttools==4.58.2
 fsspec==2025.5.1
 gast==0.6.0
 gitdb==4.0.12
 GitPython==3.1.44
 google-pasta==0.2.0
-grpcio==1.71.0
-h5py==3.13.0
+grpcio==1.72.1
+h5py==3.14.0
 identify==2.6.12
 idna==3.10
 iniconfig==2.1.0
@@ -74,9 +74,9 @@ opt_einsum==3.4.0
 optree==0.16.0
 osqp==1.0.4
 packaging==25.0
-pandas==2.2.3
+pandas==2.3.0
 pathspec==0.12.1
-PennyLane_Lightning==0.42.0.dev23
+PennyLane_Lightning==0.42.0.dev30
 pillow==11.2.1
 platformdirs==4.3.8
 pluggy==1.6.0
@@ -91,7 +91,7 @@ Pygments==2.19.1
 pylint==2.7.4
 pyparsing==3.2.3
 pyproject_hooks==1.2.0
-pytest==8.3.5
+pytest==8.4.0
 pytest-benchmark==5.1.0
 pytest-cov==6.1.1
 pytest-forked==1.6.0
@@ -106,7 +106,7 @@ requests==2.32.3
 rich==13.9.4
 rustworkx==0.16.0
 scipy==1.15.3
-scipy-openblas32==0.3.29.265.0
+scipy-openblas32==0.3.29.265.1
 scs==3.2.7.post2
 six==1.17.0
 smmap==5.0.2
@@ -121,10 +121,10 @@ tf_keras==2.18.0
 toml==0.10.2
 tomli==2.2.1
 tomli_w==1.2.0
-tomlkit==0.13.2
+tomlkit==0.13.3
 torch==2.5.1
 triton==3.1.0
-typing_extensions==4.13.2
+typing_extensions==4.14.0
 tzdata==2025.2
 urllib3==2.4.0
 virtualenv==20.31.2
diff --git a/.github/stable/core.txt b/.github/stable/core.txt
index ceb876d39..70a722bdd 100644
--- a/.github/stable/core.txt
+++ b/.github/stable/core.txt
@@ -23,7 +23,7 @@ exceptiongroup==1.3.0
 execnet==2.1.1
 filelock==3.18.0
 flaky==3.8.1
-fonttools==4.58.1
+fonttools==4.58.2
 gast==0.6.0
 gitdb==4.0.12
 GitPython==3.1.44
@@ -46,9 +46,9 @@ nodeenv==1.9.1
 numpy==2.2.6
 osqp==1.0.4
 packaging==25.0
-pandas==2.2.3
+pandas==2.3.0
 pathspec==0.12.1
-PennyLane_Lightning==0.42.0.dev23
+PennyLane_Lightning==0.42.0.dev30
 pillow==11.2.1
 platformdirs==4.3.8
 pluggy==1.6.0
@@ -62,7 +62,7 @@ Pygments==2.19.1
 pylint==2.7.4
 pyparsing==3.2.3
 pyproject_hooks==1.2.0
-pytest==8.3.5
+pytest==8.4.0
 pytest-benchmark==5.1.0
 pytest-cov==6.1.1
 pytest-forked==1.6.0
@@ -77,7 +77,7 @@ requests==2.32.3
 rich==13.9.4
 rustworkx==0.16.0
 scipy==1.15.3
-scipy-openblas32==0.3.29.265.0
+scipy-openblas32==0.3.29.265.1
 scs==3.2.7.post2
 six==1.17.0
 smmap==5.0.2
@@ -86,8 +86,8 @@ termcolor==3.1.0
 toml==0.10.2
 tomli==2.2.1
 tomli_w==1.2.0
-tomlkit==0.13.2
-typing_extensions==4.13.2
+tomlkit==0.13.3
+typing_extensions==4.14.0
 tzdata==2025.2
 urllib3==2.4.0
 virtualenv==20.31.2
diff --git a/.github/stable/external.txt b/.github/stable/external.txt
index bc83fb928..7b4bdb198 100644
--- a/.github/stable/external.txt
+++ b/.github/stable/external.txt
@@ -4,7 +4,7 @@ antlr4-python3-runtime==4.13.2
 anyio==4.9.0
 anywidget==0.9.18
 appdirs==1.4.4
-argon2-cffi==23.1.0
+argon2-cffi==25.1.0
 argon2-cffi-bindings==21.2.0
 arrow==1.3.0
 astroid==2.6.6
@@ -52,11 +52,12 @@ exceptiongroup==1.3.0
 execnet==2.1.1
 executing==2.2.0
 fastjsonschema==2.21.1
+filecheck==1.0.2
 filelock==3.18.0
 flaky==3.8.1
 Flask==3.0.3
 flatbuffers==25.2.10
-fonttools==4.58.1
+fonttools==4.58.2
 fqdn==1.5.1
 fxpmath==0.4.9
 galois==0.4.6
@@ -65,9 +66,9 @@ gitdb==4.0.12
 GitPython==3.1.44
 google-pasta==0.2.0
 graphviz==0.20.3
-grpcio==1.71.0
+grpcio==1.72.1
 h11==0.16.0
-h5py==3.13.0
+h5py==3.14.0
 httpcore==1.0.9
 httpx==0.28.1
 ibm-cloud-sdk-core==3.24.1
@@ -121,7 +122,7 @@ ml_dtypes==0.5.1
 mpmath==1.3.0
 mypy_extensions==1.1.0
 namex==0.1.0
-narwhals==1.41.0
+narwhals==1.41.1
 nbclient==0.10.2
 nbconvert==7.16.6
 nbformat==5.10.4
@@ -140,12 +141,12 @@ ordered-set==4.1.0
 osqp==1.0.4
 overrides==7.7.0
 packaging==25.0
-pandas==2.2.3
+pandas==2.3.0
 pandocfilters==1.5.1
 parso==0.8.4
 pathspec==0.12.1
 pbr==6.1.1
-PennyLane-Catalyst==0.12.0.dev45
+PennyLane-Catalyst==0.12.0.dev50
 PennyLane-qiskit @ git+https://github.com/PennyLaneAI/pennylane-qiskit.git@5544c00d2022d8ca07fca97a0d3ec8de37811a88
 PennyLane_Lightning==0.42.0.dev16
 PennyLane_Lightning_Kokkos==0.42.0.dev16
@@ -156,7 +157,7 @@ plotly==6.1.2
 pluggy==1.6.0
 ply==3.11
 pre_commit==4.2.0
-prometheus_client==0.22.0
+prometheus_client==0.22.1
 prompt_toolkit==3.0.51
 protobuf==5.29.5
 psutil==7.0.0
@@ -176,7 +177,7 @@ pyparsing==3.1.4
 pyperclip==1.9.0
 pyproject_hooks==1.2.0
 pyspnego==0.11.2
-pytest==8.3.5
+pytest==8.4.0
 pytest-benchmark==5.1.0
 pytest-cov==6.1.1
 pytest-forked==1.6.0
@@ -209,7 +210,7 @@ rich==13.9.4
 rpds-py==0.25.1
 rustworkx==0.16.0
 scipy==1.15.2
-scipy-openblas32==0.3.29.265.0
+scipy-openblas32==0.3.29.265.1
 scs==3.2.7.post2
 Send2Trash==1.8.3
 six==1.17.0
@@ -235,14 +236,14 @@ tinycss2==1.4.0
 toml==0.10.2
 tomli==2.2.1
 tomli_w==1.2.0
-tomlkit==0.13.2
+tomlkit==0.13.3
 toolz==1.0.0
 tornado==6.5.1
 tqdm==4.67.1
 traitlets==5.14.3
 types-python-dateutil==2.9.0.20250516
 typing-inspection==0.4.1
-typing_extensions==4.13.2
+typing_extensions==4.14.0
 tzdata==2025.2
 uri-template==1.3.0
 urllib3==2.4.0
@@ -256,4 +257,4 @@ Werkzeug==3.0.6
 widgetsnbextension==4.0.14
 wrapt==1.12.1
 xdsl==0.38.0
-zipp==3.22.0
+zipp==3.23.0
diff --git a/.github/stable/jax.txt b/.github/stable/jax.txt
index 6423c8bf1..1f0f81908 100644
--- a/.github/stable/jax.txt
+++ b/.github/stable/jax.txt
@@ -23,7 +23,7 @@ exceptiongroup==1.3.0
 execnet==2.1.1
 filelock==3.18.0
 flaky==3.8.1
-fonttools==4.58.1
+fonttools==4.58.2
 gast==0.6.0
 gitdb==4.0.12
 GitPython==3.1.44
@@ -50,9 +50,9 @@ numpy==2.2.6
 opt_einsum==3.4.0
 osqp==1.0.4
 packaging==25.0
-pandas==2.2.3
+pandas==2.3.0
 pathspec==0.12.1
-PennyLane_Lightning==0.42.0.dev23
+PennyLane_Lightning==0.42.0.dev30
 pillow==11.2.1
 platformdirs==4.3.8
 pluggy==1.6.0
@@ -66,7 +66,7 @@ Pygments==2.19.1
 pylint==2.7.4
 pyparsing==3.2.3
 pyproject_hooks==1.2.0
-pytest==8.3.5
+pytest==8.4.0
 pytest-benchmark==5.1.0
 pytest-cov==6.1.1
 pytest-forked==1.6.0
@@ -81,7 +81,7 @@ requests==2.32.3
 rich==13.9.4
 rustworkx==0.16.0
 scipy==1.15.3
-scipy-openblas32==0.3.29.265.0
+scipy-openblas32==0.3.29.265.1
 scs==3.2.7.post2
 six==1.17.0
 smmap==5.0.2
@@ -90,8 +90,8 @@ termcolor==3.1.0
 toml==0.10.2
 tomli==2.2.1
 tomli_w==1.2.0
-tomlkit==0.13.2
-typing_extensions==4.13.2
+tomlkit==0.13.3
+typing_extensions==4.14.0
 tzdata==2025.2
 urllib3==2.4.0
 virtualenv==20.31.2
diff --git a/.github/stable/tf.txt b/.github/stable/tf.txt
index 8ac2d062d..8eda3cb47 100644
--- a/.github/stable/tf.txt
+++ b/.github/stable/tf.txt
@@ -25,13 +25,13 @@ execnet==2.1.1
 filelock==3.18.0
 flaky==3.8.1
 flatbuffers==25.2.10
-fonttools==4.58.1
+fonttools==4.58.2
 gast==0.6.0
 gitdb==4.0.12
 GitPython==3.1.44
 google-pasta==0.2.0
-grpcio==1.71.0
-h5py==3.13.0
+grpcio==1.72.1
+h5py==3.14.0
 identify==2.6.12
 idna==3.10
 iniconfig==2.1.0
@@ -58,9 +58,9 @@ opt_einsum==3.4.0
 optree==0.16.0
 osqp==1.0.4
 packaging==25.0
-pandas==2.2.3
+pandas==2.3.0
 pathspec==0.12.1
-PennyLane_Lightning==0.42.0.dev23
+PennyLane_Lightning==0.42.0.dev30
 pillow==11.2.1
 platformdirs==4.3.8
 pluggy==1.6.0
@@ -75,7 +75,7 @@ Pygments==2.19.1
 pylint==2.7.4
 pyparsing==3.2.3
 pyproject_hooks==1.2.0
-pytest==8.3.5
+pytest==8.4.0
 pytest-benchmark==5.1.0
 pytest-cov==6.1.1
 pytest-forked==1.6.0
@@ -90,7 +90,7 @@ requests==2.32.3
 rich==13.9.4
 rustworkx==0.16.0
 scipy==1.15.3
-scipy-openblas32==0.3.29.265.0
+scipy-openblas32==0.3.29.265.1
 scs==3.2.7.post2
 six==1.17.0
 smmap==5.0.2
@@ -104,8 +104,8 @@ tf_keras==2.18.0
 toml==0.10.2
 tomli==2.2.1
 tomli_w==1.2.0
-tomlkit==0.13.2
-typing_extensions==4.13.2
+tomlkit==0.13.3
+typing_extensions==4.14.0
 tzdata==2025.2
 urllib3==2.4.0
 virtualenv==20.31.2
diff --git a/.github/stable/torch.txt b/.github/stable/torch.txt
index 4f572a663..d0c5a4e39 100644
--- a/.github/stable/torch.txt
+++ b/.github/stable/torch.txt
@@ -23,7 +23,7 @@ exceptiongroup==1.3.0
 execnet==2.1.1
 filelock==3.18.0
 flaky==3.8.1
-fonttools==4.58.1
+fonttools==4.58.2
 fsspec==2025.5.1
 gast==0.6.0
 gitdb==4.0.12
@@ -60,9 +60,9 @@ nvidia-nvjitlink-cu12==12.4.127
 nvidia-nvtx-cu12==12.4.127
 osqp==1.0.4
 packaging==25.0
-pandas==2.2.3
+pandas==2.3.0
 pathspec==0.12.1
-PennyLane_Lightning==0.42.0.dev23
+PennyLane_Lightning==0.42.0.dev30
 pillow==11.2.1
 platformdirs==4.3.8
 pluggy==1.6.0
@@ -76,7 +76,7 @@ Pygments==2.19.1
 pylint==2.7.4
 pyparsing==3.2.3
 pyproject_hooks==1.2.0
-pytest==8.3.5
+pytest==8.4.0
 pytest-benchmark==5.1.0
 pytest-cov==6.1.1
 pytest-forked==1.6.0
@@ -91,7 +91,7 @@ requests==2.32.3
 rich==13.9.4
 rustworkx==0.16.0
 scipy==1.15.3
-scipy-openblas32==0.3.29.265.0
+scipy-openblas32==0.3.29.265.1
 scs==3.2.7.post2
 six==1.17.0
 smmap==5.0.2
@@ -101,10 +101,10 @@ termcolor==3.1.0
 toml==0.10.2
 tomli==2.2.1
 tomli_w==1.2.0
-tomlkit==0.13.2
+tomlkit==0.13.3
 torch==2.5.1
 triton==3.1.0
-typing_extensions==4.13.2
+typing_extensions==4.14.0
 tzdata==2025.2
 urllib3==2.4.0
 virtualenv==20.31.2
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 1c924a7fa..7eb9678cc 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -202,34 +202,27 @@
   Here, when the Hadamard and ``CRZ`` have relatively high weights, a decomposition involving them is considered *less* 
   efficient. When they have relatively low weights, a decomposition involving them is considered *more* efficient.
 
-* The decomposition of `qml.PCPhase` is now significantly more efficient for more than 2 qubits.
-  [(#7166)](https://github.com/PennyLaneAI/pennylane/pull/7166)
+* Decomposition rules that can be accessed with the new graph-based decomposition system are
+  implemented for the following operators:
 
-* New decomposition rules comprising rotation gates and global phases have been added to `QubitUnitary` 
-  and `ControlledQubitUnitary` that can be accessed with the new graph-based decomposition system. 
-  The most efficient set of rotations to decompose into will be chosen based on the target gate set.
-  [(#7211)](https://github.com/PennyLaneAI/pennylane/pull/7211)
-  [(#7371)](https://github.com/PennyLaneAI/pennylane/pull/7371)
+  * :class:`~.QubitUnitary`
+    [(#7211)](https://github.com/PennyLaneAI/pennylane/pull/7211)
 
-  ```python
-  from functools import partial
-  import numpy as np
-  import pennylane as qml
-  
-  qml.decomposition.enable_graph()
-  
-  U = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
-  
-  @partial(qml.transforms.decompose, gate_set={"RX", "RY", "GlobalPhase"})
-  @qml.qnode(qml.device("default.qubit"))
-  def circuit():
-      qml.QubitUnitary(np.array([[1, 1], [1, -1]]) / np.sqrt(2), wires=[0])
-      return qml.expval(qml.PauliZ(0))
-  ```
-  ```pycon
-  >>> print(qml.draw(circuit)())
-  0: ──RX(0.00)──RY(1.57)──RX(3.14)──GlobalPhase(-1.57)─┤  <Z>
-  ```
+  * :class:`~.ControlledQubitUnitary`
+    [(#7371)](https://github.com/PennyLaneAI/pennylane/pull/7371)
+
+  * :class:`~.MultiControlledX`
+    [(#7405)](https://github.com/PennyLaneAI/pennylane/pull/7405)
+
+  * :class:`~pennylane.ops.Exp`. 
+    [(#7489)](https://github.com/PennyLaneAI/pennylane/pull/7489)
+
+    Specifically, the following decompositions have been added:
+    * Suzuki-Trotter decomposition when the `num_steps` keyword argument is specified.
+    * Decomposition to a :class:`~pennylane.PauliRot` when the base is a single-term Pauli word.
+
+  * :class:`~.PCPhase`
+    [(#7591)](https://github.com/PennyLaneAI/pennylane/pull/7591)
 
 * A new decomposition rule that uses a single work wire for decomposing multi-controlled operators is added.
   [(#7383)](https://github.com/PennyLaneAI/pennylane/pull/7383)
@@ -325,12 +318,6 @@
   0: ──RX(-0.50)─╭●────────────┤  <Z>
   1: ────────────╰X──RY(-0.50)─┤
   ```
-  
-* Decomposition rules compatible with the new graph-based decomposition system have been implemented
-  for :class:`~pennylane.ops.Exp`. Specifically, the following decompositions have been added:
-  * Suzuki-Trotter decomposition when the `num_steps` keyword argument is specified.
-  * Decomposition to a :class:`~pennylane.PauliRot` when the base is a single-term Pauli word.
-  [(#7489)](https://github.com/PennyLaneAI/pennylane/pull/7489)
 
 * A `work_wire_type` argument has been added to :func:`~pennylane.ctrl` and :class:`~pennylane.ControlledQubitUnitary`
   for more fine-grained control over the type of work wire used in their decompositions.
@@ -352,7 +339,10 @@
   no longer depends on the wires of `GlobalPhase` or `Identity`. Control nodes of controlled global phases
   and identities no longer receive the operator label, which is in line with other controlled operations.
   [(#7457)](https://github.com/PennyLaneAI/pennylane/pull/7457)
-  
+
+* The decomposition of `qml.PCPhase` is now significantly more efficient for more than 2 qubits.
+  [(#7166)](https://github.com/PennyLaneAI/pennylane/pull/7166)
+
 * :class:`~.QubitUnitary` now supports a decomposition that is compatible with an arbitrary number of qubits. 
   This represents a fundamental improvement over the previous implementation, which was limited to two-qubit systems.
   [(#7277)](https://github.com/PennyLaneAI/pennylane/pull/7277)
diff --git a/pennylane/ops/qubit/parametric_ops_multi_qubit.py b/pennylane/ops/qubit/parametric_ops_multi_qubit.py
index 83819b6af..de8d89a82 100644
--- a/pennylane/ops/qubit/parametric_ops_multi_qubit.py
+++ b/pennylane/ops/qubit/parametric_ops_multi_qubit.py
@@ -18,17 +18,18 @@ core parametrized gates.
 """
 # pylint: disable=arguments-differ
 import functools
+from collections import Counter
 from operator import matmul
 from typing import Optional, Union
 
 import numpy as np
 
 import pennylane as qml
-from pennylane.decomposition import add_decomps, register_resources
+from pennylane import math, queuing
+from pennylane.decomposition import add_decomps, controlled_resource_rep, register_resources
 from pennylane.decomposition.symbolic_decomposition import adjoint_rotation, pow_rotation
-from pennylane.math import expand_matrix
 from pennylane.math.decomposition import decomp_int_to_powers_of_two
-from pennylane.operation import FlatPytree, Operation
+from pennylane.operation import FlatPytree, Operation, Operator
 from pennylane.typing import TensorLike
 from pennylane.wires import Wires, WiresLike
 
@@ -110,19 +111,17 @@ class MultiRZ(Operation):
                 [0.0000+0.0000j, 0.0000+0.0000j, 0.9988+0.0500j, 0.0000+0.0000j],
                 [0.0000+0.0000j, 0.0000+0.0000j, 0.0000+0.0000j, 0.9988-0.0500j]])
         """
-        eigs = qml.math.convert_like(qml.pauli.pauli_eigs(num_wires), theta)
+        eigs = math.convert_like(qml.pauli.pauli_eigs(num_wires), theta)
 
-        if qml.math.get_interface(theta) == "tensorflow":
-            theta = qml.math.cast_like(theta, 1j)
-            eigs = qml.math.cast_like(eigs, 1j)
+        if math.get_interface(theta) == "tensorflow":
+            theta = math.cast_like(theta, 1j)
+            eigs = math.cast_like(eigs, 1j)
 
-        if qml.math.ndim(theta) == 0:
-            return qml.math.diag(qml.math.exp(-0.5j * theta * eigs))
+        if math.ndim(theta) == 0:
+            return math.diag(math.exp(-0.5j * theta * eigs))
 
-        diags = qml.math.exp(qml.math.outer(-0.5j * theta, eigs))
-        return diags[:, :, np.newaxis] * qml.math.cast_like(
-            qml.math.eye(2**num_wires, like=diags), diags
-        )
+        diags = math.exp(math.outer(-0.5j * theta, eigs))
+        return diags[:, :, np.newaxis] * math.cast_like(math.eye(2**num_wires, like=diags), diags)
 
     def generator(self) -> "qml.Hamiltonian":
         return qml.Hamiltonian([-0.5], [functools.reduce(matmul, [PauliZ(w) for w in self.wires])])
@@ -156,21 +155,21 @@ class MultiRZ(Operation):
         tensor([0.9689-0.2474j, 0.9689+0.2474j, 0.9689+0.2474j, 0.9689-0.2474j,
                 0.9689+0.2474j, 0.9689-0.2474j, 0.9689-0.2474j, 0.9689+0.2474j])
         """
-        eigs = qml.math.convert_like(qml.pauli.pauli_eigs(num_wires), theta)
+        eigs = math.convert_like(qml.pauli.pauli_eigs(num_wires), theta)
 
-        if qml.math.get_interface(theta) == "tensorflow":
-            theta = qml.math.cast_like(theta, 1j)
-            eigs = qml.math.cast_like(eigs, 1j)
+        if math.get_interface(theta) == "tensorflow":
+            theta = math.cast_like(theta, 1j)
+            eigs = math.cast_like(eigs, 1j)
 
-        if qml.math.ndim(theta) == 0:
-            return qml.math.exp(-0.5j * theta * eigs)
+        if math.ndim(theta) == 0:
+            return math.exp(-0.5j * theta * eigs)
 
-        return qml.math.exp(qml.math.outer(-0.5j * theta, eigs))
+        return math.exp(math.outer(-0.5j * theta, eigs))
 
     @staticmethod
     def compute_decomposition(  # pylint: disable=arguments-differ,unused-argument
         theta: TensorLike, wires: WiresLike, **kwargs
-    ) -> list["qml.operation.Operator"]:
+    ) -> list[Operator]:
         r"""Representation of the operator as a product of other operators (static method). :
 
         .. math:: O = O_1 O_2 \dots O_n.
@@ -204,7 +203,7 @@ class MultiRZ(Operation):
     def adjoint(self) -> "MultiRZ":
         return MultiRZ(-self.parameters[0], wires=self.wires)
 
-    def pow(self, z: Union[int, float]) -> list["qml.operation.Operator"]:
+    def pow(self, z: Union[int, float]) -> list[Operator]:
         return [MultiRZ(self.data[0] * z, wires=self.wires)]
 
     def simplify(self) -> "MultiRZ":
@@ -373,7 +372,7 @@ class PauliRot(Operation):
 
         # TODO[dwierichs]: Implement a proper label for parameter-broadcasted operators
         if decimals is not None and self.batch_size is None:
-            param_string = f"\n({qml.math.asarray(self.parameters[0]):.{decimals}f})"
+            param_string = f"\n({math.asarray(self.parameters[0]):.{decimals}f})"
             op_label += param_string
 
         return op_label
@@ -423,10 +422,10 @@ class PauliRot(Operation):
                 "Allowed characters are I, X, Y and Z"
             )
 
-        interface = qml.math.get_interface(theta)
+        interface = math.get_interface(theta)
 
         if interface == "tensorflow":
-            theta = qml.math.cast_like(theta, 1j)
+            theta = math.cast_like(theta, 1j)
 
         # Simplest case is if the Pauli is the identity matrix
         if set(pauli_word) == {"I"}:
@@ -442,18 +441,18 @@ class PauliRot(Operation):
 
         # now we conjugate with Hadamard and RX to create the Pauli string
         conjugation_matrix = functools.reduce(
-            qml.math.kron,
+            math.kron,
             [PauliRot._PAULI_CONJUGATION_MATRICES[gate] for gate in non_identity_gates],
         )
         if interface == "tensorflow":
-            conjugation_matrix = qml.math.cast_like(conjugation_matrix, 1j)
+            conjugation_matrix = math.cast_like(conjugation_matrix, 1j)
         # Note: we use einsum with reverse arguments here because it is not multi-dispatched
         # and the tensordot containing multi_Z_rot_matrix should decide about the interface
-        return expand_matrix(
-            qml.math.einsum(
+        return math.expand_matrix(
+            math.einsum(
                 "...jk,ij->...ik",
-                qml.math.tensordot(multi_Z_rot_matrix, conjugation_matrix, axes=[[-1], [0]]),
-                qml.math.conj(conjugation_matrix),
+                math.tensordot(multi_Z_rot_matrix, conjugation_matrix, axes=[[-1], [0]]),
+                math.conj(conjugation_matrix),
             ),
             non_identity_wires,
             list(range(len(pauli_word))),
@@ -491,8 +490,8 @@ class PauliRot(Operation):
         >>> qml.PauliRot.compute_eigvals(torch.tensor(0.5), "X")
         tensor([0.9689-0.2474j, 0.9689+0.2474j])
         """
-        if qml.math.get_interface(theta) == "tensorflow":
-            theta = qml.math.cast_like(theta, 1j)
+        if math.get_interface(theta) == "tensorflow":
+            theta = math.cast_like(theta, 1j)
 
         # Identity must be treated specially because its eigenvalues are all the same
         if set(pauli_word) == {"I"}:
@@ -503,7 +502,7 @@ class PauliRot(Operation):
     @staticmethod
     def compute_decomposition(
         theta: TensorLike, wires: WiresLike, pauli_word: str
-    ) -> list["qml.operation.Operator"]:
+    ) -> list[Operator]:
         r"""Representation of the operator as a product of other operators (static method). :
 
         .. math:: O = O_1 O_2 \dots O_n.
@@ -600,71 +599,6 @@ add_decomps("Adjoint(PauliRot)", adjoint_rotation)
 add_decomps("Pow(PauliRot)", pow_rotation)
 
 
-def _ctrl_phase_shift(phi, wire, subspace, control_data):
-    r"""Implement a ((multi-)controlled) phase shift on the specified subspace of a
-    target qubit/wire.
-
-    Args:
-        phi (float): Phase shift angle
-        wire (~.Wires): Target wire
-        subspace (int): which subspace of the target wire the phase shift is applied to. 0 or 1.
-        control_data (tuple): A three-tuple containing control wires, control values, and work
-            wires of the ((multi-)controlled) phase shift gate.
-
-    Returns:
-        list[~.Operation]: The operations implementing the phase shift.
-
-    The decomposition for subspace=1 always is a simple ``PhaseShift`` gate, or its controlled
-    counterpart.
-    The decomposition of a non-controlled phase shift for subspace=0 can be achieved in two ways:
-    The first is to flip the angle of the phase shift and complementing it with a global phase,
-    so that the (diagonal of the) gate is decomposed as
-
-    .. math::
-
-        (\exp(i\phi), 1) = (1, \exp(-i\phi)) (\exp(i\phi), \exp(i\phi)).
-
-    The second is to conjugate a phase shift by Pauli-X operators on the same qubit, decomposing
-    the gate matrix as
-
-    / e^(i\phi)   0 \ -- / 0   1 \/ 1     0     \/ 0   1 \
-    \     0       1 / -- \ 1   0 /\ 0 e^(i\phi) /\ 1   0 /.
-
-    Without controls, the first approach is nicer, because global phases usually are free
-    operations. With controls, however, this approach would lead to a controlled global
-    phase, which is equivalent to a phase shift on top of the controlled phase shift required
-    anyways. For the second approach, we may use
-    `the ctrl(compute-uncompute) pattern <https://iopscience.iop.org/article/10.1088/2058-9565/aaa5cc>`__
-    to avoid controlling the Pauli-X operations, yielding two non-controlled Pauli-X gates and the
-    "main" controlled phase shift operation. We deem this decomposition to be better.
-
-    Note that such considerations will be less relevant in the new decomposition framework in
-    PennyLane.
-    """
-    control_wires, control_values, work_wires = control_data
-    if subspace == 1:
-        base_op = qml.PhaseShift(phi, wires=wire)
-        # If there are no control wires, we are dealing with the very first phase shift of
-        # the decomposition, which should be adding projectors. So subspace should have been 0.
-        assert len(control_wires) > 0
-        return [qml.ctrl(base_op, control_wires, control_values, work_wires)]
-
-    if len(control_wires) == 0:
-        # Flip angle for phase_shift(subspace=0) = phase_shift(subspace=1)*global_phase
-        phase_shift = qml.PhaseShift(-phi, wires=wire)
-        return [phase_shift, qml.GlobalPhase(-phi, wires=wire)]
-    return [
-        qml.X(wire),
-        qml.ctrl(
-            qml.PhaseShift(phi, wires=wire),
-            control=control_wires,
-            control_values=control_values,
-            work_wires=work_wires,
-        ),
-        qml.X(wire),
-    ]
-
-
 class PCPhase(Operation):
     r"""PCPhase(phi, dim, wires)
     A projector-controlled phase gate.
@@ -751,6 +685,8 @@ class PCPhase(Operation):
     grad_method = "A"
     parameter_frequencies = [(2,)]
 
+    resource_keys = {"num_wires", "dim"}
+
     def generator(self) -> "qml.Hermitian":
         r"""Generator of the ``PCPhase`` operator, which is in single-parameter-form.
         The operator reads
@@ -789,15 +725,19 @@ class PCPhase(Operation):
         super().__init__(phi, wires=wires, id=id)
         self.hyperparameters["dimension"] = (dim, 2 ** len(wires))
 
+    @property
+    def resource_params(self) -> dict:
+        return {"num_wires": len(self.wires), "dim": self.hyperparameters["dimension"][0]}
+
     @staticmethod
     def compute_matrix(phi: TensorLike, dimension: tuple[int, int]) -> TensorLike:
         """Get the matrix representation of Pi-controlled phase unitary."""
         d, t = dimension
 
-        if qml.math.get_interface(phi) == "tensorflow":
-            p = qml.math.exp(1j * qml.math.cast_like(phi, 1j))
-            minus_p = qml.math.exp(-1j * qml.math.cast_like(phi, 1j))
-            zeros = qml.math.zeros_like(p)
+        if math.get_interface(phi) == "tensorflow":
+            p = math.exp(1j * math.cast_like(phi, 1j))
+            minus_p = math.exp(-1j * math.cast_like(phi, 1j))
+            zeros = math.zeros_like(p)
 
             columns = []
             for i in range(t):
@@ -806,17 +746,17 @@ class PCPhase(Operation):
                     if i < d
                     else [minus_p if j == i else zeros for j in range(t)]
                 )
-            r = qml.math.stack(columns, like="tensorflow", axis=-2)
+            r = math.stack(columns, like="tensorflow", axis=-2)
             return r
 
         arg = 1j * phi
-        prefactors = qml.math.array([1] * d + [-1] * (t - d), like=phi)
+        prefactors = math.array([1] * d + [-1] * (t - d), like=phi)
 
-        if qml.math.ndim(arg) == 0:
-            return qml.math.diag(qml.math.exp(arg * prefactors))
+        if math.ndim(arg) == 0:
+            return math.diag(math.exp(arg * prefactors))
 
-        diags = qml.math.exp(qml.math.outer(arg, prefactors))
-        return qml.math.stack([qml.math.diag(d) for d in diags])
+        diags = math.exp(math.outer(arg, prefactors))
+        return math.stack([math.diag(d) for d in diags])
 
     @staticmethod
     def compute_eigvals(*params: TensorLike, **hyperparams) -> TensorLike:
@@ -824,24 +764,24 @@ class PCPhase(Operation):
         phi = params[0]
         d, t = hyperparams["dimension"]
 
-        if qml.math.get_interface(phi) == "tensorflow":
-            phase = qml.math.exp(1j * qml.math.cast_like(phi, 1j))
-            minus_phase = qml.math.exp(-1j * qml.math.cast_like(phi, 1j))
+        if math.get_interface(phi) == "tensorflow":
+            phase = math.exp(1j * math.cast_like(phi, 1j))
+            minus_phase = math.exp(-1j * math.cast_like(phi, 1j))
             return stack_last([phase if index < d else minus_phase for index in range(t)])
 
         arg = 1j * phi
-        prefactors = qml.math.array([1] * d + [-1] * (t - d), like=phi)
+        prefactors = math.array([1] * d + [-1] * (t - d), like=phi)
 
-        if qml.math.ndim(phi) == 0:
+        if math.ndim(phi) == 0:
             product = arg * prefactors
         else:
-            product = qml.math.outer(arg, prefactors)
-        return qml.math.exp(product)
+            product = math.outer(arg, prefactors)
+        return math.exp(product)
 
     @staticmethod
     def compute_decomposition(
         *params: TensorLike, wires: WiresLike, **hyperparams
-    ) -> list["qml.operation.Operator"]:
+    ) -> list[Operator]:
         r"""Representation of the PCPhase operator as a product of other operators (static method).
 
         Args:
@@ -939,54 +879,14 @@ class PCPhase(Operation):
         3: ──GlobalPhase(1.23)─────────────────╰Rϕ(-2.46)─┤
 
         """
-        phi = params[0]
-        dim, _ = hyperparams["dimension"]
-
-        # Unused wires of PCPhase can be used as work wires for controlled phase shifts
-        work_wires = list(wires)
-        # Use one more bit than there are wires, according to flipping all relevant bits for the
-        # projector decomposition, or a global phase on the gate level. Afterwards, we have
-        # dim <= 2**len(wires)=2**(n-1), which is a requirement of decomp_int_to_powers_of_two.
-        flipped, *powers_of_two = decomp_int_to_powers_of_two(dim, n=len(wires) + 1)
-        sigma = (-1) ** flipped
-        # Overall global phase to implement I_N part of the generator
-        decomp = [qml.GlobalPhase(sigma * phi)]
-        # If in flipped (sigma=-1) mode, reverse the sign of all coefficients
-        powers_of_two = [sigma * val for val in powers_of_two]
-
-        phi = 2 * sigma * phi
-
-        control_wires = []
-        control_values = []
-        for i, (c_i, wire) in enumerate(zip(powers_of_two, wires, strict=True)):
-            # Remove newly considered wire from work_wires
-            work_wires = work_wires[1:]
-            if c_i != 0:
-                # Projector with rank 2**(n-1-i) needs to be added/subtracted
-                subspace = int(c_i < 0)  # If c_i<0, target |1> subspace, else target |0> subspace
-                if flipped:  # Flip subspace if in flipped (sigma=-1) mode
-                    subspace = 1 - subspace
-                control_data = (control_wires, control_values, work_wires)
-                # Here c_i actually controls whether we add/subtract the projector
-                decomp.extend(_ctrl_phase_shift(c_i * phi, wire, subspace, control_data))
-
-            # Now that the effect on the current wire is fixed, we will only add more fine-grained
-            # projectors (i.e. smaller rank), and they will be controlled on the current wire
-            control_wires.append(wire)
-            # The control value to be used on the current wire (it will be the same for all
-            # subsequent operations) depends on whether we add or subtract next)
-            d_i = next(iter(val for val in powers_of_two[i + 1 :] if val != 0), None)
-            # If we add next, control into the |1> subspace, otherwise into |0> subspace
-            # The control value is modified both if we are in the flipped global phase mode and if
-            # the current loop iteration did not add a bit string/projector/gate (c_i==0)
-            next_cval = d_i == 1
-            if c_i == 0:
-                next_cval = not next_cval
-            if flipped:
-                next_cval = not next_cval
-            control_values.append(next_cval)
+        with queuing.AnnotatedQueue() as q:
+            _decompose_pcphase(*params, wires=wires, **hyperparams)
+
+        if queuing.QueuingManager.recording():
+            for op in q.queue:
+                queuing.apply(op)
 
-        return decomp
+        return q.queue
 
     def adjoint(self) -> "PCPhase":
         """Computes the adjoint of the operator."""
@@ -994,7 +894,7 @@ class PCPhase(Operation):
         dim, _ = self.hyperparameters["dimension"]
         return PCPhase(-1 * phi, dim=dim, wires=self.wires)
 
-    def pow(self, z: Union[int, float]) -> list["qml.operation.Operator"]:
+    def pow(self, z: Union[int, float]) -> list[Operator]:
         """Computes the operator raised to z."""
         phi = self.parameters[0]
         dim, _ = self.hyperparameters["dimension"]
@@ -1020,6 +920,187 @@ class PCPhase(Operation):
         return super().label(decimals=decimals, base_label=base_label or "∏_ϕ", cache=cache)
 
 
+def _ctrl_phase_shift_resource(subspace, n_control_wires, n_zero_control_values, n_work_wires):
+    if n_control_wires == 0:
+        return {qml.PhaseShift: 1}
+    return {
+        controlled_resource_rep(
+            qml.PhaseShift,
+            {},
+            num_control_wires=n_control_wires,
+            num_zero_control_values=n_zero_control_values,
+            num_work_wires=n_work_wires,
+        ): 1,
+        qml.X: 2 * (1 - subspace),
+    }
+
+
+def _ctrl_phase_shift(
+    phi, target_wire, subspace, control_wires, control_values, work_wires
+):  # pylint: disable=too-many-arguments
+    r"""Implement a ((multi-)controlled) phase shift on the specified subspace of a
+    target qubit/wire.
+
+    Args:
+        phi (float): Phase shift angle
+        target_wire (~.Wires): the target wire to apply phase shift to.
+        subspace (int): which subspace of the target wire the phase shift is applied to. 0 or 1.
+        control_wires (WiresLike): the control wires
+        control_values (Iterable[bool | int]): the control values.
+        work_wires (WiresLike): the work wires
+
+    Returns:
+        float: any global phase produced in the process.
+
+    The decomposition for subspace=1 always is a simple ``PhaseShift`` gate, or its controlled
+    counterpart. The decomposition of a non-controlled phase shift for subspace=0 can be achieved
+    in two ways: The first is to flip the angle of the phase shift and complementing it with a
+    global phase, so that the (diagonal of the) gate is decomposed as
+
+    .. math::
+
+        (\exp(i\phi), 1) = (1, \exp(-i\phi)) (\exp(i\phi), \exp(i\phi)).
+
+    The second is to conjugate a phase shift by Pauli-X operators on the same qubit, decomposing
+    the gate matrix as
+
+    / e^(i\phi)   0 \ -- / 0   1 \/ 1     0     \/ 0   1 \
+    \     0       1 / -- \ 1   0 /\ 0 e^(i\phi) /\ 1   0 /.
+
+    Without controls, the first approach is nicer, because global phases usually are free
+    operations. With controls, however, this approach would lead to a controlled global
+    phase, which is equivalent to a phase shift on top of the controlled phase shift. For the
+    second approach, we may use
+    `the ctrl(compute-uncompute) pattern <https://iopscience.iop.org/article/10.1088/2058-9565/aaa5cc>`__
+    to avoid controlling the Pauli-X operations, yielding two non-controlled Pauli-X gates and the
+    "main" controlled phase shift operation. We deem this decomposition to be better.
+
+    """
+
+    if subspace == 1:
+        # If there are no control wires, we are dealing with the very first phase shift of
+        # the decomposition, which should be adding projectors. So subspace should have been 0.
+        assert len(control_wires) > 0
+        qml.ctrl(
+            qml.PhaseShift(phi, wires=target_wire),
+            control=control_wires,
+            control_values=control_values,
+            work_wires=work_wires,
+        )
+        return 0.0
+
+    if len(control_wires) == 0:
+        # Flip angle for phase_shift(subspace=0) = phase_shift(subspace=1)*global_phase
+        qml.PhaseShift(-phi, wires=target_wire)
+        return -phi
+
+    qml.X(target_wire)
+    qml.ctrl(
+        qml.PhaseShift(phi, wires=target_wire),
+        control=control_wires,
+        control_values=control_values,
+        work_wires=work_wires,
+    )
+    qml.X(target_wire)
+    return 0.0
+
+
+def _decompose_pcphase_resource(num_wires, dim):
+    """Decompose the PCPhase operation into controlled phase shifts and Pauli-X gates."""
+
+    gate_count = Counter()
+    flipped, *powers_of_two = decomp_int_to_powers_of_two(dim, num_wires + 1)
+    sigma = (-1) ** flipped
+    powers_of_two = [sigma * val for val in powers_of_two]
+
+    n_zero_control_values = 0
+    for i, c_i in enumerate(powers_of_two):
+
+        if c_i != 0:
+            subspace = int(c_i < 0)
+            if flipped:
+                subspace = 1 - subspace
+            gate_count.update(
+                _ctrl_phase_shift_resource(
+                    subspace,
+                    n_control_wires=i,
+                    n_zero_control_values=n_zero_control_values,
+                    n_work_wires=num_wires - i - 1,
+                )
+            )
+
+        d_i = next(iter(val for val in powers_of_two[i + 1 :] if val != 0), None)
+        next_cval = d_i == 1
+        if c_i == 0:
+            next_cval = not next_cval
+        if flipped:
+            next_cval = not next_cval
+        if not next_cval:
+            n_zero_control_values += 1
+
+    gate_count[qml.GlobalPhase] += 1
+    return dict(gate_count)
+
+
+@register_resources(_decompose_pcphase_resource)
+def _decompose_pcphase(phi, wires, dimension):
+    """Decompose the PCPhase operation into controlled phase shifts and Pauli-X gates."""
+
+    dim, _ = dimension
+
+    # Use one more bit than there are wires, according to flipping all relevant bits for the
+    # projector decomposition, or a global phase on the gate level. Afterwards, we have
+    # dim <= 2**len(wires)=2**(n-1), which is a requirement of decomp_int_to_powers_of_two.
+    flipped, *powers_of_two = decomp_int_to_powers_of_two(dim, len(wires) + 1)
+
+    sigma = (-1) ** flipped
+
+    # Overall global phase to implement I_N part of the generator
+    global_phase = sigma * phi
+
+    phi = 2 * sigma * phi
+
+    # If in flipped (sigma=-1) mode, reverse the sign of all coefficients
+    powers_of_two = [sigma * val for val in powers_of_two]
+
+    assert len(powers_of_two) == len(wires)
+
+    control_values = []
+    for i, c_i in enumerate(powers_of_two):
+
+        if c_i != 0:
+            # Projector with rank 2**(n-1-i) needs to be added/subtracted
+            subspace = int(c_i < 0)  # If c_i < 0, target |1> subspace, else target |0> subspace
+            if flipped:  # Flip subspace if in flipped (sigma=-1) mode
+                subspace = 1 - subspace
+            global_phase += _ctrl_phase_shift(
+                c_i * phi,
+                target_wire=wires[i],
+                subspace=subspace,
+                control_wires=wires[:i],
+                control_values=control_values,
+                work_wires=wires[i + 1 :],  # Unused wires of PCPhase can be used as work wires
+            )
+
+        # The control value to be used on the current wire (it will be the same for all
+        # subsequent operations) depends on whether we add or subtract next)
+        d_i = next(iter(val for val in powers_of_two[i + 1 :] if val != 0), None)
+        # If we add next, control into the |1> subspace, otherwise into |0> subspace
+        # The control value is modified both if we are in the flipped global phase mode and if
+        # the current loop iteration did not add a bit string/projector/gate (c_i==0)
+        next_cval = d_i == 1
+        if c_i == 0:
+            next_cval = not next_cval
+        if flipped:
+            next_cval = not next_cval
+        control_values.append(next_cval)
+
+    qml.GlobalPhase(global_phase)
+
+
+add_decomps(PCPhase, _decompose_pcphase)
+
+
 class IsingXX(Operation):
     r"""
     Ising XX coupling gate
@@ -1099,26 +1180,26 @@ class IsingXX(Operation):
                 [0.0000-0.2474j, 0.0000+0.0000j, 0.0000+0.0000j, 0.9689+0.0000j]],
                dtype=torch.complex128)
         """
-        c = qml.math.cos(phi / 2)
-        s = qml.math.sin(phi / 2)
+        c = math.cos(phi / 2)
+        s = math.sin(phi / 2)
 
-        eye = qml.math.eye(4, like=phi)
-        rev_eye = qml.math.convert_like(np.eye(4)[::-1].copy(), phi)
-        if qml.math.get_interface(phi) == "tensorflow":
-            c = qml.math.cast_like(c, 1j)
-            s = qml.math.cast_like(s, 1j)
-            eye = qml.math.cast_like(eye, 1j)
-            rev_eye = qml.math.cast_like(rev_eye, 1j)
+        eye = math.eye(4, like=phi)
+        rev_eye = math.convert_like(np.eye(4)[::-1].copy(), phi)
+        if math.get_interface(phi) == "tensorflow":
+            c = math.cast_like(c, 1j)
+            s = math.cast_like(s, 1j)
+            eye = math.cast_like(eye, 1j)
+            rev_eye = math.cast_like(rev_eye, 1j)
 
         # The following avoids casting an imaginary quantity to reals when backpropagating
         js = -1j * s
-        if qml.math.ndim(phi) == 0:
+        if math.ndim(phi) == 0:
             return c * eye + js * rev_eye
 
-        return qml.math.tensordot(c, eye, axes=0) + qml.math.tensordot(js, rev_eye, axes=0)
+        return math.tensordot(c, eye, axes=0) + math.tensordot(js, rev_eye, axes=0)
 
     @staticmethod
-    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list["qml.operation.Operator"]:
+    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list[Operator]:
         r"""Representation of the operator as a product of other operators (static method). :
 
         .. math:: O = O_1 O_2 \dots O_n.
@@ -1150,7 +1231,7 @@ class IsingXX(Operation):
         (phi,) = self.parameters
         return IsingXX(-phi, wires=self.wires)
 
-    def pow(self, z: Union[int, float]) -> list["qml.operation.Operator"]:
+    def pow(self, z: Union[int, float]) -> list[Operator]:
         return [IsingXX(self.data[0] * z, wires=self.wires)]
 
     def simplify(self) -> "IsingXX":
@@ -1234,7 +1315,7 @@ class IsingYY(Operation):
         return {}
 
     @staticmethod
-    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list["qml.operation.Operator"]:
+    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list[Operator]:
         r"""Representation of the operator as a product of other operators (static method). :
 
         .. math:: O = O_1 O_2 \dots O_n.
@@ -1285,16 +1366,16 @@ class IsingYY(Operation):
                 [0.0000+0.0000j, 0.0000-0.2474j, 0.9689+0.0000j, 0.0000+0.0000j],
                 [0.0000+0.2474j, 0.0000+0.0000j, 0.0000+0.0000j, 0.9689+0.0000j]])
         """
-        c = qml.math.cos(phi / 2)
-        s = qml.math.sin(phi / 2)
+        c = math.cos(phi / 2)
+        s = math.sin(phi / 2)
 
-        if qml.math.get_interface(phi) == "tensorflow":
-            c = qml.math.cast_like(c, 1j)
-            s = qml.math.cast_like(s, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            c = math.cast_like(c, 1j)
+            s = math.cast_like(s, 1j)
 
         js = 1j * s
-        r_term = qml.math.cast_like(
-            qml.math.array(
+        r_term = math.cast_like(
+            math.array(
                 [
                     [0.0, 0.0, 0.0, 1.0],
                     [0.0, 0.0, -1.0, 0.0],
@@ -1305,16 +1386,16 @@ class IsingYY(Operation):
             ),
             1j,
         )
-        if qml.math.ndim(phi) == 0:
-            return c * qml.math.cast_like(qml.math.eye(4, like=c), c) + js * r_term
+        if math.ndim(phi) == 0:
+            return c * math.cast_like(math.eye(4, like=c), c) + js * r_term
 
-        return qml.math.tensordot(c, np.eye(4), axes=0) + qml.math.tensordot(js, r_term, axes=0)
+        return math.tensordot(c, np.eye(4), axes=0) + math.tensordot(js, r_term, axes=0)
 
     def adjoint(self) -> "IsingYY":
         (phi,) = self.parameters
         return IsingYY(-phi, wires=self.wires)
 
-    def pow(self, z: Union[int, float]) -> list["qml.operation.Operator"]:
+    def pow(self, z: Union[int, float]) -> list[Operator]:
         return [IsingYY(self.data[0] * z, wires=self.wires)]
 
     def simplify(self) -> "IsingYY":
@@ -1450,22 +1531,22 @@ class IsingZZ(Operation):
                 [0.0000+0.0000j, 0.0000+0.0000j, 0.9689+0.2474j, 0.0000+0.0000j],
                 [0.0000+0.0000j, 0.0000+0.0000j, 0.0000+0.0000j, 0.9689-0.2474j]])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            p = qml.math.exp(-0.5j * qml.math.cast_like(phi, 1j))
-            if qml.math.ndim(p) == 0:
-                return qml.math.diag([p, qml.math.conj(p), qml.math.conj(p), p])
+        if math.get_interface(phi) == "tensorflow":
+            p = math.exp(-0.5j * math.cast_like(phi, 1j))
+            if math.ndim(p) == 0:
+                return math.diag([p, math.conj(p), math.conj(p), p])
 
-            diags = stack_last([p, qml.math.conj(p), qml.math.conj(p), p])
-            return diags[:, :, np.newaxis] * qml.math.cast_like(qml.math.eye(4, like=diags), diags)
+            diags = stack_last([p, math.conj(p), math.conj(p), p])
+            return diags[:, :, np.newaxis] * math.cast_like(math.eye(4, like=diags), diags)
 
-        signs = qml.math.array([1, -1, -1, 1], like=phi)
+        signs = math.array([1, -1, -1, 1], like=phi)
         arg = -0.5j * phi
 
-        if qml.math.ndim(arg) == 0:
-            return qml.math.diag(qml.math.exp(arg * signs))
+        if math.ndim(arg) == 0:
+            return math.diag(math.exp(arg * signs))
 
-        diags = qml.math.exp(qml.math.outer(arg, signs))
-        return diags[:, :, np.newaxis] * qml.math.cast_like(qml.math.eye(4, like=diags), diags)
+        diags = math.exp(math.outer(arg, signs))
+        return diags[:, :, np.newaxis] * math.cast_like(math.eye(4, like=diags), diags)
 
     @staticmethod
     def compute_eigvals(phi: TensorLike) -> TensorLike:  # pylint: disable=arguments-differ
@@ -1494,22 +1575,22 @@ class IsingZZ(Operation):
         >>> qml.IsingZZ.compute_eigvals(torch.tensor(0.5))
         tensor([0.9689-0.2474j, 0.9689+0.2474j, 0.9689+0.2474j, 0.9689-0.2474j])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phase = qml.math.exp(-0.5j * qml.math.cast_like(phi, 1j))
-            return stack_last([phase, qml.math.conj(phase), qml.math.conj(phase), phase])
+        if math.get_interface(phi) == "tensorflow":
+            phase = math.exp(-0.5j * math.cast_like(phi, 1j))
+            return stack_last([phase, math.conj(phase), math.conj(phase), phase])
 
-        prefactors = qml.math.array([-0.5j, 0.5j, 0.5j, -0.5j], like=phi)
-        if qml.math.ndim(phi) == 0:
+        prefactors = math.array([-0.5j, 0.5j, 0.5j, -0.5j], like=phi)
+        if math.ndim(phi) == 0:
             product = phi * prefactors
         else:
-            product = qml.math.outer(phi, prefactors)
-        return qml.math.exp(product)
+            product = math.outer(phi, prefactors)
+        return math.exp(product)
 
     def adjoint(self) -> "IsingZZ":
         (phi,) = self.parameters
         return IsingZZ(-phi, wires=self.wires)
 
-    def pow(self, z: Union[int, float]) -> list["qml.operation.Operator"]:
+    def pow(self, z: Union[int, float]) -> list[Operator]:
         return [IsingZZ(self.data[0] * z, wires=self.wires)]
 
     def simplify(self) -> "IsingZZ":
@@ -1611,7 +1692,7 @@ class IsingXY(Operation):
         return {}
 
     @staticmethod
-    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list["qml.operation.Operator"]:
+    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list[Operator]:
         r"""Representation of the operator as a product of other operators (static method). :
 
         .. math:: O = O_1 O_2 \dots O_n.
@@ -1665,16 +1746,16 @@ class IsingXY(Operation):
                [0.        +0.j        , 0.        +0.24740396j,        0.96891242+0.j        , 0.        +0.j        ],
                [0.        +0.j        , 0.        +0.j        ,        0.        +0.j        , 1.        +0.j        ]])
         """
-        c = qml.math.cos(phi / 2)
-        s = qml.math.sin(phi / 2)
+        c = math.cos(phi / 2)
+        s = math.sin(phi / 2)
 
-        if qml.math.get_interface(phi) == "tensorflow":
-            c = qml.math.cast_like(c, 1j)
-            s = qml.math.cast_like(s, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            c = math.cast_like(c, 1j)
+            s = math.cast_like(s, 1j)
 
         js = 1j * s
-        off_diag = qml.math.cast_like(
-            qml.math.array(
+        off_diag = math.cast_like(
+            math.array(
                 [
                     [0.0, 0.0, 0.0, 0.0],
                     [0.0, 0.0, 1.0, 0.0],
@@ -1685,12 +1766,12 @@ class IsingXY(Operation):
             ),
             1j,
         )
-        if qml.math.ndim(phi) == 0:
-            return qml.math.diag([1, c, c, 1]) + js * off_diag
+        if math.ndim(phi) == 0:
+            return math.diag([1, c, c, 1]) + js * off_diag
 
-        ones = qml.math.ones_like(c)
+        ones = math.ones_like(c)
         diags = stack_last([ones, c, c, ones])[:, :, np.newaxis]
-        return diags * np.eye(4) + qml.math.tensordot(js, off_diag, axes=0)
+        return diags * np.eye(4) + math.tensordot(js, off_diag, axes=0)
 
     @staticmethod
     def compute_eigvals(phi: TensorLike) -> TensorLike:  # pylint: disable=arguments-differ
@@ -1719,20 +1800,20 @@ class IsingXY(Operation):
         >>> qml.IsingXY.compute_eigvals(0.5)
         array([0.96891242+0.24740396j, 0.96891242-0.24740396j,       1.        +0.j        , 1.        +0.j        ])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            phi = math.cast_like(phi, 1j)
 
         signs = np.array([1, -1, 0, 0])
-        if qml.math.ndim(phi) == 0:
-            return qml.math.exp(0.5j * phi * signs)
+        if math.ndim(phi) == 0:
+            return math.exp(0.5j * phi * signs)
 
-        return qml.math.exp(qml.math.tensordot(0.5j * phi, signs, axes=0))
+        return math.exp(math.tensordot(0.5j * phi, signs, axes=0))
 
     def adjoint(self) -> "IsingXY":
         (phi,) = self.parameters
         return IsingXY(-phi, wires=self.wires)
 
-    def pow(self, z: Union[int, float]) -> list["qml.operation.Operator"]:
+    def pow(self, z: Union[int, float]) -> list[Operator]:
         return [IsingXY(self.data[0] * z, wires=self.wires)]
 
     def simplify(self) -> "IsingXY":
@@ -1809,7 +1890,7 @@ class PSWAP(Operation):
         return {}
 
     @staticmethod
-    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list["qml.operation.Operator"]:
+    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list[Operator]:
         r"""Representation of the operator as a product of other operators (static method). :
 
         .. math:: O = O_1 O_2 \dots O_n.
@@ -1860,14 +1941,14 @@ class PSWAP(Operation):
               [0.        +0.j, 0.87758256+0.47942554j, 0.        +0.j        , 0.        +0.j],
               [0.        +0.j, 0.        +0.j        , 0.        +0.j        , 1.        +0.j]])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            phi = math.cast_like(phi, 1j)
 
-        e = qml.math.exp(1j * phi)
-        zero = qml.math.zeros_like(phi)
-        one = qml.math.ones_like(phi)
+        e = math.exp(1j * phi)
+        zero = math.zeros_like(phi)
+        one = math.ones_like(phi)
 
-        return qml.math.stack(
+        return math.stack(
             [
                 stack_last([one, zero, zero, zero]),
                 stack_last([zero, zero, e, zero]),
@@ -1904,12 +1985,12 @@ class PSWAP(Operation):
         >>> qml.PSWAP.compute_eigvals(0.5)
         array([ 1.        +0.j        ,  1.        +0.j,       -0.87758256-0.47942554j,  0.87758256+0.47942554j])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            phi = math.cast_like(phi, 1j)
 
-        e = qml.math.exp(1j * phi)
-        one = qml.math.ones_like(phi)
-        return qml.math.transpose(qml.math.stack([one, one, -e, e]))
+        e = math.exp(1j * phi)
+        one = math.ones_like(phi)
+        return math.transpose(math.stack([one, one, -e, e]))
 
     def adjoint(self) -> "PSWAP":
         (phi,) = self.parameters
@@ -2025,14 +2106,14 @@ class CPhaseShift00(Operation):
                     [0.0000+0.0000j, 0.0+0.0j, 1.0+0.0j, 0.0+0.0j],
                     [0.0000+0.0000j, 0.0+0.0j, 0.0+0.0j, 1.0+0.0j]])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            phi = math.cast_like(phi, 1j)
 
-        exp_part = qml.math.exp(1j * phi)
+        exp_part = math.exp(1j * phi)
 
-        if qml.math.ndim(phi) > 0:
-            ones = qml.math.ones_like(exp_part)
-            zeros = qml.math.zeros_like(exp_part)
+        if math.ndim(phi) > 0:
+            ones = math.ones_like(exp_part)
+            zeros = math.zeros_like(exp_part)
             matrix = [
                 [exp_part, zeros, zeros, zeros],
                 [zeros, ones, zeros, zeros],
@@ -2040,9 +2121,9 @@ class CPhaseShift00(Operation):
                 [zeros, zeros, zeros, ones],
             ]
 
-            return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
+            return math.stack([stack_last(row) for row in matrix], axis=-2)
 
-        return qml.math.diag([exp_part, 1, 1, 1])
+        return math.diag([exp_part, 1, 1, 1])
 
     @staticmethod
     def compute_eigvals(phi: TensorLike) -> TensorLike:  # pylint: disable=arguments-differ
@@ -2071,15 +2152,15 @@ class CPhaseShift00(Operation):
         >>> qml.CPhaseShift00.compute_eigvals(torch.tensor(0.5))
         tensor([0.8776+0.4794j, 1.0000+0.0000j, 1.0000+0.0000j, 1.0000+0.0000j])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            phi = math.cast_like(phi, 1j)
 
-        exp_part = qml.math.exp(1j * phi)
-        ones = qml.math.ones_like(exp_part)
+        exp_part = math.exp(1j * phi)
+        ones = math.ones_like(exp_part)
         return stack_last([exp_part, ones, ones, ones])
 
     @staticmethod
-    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list["qml.operation.Operator"]:
+    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list[Operator]:
         r"""Representation of the operator as a product of other operators (static method). :
 
         .. math:: O = O_1 O_2 \dots O_n.
@@ -2245,14 +2326,14 @@ class CPhaseShift01(Operation):
                     [0.0+0.0j, 0.0000+0.0000j, 1.0+0.0j, 0.0+0.0j],
                     [0.0+0.0j, 0.0000+0.0000j, 0.0+0.0j, 1.0+0.0j]])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            phi = math.cast_like(phi, 1j)
 
-        exp_part = qml.math.exp(1j * phi)
+        exp_part = math.exp(1j * phi)
 
-        if qml.math.ndim(phi) > 0:
-            ones = qml.math.ones_like(exp_part)
-            zeros = qml.math.zeros_like(exp_part)
+        if math.ndim(phi) > 0:
+            ones = math.ones_like(exp_part)
+            zeros = math.zeros_like(exp_part)
             matrix = [
                 [ones, zeros, zeros, zeros],
                 [zeros, exp_part, zeros, zeros],
@@ -2260,9 +2341,9 @@ class CPhaseShift01(Operation):
                 [zeros, zeros, zeros, ones],
             ]
 
-            return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
+            return math.stack([stack_last(row) for row in matrix], axis=-2)
 
-        return qml.math.diag([1, exp_part, 1, 1])
+        return math.diag([1, exp_part, 1, 1])
 
     @staticmethod
     def compute_eigvals(phi: TensorLike) -> TensorLike:  # pylint: disable=arguments-differ
@@ -2291,15 +2372,15 @@ class CPhaseShift01(Operation):
         >>> qml.CPhaseShift01.compute_eigvals(torch.tensor(0.5))
         tensor([1.0000+0.0000j, 0.8776+0.4794j, 1.0000+0.0000j, 1.0000+0.0000j])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            phi = math.cast_like(phi, 1j)
 
-        exp_part = qml.math.exp(1j * phi)
-        ones = qml.math.ones_like(exp_part)
+        exp_part = math.exp(1j * phi)
+        ones = math.ones_like(exp_part)
         return stack_last([ones, exp_part, ones, ones])
 
     @staticmethod
-    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list["qml.operation.Operator"]:
+    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list[Operator]:
         r"""Representation of the operator as a product of other operators (static method). :
 
         .. math:: O = O_1 O_2 \dots O_n.
@@ -2455,14 +2536,14 @@ class CPhaseShift10(Operation):
                     [0.0+0.0j, 0.0+0.0j, 0.8776+0.4794j, 0.0+0.0j],
                     [0.0+0.0j, 0.0+0.0j, 0.0000+0.0000j, 1.0+0.0j]])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            phi = math.cast_like(phi, 1j)
 
-        exp_part = qml.math.exp(1j * phi)
+        exp_part = math.exp(1j * phi)
 
-        if qml.math.ndim(phi) > 0:
-            ones = qml.math.ones_like(exp_part)
-            zeros = qml.math.zeros_like(exp_part)
+        if math.ndim(phi) > 0:
+            ones = math.ones_like(exp_part)
+            zeros = math.zeros_like(exp_part)
             matrix = [
                 [ones, zeros, zeros, zeros],
                 [zeros, ones, zeros, zeros],
@@ -2470,9 +2551,9 @@ class CPhaseShift10(Operation):
                 [zeros, zeros, zeros, ones],
             ]
 
-            return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
+            return math.stack([stack_last(row) for row in matrix], axis=-2)
 
-        return qml.math.diag([1, 1, exp_part, 1])
+        return math.diag([1, 1, exp_part, 1])
 
     @staticmethod
     def compute_eigvals(phi: TensorLike) -> TensorLike:  # pylint: disable=arguments-differ
@@ -2501,15 +2582,15 @@ class CPhaseShift10(Operation):
         >>> qml.CPhaseShift10.compute_eigvals(torch.tensor(0.5))
         tensor([1.0000+0.0000j, 1.0000+0.0000j, 0.8776+0.4794j, 1.0000+0.0000j])
         """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
+        if math.get_interface(phi) == "tensorflow":
+            phi = math.cast_like(phi, 1j)
 
-        exp_part = qml.math.exp(1j * phi)
-        ones = qml.math.ones_like(exp_part)
+        exp_part = math.exp(1j * phi)
+        ones = math.ones_like(exp_part)
         return stack_last([ones, ones, exp_part, ones])
 
     @staticmethod
-    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list["qml.operation.Operator"]:
+    def compute_decomposition(phi: TensorLike, wires: WiresLike) -> list[Operator]:
         r"""Representation of the operator as a product of other operators (static method). :
 
         .. math:: O = O_1 O_2 \dots O_n.


Test patch:
diff --git a/tests/ops/qubit/test_parametric_ops.py b/tests/ops/qubit/test_parametric_ops.py
index 1913239ba..588d7801a 100644
--- a/tests/ops/qubit/test_parametric_ops.py
+++ b/tests/ops/qubit/test_parametric_ops.py
@@ -678,6 +678,31 @@ class TestDecompositions:
         for expected_mat, decomp_mat in zip(expected_mats, decomp_mats):
             assert np.allclose(expected_mat, decomp_mat)
 
+    @pytest.mark.unit
+    @pytest.mark.parametrize("dim, wires", two_wire_pcphases + five_wire_pcphases + other_pcphases)
+    def test_pcphase_decomposition_new(self, dim, wires):
+        """Tests the PCPhase decomposition rules."""
+
+        op = qml.PCPhase(np.random.random(), dim, wires=wires)
+        for rule in qml.list_decomps(qml.PCPhase):
+            _test_decomposition_rule(op, rule)
+
+    @pytest.mark.integration
+    @pytest.mark.usefixtures("enable_graph_decomposition")
+    def test_pcphase_decomposition_graph(self):
+        """Tests that the PCPhase can be decomposed all the way down."""
+
+        tape = qml.tape.QuantumScript(
+            [qml.PCPhase(0.123, 12, wires=[0, 1, 2, 3, 4, 5, 6, 7, 8])], []
+        )
+        expected_matrix = qml.matrix(tape, wire_order=[0, 1, 2, 3, 4, 5, 6, 7, 8])
+
+        [decomp], _ = qml.transforms.decompose(
+            tape, gate_set={qml.RX, qml.RY, qml.RZ, qml.CNOT, qml.X, qml.Toffoli, qml.GlobalPhase}
+        )
+        mat = qml.matrix(decomp, wire_order=[0, 1, 2, 3, 4, 5, 6, 7, 8])
+        assert qml.math.allclose(mat, expected_matrix)
+
     @pytest.mark.parametrize("phi", [-0.1, 0.2, 0.5])
     @pytest.mark.parametrize(
         "cphase_op,lam_pos",