Repository: PennyLaneAI/pennylane
Pull request #7070 (2025-03-11 14:37:49)

[Decomposition] decomposition rules for the `PhaseAdder` template
**Taken over by @astralcai**

**Context:**

**Description of the Change:**

Translate the decomposition rule of `PhaseAdder` to a quantum function compatible with the new system.

**Benefits:**

**Possible Drawbacks:**

**Related GitHub Issues:**

[sc-85617]
[sc-92927]


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7070?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 99.68%. Comparing base [(`480e124`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/480e1240d6b687233b12e3e07efea04f2797dcc8?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`835b093`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/835b093f1747dfc86f41f3cb312dc078d37eec6b?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

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


```diff
@@           Coverage Diff           @@
##           master    #7070   +/-   ##
=======================================
  Coverage   99.68%   99.68%           
=======================================
  Files         530      530           
  Lines       52140    52173   +33     
=======================================
+ Hits        51977    52010   +33     
  Misses        163      163           
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7070?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 1b08e2490..9a1fe7e37 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -4,6 +4,14 @@
 
 <h3>New features since last release</h3>
 
+* A new decomposition based on *unary iteration* has been added to :class:`qml.Select`.
+  This decomposition reduces the :class:`T` count significantly, and uses :math:`c-1`
+  auxiliary wires for a :class:`qml.Select` operation with :math:`c` control wires.
+  Unary iteration leverages these auxiliary wires to store intermediate values for reuse
+  among the different multi-controlled operators, avoiding unnecessary recomputation.
+  Check out the documentation for a thorough explanation.
+  [(#7623)](https://github.com/PennyLaneAI/pennylane/pull/7623)
+
 * A new function called :func:`qml.from_qasm3` has been added, which converts OpenQASM 3.0 circuits into quantum functions
   that can be subsequently loaded into QNodes and executed. 
   [(#7432)](https://github.com/PennyLaneAI/pennylane/pull/7432)
@@ -140,6 +148,7 @@
 * Two new functions called :func:`~.math.convert_to_su2` and :func:`~.math.convert_to_su4` have been added to `qml.math`, which convert unitary matrices to SU(2) or SU(4), respectively, and optionally a global phase.
   [(#7211)](https://github.com/PennyLaneAI/pennylane/pull/7211)
 
+
 * A new template :class:`~.TemporaryAND` has been added. The  :class:`~.TemporaryAND` (a.k.a.  :class:`~.Elbow`)
   operation is a three-qubit gate equivalent to an ``AND``, or reversible :class:`~pennylane.Toffoli`, gate
   that leverages extra information about the target wire to enable more efficient circuit decompositions.
@@ -258,13 +267,14 @@
   * :class:`~.PCPhase`
     [(#7591)](https://github.com/PennyLaneAI/pennylane/pull/7591)
 
-<<<<<<< qpe-decomp
   * :class:`~.QuantumPhaseEstimation`
     [(#7637)](https://github.com/PennyLaneAI/pennylane/pull/7637)
-=======
+
   * :class:`~.BasisRotation`
     [(#7074)](https://github.com/PennyLaneAI/pennylane/pull/7074)
->>>>>>> master
+
+  * :class:`~.PhaseAdder`
+    [(#7070)](https://github.com/PennyLaneAI/pennylane/pull/7070)
 
   * :class:`~.IntegerComparator`
     [(#7636)](https://github.com/PennyLaneAI/pennylane/pull/7636)
@@ -535,7 +545,7 @@
   these dependencies only have to be installed for the functions that use them, not to use
   ``labs`` functionalities in general. This decouples the various submodules, and even functions
   within the same submodule, from each other.
-  [(#7xxx)](https://github.com/PennyLaneAI/pennylane/pull/7xxx)
+  [(#7650)](https://github.com/PennyLaneAI/pennylane/pull/7650)
 
 * A new module :mod:`pennylane.labs.intermediate_reps <pennylane.labs.intermediate_reps>`
   provides functionality to compute intermediate representations for particular circuits.
@@ -737,8 +747,11 @@ Here's a list of deprecations made this release. For a more detailed breakdown o
 
 <h3>Documentation 📝</h3>
 
+* The usage examples for `qml.decomposition.DecompositionGraph` have been updated.
+  [(#7692)](https://github.com/PennyLaneAI/pennylane/pull/7692)
+
 * The entry in the :doc:`/news/program_capture_sharp_bits` has been updated to include
-  additional supported lightning devices: ``lightning.kokkos`` and ``lightning.gpu``.
+  additional supported lightning devices: `lightning.kokkos` and `lightning.gpu`.
   [(#7674)](https://github.com/PennyLaneAI/pennylane/pull/7674)
 
 * Updated the circuit drawing for `qml.Select` to include two commonly used symbols for 
diff --git a/pennylane/_version.py b/pennylane/_version.py
index 25bf54eac..e762c8eed 100644
--- a/pennylane/_version.py
+++ b/pennylane/_version.py
@@ -16,4 +16,4 @@
 Version number (major.minor.patch[-label])
 """
 
-__version__ = "0.42.0-dev53"
+__version__ = "0.42.0-dev54"
diff --git a/pennylane/decomposition/decomposition_graph.py b/pennylane/decomposition/decomposition_graph.py
index 81a448075..a160c1591 100644
--- a/pennylane/decomposition/decomposition_graph.py
+++ b/pennylane/decomposition/decomposition_graph.py
@@ -109,6 +109,8 @@ class DecompositionGraph:  # pylint: disable=too-many-instance-attributes
 
     .. code-block:: python
 
+        from pennylane.decomposition import DecompositionGraph
+
         op = qml.CRX(0.5, wires=[0, 1])
         graph = DecompositionGraph(
             operations=[op],
@@ -126,7 +128,7 @@ class DecompositionGraph:  # pylint: disable=too-many-instance-attributes
      CNOT(wires=[0, 1]),
      RZ(-1.5707963267948966, wires=[1])]
     >>> graph.resource_estimate(op)
-    <num_gates=10, gate_counts={RZ: 6, CNOT: 2, RX: 2}>
+    <num_gates=10, gate_counts={RZ: 6, CNOT: 2, RX: 2}, weighted_cost=10.0>
 
     """
 
@@ -381,7 +383,7 @@ class DecompositionGraph:  # pylint: disable=too-many-instance-attributes
          CNOT(wires=[0, 1]),
          RZ(-1.5707963267948966, wires=[1])]
         >>> graph.resource_estimate(op)
-        <num_gates=10, gate_counts={RZ: 6, CNOT: 2, RX: 2}>
+        <num_gates=10, gate_counts={RZ: 6, CNOT: 2, RX: 2}, weighted_cost=10.0>
 
         """
         if not self.is_solved_for(op):
@@ -407,7 +409,7 @@ class DecompositionGraph:  # pylint: disable=too-many-instance-attributes
 
         .. code-block:: python
 
-            op = qml.CRX(0.5, wires=[0, 1])
+            op = qml.CRY(0.2, wires=[0, 2])
             graph = DecompositionGraph(
                 operations=[op],
                 gate_set={"RZ", "RX", "CNOT", "GlobalPhase"},
@@ -418,12 +420,10 @@ class DecompositionGraph:  # pylint: disable=too-many-instance-attributes
         >>> with qml.queuing.AnnotatedQueue() as q:
         ...     rule(*op.parameters, wires=op.wires, **op.hyperparameters)
         >>> q.queue
-        [RZ(1.5707963267948966, wires=[1]),
-         RY(0.25, wires=[1]),
-         CNOT(wires=[0, 1]),
-         RY(-0.25, wires=[1]),
-         CNOT(wires=[0, 1]),
-         RZ(-1.5707963267948966, wires=[1])]
+        [RY(0.1, wires=[2]),
+         CNOT(wires=[0, 2]),
+         RY(-0.1, wires=[2]),
+         CNOT(wires=[0, 2])]
 
         """
         if not self.is_solved_for(op):
diff --git a/pennylane/decomposition/decomposition_rule.py b/pennylane/decomposition/decomposition_rule.py
index 4a01cccd0..13a3f6e61 100644
--- a/pennylane/decomposition/decomposition_rule.py
+++ b/pennylane/decomposition/decomposition_rule.py
@@ -17,7 +17,7 @@
 from __future__ import annotations
 
 import inspect
-from collections import defaultdict
+from collections import Counter, defaultdict
 from textwrap import dedent
 from typing import Callable, Optional, Type, overload
 
@@ -294,10 +294,13 @@ class DecompositionRule:
         """Computes the resources required to implement this decomposition rule."""
         if self._compute_resources is None:
             raise NotImplementedError("No resource estimation found for this decomposition rule.")
-        gate_counts = self._compute_resources(*args, **kwargs)
-        assert isinstance(gate_counts, dict), "Resource function must return a dictionary."
-        gate_counts = {_auto_wrap(op): count for op, count in gate_counts.items() if count > 0}
-        return Resources(gate_counts)
+        raw_gate_counts = self._compute_resources(*args, **kwargs)
+        assert isinstance(raw_gate_counts, dict), "Resource function must return a dictionary."
+        gate_counter = Counter()
+        for op, count in raw_gate_counts.items():
+            if count > 0:
+                gate_counter.update({_auto_wrap(op): count})
+        return Resources(dict(gate_counter))
 
     def is_applicable(self, *args, **kwargs) -> bool:
         """Checks whether this decomposition rule is applicable."""
diff --git a/pennylane/templates/subroutines/phase_adder.py b/pennylane/templates/subroutines/phase_adder.py
index 0b4bfd887..8ed5f6fef 100644
--- a/pennylane/templates/subroutines/phase_adder.py
+++ b/pennylane/templates/subroutines/phase_adder.py
@@ -18,6 +18,12 @@ Contains the PhaseAdder template.
 import numpy as np
 
 import pennylane as qml
+from pennylane.decomposition import (
+    add_decomps,
+    adjoint_resource_rep,
+    register_resources,
+    resource_rep,
+)
 from pennylane.operation import Operation
 from pennylane.wires import Wires, WiresLike
 
@@ -62,10 +68,11 @@ class PhaseAdder(Operation):
         x_wires (Sequence[int]): the wires the operation acts on. The number of wires must be enough
             for a binary representation of the value being targeted, :math:`x`. In some cases an additional
             wire is needed, see usage details below. The number of wires also limits the maximum
-            value for `mod`.
+            value for ``mod``.
         mod (int): the modulo for performing the addition. If not provided, it will be set to its maximum value, :math:`2^{\text{len(x_wires)}}`.
         work_wire (Sequence[int] or int): the auxiliary wire to use for the addition. Optional
-            when `mod` is :math:`2^{len(x\_wires)}`. Defaults to empty tuple.
+            when ``mod`` is :math:`2^{\text{len(x_wires)}}`. Defaults to empty tuple. The work wire must
+            be in the :math:`|0\rangle` state.
 
     **Example**
 
@@ -124,9 +131,10 @@ class PhaseAdder(Operation):
 
     grad_method = None
 
-    def __init__(
-        self, k, x_wires: WiresLike, mod=None, work_wire: WiresLike = (), id=None
-    ):  # pylint: disable=too-many-arguments
+    resource_keys = {"num_x_wires", "mod"}
+
+    # pylint: disable=too-many-arguments
+    def __init__(self, k, x_wires: WiresLike, mod=None, work_wire: WiresLike = (), id=None):
 
         work_wire = Wires(() if work_wire is None else work_wire)
         x_wires = Wires(x_wires)
@@ -153,13 +161,18 @@ class PhaseAdder(Operation):
                         "None of the wires in work_wire should be included in x_wires."
                     )
 
-        all_wires = x_wires + work_wire
-
         self.hyperparameters["k"] = k % mod
         self.hyperparameters["mod"] = mod
         self.hyperparameters["work_wire"] = work_wire
         self.hyperparameters["x_wires"] = x_wires
-        super().__init__(wires=all_wires, id=id)
+        super().__init__(wires=x_wires, id=id)
+
+    @property
+    def resource_params(self) -> dict:
+        return {
+            "num_x_wires": len(self.hyperparameters["x_wires"]),
+            "mod": self.hyperparameters["mod"],
+        }
 
     @property
     def num_params(self):
@@ -244,3 +257,53 @@ class PhaseAdder(Operation):
             op_list.extend(_add_k_fourier(k, x_wires))
 
         return op_list
+
+
+def _phase_adder_decomposition_resources(num_x_wires, mod) -> dict:
+
+    if mod == 2**num_x_wires:
+        return {qml.PhaseShift: num_x_wires}
+
+    return {
+        qml.PhaseShift: 2 * num_x_wires,
+        adjoint_resource_rep(qml.PhaseShift, {}): 2 * num_x_wires,
+        resource_rep(qml.QFT, num_wires=num_x_wires): 2,
+        adjoint_resource_rep(qml.QFT, {"num_wires": num_x_wires}): 2,
+        qml.CNOT: 2,
+        qml.X: 2,
+        qml.ControlledPhaseShift: num_x_wires,
+    }
+
+
+# pylint: disable=no-value-for-parameter
+@register_resources(_phase_adder_decomposition_resources)
+def _phase_adder_decomposition(k, x_wires: WiresLike, mod, work_wire, **__):
+
+    n_wires = len(x_wires)
+
+    @qml.for_loop(n_wires)
+    def _add_k_fourier_loop(i, _k):
+        qml.PhaseShift(_k * np.pi / (2**i), wires=x_wires[i])
+        return _k
+
+    if mod == 2**n_wires:
+        _add_k_fourier_loop(k)
+        return
+
+    aux_k = x_wires[0]
+    _add_k_fourier_loop(k)
+    qml.adjoint(_add_k_fourier_loop)(mod)
+    qml.adjoint(qml.QFT)(wires=x_wires)
+    qml.CNOT(wires=[aux_k, work_wire[0]])
+    qml.QFT(wires=x_wires)
+    qml.ctrl(_add_k_fourier_loop, control=work_wire)(mod)
+    qml.adjoint(_add_k_fourier_loop)(k)
+    qml.adjoint(qml.QFT)(wires=x_wires)
+    qml.X(aux_k)
+    qml.CNOT(wires=[aux_k, work_wire[0]])
+    qml.X(aux_k)
+    qml.QFT(wires=x_wires)
+    _add_k_fourier_loop(k)
+
+
+add_decomps(PhaseAdder, _phase_adder_decomposition)
diff --git a/pennylane/templates/subroutines/select.py b/pennylane/templates/subroutines/select.py
index 7d14dd5c6..392d32b7c 100644
--- a/pennylane/templates/subroutines/select.py
+++ b/pennylane/templates/subroutines/select.py
@@ -15,18 +15,28 @@
 Contains the Select template.
 """
 
-
 import copy
 import itertools
+from collections import Counter
 
 import pennylane as qml
 from pennylane import math
+from pennylane.decomposition import (
+    add_decomps,
+    adjoint_resource_rep,
+    controlled_resource_rep,
+    register_resources,
+)
+from pennylane.decomposition.resources import resource_rep
 from pennylane.operation import Operation
+from pennylane.ops import CNOT, X, adjoint, ctrl
+from pennylane.templates.subroutines.temporary_and import TemporaryAND
+from pennylane.wires import Wires
 
 
 class Select(Operation):
-    r"""Applies different operations depending on the state of
-    designated control qubits.
+    r"""The ``Select`` operator, also known as multiplexer or multiplexed operation,
+    applies different operations depending on the state of designated control wires.
 
     .. math:: Select|i\rangle \otimes |\psi\rangle = |i\rangle \otimes U_i |\psi\rangle
 
@@ -35,16 +45,16 @@ class Select(Operation):
                     :width: 70%
                     :target: javascript:void(0);
 
-    This operator is also known as **multiplexer**, or multiplexed operation.
     If the applied operations :math:`\{U_i\}` are all single-qubit Pauli rotations about the
-    same axis, with the angle determined by the control qubits, this is also called a
+    same axis, with the angle determined by the control wires, this is also called a
     **uniformly controlled rotation** gate.
 
     .. seealso:: :class:`~.SelectPauliRot`
 
     Args:
         ops (list[Operator]): operations to apply
-        control (Sequence[int]): the wires controlling which operation is applied
+        control (Sequence[int]): the wires controlling which operation is applied.
+            At least :math:`\lceil \log_2 K\rceil` wires are required for :math:`K` operations.
         id (str or None): String representing the operation (optional)
 
     .. note::
@@ -70,10 +80,202 @@ class Select(Operation):
     2: ─╰X─│──╰Y─├SWAP─┤  State
     3: ────╰X────╰SWAP─┤  State
 
+    .. details::
+        :title: Unary iterator decomposition
+
+        Generically, ``Select`` is decomposed into one multi-controlled operator for each target
+        operator. However, if auxiliary wires are available, a decomposition using a
+        "unary iterator" can be applied. It was introduced by
+        `Babbush et al. (2018) <https://arxiv.org/abs/1805.03662>`__.
+
+        **Principle**
+
+        Unary iteration leverages auxiliary wires to store intermediate values for reuse between
+        the different multi-controlled operators, avoiding unnecessary recomputation.
+        In addition to this caching functionality, unary iteration reduces the cost of the
+        computation directly, because the involved reversible AND (or Toffoli) gates can be
+        implemented at lower cost if the target is known to be in the :math:`|0\rangle` state
+        (see :class:`~TemporaryAND`).
+
+        For :math:`K` operators to be Select-applied, :math:`c=\lceil\log_2 K\rceil` control
+        wires are required. Unary iteration demands an additional :math:`c-1` auxiliary wires.
+        Below we first show an example for :math:`K` being a power of two, i.e., :math:`K=2^c`.
+        Then we elaborate on implementation details for the case :math:`K<2^c`, which we call
+        a *partial Select* operator.
+
+        **Example**
+
+        Assume that we want to Select-apply :math:`K=8=2^3` operators to two target wires,
+        which requires :math:`c=\lceil \log_2 K\rceil=3` control wires. The generic
+        decomposition for this takes the form
+
+        .. code-block::
+
+            0: ─╭○─────╭○─────╭○─────╭○─────╭●─────╭●─────╭●─────╭●─────┤
+            1: ─├○─────├○─────├●─────├●─────├○─────├○─────├●─────├●─────┤
+            2: ─├○─────├●─────├○─────├●─────├○─────├●─────├○─────├●─────┤
+            3: ─├U(M0)─├U(M1)─├U(M2)─├U(M3)─├U(M4)─├U(M5)─├U(M6)─├U(M7)─┤
+            4: ─╰U(M0)─╰U(M1)─╰U(M2)─╰U(M3)─╰U(M4)─╰U(M5)─╰U(M6)─╰U(M7)─┤.
+
+        Unary iteration then uses :math:`c-1=2` auxiliary wires, denoted ``aux0`` and ``aux1``
+        below, to first rewrite the control structure:
+
+        .. code-block::
+
+            0:    ─╭○───────○╮─╭○───────○╮─╭○───────○╮─╭○───────○╮─╭●───────●╮─╭●───────●╮─╭●───────●╮─╭●───────●╮─┤
+            1:    ─├○───────○┤─├○───────○┤─├●───────●┤─├●───────●┤─├○───────○┤─├○───────○┤─├●───────●┤─├●───────●┤─┤
+            aux0:  ╰─╭●───●╮─╯ ╰─╭●───●╮─╯ ╰─╭●───●╮─╯ ╰─╭●───●╮─╯ ╰─╭●───●╮─╯ ╰─╭●───●╮─╯ ╰─╭●───●╮─╯ ╰─╭●───●╮─╯ │
+            2:    ───├○───○┤─────├●───●┤─────├○───○┤─────├●───●┤─────├○───○┤─────├●───●┤─────├○───○┤─────├●───●┤───┤
+            aux1:    ╰─╭●──╯     ╰─╭●──╯     ╰─╭●──╯     ╰─╭●──╯     ╰─╭●──╯     ╰─╭●──╯     ╰─╭●──╯     ╰─╭●──╯   │
+            3:    ─────├U(M0)──────├U(M1)──────├U(M2)──────├U(M3)──────├U(M4)──────├U(M5)──────├U(M6)──────├U(M7)──┤
+            4:    ─────╰U(M0)──────╰U(M1)──────╰U(M2)──────╰U(M3)──────╰U(M4)──────╰U(M5)──────╰U(M6)──────╰U(M7)──┤
+
+        Here, we used the symbols
+
+        .. code-block::
+
+            0: ─╭●──       ─●─╮─
+            1: ─├●──  and  ─●─┤─
+            2:  ╰───       ───╯
+
+        for :class:`~.TemporaryAND` and its adjoint, respectively, and skipped drawing the
+        auxiliary wires in areas where they are guaranteed to be in the state :math:`|0\rangle`.
+        We will need three simplification rules for pairs of ``TemporaryAND`` gates:
+
+        .. code-block::
+
+            ─○─╮─╭○──   ──     ─○─╮─╭○──   ─╭○─       ─○─╮─╭●──   ─╭●────
+            ─○─┤─├○── = ──,    ─○─┤─├●── = ─│──, and  ─●─┤─├○── = ─│──╭●─.
+            ───╯ ╰───   ──     ───╯ ╰───   ─╰X─       ───╯ ╰───   ─╰X─╰X─
+
+        Applying these simplifications reduces the computational cost of the ``Select``
+        template:
+
+        .. code-block::
+
+            0:    ─╭○────────────────╭○──────────────────╭●─────────────────────╭●─────────────────●╮─┤
+            1:    ─├○────────────────│───────────────────│──╭●──────────────────│──────────────────●┤─┤
+            aux0:  ╰─╭●─────╭●────●╮─╰X─╭●─────╭●─────●╮─╰X─╰X─╭●─────╭●─────●╮─╰X─╭●─────╭●─────●╮─╯ │
+            2:    ───├○─────│─────●┤────├○─────│──────●┤───────├○─────│──────●┤────├○─────│──────●┤───┤
+            aux1:    ╰─╭●───╰X─╭●──╯    ╰─╭●───╰X──╭●──╯       ╰─╭●───╰X──╭●──╯    ╰─╭●───╰X──╭●──╯   │
+            3:    ─────├U(M0)──├U(M1)─────├U(M2)───├U(M3)────────├U(M4)───├U(M5)─────├U(M6)───├U(M7)──┤
+            4:    ─────╰U(M0)──╰U(M1)─────╰U(M2)───╰U(M3)────────╰U(M4)───╰U(M5)─────╰U(M6)───╰U(M7)──┤
+
+        An additional cost reduction then results from the fact that the ``TemporaryAND``
+        gate and its adjoint require four and zero :class:`~T` gates, respectively,
+        in contrast to the seven ``T`` gates required by a decomposition of :class:`~Toffoli`.
+
+        For general :math:`c` and :math:`K=2^c`, the decomposition takes a similar form, with
+        alternating control and auxiliary wires.
+
+        An implementation of the unary iterator is achieved in the following steps:
+        We first define a recursive sub-circuit ``R``;
+        given :math:`L` operators and :math:`2 \lceil\log_2(L)\rceil + 1` control and
+        auxiliary wires, there are three cases that ``R`` distinguishes. First, if ``L>1``,
+        it applies the circuit
+
+        .. code-block::
+
+            aux_j:   ╭R   ─╭●────╭●────●─╮─
+            j+1:     ├R = ─├○────│─────●─┤─
+            aux_j+1: ╰R    ╰──R──╰X─R────╯ ,
+
+        where each label ``R`` symbolizes a call to ``R`` itself, on the next recursion level.
+        These next-level calls use
+        :math:`L' = 2^{\lceil\log_2(L)\rceil-1}` (i.e. half of :math:`L`, rounded up to the next
+        power of two) and :math:`L-L'` (i.e. the rest) operators, respectively.
+
+        Second, if ``L=1``, the single operator is applied, controlled on the first control wire.
+        Finally, if ``L=0``, ``R`` does not apply any operators.
+
+        With ``R`` defined, we are ready to outline the main circuit structure:
+
+        #. Apply the left-most ``TemporaryAND`` controlled on qubits ``0`` and ``1``.
+        #. Split the target operators into four "quarters" (often with varying sizes)
+           and apply the first quarter using ``R``.
+        #. Apply ``[X(0), CNOT([0, "aux0"]), X(0)]``.
+        #. Apply the second quarter using ``R``.
+        #. Apply ``[CNOT([0, "aux0"]), CNOT([1, "aux0"])]``.
+        #. Apply the third quarter using ``R``.
+        #. Apply ``[CNOT([0, "aux0"])]``.
+        #. Apply the last quarter using ``R``.
+        #. Apply the right-most ``adjoint(TemporaryAND)`` controlled on qubits ``0`` and ``1``.
+
+        **Partial Select decomposition**
+
+        The unary iterator decomposition of the ``Select`` template can be
+        simplified further if both of the following criteria are met:
+
+        #. There are fewer target operators than would maximally be possible for the given
+           number of control wires, i.e. :math:`K<2^c`.
+
+        #. The state :math:`|\psi\rangle` of the control wires satisfies
+           :math:`\langle j | \psi\rangle=0` for all computational basis states with :math:`j\geq K`.
+
+        We do not derive this reduction here but discuss the modifications to the implementation
+        above that result from it.
+
+        Given :math:`K=2^c-b` operators, where :math:`c` is defined as above and we
+        have :math:`0\leq b<2^{c-1}`, the nine steps above are modified into one of three variants.
+        In each variant, the first :math:`2^{c-1}` operators are applied in two equal portions,
+        containing :math:`2^{c-2}` operators each.
+        After this, :math:`\ell=2^{c-1} -b` operators remain and the three circuit variants are
+        distinguished, based on :math:`\ell`:
+
+        - if :math:`\ell \geq 2^{c-2}`, the following, rather generic, circuit is applied:
+
+          .. code-block::
+
+              0:    ─╭○─────╭○─────╭●────────╭●─────●─╮─
+              1:    ─├○─────│──────│──╭●─────│──────●─┤─
+              aux0:  ╰──╭R──╰X─╭R──╰X─╰X─╭R──╰X─╭R────╯
+              2:    ────├R─────├R────────├R─────├R──────
+              aux1:     ╰R     ╰R        ╰R     ╰R      .
+
+          Here, each operator with three ``R`` labels symbolizes a call to ``R``. The first
+          call in the second half applies :math:`2^{\lceil\log_2(\ell)\rceil-1}` operators.
+          Note that this case is triggered if :math:`K` is larger than or equal to
+          :math:`\tfrac{3}{4}` of the maximal capacity for :math:`c` control wires.
+          Also note how the two middle ``TemporaryAND`` gates were merged into two CNOTs,
+          like for the non-partial Select operator.
+
+        - if :math:`1<\ell < 2^{c-2}`, the following circuit is applied:
+
+          .. code-block::
+
+              0:    ─╭○─────╭○─────○─╮╭●─────╭●─────●─╮─
+              1:    ─├○─────│──────●─┤│──────│────────│─
+              aux0:  ╰──╭R──╰X─╭R────╯│      │        │
+              2:    ────├R─────├R─────├○─────│──────●─┤─
+              aux1:     ╰R     ╰R     ╰───R──╰X──R────╯
+
+          where the second half may skip more than one control and auxiliary wire each.
+          In this diagram, both the operators with three and one ``R`` labels represent calls to
+          ``R``, with single-label instances applying fewer operators.
+          The first call to ``R`` in the second half applies :math:`2^{\lceil\log_2(\ell)\rceil-1}`
+          operators. The middle elbows act on distinct wire triples and can not be merged as
+          above.
+
+        - if :math:`\ell=1`, the following circuit is applied:
+
+          .. code-block::
+
+              0:    ─╭○─────╭○─────○─╮╭●──
+              1:    ─├○─────│──────●─┤│───
+              aux0:  ╰──╭R──╰X─╭R────╯│───
+              2:    ────├R─────├R─────│───
+              aux1:     ╰R     ╰R     ╰U  .
+
+          Here, the three connected ``R`` labels symbolize a call to ``R`` and
+          apply :math:`2^{c-2}` operators each.
+          The controlled gate on the right applies the single remaining operator.
+
     """
 
+    resource_keys = set(["ops"])
+
     def _flatten(self):
-        return (self.ops), (self.control)
+        return (self.ops), (self.control, self.hyperparameters["work_wires"])
 
     @classmethod
     def _primitive_bind_call(cls, *args, **kwargs):
@@ -81,33 +283,34 @@ class Select(Operation):
 
     @classmethod
     def _unflatten(cls, data, metadata) -> "Select":
-        return cls(data, metadata)
+        return cls(data, control=metadata[0], work_wires=metadata[1])
 
     def __repr__(self):
         return f"Select(ops={self.ops}, control={self.control})"
 
-    def __init__(self, ops, control, id=None):
-        control = qml.wires.Wires(control)
+    def __init__(self, ops, control, work_wires=None, id=None):
+        control = Wires(control)
+        work_wires = Wires(() if work_wires is None else work_wires)
         self.hyperparameters["ops"] = tuple(ops)
         self.hyperparameters["control"] = control
+        self.hyperparameters["work_wires"] = work_wires
 
         if 2 ** len(control) < len(ops):
             raise ValueError(
                 f"Not enough control wires ({len(control)}) for the desired number of "
-                + f"operations ({len(ops)}). At least {int(math.ceil(math.log2(len(ops))))} control "
+                + f"operations ({len(ops)}). At least {_ceil_log(len(ops))} control "
                 + "wires required."
             )
 
         if any(
-            control_wire in qml.wires.Wires.all_wires([op.wires for op in ops])
-            for control_wire in control
+            control_wire in Wires.all_wires([op.wires for op in ops]) for control_wire in control
         ):
             raise ValueError("Control wires should be different from operation wires.")
 
         for op in ops:
             qml.QueuingManager.remove(op)
 
-        target_wires = qml.wires.Wires.all_wires([op.wires for op in ops])
+        target_wires = Wires.all_wires([op.wires for op in ops])
         self.hyperparameters["target_wires"] = target_wires
 
         all_wires = target_wires + control
@@ -205,7 +408,7 @@ class Select(Operation):
         """
         states = list(itertools.product([0, 1], repeat=len(control)))
         decomp_ops = [
-            qml.ctrl(op, control, control_values=states[index]) for index, op in enumerate(ops)
+            ctrl(op, control, control_values=states[index]) for index, op in enumerate(ops)
         ]
         return decomp_ops
 
@@ -221,10 +424,227 @@ class Select(Operation):
 
     @property
     def target_wires(self):
-        """The wires of the input operators."""
+        """The wires of the target operators."""
         return self.hyperparameters["target_wires"]
 
     @property
     def wires(self):
         """All wires involved in the operation."""
         return self.hyperparameters["control"] + self.hyperparameters["target_wires"]
+
+    @property
+    def resource_params(self) -> dict:
+        return {
+            "ops": [
+                resource_rep(type(op), **op.resource_params) for op in self.hyperparameters["ops"]
+            ]
+        }
+
+
+def _ceil(a):
+    return int(math.ceil(a))
+
+
+def _ceil_log(a):
+    return _ceil(math.log2(a))
+
+
+def _add_first_k_units(ops, controls, work_wires, k):
+    """Add all controlled-applied operators within the unary iterator scheme.
+
+    This function is used for the outer-most recursion level, and then calls _add_k_units
+    for the inner recursion levels.
+    In the documentation, this is the outer level that is being described in detail, and
+    ``_add_k_units`` corresponds to the subroutine ``R``.
+
+    """
+    assert k == len(ops) > 2
+
+    needed_controls = 2 * _ceil_log(k) - 1
+    assert len(controls) >= needed_controls, f"{len(controls)=}, {needed_controls=}"
+    controls = controls[:needed_controls]
+
+    and_wires = controls[:3]
+    new_work_wires = work_wires + controls[:2]
+    new_controls = controls[2:]
+
+    a = _ceil_log(k)  # a >= 2 because k>2 by assertion above
+    k01 = 2 ** (a - 1)  # First half of circuit will implement 2^(a-1)>=2 operators
+    k0 = k1 = 2 ** (a - 2)  # First two quarters of circuit each implement 2^(a-2)>=1 operator(s).
+    l = k - k01
+    k2 = _ceil(2 ** (_ceil_log(l) - 1))
+    k3 = k - k01 - k2
+
+    # Open TemporaryAND (controlled on |00>) + first quarter + CX (controlled on |0>) + second quarter
+    first_half = (
+        [TemporaryAND(and_wires, control_values=(0, 0))]
+        + _add_k_units(ops[:k0], new_controls, new_work_wires, k0)
+        + [X(controls[0])]
+        + [CNOT([controls[0], controls[2]])]
+        + [X(controls[0])]
+        + _add_k_units(ops[k0:k01], new_controls, new_work_wires, k1)
+    )
+
+    if l == 1:  # first variant
+
+        # Single operation left to apply: Only the third quarter will be needed, and it will not need
+        # TemporaryAND gates at all
+        and_wires_sec_half = []
+        new_controls_sec_half = controls
+        new_work_wires_sec_half = work_wires
+        # Closing TemporaryAND for first half
+        middle_part = [adjoint(TemporaryAND)(and_wires, control_values=(0, 1))]
+
+    else:
+        c_bar = 2 * (_ceil_log(k) - _ceil_log(k - k01) - 1)
+        and_wires_sec_half = [controls[0], controls[c_bar + 1], controls[c_bar + 2]]
+        new_controls_sec_half = controls[c_bar + 2 :]
+        new_work_wires_sec_half = work_wires + controls[: c_bar + 2]
+
+        if c_bar > 0:  # second variant
+            # Closing TemporaryAND for first half, opening TemporaryAND for second half
+            middle_part = [
+                adjoint(TemporaryAND)(and_wires, control_values=(0, 1)),
+                TemporaryAND(and_wires_sec_half, control_values=(1, 0)),
+            ]
+        else:  # third variant
+            middle_part = [CNOT(and_wires[::2]), CNOT(and_wires[1:])]
+
+    second_half = _add_k_units(
+        ops[k01 : k01 + k2], new_controls_sec_half, new_work_wires_sec_half, k2
+    )
+    if and_wires_sec_half:
+        second_half += (
+            [CNOT(and_wires_sec_half[::2])]
+            + _add_k_units(ops[k0 + k1 + k2 :], new_controls_sec_half, new_work_wires_sec_half, k3)
+            + [adjoint(TemporaryAND)(and_wires_sec_half)]
+        )
+
+    return first_half + middle_part + second_half
+
+
+def _add_k_units(ops, controls, work_wires, k):
+    """Add k controlled-applied operators within the unary iterator scheme, in a recursive
+    manner.
+
+    In the documentation, this subroutine is called ``R``.
+    This is _not_ used for the outer-most recursion level, see _add_first_k_units instead.
+
+    We are given ``K=len(ops)`` operators and ``2 * ⌈log_2(K)⌉ + 1`` control and auxiliary wires.
+    If ``K=0``, nothing is applied.
+    If ``K=1``, the single operator is applied, controlled on the first control wire.
+
+    In all other cases, this function applies the circuit
+
+    .. code-block::
+
+        ─╭●────╭●────●─╮─
+        ─├○────│─────●─┤─
+         ╰──■──╰X─■────╯
+
+    where each box symbolizes calls to itself on the next recursion level.
+    The next-level calls to ``_add_k_units` use
+
+    ``k_first = 2 ** (⌈log_2(k)⌉-1)`` (i.e. half of ``k``, rounded up to the next power of two)
+    and
+    ``k_second = k-k_first`` (i.e. the rest)
+
+    operators, respectively. Accordingly, two fewer control wires are used.
+
+    """
+    assert k == len(ops) > 0
+    num_bits = _ceil_log(k)
+    needed_controls = 2 * num_bits + 1
+    assert len(controls) >= needed_controls, f"{len(controls)=}, {needed_controls=}"
+
+    if k == 1:
+        assert num_bits == 0
+        return [ctrl(ops[0], control=controls[0], control_values=[1], work_wires=work_wires)]
+
+    assert num_bits != 0
+    controls = controls[:1] + controls[-(needed_controls - 1) :]
+
+    and_wires = controls[:3]
+    new_work_wires = work_wires + controls[:2]
+    new_controls = controls[2:]
+    k_first = 2 ** (num_bits - 1)
+    return (
+        [TemporaryAND(and_wires, control_values=(1, 0))]
+        + _add_k_units(ops[:k_first], new_controls, new_work_wires, k_first)
+        + [CNOT(and_wires[::2])]
+        + _add_k_units(ops[k_first:], new_controls, new_work_wires, k - k_first)
+        + [adjoint(TemporaryAND)(and_wires)]
+    )
+
+
+def _unary_select_resources(ops):
+    num_ops = len(ops)
+    cnt = Counter()
+
+    if num_ops / 2 ** _ceil_log(num_ops) > 3 / 4:
+        cnt.update(
+            {
+                resource_rep(TemporaryAND): num_ops - 3,
+                adjoint_resource_rep(TemporaryAND): num_ops - 3,
+                CNOT: num_ops,
+                X: 2,
+            }
+        )
+    else:
+        cnt.update(
+            {
+                resource_rep(TemporaryAND): num_ops - 2,
+                adjoint_resource_rep(TemporaryAND): num_ops - 2,
+                CNOT: num_ops - 2,
+                X: 2,
+            }
+        )
+
+    for op in ops:
+        key = controlled_resource_rep(op.op_type, op.params, num_control_wires=1)
+        print(key)
+        cnt[key] += 1
+
+    return dict(cnt)
+
+
+@register_resources(_unary_select_resources)
+def _unary_select(ops, control, work_wires, **_):
+    r"""This function reproduces the unary iterator behaviour in https://arxiv.org/abs/1805.03662.
+    For :math:`K` operators this decomposition requires at least :math:`c=\lceil\log_2 K\rceil`
+    control wires (as usual for Select), and :math:`c-1` additional work wires.
+    See the documentation of ``Select`` for details.
+
+    .. note::
+
+        This decomposition assumes that the state on the control wires does not have any overlap
+        with :math:`|i\rangle` for :math:`i\geq K`.
+    """
+    # Note that there might be a slightly improved implementation if a modified embedding is used
+    # during state preparation in a PrepSelPrep pattern.
+
+    if len(ops) == 0:
+        return []
+
+    min_num_controls = max(_ceil_log(len(ops)), 1)
+    assert len(control) >= min_num_controls
+    control = control[-min_num_controls:]
+    if len(work_wires) < len(control) - 1:
+        raise ValueError(
+            f"Can't use this decomposition with less than {len(control) - 1} work wires for {len(control)} controls."
+        )
+    if 1 <= len(ops) <= 2:
+        # Don't need unary iterator, just control-apply the one/two operator(s) directly.
+        return [
+            ctrl(op, control=control[0], control_values=[i], work_wires=work_wires)
+            for i, op in enumerate(ops)
+        ]
+    aux_control = [control[0]]
+    for i in range(min_num_controls - 1):
+        aux_control.append(control[i + 1])
+        aux_control.append(work_wires[i])
+    work_wires = work_wires[min_num_controls - 1 :]
+    return _add_first_k_units(ops, aux_control, work_wires, len(ops))
+
+
+add_decomps(Select, _unary_select)


Test patch:
diff --git a/tests/templates/test_subroutines/test_phase_adder.py b/tests/templates/test_subroutines/test_phase_adder.py
index ee30e530c..8c9eded78 100644
--- a/tests/templates/test_subroutines/test_phase_adder.py
+++ b/tests/templates/test_subroutines/test_phase_adder.py
@@ -19,6 +19,7 @@ import pytest
 
 import pennylane as qml
 from pennylane import numpy as np
+from pennylane.ops.functions.assert_valid import _test_decomposition_rule
 from pennylane.templates.subroutines.phase_adder import _add_k_fourier
 
 
@@ -227,6 +228,7 @@ class TestPhaseAdder:
 
     def test_decomposition(self):
         """Test that compute_decomposition and decomposition work as expected."""
+
         k = 4
         x_wires = [1, 2, 3]
         mod = 7
@@ -256,6 +258,17 @@ class TestPhaseAdder:
         for op1, op2 in zip(phase_adder_decomposition, op_list):
             qml.assert_equal(op1, op2)
 
+    @pytest.mark.parametrize("mod", [7, 8])
+    def test_decomposition_new(self, mod):
+        """Tests the decomposition rule implemented with the new system."""
+
+        k = 4
+        x_wires = [1, 2, 3]
+        work_wire = [0]
+        op = qml.PhaseAdder(k, x_wires, mod, work_wire)
+        for rule in qml.list_decomps(qml.PhaseAdder):
+            _test_decomposition_rule(op, rule)
+
     def test_work_wires_added_correctly(self):
         """Test that no work wires are added if work_wire = None"""
         wires = qml.PhaseAdder(1, x_wires=[1, 2]).wires
diff --git a/tests/templates/test_subroutines/test_select.py b/tests/templates/test_subroutines/test_select.py
index eaa990d78..12b3fc61b 100644
--- a/tests/templates/test_subroutines/test_select.py
+++ b/tests/templates/test_subroutines/test_select.py
@@ -19,17 +19,21 @@ import copy
 
 import numpy as np
 import pytest
+from scipy.stats import unitary_group
 
 import pennylane as qml
 from pennylane import numpy as pnp
+from pennylane.templates.subroutines.select import _unary_select
 
 
-def test_standard_checks():
+@pytest.mark.parametrize("num_ops", [3, 10, 15, 16])
+def test_standard_checks(num_ops):
     """Run standard validity tests."""
-    ops = [qml.PauliX(0), qml.PauliY(0)]
-    control = [1]
+    ops = [qml.PauliX(0) for _ in range(num_ops)]
+    control = [1, 2, 3, 4]
+    work_wires = [5, 6, 7]
 
-    op = qml.Select(ops, control)
+    op = qml.Select(ops, control, work_wires)
     assert op.target_wires == qml.wires.Wires(0)
     qml.ops.functions.assert_valid(op)
 
@@ -425,3 +429,98 @@ class TestInterfaces:
         jit_circuit = jax.jit(circuit)
 
         assert qml.math.allclose(circuit(), jit_circuit())
+
+
+num_controls_and_num_ops = (
+    [(nc, i) for nc in range(1, 5) for i in range(1, 2**nc + 1)]
+    + [(5, 1), (5, 2), (5, 17), (5, 24), (5, 31)]
+    + [(6, 3), (6, 8), (6, 27)]
+)
+
+
+class TestUnaryIterator:
+    """Tests for the auxiliary qubit-based unary iterator decomposition of Select."""
+
+    @pytest.mark.parametrize("num_controls", [0, 1, 2, 3])
+    def test_no_ops(self, num_controls):
+        """Test that the unary iterator does not return any operators for an empty list
+        of target operators."""
+
+        control = list(range(num_controls))
+        work = list(range(num_controls, 2 * num_controls - 1))
+
+        decomp = _unary_select([], control=control, work_wires=work)
+        assert decomp == []
+
+    @pytest.mark.parametrize("num_controls, num_ops", num_controls_and_num_ops)
+    def test_identity_with_basis_states(self, num_controls, num_ops):
+        """Test that the unary iterator is correct by asserting that the identity
+        matrix is created by preparing the i-th computational basis state conditioned on the
+        i-th basis state in the control qubits."""
+
+        dev = qml.device("default.qubit")
+
+        # Create angle set so that feeding angles[i] into RX on the i-th control wire will
+        # yield broadcasted BasisEmbedding (which does not support broadcasting atm)
+        angles = [list(map(int, np.binary_repr(i, width=num_controls))) for i in range(num_ops)]
+        angles = np.pi * np.array(angles).T
+        control = list(range(num_controls))
+        work = list(range(num_controls, 2 * num_controls - 1))
+        target = list(range(2 * num_controls - 1, 3 * num_controls - 1))
+
+        ops = [qml.BasisEmbedding(i, wires=target) for i in range(num_ops)]
+
+        @qml.qnode(dev)
+        def circuit():
+            for w, angle in zip(control, angles, strict=True):
+                qml.RX(angle, w)
+            _unary_select(ops, control=control, work_wires=work)
+            return qml.probs(target)
+
+        probs = circuit()
+        assert np.allclose(probs, np.eye(2**num_controls)[:num_ops])
+
+    @pytest.mark.parametrize(
+        ("num_ops", "control", "work", "msg_match"),
+        [(9, 4, 1, "Can't use this decomposition")],
+    )
+    def test_operation_and_test_wires_error(
+        self, num_ops, control, work, msg_match
+    ):  # pylint: disable=too-many-arguments
+        """Test that proper errors are raised"""
+
+        wires = qml.registers({"target": num_ops, "control": control, "work": work})
+        ops = [qml.BasisEmbedding(i, wires=wires["target"]) for i in range(num_ops)]
+
+        with pytest.raises(ValueError, match=msg_match):
+            _unary_select(ops, control=wires["control"], work_wires=wires["work"])
+
+    @pytest.mark.parametrize("num_controls, num_ops", num_controls_and_num_ops)
+    def test_comparison_with_select(self, num_controls, num_ops, seed):
+        """Test that the unary iterator is correct by comparing it to the standard Select
+        decomposition."""
+
+        angles = [list(map(int, np.binary_repr(i, width=num_controls))) for i in range(num_ops)]
+        angles = np.pi * np.array(angles).T
+        control = list(range(num_controls))
+        work = list(range(num_controls, 2 * num_controls - 1))
+        target = [2 * num_controls - 1, 2 * num_controls]
+
+        dev = qml.device("default.qubit", wires=control + work + target)
+        unitaries = unitary_group.rvs(4, size=num_ops, random_state=seed)
+        if num_ops == 1:
+            unitaries = np.array([unitaries])
+        ops = [qml.QubitUnitary(U, wires=target) for U in unitaries]
+        adj_ops = [qml.QubitUnitary(U.conj().T, wires=target) for U in unitaries]
+
+        @qml.qnode(dev)
+        def circuit():
+            for w, angle in zip(control, angles, strict=True):
+                qml.RX(angle, w)
+            _unary_select(ops, control=control, work_wires=work)
+            qml.Select(adj_ops, control=control, work_wires=None)
+            return qml.probs(target)
+
+        probs = circuit()
+        exp = np.eye(4)[0]
+        assert np.allclose(probs, exp)