Repository: PennyLaneAI/pennylane
Pull request #7333 (2025-04-24 13:28:44)

[WIP] [Labs] phase polynomial intermediate representation
Building on https://github.com/PennyLaneAI/pennylane/pull/7229
Implementing https://www.notion.so/xanaduai/25-Q1-Intermediate-Representation-zoo-177bc6bd176481d6b67ad36143e61e84?pvs=4#177bc6bd176480c0a9d5ce8c4ef01a26

Create Parity table IR prototype [sc-88507]

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/7333?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 99.64%. Comparing base [(`2469b32`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/2469b321c3f601e8a4640f18644304b1e9d1d8d1?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`cb35357`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/cb35357e8a4362b867309b3b78442780847e8cae?dropdown=coverage&el=desc&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                 @@
##           phase-polynomials    #7333   +/-   ##
==================================================
  Coverage              99.64%   99.64%           
==================================================
  Files                    514      514           
  Lines                  50088    50088           
==================================================
  Hits                   49912    49912           
  Misses                   176      176           
```

</details>

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

<details><summary> :rocket: New features to boost your workflow: </summary>

- :snowflake: [Test Analytics](https://docs.codecov.com/docs/test-analytics): Detect flaky tests, report on failures, and find test suite problems.
</details>
Qottmann: > Are we adding a page about this IR to the Compilation hub?

@AntonNI8 yes! :) 

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 83a01c32a..6689a7702 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -150,6 +150,12 @@
 
 <h3>Improvements 🛠</h3>
 
+* The decomposition of `DiagonalQubitUnitary` has been updated to a recursive decomposition
+  into a smaller `DiagonalQubitUnitary` and a `SelectPauliRot` operation. This is a known
+  decomposition [Theorem 7 in Shende et al.](https://arxiv.org/abs/quant-ph/0406176)
+  that contains fewer gates than the previous decomposition.
+  [(#7370)](https://github.com/PennyLaneAI/pennylane/pull/7370)
+
 * PennyLane supports `JAX` version 0.6.0.
   [(#7299)](https://github.com/PennyLaneAI/pennylane/pull/7299)
 
@@ -235,11 +241,17 @@
 
 <h3>Labs: a place for unified and rapid prototyping of research software 🧪</h3>
 
-* A :func:`parity_matrix <pennylane.labs.intermediate_reps.parity_matrix>` function is now available
-  in :mod:`pennylane.labs.intermediate_reps <pennylane.labs.intermediate_reps>`.
-  It allows computation of the parity matrix of a CNOT circuit; an efficient intermediate representation.
-  It is important for CNOT routing algorithms and other quantum compilation routines.
+
+* A new module :mod:`pennylane.labs.intermediate_reps <pennylane.labs.intermediate_reps>`
+  provides functionality to compute intermediate representations for particular circuits.
+  :func:`parity_matrix <pennylane.labs.intermediate_reps.parity_matrix>` computes
+  the parity matrix intermediate representation for CNOT circuits.
+  :func:`phase_polynomial <pennylane.labs.intermediate_reps.phase_polynomial>` computes
+  the phase polynomial intermediate representation for {CNOT, RZ} circuits.
+  These efficient intermediate representations are important
+  for CNOT routing algorithms and other quantum compilation routines.
   [(#7229)](https://github.com/PennyLaneAI/pennylane/pull/7229)
+  [(#7333)](https://github.com/PennyLaneAI/pennylane/pull/7333)
 
 
 <h3>Breaking changes 💔</h3>
diff --git a/pennylane/labs/intermediate_reps/__init__.py b/pennylane/labs/intermediate_reps/__init__.py
index e7042ac31..223cb5947 100644
--- a/pennylane/labs/intermediate_reps/__init__.py
+++ b/pennylane/labs/intermediate_reps/__init__.py
@@ -23,8 +23,10 @@ Intermediate representations (IRs) are alternative representations of quantum ci
     :toctree: api
 
     ~parity_matrix
+    ~phase_polynomial
 
 
 """
 
 from .parity_matrix import parity_matrix
+from .phase_polynomial import phase_polynomial
diff --git a/pennylane/labs/intermediate_reps/phase_polynomial.py b/pennylane/labs/intermediate_reps/phase_polynomial.py
new file mode 100644
index 000000000..14b2ed170
--- /dev/null
+++ b/pennylane/labs/intermediate_reps/phase_polynomial.py
@@ -0,0 +1,194 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Phase polynomial intermediate representation"""
+from typing import Sequence
+
+import numpy as np
+
+import pennylane as qml
+
+
+def phase_polynomial(
+    circ: qml.tape.QuantumScript, wire_order: Sequence = None, verbose: bool = False
+):
+    r"""
+    Phase polynomial intermediate representation for circuits consisting of CNOT and RZ gates.
+
+    The action of such circuits can be described by a phase polynomial :math:`p(\boldsymbol{x})` and a :func:`~parity_matrix` :math:`P` acting on a computational basis state :math:`|\boldsymbol{x}\rangle = |x_1, x_2, .., x_n\rangle` in the following way:
+
+    .. math:: U |\boldsymbol{x}\rangle = e^{i p(\boldsymbol{x})} |P \boldsymbol{x}\rangle.
+
+    Since the parity matrix :math:`P` is part of this description, :math:`p` and :math:`P` in conjunction are sometimes referred to as the phase polynomial intermediate representation (IR).
+
+    The phase polynomial :math:`p(\boldsymbol{x})` is described in terms of its parity table :math:`P_T` and associated angles. For this, note that
+    the action of a :class:`~RZ` gate onto a computational basis state :math:`|x\rangle` is given by
+
+    .. math:: R_Z(\theta) |x\rangle = e^{-i \frac{\theta}{2} (1 - 2x)} |x\rangle.
+
+    The parity table :math:`P_T` is made up of the `parities` :math:`\boldsymbol{x}` at the point in the circuit where the associated :class:`~RZ` gate is acting.
+    To track the impact of the gate, we thus simply collect the current parity and remember the angle.
+    Take for example the circuit ``[CNOT((0, 1)), RZ(theta, 1), CNOT((0, 1))]`` (read from left to right like a circuit diagram). We start in some arbitrary computational basis state
+    ``x = [x1, x2]``. The first CNOT is transforming the input state to ``[x1, x1 ⊕ x2]``.
+    For the action of ``RZ`` we remember the angle ``theta`` as well as the current parity ``x1 ⊕ x2`` on that wire.
+    The second CNOT gate undoes the parity change and restores the original computational basis state ``[x1, x2]``.
+
+    Hence, the parity matrix is simply the identity, but the parity table for the phase polynomial is ``P_T = [[x1 ⊕ x2]]`` (or ``[[1, 1]]``) together with the angle ``theta`` in the list of angles ``[theta]``.
+    The computation of the circuit is thus simply
+
+    .. math:: U |x_1, x_2\rangle = e^{-i \frac{\theta}{2} \left(1 - 2(x_1 \oplus x_2) \right)} |x_1, x_2\rangle
+
+    The semantics of this function is roughly given by the following implementation:
+
+    .. code-block:: python
+
+        def compute_phase_polynomial(circ, verbose=False):
+            wires = circ.wires
+            parity_matrix = np.eye(len(wires), dtype=int)
+            parity_table = []
+            angles = []
+
+            for op in circ.operations:
+
+                if op.name == "CNOT":
+                    control, target = op.wires
+                    parity_matrix[target] = (parity_matrix[target] + parity_matrix[control]) % 2
+
+                elif op.name == "RZ":
+                    angles.append(op.data[0]) # append theta_i
+                    parity_table.append(parity_matrix[op.wires[0]].copy()) # append _current_ parity (hence the copy)
+
+            return parity_matrix, np.array(parity_table).T, angles
+
+    Args:
+        circ (qml.tape.QuantumScript): Quantum circuit containing only CNOT gates.
+        wire_order (Iterable): ``wire_order`` indicating how rows and columns should be ordered. If ``None`` is provided, we take the wires of the input circuit (``circ.wires``).
+        verbose (bool): Whether or not progress should be printed during computation.
+
+    Returns:
+        tuple(np.ndarray, np.ndarray, np.ndarray): A tuple consisting of the :func:`~parity_matrix`, parity table and corresponding angles for each parity.
+
+    **Example**
+
+    We look at the circuit in Figure 1 in `arXiv:2104.00934 <https://arxiv.org/abs/2104.00934>`__.
+
+    >>> circ = qml.tape.QuantumScript([
+    ...     qml.CNOT((1, 0)),
+    ...     qml.RZ(1, 0),
+    ...     qml.CNOT((2, 0)),
+    ...     qml.RZ(2, 0),
+    ...     qml.CNOT((0, 1)),
+    ...     qml.CNOT((3, 1)),
+    ...     qml.RZ(3, 1)
+    ... ])
+    >>> print(qml.drawer.tape_text(circ, decimals=0, wire_order=range(4)))
+    0: ─╭X──RZ(1)─╭X──RZ(2)─╭●───────────┤
+    1: ─╰●────────│─────────╰X─╭X──RZ(3)─┤
+    2: ───────────╰●───────────│─────────┤
+    3: ────────────────────────╰●────────┤
+
+    The phase polynomial representation consisting of the parity matrix, parity table and associated
+    angles are computed by ``phase_polynomial``.
+
+    >>> pmat, ptab, angles = phase_polynomial(circ, wire_order=range(4))
+    >>> pmat
+    array([[1, 1, 1, 0],
+           [1, 0, 1, 1],
+           [0, 0, 1, 0],
+           [0, 0, 0, 1]])
+    >>> ptab
+    array([[1, 1, 1],
+           [1, 1, 0],
+           [0, 1, 1],
+           [0, 0, 1]])
+    >>> angles
+    array([1, 2, 3])
+
+    .. details::
+        :title: Details
+
+        We can go through explicitly reconstructing the output wavefunction.
+        First, let us compute the exact wavefunction from the circuit.
+
+        .. code-block:: python
+
+            input = np.array([1, 1, 1, 1]) # computational basis state
+
+            def comp_basis_to_wf(basis_state):
+                return qml.BasisState(np.array(basis_state), range(4)).state_vector().reshape(-1)
+
+            input_wf = comp_basis_to_wf(input)
+            output_wf = qml.matrix(circ, wire_order=range(4)) @ input_wf
+
+        The output wavefunction is given by :math:`e^{2i} * |1 1 1 1\rangle`, which we can confirm:
+
+        >>> np.allclose(output_wf, np.exp(2j) * input_wf)
+        True
+
+        Note that the action of an :class:`~RZ` gate is given by
+
+        .. math:: R_Z(\theta) |x\rangle = e^{-i \frac{\theta}{2} Z} |x\rangle = e^{-i \frac{\theta}{2} (1 - 2x)} |x\rangle
+
+        Hence, we need to convert the collected parities :math:`\boldsymbol{x}` as :math:`-(1 - 2\boldsymbol{x})/2`, accordingly. In particular, the collected phase :math:`p(x)` is given by
+
+        >>> output_phase = -(1 - 2 * ((input @ ptab) % 2))/2
+        >>> output_phase = output_phase @ angles
+
+        The final output wavefunction from the phase polynomial description is then given by the following.
+
+        >>> output_wf_re = np.exp(1j * output_phase) * comp_basis_to_wf(pmat @ input % 2)
+
+        We can compare it to the exact output wavefunction and see that they match:
+
+        >>> np.allclose(output_wf_re, output_wf)
+        True
+
+    """
+    wires = circ.wires
+
+    if wire_order is None:
+        wire_order = wires
+
+    wire_map = {wire: idx for idx, wire in enumerate(wire_order)}
+
+    parity_matrix = np.eye(len(wire_map), dtype=int)
+    parity_table = []
+    angles = []
+    i = 0
+    for op in circ.operations:
+
+        if op.name == "CNOT":
+            control, target = op.wires
+            parity_matrix[wire_map[target]] = (
+                parity_matrix[wire_map[target]] + parity_matrix[wire_map[control]]
+            ) % 2
+        elif op.name == "RZ":
+            angles.append(op.data[0])  # append theta_i
+            RZ_wire = wire_map[op.wires[0]]
+
+            # append _current_ parity (hence the copy)
+            parity_table.append(parity_matrix[RZ_wire].copy())
+        else:
+            raise TypeError(
+                f"phase_polynomial can only handle CNOT and RZ operators, received {op}"
+            )
+
+        if verbose:
+            print(f"Operator {op.name} - #{i}")
+            print(
+                f"parity matrix: \n{parity_matrix}\n parity table: \n{np.array(parity_table).T}\n angles: {angles}"
+            )
+
+        i += 1
+
+    return parity_matrix, np.array(parity_table).T, np.array(angles)
diff --git a/pennylane/ops/qubit/matrix_ops.py b/pennylane/ops/qubit/matrix_ops.py
index 8c3600ef2..23961299c 100644
--- a/pennylane/ops/qubit/matrix_ops.py
+++ b/pennylane/ops/qubit/matrix_ops.py
@@ -17,7 +17,6 @@ accept a hermitian or an unitary matrix as a parameter.
 """
 # pylint:disable=arguments-differ
 import warnings
-from itertools import product
 from typing import Optional, Union
 
 import numpy as np
@@ -480,49 +479,86 @@ class DiagonalQubitUnitary(Operation):
         Returns:
             list[Operator]: decomposition into lower level operations
 
+        Implements Theorem 7 of `Shende et al. <https://arxiv.org/abs/quant-ph/0406176>`__.
+        Decomposing a ``DiagonalQubitUnitary`` on :math:`n` wires (:math:`n>1`) yields a
+        uniformly-controlled :math:`R_Z` gate, or :class:`~.SelectPauliRot` gate, as well as a
+        ``DiagonalQubitUnitary`` on :math:`n-1` wires. For :math:`n=1` wires, the decomposition
+        yields a :class:`~.RZ` gate and a :class:`~.GlobalPhase`.
+        Resolving this recursion relationship, one would obtain :math:`n-1` ``SelectPauliRot``
+        gates with :math:`n, n-1, \dots, 1` controls each, a single ``RZ`` gate, and
+        a ``GlobalPhase``.
+
         **Example:**
 
         >>> diag = np.exp(1j * np.array([0.4, 2.1, 0.5, 1.8]))
         >>> qml.DiagonalQubitUnitary.compute_decomposition(diag, wires=[0, 1])
-        [QubitUnitary(array([[0.36235775+0.93203909j, 0.        +0.j        ],
-         [0.        +0.j        , 0.36235775+0.93203909j]]), wires=[0]),
-         RZ(1.5000000000000002, wires=[1]),
-         RZ(-0.10000000000000003, wires=[0]),
-         IsingZZ(0.2, wires=[0, 1])]
+        [SelectPauliRot(array([1.7, 1.3]), wires=[0, 1]),
+         DiagonalQubitUnitary(array([0.31532236+0.94898462j, 0.40848744+0.91276394j]), wires=[0])]
 
-        """
-        n = len(wires)
-
-        # Cast the diagonal into a complex dtype so that the logarithm works as expected
-        D_casted = qml.math.cast(D, "complex128")
-
-        phases = qml.math.real(qml.math.log(D_casted) * (-1j))
-        coeffs = _walsh_hadamard_transform(phases, n).T
-        global_phase = qml.math.exp(1j * coeffs[0])
-        # For all other gates, there is a prefactor -1/2 to be compensated.
-        coeffs = coeffs * (-2.0)
-
-        # TODO: Replace the following by a GlobalPhase gate.
-        ops = [QubitUnitary(qml.math.tensordot(global_phase, qml.math.eye(2), axes=0), wires[0])]
-        for wire0 in range(n):
-            # Single PauliZ generators correspond to the coeffs at powers of two
-            ops.append(qml.RZ(coeffs[1 << wire0], wires[n - 1 - wire0]))
-            # Double PauliZ generators correspond to the coeffs at the sum of two powers of two
-            ops.extend(
-                qml.IsingZZ(
-                    coeffs[(1 << wire0) + (1 << wire1)],
-                    [wires[n - 1 - wire0], wires[n - 1 - wire1]],
-                )
-                for wire1 in range(wire0)
-            )
+        .. details::
 
-        # Add all multi RZ gates that are not generated by single or double PauliZ generators
-        ops.extend(
-            qml.MultiRZ(c, [wires[k] for k in np.where(term)[0]])
-            for c, term in zip(coeffs, product((0, 1), repeat=n))
-            if sum(term) > 2
-        )
-        return ops
+            :title: Finding the parameters
+
+            Theorem 7 referenced above only tells us the structure of the circuit, but not the
+            parameters for the ``SelectPauliRot`` and ``DiagonalQubitUnitary`` in the decomposition.
+            In the following, we will only write out the diagonals of all gates.
+            Consider a ``DiagonalQubitUnitary`` on :math:`n` qubits that we want to decompose:
+
+            .. math::
+
+                D(\theta) = (\exp(i\theta_0), \exp(i\theta_1), \dots,
+                \exp(i\theta_{N-2}), \exp(i\theta_{N-1})).
+
+            Here, :math:`N=2^n` is the Hilbert space dimension for :math:`n` qubits, which is
+            the same as the number of parameters in :math:`D`.
+
+            A ``SelectPauliRot`` gate using ``RZ`` rotations, or multiplexed ``RZ`` rotation, using
+            the first :math:`n-1` qubits as controls and the last qubit as target, takes the form
+
+            .. math::
+
+                UCR_Z(\phi) = (\exp(-\frac{i}{2}\phi_0), \exp(\frac{i}{2}\phi_0), \dots,
+                \exp(-\frac{i}{2}\phi_{N/2-1}), \exp(\frac{i}{2}\phi_{N/2-1})),
+
+            i.e., it moves the phase of neighbouring pairs of computational basis states by
+            the same amount, but in opposite direction. There are :math:`N/2` parameters
+            in this gate.
+            Similarly, a ``DiagonalQubitUnitary`` acting on the first :math:`n-1` qubits only (the
+            ones that were controls for ``SelectPauliRot``) takes the form
+
+            .. math::
+
+                D'(\theta') = (\exp(i\theta'_0), \exp(i\theta'_0), \dots,
+                \exp(i\theta'_{N/2-1}), \exp(i\theta'_{N/2-1})).
+
+            That is, :math:`D'` moves the phase of neighbouring pairs of basis states by the same
+            amount and in the same direction. It, too, has :math:`N/2` parameters.
+            Now, we see that we can compute the rotation angles, or phases, :math:`\phi` and
+            :math:`\theta'` quite easily from the original :math:`\theta`:
+
+            .. math::
+
+                (\exp(i\theta_{2i}), \exp(i\theta_{2i+1})) &=
+                (\exp(-\frac{i}{2}\phi_i)\exp(i\theta'_i), \exp(\frac{i}{2}\phi_i)\exp(i\theta'_i))\\
+                \Rightarrow \qquad \theta'_i &=\frac{1}{2}(\theta_{2i}+\theta_{2i+1})\\
+                \phi_i &=\theta_{2i+1}-\theta_{2i}.
+
+            So the phases for the new gates arise simply as difference and average of the
+            odd-indexed and even-indexed phases.
+
+        """
+        angles = qml.math.angle(D)
+        diff = angles[..., 1::2] - angles[..., ::2]
+        mean = (angles[..., ::2] + angles[..., 1::2]) / 2
+        if len(wires) == 1:
+            return [  # Squeeze away non-broadcasting axis (there is just one angle for RZ/GPhase
+                qml.GlobalPhase(-qml.math.squeeze(mean, axis=-1), wires=wires),
+                qml.RZ(qml.math.squeeze(diff, axis=-1), wires=wires),
+            ]
+        return [  # Note that we use the first qubits as control, the reference uses the last qubits
+            qml.DiagonalQubitUnitary(np.exp(1j * mean), wires=wires[:-1]),
+            qml.SelectPauliRot(diff, control_wires=wires[:-1], target_wire=wires[-1]),
+        ]
 
     def adjoint(self) -> "DiagonalQubitUnitary":
         return DiagonalQubitUnitary(qml.math.conj(self.parameters[0]), wires=self.wires)


Test patch:
diff --git a/pennylane/labs/tests/intermediate_reps/test_phase_polynomial.py b/pennylane/labs/tests/intermediate_reps/test_phase_polynomial.py
new file mode 100644
index 000000000..287b4f69f
--- /dev/null
+++ b/pennylane/labs/tests/intermediate_reps/test_phase_polynomial.py
@@ -0,0 +1,85 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for pennylane/labs/phase_polynomials/phase_polynomial.py"""
+# pylint: disable = no-self-use
+import numpy as np
+import pytest
+
+import pennylane as qml
+from pennylane.labs.intermediate_reps import phase_polynomial
+
+# trivial example
+circ0 = qml.tape.QuantumScript([])
+pmat0 = np.eye(4)
+ptab0 = np.array([])
+angles0 = np.array([])
+
+# example from fig. 1 in https://arxiv.org/abs/2104.00934
+circ1 = qml.tape.QuantumScript(
+    [
+        qml.CNOT((1, 0)),
+        qml.RZ(1, 0),
+        qml.CNOT((2, 0)),
+        qml.RZ(2, 0),
+        qml.CNOT((0, 1)),
+        qml.CNOT((3, 1)),
+        qml.RZ(3, 1),
+    ],
+    [],
+)
+
+pmat1 = np.array([[1, 1, 1, 0], [1, 0, 1, 1], [0, 0, 1, 0], [0, 0, 0, 1]])
+ptab1 = np.array([[1, 1, 1], [1, 1, 0], [0, 1, 1], [0, 0, 1]])
+angles1 = np.array([1, 2, 3])
+
+
+class TestPhasePolynomial:
+    """Tests for qml.labs.dla.phase_polynomials.phase_polynomial"""
+
+    @pytest.mark.parametrize(
+        "circ, res",
+        (
+            (circ0, (pmat0, ptab0, angles0)),
+            (circ1, (pmat1, ptab1, angles1)),
+        ),
+    )
+    def test_computation(self, circ, res):
+        """Test phase_polynomial computes the correct parity matrix, parity table and angles"""
+
+        pmat, ptab, angles = phase_polynomial(circ, wire_order=range(4))
+        pmat_true, ptab_true, angles_true = res
+
+        assert np.allclose(pmat, pmat_true)
+        assert np.allclose(ptab, ptab_true)
+        assert np.allclose(angles, angles_true)
+
+    def test_verbose(self, capsys):
+        """Test the verbose output of phase_polynomial"""
+        _ = phase_polynomial(circ1, verbose=True)
+        captured = capsys.readouterr()
+        assert "Operator CNOT - #0" in captured.out
+        assert "Operator RZ - #1" in captured.out
+        assert "Operator CNOT - #2" in captured.out
+
+    def test_wire_order_string_wires(self):
+        """Test wire_order with string wires"""
+        wire_order_abcd = ["a", "b", "c", "d"]
+
+        (circ1_abcd,), _ = qml.map_wires(circ1, dict(enumerate(wire_order_abcd)))
+
+        pmat, ptab, angles = phase_polynomial(circ1_abcd, wire_order=["a", "b", "c", "d"])
+
+        assert np.allclose(pmat, pmat1)
+        assert np.allclose(ptab, ptab1)
+        assert np.allclose(angles, angles1)
diff --git a/tests/ops/qubit/test_matrix_ops.py b/tests/ops/qubit/test_matrix_ops.py
index c8379033b..81104fa41 100644
--- a/tests/ops/qubit/test_matrix_ops.py
+++ b/tests/ops/qubit/test_matrix_ops.py
@@ -641,26 +641,23 @@ class TestDiagonalQubitUnitary:
         decomp = qml.DiagonalQubitUnitary.compute_decomposition(D, [0])
         decomp2 = qml.DiagonalQubitUnitary(D, wires=[0]).decomposition()
 
-        ph = np.exp(3j * np.pi / 4)
         for dec in (decomp, decomp2):
             assert len(dec) == 2
-            qml.assert_equal(decomp[0], qml.QubitUnitary(np.eye(2) * ph, 0))
+            qml.assert_equal(decomp[0], qml.GlobalPhase(-3 * np.pi / 4, 0))
             qml.assert_equal(decomp[1], qml.RZ(np.pi / 2, 0))
 
     def test_decomposition_single_qubit_broadcasted(self):
         """Test that a broadcasted single-qubit DiagonalQubitUnitary is decomposed correctly."""
-        D = np.stack(
-            [[1j, -1], np.exp(1j * np.array([np.pi / 8, -np.pi / 8])), [1j, -1j], [-1, -1]]
-        )
-        angles = np.array([np.pi / 2, -np.pi / 4, -np.pi, 0])
+        D = np.exp(1j * np.pi * np.array([[1 / 2, 1], [1 / 8, -1 / 8], [1 / 2, -1 / 2], [1, 1]]))
 
         decomp = qml.DiagonalQubitUnitary.compute_decomposition(D, [0])
         decomp2 = qml.DiagonalQubitUnitary(D, wires=[0]).decomposition()
 
-        ph = [np.exp(3j * np.pi / 4), 1, 1, -1]
+        angles = np.array([1 / 2, -1 / 4, -1, 0]) * np.pi
+        global_angles = np.array([3 / 4, 0, 0, 1]) * np.pi
         for dec in (decomp, decomp2):
             assert len(dec) == 2
-            qml.assert_equal(decomp[0], qml.QubitUnitary(np.array([np.eye(2) * p for p in ph]), 0))
+            qml.assert_equal(decomp[0], qml.GlobalPhase(-global_angles, 0))
             qml.assert_equal(decomp[1], qml.RZ(angles, 0))
 
     def test_decomposition_two_qubits(self):
@@ -670,12 +667,13 @@ class TestDiagonalQubitUnitary:
         decomp = qml.DiagonalQubitUnitary.compute_decomposition(D, [0, 1])
         decomp2 = qml.DiagonalQubitUnitary(D, wires=[0, 1]).decomposition()
 
+        angles = np.array([-2, 0.5])
+        new_D = np.exp(1j * np.array([0, 3 / 4]))
+
         for dec in (decomp, decomp2):
-            assert len(dec) == 4
-            qml.assert_equal(decomp[0], qml.QubitUnitary(np.eye(2) * np.exp(0.375j), 0))
-            qml.assert_equal(decomp[1], qml.RZ(-0.75, 1))
-            qml.assert_equal(decomp[2], qml.RZ(0.75, 0))
-            qml.assert_equal(decomp[3], qml.IsingZZ(-1.25, [0, 1]))
+            assert len(dec) == 2
+            qml.assert_equal(decomp[0], qml.DiagonalQubitUnitary(new_D, wires=[0]))
+            qml.assert_equal(decomp[1], qml.SelectPauliRot(angles, [0], target_wire=1))
 
     def test_decomposition_two_qubits_broadcasted(self):
         """Test that a broadcasted two-qubit DiagonalQubitUnitary is decomposed correctly."""
@@ -684,14 +682,13 @@ class TestDiagonalQubitUnitary:
         decomp = qml.DiagonalQubitUnitary.compute_decomposition(D, [0, 1])
         decomp2 = qml.DiagonalQubitUnitary(D, wires=[0, 1]).decomposition()
 
-        angles = [[-0.75, -0.8, 0.65], [0.75, -2.4, -0.55], [-1.25, 0.4, -1.35]]
-        ph = [np.exp(1j * 0.375), np.exp(1j * 0.9), np.exp(1j * 0.475)]
+        angles = np.array([[-2, 0.5], [-0.4, -1.2], [-0.7, 2.0]])
+        new_D = np.exp(1j * np.array([[0, 3 / 4], [2.1, -0.3], [0.75, 0.2]]))
+
         for dec in (decomp, decomp2):
-            assert len(dec) == 4
-            qml.assert_equal(decomp[0], qml.QubitUnitary(np.array([np.eye(2) * p for p in ph]), 0))
-            qml.assert_equal(decomp[1], qml.RZ(angles[0], 1))
-            qml.assert_equal(decomp[2], qml.RZ(angles[1], 0))
-            qml.assert_equal(decomp[3], qml.IsingZZ(angles[2], [0, 1]))
+            assert len(dec) == 2
+            qml.assert_equal(decomp[0], qml.DiagonalQubitUnitary(new_D, wires=[0]))
+            qml.assert_equal(decomp[1], qml.SelectPauliRot(angles, [0], target_wire=1))
 
     def test_decomposition_three_qubits(self):
         """Test that a three-qubit DiagonalQubitUnitary is decomposed correctly."""
@@ -700,16 +697,12 @@ class TestDiagonalQubitUnitary:
         decomp = qml.DiagonalQubitUnitary.compute_decomposition(D, [0, 1, 2])
         decomp2 = qml.DiagonalQubitUnitary(D, wires=[0, 1, 2]).decomposition()
 
+        angles = np.array([-2, 0.5, -0.1, 1.7])
+        new_D = np.exp(1j * np.array([0, 3 / 4, 0.15, 1.45]))
         for dec in (decomp, decomp2):
-            assert len(dec) == 8
-            qml.assert_equal(decomp[0], qml.QubitUnitary(np.eye(2) * np.exp(0.5875j), 0))
-            qml.assert_equal(decomp[1], qml.RZ(0.025, 2))
-            qml.assert_equal(decomp[2], qml.RZ(1.025, 1))
-            qml.assert_equal(decomp[3], qml.IsingZZ(-1.075, [1, 2]))
-            qml.assert_equal(decomp[4], qml.RZ(0.425, 0))
-            qml.assert_equal(decomp[5], qml.IsingZZ(-0.775, [0, 2]))
-            qml.assert_equal(decomp[6], qml.IsingZZ(-0.275, [0, 1]))
-            qml.assert_equal(decomp[7], qml.MultiRZ(-0.175, [0, 1, 2]))
+            assert len(dec) == 2
+            qml.assert_equal(decomp[0], qml.DiagonalQubitUnitary(new_D, wires=[0, 1]))
+            qml.assert_equal(decomp[1], qml.SelectPauliRot(angles, [0, 1], target_wire=2))
 
     def test_decomposition_three_qubits_broadcasted(self):
         """Test that a broadcasted three-qubit DiagonalQubitUnitary is decomposed correctly."""
@@ -723,26 +716,12 @@ class TestDiagonalQubitUnitary:
         decomp = qml.DiagonalQubitUnitary.compute_decomposition(D, [0, 1, 2])
         decomp2 = qml.DiagonalQubitUnitary(D, wires=[0, 1, 2]).decomposition()
 
-        angles = [
-            [0.025, -0.55],
-            [1.025, -0.75],
-            [-1.075, -0.75],
-            [0.425, 0.3],
-            [-0.775, 0.4],
-            [-0.275, 0.5],
-            [-0.175, 0.1],
-        ]
-        ph = [np.exp(0.5875j), np.exp(0.625j)]
+        angles = np.array([[-2, 0.5, -0.1, 1.7], [-0.8, 0.5, -1.8, -0.1]])
+        new_D = np.exp(1j * np.array([[0, 3 / 4, 0.15, 1.45], [0.6, 0.35, 1.4, 0.15]]))
         for dec in (decomp, decomp2):
-            assert len(dec) == 8
-            qml.assert_equal(decomp[0], qml.QubitUnitary(np.array([np.eye(2) * p for p in ph]), 0))
-            qml.assert_equal(decomp[1], qml.RZ(angles[0], 2))
-            qml.assert_equal(decomp[2], qml.RZ(angles[1], 1))
-            qml.assert_equal(decomp[3], qml.IsingZZ(angles[2], [1, 2]))
-            qml.assert_equal(decomp[4], qml.RZ(angles[3], 0))
-            qml.assert_equal(decomp[5], qml.IsingZZ(angles[4], [0, 2]))
-            qml.assert_equal(decomp[6], qml.IsingZZ(angles[5], [0, 1]))
-            qml.assert_equal(decomp[7], qml.MultiRZ(angles[6], [0, 1, 2]))
+            assert len(dec) == 2
+            qml.assert_equal(decomp[0], qml.DiagonalQubitUnitary(new_D, wires=[0, 1]))
+            qml.assert_equal(decomp[1], qml.SelectPauliRot(angles, [0, 1], target_wire=2))
 
     @pytest.mark.parametrize("n", [1, 2, 3])
     def test_decomposition_matrix_match(self, n, seed):
@@ -775,7 +754,7 @@ class TestDiagonalQubitUnitary:
         assert qml.math.allclose(orig_mat, decomp_mat2)
 
     @pytest.mark.parametrize(
-        "dtype", [np.float64, np.float32, np.int64, np.int32, np.complex128, np.complex64]
+        "dtype", [np.float64, np.float32, np.int64, np.int32, np.int16, np.complex128, np.complex64]
     )
     def test_decomposition_cast_to_complex128(self, dtype):
         """Test that the parameters of decomposed operations are of the correct dtype."""
@@ -784,10 +763,18 @@ class TestDiagonalQubitUnitary:
         decomp1 = qml.DiagonalQubitUnitary(D, wires).decomposition()
         decomp2 = qml.DiagonalQubitUnitary.compute_decomposition(D, wires)
 
-        assert decomp1[0].data[0].dtype == np.complex128
-        assert decomp2[0].data[0].dtype == np.complex128
-        assert all(op.data[0].dtype == np.float64 for op in decomp1[1:])
-        assert all(op.data[0].dtype == np.float64 for op in decomp2[1:])
+        r_dtype = (
+            np.float64 if dtype in [np.float64, np.int64, np.int32, np.complex128] else np.float32
+        )
+        c_dtype = (
+            np.complex128
+            if dtype in [np.float64, np.int64, np.int32, np.complex128]
+            else np.complex64
+        )
+        assert decomp1[0].data[0].dtype == c_dtype
+        assert decomp2[0].data[0].dtype == c_dtype
+        assert decomp1[1].data[0].dtype == r_dtype
+        assert decomp2[1].data[0].dtype == r_dtype
 
     def test_controlled(self):
         """Test that the correct controlled operation is created when controlling a qml.DiagonalQubitUnitary."""