Repository: PennyLaneAI/pennylane
Pull request #7436 (2025-05-15 08:18:52)

Reduce some `templates/subroutines` tests
**Context:**
Templates quickly become large operations, which can cause template tests to become very expensive.

- a "simple" `assert_valid` check on `qml.kUpCCGSD` for 8 wires recently took >8 minutes in CI.
- there are other tests that do not really add test coverage but come with significant cost for the test suite.

**Description of the Change:**
Mildly reduce some test cases, with the constraint that no coverage/certainty may be lost.

**Benefits:**
Lighter test suite

**Possible Drawbacks:**
None identified

**Related GitHub Issues:**

[sc-91348]

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.

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 076400a4d..f3968b71f 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -73,6 +73,9 @@
 
 <h4>Resource-efficient Decompositions 🔎</h4>
 
+* The decomposition of `qml.PCPhase` is now significantly more efficient for more than 2 qubits.
+  [(#7166)](https://github.com/PennyLaneAI/pennylane/pull/7166)
+
 * New decomposition rules comprising rotation gates and global phases have been added to `QubitUnitary` 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.
@@ -315,6 +318,9 @@ Here's a list of deprecations made this release. For a more detailed breakdown o
 
 <h3>Internal changes ⚙️</h3>
 
+* Remove and reduce excessively expensive test cases in `tests/templates/test_subroutines/` that do not add value.
+  [(#7436)](https://github.com/PennyLaneAI/pennylane/pull/7436)
+
 * Stop using `pytest-timeout` in the PennyLane CI/CD pipeline.
   [(#7451)](https://github.com/PennyLaneAI/pennylane/pull/7451)
 
diff --git a/pennylane/math/decomp_int_to_powers_of_two.md b/pennylane/math/decomp_int_to_powers_of_two.md
new file mode 100644
index 000000000..d87cfbea4
--- /dev/null
+++ b/pennylane/math/decomp_int_to_powers_of_two.md
@@ -0,0 +1,100 @@
+## Decomposing an integer into powers of two
+
+Here we decompose an integer $d<=2^{n-1}$ into additions and subtractions of the
+smallest-possible number of powers of two.
+This is equivalent to computing the minimal bit operations required to produce
+a bit string by subtraction and addition of weight-one bit strings:
+
+
+$$
+d = \sum_{i=0}^{n-1} c_i 2^{n-1-i}, \quad c_i\in\{-1, 0, 1\}.
+$$
+
+The number of powers of two is given by the OEIS sequence
+[A007302](https://oeis.org/A007302) and can be computed as
+`(d ^ (3 * d)).bit_count()`. Note that we use zero-based indexing.
+
+### Algorithm
+
+The algorithm works in the following steps, exploiting the insight from
+[this post](https://stackoverflow.com/questions/57797157/minimum-number-of-powers-of-2-to-get-an-integer/57805889#57805889>) on stack overflow that two bits of $d$ can be set
+correctly at a time. That is, given two bits of the target $d$ and the
+corresponding two bits of some intermediate working bit string $s$ that we try
+to set to $d$, we have 16 scenarios. We will only need to think about those
+scenarios where the less significant bit already differs between $d$ and
+$s$. In each such scenario, a single addition or subtraction to $s$
+suffices to set the bits of $s$ to those of $d$:
+
+| Bits | d=00 | d=01 | d=10 | d=11 |
+|------|------|------|------|------|
+| s=00 | 0    | +1   | 0    | -1   |
+| s=01 | -1   | 0    | +1   | 0    |
+| s=10 | 0    | -1   | 0    | +1   |
+| s=11 | +1   | 0    | -1   | 0    |
+
+Now let's lay out the algorithm itself. We will use this operation in step 3 below.
+
+1. Set $s=0=0\dots 0_2$ to be the zero bit string of length $n$.
+
+2. Initialize $i=n-1$ and an empty list $R=[]$ that will hold 2-tuples, each
+   containing an exponent and a factor.
+
+3. If $d_i=s_i$, insert $(i, 0)$ at the beginning of $R$.
+   Else, i.e., if $d_i=1-s_i$, set $\sigma=1$ if $i<2$ and otherwise
+   look at $d_{i-1}$ and select the sign $\sigma$ from the table above that
+   will make $s_{i-1}s_i$ match $d_{i-1}d_i$.
+   Insert $(i, \sigma)$ at the beginning of $R$ and add $\sigma 2^{n-i-1}$
+   to $s$.
+
+4. If $i=0$, terminate. Otherwise, update $i=i-1$ and go to step 3.
+
+### Calculation example
+
+Take $n=5$ and $d=7=00111_2$. If we only allowed for addition of powers of two,
+we would simply aggregate the single non-zero bits $00001_2$, $00010_2$ and
+$00100_2$, requiring three powers of two to be added. However, including substraction
+allows us to decompose $d$ as $d=8-1=01000_2 - 00001_2$ instead, only
+requiring two powers of two.
+To compute this result, we go through the steps from above:
+
+1. Set $s=00000_2$.
+
+2. Initialize $i=n-1=4$ and $R=[]$.
+
+3. $d_4=1=1-s_4$, so look at $d_3=1$ as well. The sign at $(00, 11)$
+   in the table above is $-1$, so we insert $(4, -1)$
+   and compute $s=00000_2 - 00001_2=11111_2$.
+
+4. Update $i=4-1=3$.
+
+5. $d_3=1=s_3$, so we insert $(3, 0)$ and move on.
+
+6. Update $i=3-1=2$.
+
+7. $d_2=1=s_2$, so we insert $(2, 0)$ and move on.
+
+8. Update $i=2-1=1$.
+
+9. $d_1=0=1-s_1$, and $i=1<2$, so we insert $(1, +1)$
+   and compute $s=11111_2 + 01000_2=00111_2=d$.
+
+10. Update $i=1-1=0$.
+
+11. $d_0=0=s_0$, so we insert $(0, 0)$ and move on.
+
+12. we have $i=0$, so we terminate.
+
+Overall, this creates the list $R=[(0,0), (1,+1), (2, 0), (3, 0), (4, -1)]$,
+corresponding to the decomposition we were after,
+$d=2^{n-1-1} - 2^{n-4-1}=2^3-2^0=8-1$.
+
+In the code below, we will actually *skip* the first entry of the 2-tuples in
+$R$, because they always equal the position of the tuple within $R$
+(after completing the algorithm).
+
+As a second example, consider the number
+$k=121_{10}=01111001_2$, which can be (trivially) decomposed into a sum of
+five powers of two by reading off the bits:
+$k = 2^6 + 2^5 + 2^4 + 2^3 + 2^0 = 64 + 32 + 16 + 8 + 1$.
+The decomposition here, however, allows for minus signs and achieves the decomposition
+$k = 2^7 - 2^3 + 2^0 = 128 - 8 + 1$, which only requires three powers of two.
diff --git a/pennylane/math/decomposition.py b/pennylane/math/decomposition.py
index cb9e0e5d7..5fc165370 100644
--- a/pennylane/math/decomposition.py
+++ b/pennylane/math/decomposition.py
@@ -218,3 +218,51 @@ def su2su2_to_tensor_products(U):
     )
 
     return math.convert_like(A, U), math.convert_like(B, U)
+
+
+def decomp_int_to_powers_of_two(k: int, n: int) -> list[int]:
+    r"""Decompose an integer :math:`k<=2^{n-1}` into additions and subtractions of the
+    smallest-possible number of powers of two.
+
+    Args:
+        k (int): Integer to be decomposed
+        n (int): Number of bits to consider
+
+    Returns:
+        list[int]: A list with length ``n``, with entry :math:`c_i` at position :math:`i`.
+
+    This function is documented in ``pennylane/ops/qubit/pcphase_decomposition.md``.
+
+    As an example, consider the number
+    :math:`k=121_{10}=01111001_2`, which can be (trivially) decomposed into a sum of
+    five powers of two by reading off the bits:
+    :math:`k = 2^6 + 2^5 + 2^4 + 2^3 + 2^0 = 64 + 32 + 16 + 8 + 1`.
+    The decomposition here, however, allows for minus signs and achieves the decomposition
+    :math:`k = 2^7 - 2^3 + 2^0 = 128 - 8 + 1`, which only requires three powers of two.
+    """
+    R = []
+    assert k <= 2 ** (n - 1)
+    s = 0
+    powers = 2 ** np.arange(n)
+    for p in powers:  # p = 2**(n-1-i)
+        if s & p == k & p:
+            # Equal bit, move on
+            factor = 0
+        else:
+            # Differing bit, consider pairs of bits
+            if p >= 2 ** (n - 2):
+                # 2**(n-1-i) >= 2**(n-2) is the same condition as i < 2
+                factor = 1
+            else:
+                # Table entry from documentation
+                in_middle_rows = (s & (p + 2 * p)).bit_count() == 1  # two bits of s are 01 or 10
+                in_last_cols = bool(k & (2 * p))  # latter bit of k is 1
+                if in_middle_rows != in_last_cols:  # xor between in_middle_rows and in_last_cols
+                    factor = -1
+                else:
+                    factor = 1
+
+            s += factor * p
+        R.insert(0, factor)
+
+    return R
diff --git a/pennylane/operation.py b/pennylane/operation.py
index 1e54d2b37..5a602cb7b 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -1523,7 +1523,7 @@ class Operator(abc.ABC, metaclass=ABCCaptureMeta):
         we get the generator
 
         >>> U.generator()
-          0.5 * Y(0) + Z(0) @ X(1)
+        0.5 * Y(0) + Z(0) @ X(1)
 
         The generator may also be provided in the form of a dense or sparse Hamiltonian
         (using :class:`.LinearCombination` and :class:`.SparseHamiltonian` respectively).
diff --git a/pennylane/ops/qubit/parametric_ops_multi_qubit.py b/pennylane/ops/qubit/parametric_ops_multi_qubit.py
index 1a186b203..602eda8b5 100644
--- a/pennylane/ops/qubit/parametric_ops_multi_qubit.py
+++ b/pennylane/ops/qubit/parametric_ops_multi_qubit.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2023 Xanadu Quantum Technologies Inc.
+# Copyright 2018-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.
@@ -26,6 +26,7 @@ import numpy as np
 import pennylane as qml
 from pennylane.decomposition import add_decomps, register_resources
 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.typing import TensorLike
 from pennylane.wires import Wires, WiresLike
@@ -602,13 +603,78 @@ pauli_rot_decomposition = qml.register_resources(_pauli_rot_resources, _pauli_ro
 add_decomps(PauliRot, pauli_rot_decomposition)
 
 
+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.
 
     This gate applies a complex phase :math:`e^{i\phi}` to the first :math:`dim`
     basis vectors of the input state while applying a complex phase :math:`e^{-i \phi}`
-    to the remaining basis vectors. For example, consider the 2-qubit case where :math:`dim = 3`:
+    to the remaining basis vectors. For example, consider the 2-qubit case where ``dim = 3``:
 
     .. math:: \Pi(\phi) = \begin{bmatrix}
                 e^{i\phi} & 0 & 0 & 0 \\
@@ -617,6 +683,10 @@ class PCPhase(Operation):
                 0 & 0 & 0 & e^{-i\phi}
             \end{bmatrix}.
 
+    This can also be written as :math:`\Pi(\phi) = \exp(i\phi(2\Pi-\mathbb{I}_N))`, where
+    :math:`N=2^n` is the Hilbert space dimension for :math:`n` qubits and :math:`\Pi` is
+    the diagonal projector with ``dim`` ones and ``N-dim`` zeros.
+
     **Details:**
 
     * Number of wires: Any (the operation can act on any number of wires)
@@ -633,30 +703,46 @@ class PCPhase(Operation):
 
     We can define a circuit using :class:`~.PCPhase` as follows:
 
-    >>> dev = qml.device('default.qubit', wires=2)
-    >>> @qml.qnode(dev)
-    >>> def example_circuit():
-    ...     qml.PCPhase(0.27, dim = 2, wires=range(2))
-    ...     return qml.state()
+    >>> op_3 = qml.PCPhase(0.27, dim = 3, wires=range(3))
 
-    The resulting operation applies a complex phase :math:`e^{0.27i}` to the first :math:`dim = 2`
-    basis vectors and :math:`e^{-0.27i}` to the remaining basis vectors.
+    The resulting operation applies a complex phase :math:`e^{0.27i}` to the first :math:`dim = 3`
+    basis vectors and :math:`e^{-0.27i}` to the remaining basis vectors, as we can see from
+    the diagonal of the matrix for this circuit.
 
-    >>> print(np.round(qml.matrix(example_circuit)(),2))
-    [[0.96+0.27j 0.  +0.j   0.  +0.j   0.  +0.j  ]
-     [0.  +0.j   0.96+0.27j 0.  +0.j   0.  +0.j  ]
-     [0.  +0.j   0.  +0.j   0.96-0.27j 0.  +0.j  ]
-     [0.  +0.j   0.  +0.j   0.  +0.j   0.96-0.27j]]
+    >>> print(np.round(np.diag(qml.matrix(op_3)),2))
+    [0.96+0.27j 0.96+0.27j 0.96+0.27j 0.96-0.27j 0.96-0.27j 0.96-0.27j
+     0.96-0.27j 0.96-0.27j]
 
-    We can also choose a different :math:`dim` value to apply the phase shift to a different set of
+    We can also choose a different ``dim`` value to apply the phase shift to a different set of
     basis vectors as follows:
 
-    >>> pc_op = qml.PCPhase(1.23, dim=3, wires=[1, 2])
-    >>> print(np.round(qml.matrix(pc_op),2))
-    [[0.33+0.94j 0.  +0.j   0.  +0.j   0.  +0.j  ]
-     [0.  +0.j   0.33+0.94j 0.  +0.j   0.  +0.j  ]
-     [0.  +0.j   0.  +0.j   0.33+0.94j 0.  +0.j  ]
-     [0.  +0.j   0.  +0.j   0.  +0.j   0.33-0.94j]]
+    >>> op_7 = qml.PCPhase(1.23, dim=7, wires=[1, 2, 3])
+    >>> print(np.round(np.diag(qml.matrix(op_7)),2))
+    [0.33+0.94j 0.33+0.94j 0.33+0.94j 0.33+0.94j 0.33+0.94j 0.33+0.94j
+     0.33+0.94j 0.33-0.94j]
+
+    ``PCPhase`` operations are decomposed into (multi-)controlled :class:`~.PhaseShift`
+    operations which share the same control values on common control wires, and Pauli-X operations,
+    possibly complemented by a global phase.
+
+    >>> op_13 = qml.PCPhase(1.23, dim=13, wires=[1, 2, 3, 4])
+    >>> print(qml.draw(op_13.decomposition)())
+    1: ──GlobalPhase(-1.23)─╭●─────────╭●───────────┤
+    2: ──GlobalPhase(-1.23)─╰Rϕ(-2.46)─├●───────────┤
+    3: ──GlobalPhase(-1.23)────────────├○───────────┤
+    4: ──GlobalPhase(-1.23)──X─────────╰Rϕ(2.46)──X─┤
+
+    If ``dim`` is a power of two, a single (multi-controlled) ``PhaseShift`` gate is sufficient:
+
+    >>> op_16 = qml.PCPhase(1.23, dim=16, wires=range(6))
+    >>> print(qml.draw(op_16.decomposition, wire_order=range(6), show_all_wires=True)())
+    0: ──GlobalPhase(1.23)────╭○───────────┤
+    1: ──GlobalPhase(1.23)──X─╰Rϕ(2.46)──X─┤
+    2: ──GlobalPhase(1.23)─────────────────┤
+    3: ──GlobalPhase(1.23)─────────────────┤
+    4: ──GlobalPhase(1.23)─────────────────┤
+    5: ──GlobalPhase(1.23)─────────────────┤
+
     """
 
     num_params = 1
@@ -669,8 +755,25 @@ class PCPhase(Operation):
     parameter_frequencies = [(2,)]
 
     def generator(self) -> "qml.Hermitian":
-        dim, shape = self.hyperparameters["dimension"]
-        mat = np.diag([1 if index < dim else -1 for index in range(shape)])
+        r"""Generator of the ``PCPhase`` operator, which is in single-parameter-form.
+        The operator reads
+
+        .. math:: \Pi(\phi) = e^{i\phi (2\Pi - \mathbb{I}_N)},
+
+        where :math:`\Pi` is the projector onto the first :math`d` (``dim``) computational basis
+        states and :math:`N=2^n` is the Hilbert space dimension for :math:`n` qubits.
+
+        Correspondingly, the generator is
+        :math:`2\Pi - \mathbb{I}_N=\text{diag}(\underset{d\text{ times}}{\underbrace{1, \dots, 1}},\underset{(N-d)\text{ times}}{\underbrace{-1, \dots, -1}})`:
+
+        >>> qml.PCPhase(0.5, dim=3, wires=[0, 1]).generator()
+        Hermitian(array([[ 1,  0,  0,  0],
+           [ 0,  1,  0,  0],
+           [ 0,  0,  1,  0],
+           [ 0,  0,  0, -1]]), wires=[0, 1])
+        """
+        dim, N = self.hyperparameters["dimension"]
+        mat = np.diag([1] * dim + [-1] * (N - dim))
         return qml.Hermitian(mat, wires=self.wires)
 
     def _flatten(self) -> FlatPytree:
@@ -710,7 +813,7 @@ class PCPhase(Operation):
             return r
 
         arg = 1j * phi
-        prefactors = qml.math.array([1 if index < d else -1 for index in range(t)], like=phi)
+        prefactors = qml.math.array([1] * d + [-1] * (t - d), like=phi)
 
         if qml.math.ndim(arg) == 0:
             return qml.math.diag(qml.math.exp(arg * prefactors))
@@ -730,7 +833,7 @@ class PCPhase(Operation):
             return stack_last([phase if index < d else minus_phase for index in range(t)])
 
         arg = 1j * phi
-        prefactors = qml.math.array([1 if index < d else -1 for index in range(t)], like=phi)
+        prefactors = qml.math.array([1] * d + [-1] * (t - d), like=phi)
 
         if qml.math.ndim(phi) == 0:
             product = arg * prefactors
@@ -742,58 +845,151 @@ class PCPhase(Operation):
     def compute_decomposition(
         *params: TensorLike, wires: WiresLike, **hyperparams
     ) -> list["qml.operation.Operator"]:
-        r"""Representation of the operator as a product of other operators (static method).
-
-        .. math:: O = O_1 O_2 \dots O_n.
-
-        .. note::
-
-            Operations making up the decomposition should be queued within the
-            ``compute_decomposition`` method.
-
-        .. seealso:: :meth:`~.Operator.decomposition`.
+        r"""Representation of the PCPhase operator as a product of other operators (static method).
 
         Args:
-            *params (list): trainable parameters of the operator, as stored in the ``parameters`` attribute
+            *params (list): trainable parameters of the operator, as stored in the
+                ``parameters`` attribute
             wires (Iterable[Any], Wires): wires that the operator acts on
-            **hyperparams (dict): non-trainable hyper-parameters of the operator, as stored in the ``hyperparameters`` attribute
+            **hyperparams (dict): non-trainable hyper-parameters of the operator,
+                as stored in the ``hyperparameters`` attribute
 
         Returns:
             list[Operator]: decomposition of the operator
-        """
-        phi = params[0]
-        k, n = hyperparams["dimension"]
-
-        def _get_op_from_binary_rep(binary_rep, theta, wires):
-            if len(binary_rep) == 1:
-                op = (
-                    PhaseShift(theta, wires[0])
-                    if int(binary_rep)
-                    else PauliX(wires[0]) @ PhaseShift(theta, wires[0]) @ PauliX(wires[0])
-                )
-            else:
-                base_op = (
-                    PhaseShift(theta, wires[-1])
-                    if int(binary_rep[-1])
-                    else PauliX(wires[-1]) @ PhaseShift(theta, wires[-1]) @ PauliX(wires[-1])
-                )
-                op = qml.ctrl(
-                    base_op, control=wires[:-1], control_values=[int(i) for i in binary_rep[:-1]]
-                )
-            return op
 
-        n_log2 = int(np.log2(n))
-        positive_binary_reps = [bin(_k)[2:].zfill(n_log2) for _k in range(k)]
-        negative_binary_reps = [bin(_k)[2:].zfill(n_log2) for _k in range(k, n)]
+        In short, this decomposition relies on decomposing the generator (see :meth:`~.generator`)
+        of the ``PCPhase`` gate into generators of multicontrolled :class:`~.PhaseShift` gates,
+        potentially complemented with (non-controlled) Pauli-X gates and/or a global phase.
+        For example, for ``dim=13`` on four qubits:
+
+        >>> op_13 = qml.PCPhase(1.23, dim=13, wires=[1, 2, 3, 4])
+        >>> print(qml.draw(op_13.decomposition)())
+        1: ──GlobalPhase(-1.23)─╭●─────────╭●───────────┤
+        2: ──GlobalPhase(-1.23)─╰Rϕ(-2.46)─├●───────────┤
+        3: ──GlobalPhase(-1.23)────────────├○───────────┤
+        4: ──GlobalPhase(-1.23)──X─────────╰Rϕ(2.46)──X─┤
+
+        A detailed description of the algorithm and integer decomposition subroutine
+        is provided in ``pennylane/ops/qubit/pcphase_decomposition.md`` and
+        ``pennylane/math/decomp_int_to_powers_of_two.md``. In the following
+        we provide a detailed example for illustration purposes.
+
+        **Detailed example**
+
+        Consider the projector-controlled phase gate on :math:`n=4` qubits and with
+        :math:`d=\texttt{dim}=3`, i.e,
+
+        >>> op_3 = qml.PCPhase(1.23, dim=3, wires=[0, 1, 2, 3])
+
+        It acts on :math:`N=2^n=16`-dimensional vectors and is described by
+
+        .. math:: \Pi(\phi) = \exp(i\phi G) = \exp(i\phi(2\Pi-\mathbb{I}_N)),
+
+        where :math:`G` is a diagonal matrix with :math:`d=3` ones, followed by
+        :math:`2^n-d = 16 - 3=13` negative ones. Accordingly, :math:`\Pi` is diagonal with
+        :math:`3` ones and :math:`13` zeros.
+
+        First, we implement the global phase generated by :math:`\mathbb{I}_N` with
+        a :class:`~.GlobalPhase` gate with angle :math:`-\phi`.
+        Then we decompose :math:`d` into powers of two with positive or negative sign, via
+        :math:`d=3=4-1 = 2^2-2^0`. This decomposition tells us that we can write the
+        target gate with two (multi-)controlled phase shift gates. For this, we rewrite
+        the projector :math:`\Pi` according to the decomposition as
+
+        .. math::
+
+            \Pi &= \text{diag}(1, 1, 1, 0, 0, \dots, 0)\\
+            &=\text{diag}(1, 1, 1, 1, 0, \dots, 0)
+            -\text{diag}(0, 0, 0, 1, 0, \dots, 0)
+
+        where :math:`0,\dots, 0` indicates :math:`12` zeros each time.
+        How do we realize this projector decomposition on the gate level?
+
+        A singly-controlled phase shift gate applies a phase to a quarter of all computational
+        basis states (the control filters by the state of one qubit, and the phase shift gate
+        itself filters by the :math:`|1\rangle` state of the target qubit, cutting the number
+        of states we are acting on in half each time).
+        For :math:`n=4`, this amounts to :math:`2^4/4=4` states, which is exactly
+        what we need for the first term above. To apply the phase to the *first* four states,
+        :math:`|0000\rangle`, :math:`|0001\rangle`, :math:`|0010\rangle`, and :math:`|0011\rangle`,
+        we want to "filter by" the first two qubits being in the :math:`|0\rangle` state.
+        For qubit :math:`0`, we do this by controlling on the :math:`|0\rangle` state.
+        For qubit :math:`1`, we pick it as the target of the controlled phase shift operation.
+        Generically, this would make it act on the :math:`|1\rangle` state, so we simply flip
+        qubit :math:`1` before and after the operation to apply the phase to the :math:`|0\rangle`
+        state instead.
+        Thus, we conclude this first step by applying the gates
+        ``qml.X(1)``, ``qml.ctrl(qml.PhaseShift(2 * phi, 1), control=[0], control_values=[0])``,
+        and ``qml.X(1)``.
+
+        Next, we implement the second term in the projector decomposition, applying a phase
+        to a single computational basis state. This requires us to fully control a phase shift
+        gate, i.e., we use the last qubit as target and the other three as controls (there is
+        some freedom of choice here, but this is a convenient choice).
+        We want to apply the phase to the state :math:`|3\rangle=|0011\rangle`. So the controls
+        :math:`0` and :math:`1` are set to zero and the control :math:`2` is set to one.
+        As we want to effect the phase onto the :math:`|1\rangle` state of qubit :math:`3`,
+        we don't need to flip the target bit as we did before. However, given the negative sign
+        in the projector decomposition, we need to multiply the phase with :math:`-1`.
+        Overall, we apply the gate
+        ``qml.ctrl(qml.PhaseShift(-2 * phi, 3), control=[0, 1, 2], control_values=[0, 0, 1])``,
+        which concludes the decomposition, now reading:
+
+        >>> print(qml.draw(op_3.decomposition)())
+        0: ──GlobalPhase(1.23)────╭○───────────╭○─────────┤
+        1: ──GlobalPhase(1.23)──X─╰Rϕ(2.46)──X─├○─────────┤
+        2: ──GlobalPhase(1.23)─────────────────├●─────────┤
+        3: ──GlobalPhase(1.23)─────────────────╰Rϕ(-2.46)─┤
 
-        positive_ops = [
-            _get_op_from_binary_rep(br, phi, wires=wires) for br in positive_binary_reps
-        ]
-        negative_ops = [
-            _get_op_from_binary_rep(br, -1 * phi, wires=wires) for br in negative_binary_reps
-        ]
-
-        return positive_ops + negative_ops
+        """
+        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)
+
+        return decomp
 
     def adjoint(self) -> "PCPhase":
         """Computes the adjoint of the operator."""
diff --git a/pennylane/ops/qubit/pcphase_decomposition.md b/pennylane/ops/qubit/pcphase_decomposition.md
new file mode 100644
index 000000000..3d05322cb
--- /dev/null
+++ b/pennylane/ops/qubit/pcphase_decomposition.md
@@ -0,0 +1,115 @@
+This document explains the algorithm behind the decomposition of `qml.PCPhase` gates.
+
+# Algorithm description
+
+For this decomposition, we will decompose the generator of the `PCPhase` gate,
+which we denote as $G$, into multiple projectors, which generate (multi-controlled)
+phase shift gates. The generator of `PCPhase` is given as 
+
+
+$$
+G = 2 \left(\sum_{j=0}^{d-1} |j\rangle\langle j|\right) - \mathbb{I}_N
+= \text{diag}(
+\underset{d\text{ times}}{\underbrace{1, \dots, 1}},
+\underset{(N-d)\text{ times}}{\underbrace{-1, \dots, -1}}
+),
+$$
+
+where we denoted the overall dimension as $N=2^n$ for $n$ qubits and the
+subspace dimension as $d$, which is called `dim` in code.
+Our first step is to move to a slightly different convention in order to arrive at a
+*projector* that has ones and zeros, instead of ones and
+minus ones (note how $G$ is not a projector because $G^2=\mathbb{I}_N\neq G$).
+To do this, we define
+
+$$
+\Pi = \begin{cases}
+    \sum_{j=0}^{d-1} |j\rangle\langle j| & \text{for }\ d\leq \tfrac{N}{2}\\
+    \sum_{j=d}^{N-1} |j\rangle\langle j| & \text{for }\ d> \tfrac{N}{2}
+\end{cases},
+$$
+
+so that $G=2\Pi -\mathbb{I}_N$ for $d\leq \tfrac{N}{2}$ and
+$G=\mathbb{I}_N - 2 \Pi$ for $d>\tfrac{N}{2}$. We may summarize
+this as $G=\sigma (2\Pi -\mathbb{I}_N)$ with the sign $\sigma=\pm 1$
+defined accordingly for the two cases.
+
+Now we want to decompose $\Pi$, a sum of the computational basis state projectors
+$|j\rangle\langle j|$ onto the $d$ consecutive smallest (largest, for
+$\sigma=-1$) basis states,
+into projectors that contain a power of two basis state projectors. We will see below that
+these are exactly the generators of (multi-controlled) phase shift gates, allowing us
+to implement the decomposition. As we may flip the sign of these individual phase shift
+gates, we may also use subtraction when decomposing $\Pi$.
+
+This decomposition of $\Pi$ is based on decomposing the integer $d$ into sums
+and differences of powers of two, which is detailed in 
+`pennylane/math/decomp_int_to_powers_of_two.md`. Such a decomposition reads
+
+$$
+d = \sum_{i=0}^{n-1} c_i 2^{n-1-i}, \quad c_i\in\{-1, 0, 1\}.
+$$
+
+Now we need to use this decomposition to combine projectors of $2^{n-1-i}$
+computational basis state projectors into the desired $\Pi$.
+We start with the scenario $\sigma=+1$, i.e., $d\leq \tfrac{N}{2}$, in which
+case this decomposition has the form
+
+$$
+\Pi = \sum_{i=0}^{n-1} c_i \sum_{j=s_i}^{e_i} |j\rangle\langle j|,
+$$
+
+where $s_i$ and $e_i$ are the start and end positions
+for the $i$\ th contribution to $\Pi$.
+Besides computing the integer decomposition itself, this decomposition function is mostly
+about figuring out the right control structure to realize the correct positioning
+according to $s_i$ and $e_i$.
+
+To begin, initialize empty lists of control wires and control values, corresponding
+to $s_0=0$. Then, starting at $i=0$, which also marks the current wire,
+go through the integer decomposition $\{c_i\}$ from above and perform the
+following steps:
+
+1. If $c_i=0$, go to step 2 immediately. Else, we want to
+   add/subtract computational basis state projectors to the decomposition. The number of
+   projectors $2^{n-1-i}$ is stored in the number of control wires that we already
+   aggregated. In addition, $s_i$ (and thus $e_i$) is encoded in the control
+   values, up to whether the projectors should be located in the $|0\rangle$
+   or $|1\rangle$ subspace of the current target wire $i$. If we want to
+   add projectors ($c_i=+1$), we append at the beginning of the current subspace,
+   so in the $|0\rangle$ subspace. If we want to remove projectors ($c_i=-1$),
+   we need to remove them from the end of the subspace, so from $|1\rangle$.
+   Having figured out the subspace, we add the corresponding projector to the
+   decomposition. In code, we immediately apply the corresponding phase gate, see below.
+
+2. Add the current wire $i$ to the control wires. The control value to be added
+   depends on the current $c_i$ and on the next non-zero coefficient, $d_i$.
+   For $d_i=+1$, we are going to add more projectors, so we should control into
+   the $|1\rangle$ subspace (add onto just added projectors for $c_i=+1$,
+   or re-add just removed projectors for $c_i=-1$) unless we have $c_i=0$,
+   in which case we
+   did not just modify the decomposition and need to control into the $|0\rangle$
+   subspace instead. In either case, this sets $s_{i+1}$ to $e_i+1$.
+   For $d_i=-1$, we will subtract projectors. If we just subtracted projectors (from
+   the $|1\rangle$ subspace, $c_i=-1$), we need to remove further from the
+   $|0\rangle$ subspace and control into it. If we just added projectors (to the
+   $|0\rangle$ subspace), we remove some of those again, and need to control into
+   the $|0\rangle$ subspace as well. For $c_i=0$, we need to remove from
+   the upper end, i.e., the $|1\rangle$ subspace.
+   All of this is easily summarized: If the next non-zero coefficient, $d_i$,
+   is positive (negative), set the next control value to $1$ ($0$). If
+   the current coefficient is $c_i=0$, flip the control value.
+
+3. Increment $i$ by one. If $i=n$, terminate, else go to step 1.
+
+If $\sigma=-1$ instead, i.e., $d>\tfrac{N}{2}$, the start and end points
+$s_i$ and $e_i$ need to be mapped via $x\mapsto N-1-x$. This
+implies a change of the subspace acted on in step 1 above, and a change of the control
+value in step 2 above, both via Boolean negation, i.e., $x\mapsto 1-x$.
+
+To conclude, note that the control structure above, together with the information of
+the subspace we want to add projectors to/subtract projectors from, allows us to
+collect (multi-controlled) phase shift gates as the gate decomposition.
+The sign of all phases is multiplied by $2\sigma$ to realize the first term of
+$G=\sigma(2\Pi -\mathbb{I}_N)$, and a global phase with angle multiplied by
+$\sigma$ is added to realize the second term.
diff --git a/pennylane/templates/subroutines/fable.py b/pennylane/templates/subroutines/fable.py
index 490b1d1c9..7846cecf8 100644
--- a/pennylane/templates/subroutines/fable.py
+++ b/pennylane/templates/subroutines/fable.py
@@ -137,7 +137,6 @@ class FABLE(Operation):
         Returns:
             list[.Operator]: list of gates for efficient circuit
         """
-        op_list = []
         alphas = qml.math.arccos(input_matrix).flatten()
         thetas = compute_theta(alphas)
 
@@ -150,8 +149,7 @@ class FABLE(Operation):
 
         wire_map = dict(enumerate(wires_j + wires_i))
 
-        for w in wires_i:
-            op_list.append(qml.Hadamard(w))
+        op_list = [qml.Hadamard(w) for w in wires_i]
 
         nots = {}
         for theta, control_index in zip(thetas, control_wires):
diff --git a/pennylane/templates/subroutines/kupccgsd.py b/pennylane/templates/subroutines/kupccgsd.py
index 0095ac0a9..63fb2a425 100644
--- a/pennylane/templates/subroutines/kupccgsd.py
+++ b/pennylane/templates/subroutines/kupccgsd.py
@@ -16,6 +16,7 @@ Contains the k-UpCCGSD template.
 """
 # pylint: disable-msg=too-many-branches,too-many-arguments,protected-access
 import copy
+from itertools import product
 
 import numpy as np
 
@@ -31,17 +32,15 @@ def generalized_singles(wires, delta_sz):
         \hat{T_1} = \sum_{pq} t_{p}^{q} \hat{c}^{\dagger}_{q} \hat{c}_{p}
 
     """
-    sz = np.array(
-        [0.5 if (i % 2 == 0) else -0.5 for i in range(len(wires))]
-    )  # alpha-beta electrons
+    n = len(wires)
+    sz = 0.5 * (-1) ** np.arange(n)  # alpha-beta electrons
     gen_singles_wires = []
-    for r in range(len(wires)):
-        for p in range(len(wires)):
-            if sz[p] - sz[r] == delta_sz and p != r:
-                if r < p:
-                    gen_singles_wires.append(wires[r : p + 1])
-                else:
-                    gen_singles_wires.append(wires[p : r + 1][::-1])
+    for r, p in product(range(n), repeat=2):
+        if sz[p] - sz[r] == delta_sz and p != r:
+            if r < p:
+                gen_singles_wires.append(wires[r : p + 1])
+            else:
+                gen_singles_wires.append(wires[p : r + 1][::-1])
     return gen_singles_wires
 
 
@@ -54,12 +53,9 @@ def generalized_pair_doubles(wires):
 
     """
     pair_gen_doubles_wires = [
-        [
-            wires[r : r + 2],
-            wires[p : p + 2],
-        ]  # wires for [wires[r], wires[r+1], wires[p], wires[p+1]] terms
-        for r in range(0, len(wires) - 1, 2)
-        for p in range(0, len(wires) - 1, 2)
+        # wires for [wires[r], wires[r+1], wires[p], wires[p+1]] terms
+        [wires[r : r + 2], wires[p : p + 2]]
+        for r, p in product(range(0, len(wires) - 1, 2), repeat=2)
         if p != r  # remove redundant terms
     ]
     return pair_gen_doubles_wires
@@ -125,14 +121,14 @@ class kUpCCGSD(Operation):
             # Build the electronic Hamiltonian
             symbols = ["H", "H"]
             coordinates = np.array([0.0, 0.0, -0.6614, 0.0, 0.0, 0.6614])
-            H, qubits = qml.qchem.molecular_hamiltonian(symbols, coordinates)
+            H, wires = qml.qchem.molecular_hamiltonian(symbols, coordinates)
 
             # Define the Hartree-Fock state
             electrons = 2
-            ref_state = qml.qchem.hf_state(electrons, qubits)
+            ref_state = qml.qchem.hf_state(electrons, wires)
 
             # Define the device
-            dev = qml.device('default.qubit', wires=qubits)
+            dev = qml.device('default.qubit', wires=wires)
 
             # Define the ansatz
             @qml.qnode(dev)
@@ -144,7 +140,7 @@ class kUpCCGSD(Operation):
             # Get the shape of the weights for this template
             layers = 1
             shape = qml.kUpCCGSD.shape(k=layers,
-                                n_wires=qubits, delta_sz=0)
+                                n_wires=wires, delta_sz=0)
 
             # Initialize the weight tensors
             np.random.seed(24)
@@ -230,12 +226,9 @@ class kUpCCGSD(Operation):
         d_wires = generalized_pair_doubles(list(wires))
 
         shape = qml.math.shape(weights)
-        if shape != (
-            k,
-            len(s_wires) + len(d_wires),
-        ):
+        if shape != (k, len(s_wires) + len(d_wires)):
             raise ValueError(
-                f"Weights tensor must be of shape {(k, len(s_wires) + len(d_wires),)}; got {shape}."
+                f"Weights tensor must be of shape {(k, len(s_wires) + len(d_wires))}; got {shape}."
             )
 
         init_state = qml.math.toarray(init_state)
@@ -319,7 +312,7 @@ class kUpCCGSD(Operation):
         r"""Returns the shape of the weight tensor required for this template.
         Args:
             k (int): Number of layers
-            n_wires (int): Number of qubits
+            n_wires (int): Number of wires
             delta_sz (int): Specifies the selection rules ``sz[p] - sz[r] = delta_sz``
             for the spin-projection ``sz`` of the orbitals involved in the single excitations.
             ``delta_sz`` can take the values :math:`0` and :math:`\pm 1`.
@@ -329,12 +322,12 @@ class kUpCCGSD(Operation):
 
         if n_wires < 4:
             raise ValueError(
-                f"This template requires the number of qubits to be greater than four; got 'n_wires' = {n_wires}"
+                f"This template requires the number of wires to be greater than four; got 'n_wires' = {n_wires}"
             )
 
         if n_wires % 2:
             raise ValueError(
-                f"This template requires an even number of qubits; got 'n_wires' = {n_wires}"
+                f"This template requires an even number of wires; got 'n_wires' = {n_wires}"
             )
 
         s_wires = generalized_singles(range(n_wires), delta_sz)


Test patch:
diff --git a/tests/math/test_decomposition_math.py b/tests/math/test_decomposition_math.py
new file mode 100644
index 000000000..59d0faf28
--- /dev/null
+++ b/tests/math/test_decomposition_math.py
@@ -0,0 +1,53 @@
+# 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.
+"""
+This submodule tests mathematical decomposition functions, for example of matrices or integers.
+"""
+import pytest
+
+from pennylane.math.decomposition import decomp_int_to_powers_of_two
+
+
+@pytest.mark.parametrize(
+    "k, n, exp_R",
+    [
+        (0, 1, [0]),
+        (0, 2, [0, 0]),
+        (1, 2, [0, 1]),
+        (0, 3, [0, 0, 0]),
+        (1, 3, [0, 0, 1]),
+        (2, 3, [0, 1, 0]),
+        (3, 3, [1, 0, -1]),
+        (4, 3, [1, 0, 0]),
+        (0, 4, [0, 0, 0, 0]),
+        (1, 4, [0, 0, 0, 1]),
+        (2, 4, [0, 0, 1, 0]),
+        (3, 4, [0, 1, 0, -1]),
+        (4, 4, [0, 1, 0, 0]),
+        (5, 4, [0, 1, 0, 1]),
+        (6, 4, [1, 0, -1, 0]),
+        (7, 4, [1, 0, 0, -1]),
+        (8, 4, [1, 0, 0, 0]),
+        (121, 8, [1, 0, 0, 0, -1, 0, 0, 1]),  # 121=01111001_2
+        (245, 10, [0, 1, 0, 0, 0, -1, 0, 1, 0, 1]),  # 245=0011110101_2
+        (8716, 15, [0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, -1, 0, 0]),  # 8716=010001000001100_2
+    ],
+)
+def test_decomp_int_to_powers_of_two(k, n, exp_R):
+    """Tests for ``decomp_int_to_powers_of_two``, which is used
+    to decompose, e.g., ``PCPhase`` operations."""
+
+    R = decomp_int_to_powers_of_two(k, n)
+    assert R == exp_R, f"\n{R}\n{exp_R}"
+    assert sum(val != 0 for val in R) == (k ^ (3 * k)).bit_count()
diff --git a/tests/ops/qubit/test_parametric_ops.py b/tests/ops/qubit/test_parametric_ops.py
index 7d6f0d59e..df7fa59d3 100644
--- a/tests/ops/qubit/test_parametric_ops.py
+++ b/tests/ops/qubit/test_parametric_ops.py
@@ -220,13 +220,8 @@ class TestOperations:
             return qml.state()
 
         def _get_expected_state(phase, dim, size):
-            return (
-                np.array(
-                    [np.exp(1j * phase) if i < dim else np.exp(-1j * phase) for i in range(size)]
-                )
-                * 1
-                / 2
-            )
+            phase = np.exp(1j * phase)
+            return np.concatenate([phase * np.ones(dim), phase.conj() * np.ones(size - dim)]) / 2
 
         assert np.allclose(circuit(theta, d), _get_expected_state(theta, d, 4))
 
@@ -654,26 +649,30 @@ class TestDecompositions:
 
         assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
 
-    @pytest.mark.parametrize(
-        "op", (qml.PCPhase(1.23, dim=1, wires=[0]), qml.PCPhase(1.23, dim=5, wires=[0, 1, 2]))
-    )
-    def test_pc_phase_decomposition(self, op):
+    two_wire_pcphases = [(0, [0, 1]), (1, [1, 0]), (2, ["a", 2]), (3, [1, 3]), (4, [9, 0])]
+    five_wire_pcphases = [(i, [0, 1, 3, 2, 7]) for i in range(2**5)]
+    other_pcphases = [(1, [0]), (2, [1]), (17, ["a", 2, "c", 4, 3, 0]), (3, list(range(5)))]
+
+    @pytest.mark.parametrize("dim, wires", two_wire_pcphases + five_wire_pcphases + other_pcphases)
+    def test_pcphase_decomposition(self, dim, wires):
         """Test that the PCPhase decomposition produces the same unitary"""
+        op = qml.PCPhase(np.random.random(), dim, wires=wires)
         decomp_ops = op.decomposition()
         decomp_op = qml.prod(*decomp_ops) if len(decomp_ops) > 1 else decomp_ops[0]
 
-        expected_mat = qml.matrix(op)
-        decomp_mat = qml.matrix(decomp_op)
+        expected_mat = qml.matrix(op, wire_order=wires)
+        decomp_mat = qml.matrix(decomp_op, wire_order=wires)
         assert np.allclose(expected_mat, decomp_mat)
 
-    def test_pc_phase_decomposition_broadcasted(self):
+    @pytest.mark.parametrize("dim, wires", two_wire_pcphases + other_pcphases)
+    def test_pcphase_decomposition_broadcasted(self, dim, wires):
         """Test that the broadcasted PCPhase decomposition produces the same unitary"""
-        op = qml.PCPhase([1.23, 4.56, 7.89], dim=5, wires=[0, 1, 2])
+        op = qml.PCPhase(np.random.random(3), dim, wires=wires)
         decomp_ops = op.decomposition()
         decomp_op = qml.prod(*decomp_ops) if len(decomp_ops) >= 1 else decomp_ops[0]
 
-        expected_mats = qml.matrix(op)
-        decomp_mats = qml.matrix(decomp_op)
+        expected_mats = qml.matrix(op, wire_order=wires)
+        decomp_mats = qml.matrix(decomp_op, wire_order=wires)
 
         for expected_mat, decomp_mat in zip(expected_mats, decomp_mats):
             assert np.allclose(expected_mat, decomp_mat)
diff --git a/tests/templates/test_subroutines/test_fable.py b/tests/templates/test_subroutines/test_fable.py
index d2ba5f249..6dc85a6a8 100644
--- a/tests/templates/test_subroutines/test_fable.py
+++ b/tests/templates/test_subroutines/test_fable.py
@@ -44,10 +44,7 @@ class TestFable:
     @pytest.mark.parametrize(
         ("input", "wires"),
         [
-            (
-                np.random.random((2, 2)),
-                3,
-            ),
+            (np.random.random((2, 2)), 3),
             (
                 np.array(
                     [
@@ -59,17 +56,10 @@ class TestFable:
                 ),
                 5,
             ),
-            (
-                np.random.random((8, 8)),
-                7,
-            ),
-            (
-                np.random.random((16, 16)),
-                9,
-            ),
+            (np.random.random((8, 8)), 7),
         ],
     )
-    def test_fable_real_for_variety_of_input_matrices(self, input, wires):
+    def test_fable_real_input_matrices(self, input, wires):
         """Test that FABLE produces the right circuit given a square, real-valued matrix"""
         dev = qml.device("default.qubit")
 
@@ -78,10 +68,8 @@ class TestFable:
             qml.FABLE(input_matrix=input, wires=range(wires), tol=0)
             return qml.state()
 
-        expected = (
-            len(input)
-            * qml.matrix(circuit, wire_order=range(wires))().real[0 : len(input), 0 : len(input)]
-        )
+        dim = len(input)
+        expected = dim * qml.matrix(circuit, wire_order=range(wires))().real[:dim, :dim]
         assert np.allclose(input, expected)
 
     def test_fable_imaginary_error(self, input_matrix):
@@ -391,89 +379,36 @@ class TestFable:
     @pytest.mark.parametrize(
         ("input", "wires"),
         [
-            (
-                np.random.random((1, 2)),
-                3,
-            ),
-            (
-                np.random.random((1, 1)),
-                3,
-            ),
-            (
-                np.random.random((2, 1)),
-                3,
-            ),
-            (
-                np.random.random((3, 2)),
-                5,
-            ),
-            (
-                np.random.random((4, 2)),
-                5,
-            ),
-            (
-                np.random.random((2, 3)),
-                5,
-            ),
-            (
-                np.random.random((2, 4)),
-                5,
-            ),
-            (
-                np.random.random((3, 4)),
-                5,
-            ),
-            (
-                np.random.random((4, 3)),
-                5,
-            ),
-            (
-                np.random.random((3, 5)),
-                7,
-            ),
-            (
-                np.random.random((3, 6)),
-                7,
-            ),
-            (
-                np.random.random((3, 7)),
-                7,
-            ),
-            (
-                np.random.random((4, 5)),
-                7,
-            ),
-            (
-                np.random.random((5, 5)),
-                7,
-            ),
-            (
-                np.random.random((6, 5)),
-                7,
-            ),
-            (
-                np.random.random((2, 9)),
-                9,
-            ),
-            (
-                np.random.random((9, 3)),
-                9,
-            ),
+            (np.random.random((1, 2)), 3),
+            (np.random.random((1, 1)), 3),
+            (np.random.random((2, 1)), 3),
+            (np.random.random((3, 2)), 5),
+            (np.random.random((4, 2)), 5),
+            (np.random.random((2, 3)), 5),
+            (np.random.random((2, 4)), 5),
+            (np.random.random((3, 4)), 5),
+            (np.random.random((4, 3)), 5),
+            (np.random.random((3, 5)), 7),
+            (np.random.random((3, 6)), 7),
+            (np.random.random((3, 7)), 7),
+            (np.random.random((4, 5)), 7),
+            (np.random.random((5, 5)), 7),
+            (np.random.random((6, 5)), 7),
         ],
     )
     def test_variety_of_matrix_shapes(self, input, wires):
-        """Test that FABLE runs without error for variety of cases."""
+        """Test that FABLE runs without error for a variety of input shapes."""
         dev = qml.device("default.qubit")
-        s = int(qml.math.ceil(qml.math.log2(max(len(input), len(input[0])))))
-        dimension = s**2
+        s = int(qml.math.ceil(qml.math.log2(max(input.shape))))
+        s = max(s, 1)
+        dim = 2**s
 
         @qml.qnode(dev)
         def circuit():
             qml.FABLE(input_matrix=input, wires=range(wires), tol=0)
             return qml.state()
 
-        expected = (
-            dimension
-            * qml.matrix(circuit, wire_order=range(wires))().real[0:dimension, 0:dimension]
-        )
-        assert qml.math.shape(expected) == (dimension, dimension)
+        circuit_mat = dim * qml.matrix(circuit, wire_order=range(wires))().real[:dim, :dim]
+        # Test that the matrix was encoded up to a constant
+        submat = circuit_mat[: input.shape[0], : input.shape[1]]
+        assert np.allclose(input, submat)
diff --git a/tests/templates/test_subroutines/test_kupccgsd.py b/tests/templates/test_subroutines/test_kupccgsd.py
index 98599f443..a83e3975f 100644
--- a/tests/templates/test_subroutines/test_kupccgsd.py
+++ b/tests/templates/test_subroutines/test_kupccgsd.py
@@ -54,7 +54,8 @@ def test_standard_validity(k, delta_sz, init_state, wires):
 
     op = qml.kUpCCGSD(weights, wires=wires, k=k, delta_sz=delta_sz, init_state=init_state)
 
-    qml.ops.functions.assert_valid(op)
+    skip_diff = len(wires) > 4
+    qml.ops.functions.assert_valid(op, skip_differentiation=skip_diff)
 
 
 class TestDecomposition:
@@ -160,7 +161,7 @@ class TestDecomposition:
         assert np.allclose(state1, state2, atol=tol, rtol=0)
 
     @pytest.mark.parametrize(
-        ("num_qubits", "k", "exp_state"),
+        ("num_wires", "k", "exp_state"),
         [
             (
                 4,
@@ -243,17 +244,17 @@ class TestDecomposition:
             ),
         ],
     )
-    def test_k_layers_upccgsd(self, num_qubits, k, exp_state, tol):
+    def test_k_layers_upccgsd(self, num_wires, k, exp_state, tol):
         """Test that the k-UpCCGSD template with multiple layers works correctly asserting the prepared state."""
 
-        wires = range(num_qubits)
+        wires = range(num_wires)
 
-        shape = qml.kUpCCGSD.shape(k=k, n_wires=num_qubits, delta_sz=0)
+        shape = qml.kUpCCGSD.shape(k=k, n_wires=num_wires, delta_sz=0)
         weight = np.pi / 2 * qml.math.ones(shape)
 
         dev = qml.device("default.qubit", wires=wires)
 
-        init_state = qml.math.array([1 if x < num_qubits // 2 else 0 for x in wires])
+        init_state = qml.math.array([1 if x < num_wires // 2 else 0 for x in wires])
 
         @qml.qnode(dev)
         def circuit(weight):
@@ -465,18 +466,18 @@ class TestAttributes:
         shape = qml.kUpCCGSD.shape(k, n_wires, delta_sz)
         assert shape == expected_shape
 
-    def test_shape_exception_not_enough_qubits(self):
-        """Test that the shape function warns if there are not enough qubits."""
+    def test_shape_exception_not_enough_wires(self):
+        """Test that the shape function warns if there are not enough wires."""
 
         with pytest.raises(
-            ValueError, match="This template requires the number of qubits to be greater than four"
+            ValueError, match="This template requires the number of wires to be greater than four"
         ):
             qml.kUpCCGSD.shape(k=2, n_wires=1, delta_sz=0)
 
-    def test_shape_exception_not_even_qubits(self):
-        """Test that the shape function warns if the number of qubits are not even."""
+    def test_shape_exception_not_even_wires(self):
+        """Test that the shape function warns if the number of wires are not even."""
 
-        with pytest.raises(ValueError, match="This template requires an even number of qubits"):
+        with pytest.raises(ValueError, match="This template requires an even number of wires"):
             qml.kUpCCGSD.shape(k=2, n_wires=5, delta_sz=0)
 
 
diff --git a/tests/templates/test_subroutines/test_trotter.py b/tests/templates/test_subroutines/test_trotter.py
index 4787c974a..82143dfad 100644
--- a/tests/templates/test_subroutines/test_trotter.py
+++ b/tests/templates/test_subroutines/test_trotter.py
@@ -1726,8 +1726,6 @@ class TestTrotterizedQfuncIntegration:
                 qml.ControlledPhaseShift(time * arg2, wires=[wires[1], wires[0]]),
                 qml.CNOT([wires[0], wires[1]]),
                 qml.QFT(wires=wires[1:-1]),
-                qml.QFT(wires=wires[1:-1]),
-                qml.QFT(wires=wires[1:-1]),
             ]
 
             if reverse:
@@ -1742,8 +1740,6 @@ class TestTrotterizedQfuncIntegration:
                 qml.ControlledPhaseShift((time / 2) * arg2, wires=[wires[1], wires[0]]),
                 qml.CNOT([wires[0], wires[1]]),
                 qml.QFT(wires=wires[1:-1]),
-                qml.QFT(wires=wires[1:-1]),
-                qml.QFT(wires=wires[1:-1]),
             ]
 
             if reverse:
@@ -1760,8 +1756,6 @@ class TestTrotterizedQfuncIntegration:
                 qml.ControlledPhaseShift((p_4 * time / 2) * arg2, wires=[wires[1], wires[0]]),
                 qml.CNOT([wires[0], wires[1]]),
                 qml.QFT(wires=wires[1:-1]),
-                qml.QFT(wires=wires[1:-1]),
-                qml.QFT(wires=wires[1:-1]),
             ]
 
             expected_decomp1 = (
@@ -1778,8 +1772,6 @@ class TestTrotterizedQfuncIntegration:
                 qml.ControlledPhaseShift((p4_comp * time / 2) * arg2, wires=[wires[1], wires[0]]),
                 qml.CNOT([wires[0], wires[1]]),
                 qml.QFT(wires=wires[1:-1]),
-                qml.QFT(wires=wires[1:-1]),
-                qml.QFT(wires=wires[1:-1]),
             ]
 
             expected_decomp2 = (
@@ -1803,7 +1795,7 @@ class TestTrotterizedQfuncIntegration:
         arg1 = 2.34
         arg2 = -6.78
         args = (arg1, arg2)
-        kwargs = {"kwarg1": True, "kwarg2": 3}
+        kwargs = {"kwarg1": True, "kwarg2": 1}
 
         def my_qfunc(time, arg1, arg2, wires, kwarg1=False, kwarg2=None):
             """Arbitrarily complex qfunc"""
@@ -1867,7 +1859,7 @@ class TestTrotterizedQfuncIntegration:
         arg1 = jnp.array(2.34)
         arg2 = jnp.array(-6.78)
         args = (arg1, arg2)
-        kwargs = {"kwarg1": True, "kwarg2": 3}
+        kwargs = {"kwarg1": True, "kwarg2": 1}
 
         def my_qfunc(time, arg1, arg2, wires, kwarg1=False, kwarg2=None):
             """Arbitrarily complex qfunc"""
@@ -1935,7 +1927,7 @@ class TestTrotterizedQfuncIntegration:
         arg1 = jnp.array(2.34)
         arg2 = jnp.array(-6.78)
         args = (arg1, arg2)
-        kwargs = {"kwarg1": True, "kwarg2": 3}
+        kwargs = {"kwarg1": True, "kwarg2": 1}
 
         def my_qfunc(time, arg1, arg2, wires, kwarg1=False, kwarg2=None):
             """Arbitrarily complex qfunc"""
@@ -2000,7 +1992,7 @@ class TestTrotterizedQfuncIntegration:
         wires = ["aux1", "aux2", 0, 1, "target"]
         arg1 = qnp.array(2.34)
         arg2 = qnp.array(-6.78)
-        kwargs = {"kwarg1": True, "kwarg2": 3}
+        kwargs = {"kwarg1": True, "kwarg2": 1}
 
         def my_qfunc(time, arg1, arg2, wires, kwarg1=False, kwarg2=None):
             """Arbitrarily complex qfunc"""
@@ -2070,7 +2062,7 @@ class TestTrotterizedQfuncIntegration:
         arg1 = jnp.array(2.34)
         arg2 = jnp.array(-6.78)
         args = (arg1, arg2)
-        kwargs = {"kwarg1": True, "kwarg2": 3}
+        kwargs = {"kwarg1": True, "kwarg2": 1}
 
         def my_qfunc(time, arg1, arg2, wires, kwarg1=False, kwarg2=None):
             """Arbitrarily complex qfunc"""