Repository: PennyLaneAI/pennylane
Pull request #7166 (2025-03-31 09:14:34)

Add new decomposition to `PCPhase`
**Context:**
Currently, PennyLane decomposes `PCPhase` in a doubly-unfavourable manner, by addressing each phase individually.
This not only requires exponentially many operations, but also gives each operation the maximal number of control wires, making them more expensive than less-controlled operations.

**Description of the Change:**
This PR implements a [new decomposition](https://www.notion.so/xanaduai/25-Q1-QSVT-compilation-177bc6bd1764818ca89dc3b0c74d36f7?pvs=4#190bc6bd1764809e8a4bf0ac69eb0fd0), based on an algorithm that decomposes the integer `dim` into sums and differences of powers of two. This corresponds to decomposing the projector $\Pi$ with rank `dim` of the PCPhase gate into projectors with ranks that are powers of two. The gates generated by the latter projectors are then simple (controlled) phases shifts, potentially combined with a (controlled) global phase that allows us to shift the phase of the state $|0\rangle$ instead of $|1\rangle$.

**Benefits:**
Much cheaper decomposition of `PCPhase`.

**Possible Drawbacks:**
N/A

**Related GitHub Issues:**

[sc-87926]

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7166?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
Attention: Patch coverage is `97.67442%` with `1 line` in your changes missing coverage. Please review.
> Project coverage is 99.62%. Comparing base [(`eda31c9`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/eda31c9b2518123a52238216fab4ee7fb49be347?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`d6545e8`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/d6545e8e001df39fa129d9552cf93824381b3abe?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

| [Files with missing lines](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7166?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Patch % | Lines |
|---|---|---|
| [pennylane/ops/qubit/parametric\_ops\_multi\_qubit.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7166?src=pr&el=tree&filepath=pennylane%2Fops%2Fqubit%2Fparametric_ops_multi_qubit.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9wYXJhbWV0cmljX29wc19tdWx0aV9xdWJpdC5weQ==) | 97.67% | [1 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7166?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) |

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


```diff
@@            Coverage Diff             @@
##           master    #7166      +/-   ##
==========================================
- Coverage   99.62%   99.62%   -0.01%     
==========================================
  Files         509      509              
  Lines       49602    49632      +30     
==========================================
+ Hits        49417    49446      +29     
- Misses        185      186       +1     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7166?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>

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 076400a4d..dfe0d5bcb 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.
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.


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)