Repository: PennyLaneAI/pennylane
Pull request #7793 (2025-07-05 01:36:10)

Global phase improvements to RS decomposition
**Context:** This PR fixes the missing global phase information in the `rs_decomposition`.

**Description of the Change:** This essentially does the following two key changes:
1. More robust global phase tracking via evolution of elements for the  `SU2` description of the decomposed operator. 
2. More robust custom-handling of $\pm(2k+1)\frac{\pi}{4}$ cases. This comes from the fact that $\exp{i\pi/8}$ cannot be described in either of the rings $\mathbb{D}[\omega]$ and $\mathbb{Z}[\omega]$, leading to a non-trivial decomposition attempt which was unnecessary. 

**Benefits:** Exact decomposition with correct global phase and single T-gate decomposition for $\pm(2k+1)\frac{\pi}{4}$ case.

**Possible Drawbacks:** None.

**Related GitHub Issues:** [sc-94739, sc-94740]



Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 822c19e41..e590c77e9 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -22,6 +22,9 @@
   measurement-based quantum-computing instructions in the xDSL framework.
   [(#7815)](https://github.com/PennyLaneAI/pennylane/pull/7815)
 
+* The :func:`pennylane.ops.rs_decomposition` method now performs exact decomposition and returns
+  complete global phase information when used for decomposing a phase gate to Clifford+T basis.
+  [(#7793)](https://github.com/PennyLaneAI/pennylane/pull/7793)
 
 <h3>Labs: a place for unified and rapid prototyping of research software 🧪</h3>
 
@@ -54,7 +57,8 @@
 
 This release contains contributions from (in alphabetical order):
 
-Erick Ochoa,
+Utkarsh Azad,
 Joey Carter,
+Erick Ochoa,
 Andrija Paurevic,
 Jay Soni,
diff --git a/pennylane/ops/op_math/decompositions/grid_problems.py b/pennylane/ops/op_math/decompositions/grid_problems.py
index bb5281246..e0d5bd39d 100644
--- a/pennylane/ops/op_math/decompositions/grid_problems.py
+++ b/pennylane/ops/op_math/decompositions/grid_problems.py
@@ -50,6 +50,12 @@ class Ellipse:
     Args:
         D (list[float, float, float]): The elements of the positive definite matrix.
         p (tuple[float, float]): The center of the ellipse.
+
+    .. note::
+      We obtain corresponding matrix :math:`D^{\prime}` with determinant :math:`1` from :math:`D`.
+      These are useful for computing key properties of ``EllipseState``, i.e., a pair of ellipses:
+        - a = eλ^{-z}, d = eλ^z, b^2 = e^2 - 1; (Eq. 31, arXiv:1403.2975)
+        - 2z * log(λ) = log(d / a) => z = 0.5 * log(d / a) / log(λ)
     """
 
     def __init__(
@@ -59,8 +65,6 @@ class Ellipse:
     ):
         self.a, self.b, self.d = D
         self.p = p
-        # a = eλ^{-z}, d = eλ^z, b^2 = e^2 - 1; (Eq. 31, arXiv:1403.2975)
-        # 2z * log(λ) = log(d / a) => z = 0.5 * log(d / a) / log(λ)
         self.z = 0.5 * math.log2(self.d / self.a) / math.log2(_LAMBDA)
         self.e = math.sqrt(self.a * self.d)
 
@@ -554,18 +558,16 @@ def _useful_grid_ops() -> dict["str", "GridOp"]:
 class GridIterator:
     r"""Iterate over the solutions to the scaled grid problem.
 
-    This is based on the Section 5 of arXiv:1403.2975 and
-    implements Proposition 5.22, to enumerate all solutions
-    to the scaled grid problem over the :math:`\epsilon`-region,
+    This is based on the Section 5 of arXiv:1403.2975 and implements Proposition 5.22,
+    to enumerate all solutions to the scaled grid problem over the :math:`\epsilon`-region
     and a unit disk.
 
-    It implements an ``__iter__`` method to iterate over the
-    solutions efficiently as a generator.
+    It implements an ``__iter__`` method to iterate over the solutions efficiently as a generator.
 
     Args:
         theta (float): The angle of the grid problem.
         epsilon (float): The epsilon of the grid problem.
-        max_trials (int): The maximum number of iterations. Default is 20.
+        max_trials (int): The maximum number of iterations. Default is ``20``.
     """
 
     def __init__(self, theta: float = 0.0, epsilon: float = 1e-3, max_trials: int = 20):
diff --git a/pennylane/ops/op_math/decompositions/normal_forms.py b/pennylane/ops/op_math/decompositions/normal_forms.py
index 0e0f7b137..8c868bcd4 100644
--- a/pennylane/ops/op_math/decompositions/normal_forms.py
+++ b/pennylane/ops/op_math/decompositions/normal_forms.py
@@ -15,45 +15,61 @@
 
 from copy import deepcopy
 from functools import lru_cache
+from math import pi as PI
 
 import pennylane as qml
-from pennylane.ops.op_math.decompositions.rings import DyadicMatrix, SO3Matrix, ZOmega
+from pennylane.ops.op_math.decompositions.rings import _SQRT2, DyadicMatrix, SO3Matrix, ZOmega
 
 
 @lru_cache
 def _clifford_group_to_SO3() -> dict:
-    r"""Return a dictionary mapping Clifford group elements to their corresponding SO(3) matrices
-    and global phase scaled by :math:`\pi^{-1}`."""
+    r"""Return a dictionary mapping Clifford group elements to their corresponding SO(3) matrices."""
     I, X, Y, Z = qml.I(0), qml.X(0), qml.Y(0), qml.Z(0)
     H, S, Sd = qml.H(0), qml.S(0), qml.adjoint(qml.S(0))
     # These are the Clifford group elements with :math:`\{−1, 0, 1\}` as their matrix entries.
     clifford_elems = {
-        (I,): (DyadicMatrix(ZOmega(d=1), ZOmega(), ZOmega(), ZOmega(d=1)), 0.0),
-        (H,): (-DyadicMatrix(ZOmega(b=1), ZOmega(b=1), ZOmega(b=1), ZOmega(b=-1), k=1), 0.5),
-        (S,): (DyadicMatrix(ZOmega(a=-1), ZOmega(), ZOmega(), ZOmega(c=1)), 0.25),
-        (X,): (DyadicMatrix(ZOmega(), ZOmega(b=-1), ZOmega(b=-1), ZOmega()), 0.5),
-        (Y,): (DyadicMatrix(ZOmega(), ZOmega(d=-1), ZOmega(d=1), ZOmega()), 0.5),
-        (Z,): (DyadicMatrix(ZOmega(b=-1), ZOmega(), ZOmega(), ZOmega(b=1)), 0.5),
-        (Sd,): (DyadicMatrix(ZOmega(c=-1), ZOmega(), ZOmega(), ZOmega(a=1)), 0.75),
-        (H, S): (DyadicMatrix(ZOmega(c=-1), ZOmega(a=-1), ZOmega(c=-1), ZOmega(a=1), k=1), 0.75),
-        (H, Z): (DyadicMatrix(ZOmega(d=1), ZOmega(d=-1), ZOmega(d=1), ZOmega(d=1), k=1), 0.0),
-        (H, Sd): (DyadicMatrix(ZOmega(a=-1), ZOmega(c=-1), ZOmega(a=-1), ZOmega(c=1), k=1), 0.25),
-        (S, H): (DyadicMatrix(ZOmega(c=-1), ZOmega(c=-1), ZOmega(a=-1), ZOmega(a=1), k=1), 0.75),
-        (S, X): (DyadicMatrix(ZOmega(), ZOmega(c=-1), ZOmega(a=-1), ZOmega()), 0.75),
-        (S, Y): (DyadicMatrix(ZOmega(), ZOmega(a=1), ZOmega(c=1), ZOmega()), 0.75),
-        (Z, H): (DyadicMatrix(ZOmega(d=1), ZOmega(d=1), ZOmega(d=-1), ZOmega(d=1), k=1), 0.0),
-        (Sd, H): (DyadicMatrix(ZOmega(a=-1), ZOmega(a=-1), ZOmega(c=-1), ZOmega(c=1), k=1), 0.25),
-        (S, H, S): (DyadicMatrix(ZOmega(d=1), ZOmega(b=1), ZOmega(b=1), ZOmega(d=1), k=1), 0.0),
-        (S, H, Z): (DyadicMatrix(ZOmega(a=-1), ZOmega(a=1), ZOmega(c=1), ZOmega(c=1), k=1), 0.25),
-        (S, H, Sd): (DyadicMatrix(ZOmega(b=-1), ZOmega(d=-1), ZOmega(d=1), ZOmega(b=1), k=1), 0.5),
-        (Z, H, S): (DyadicMatrix(ZOmega(a=-1), ZOmega(c=1), ZOmega(a=1), ZOmega(c=1), k=1), 0.25),
-        (Z, H, Z): (DyadicMatrix(ZOmega(b=-1), ZOmega(b=1), ZOmega(b=1), ZOmega(b=1), k=1), 0.5),
-        (Z, H, Sd): (DyadicMatrix(ZOmega(c=-1), ZOmega(a=1), ZOmega(c=1), ZOmega(a=1), k=1), 0.75),
-        (Sd, H, S): (DyadicMatrix(ZOmega(b=-1), ZOmega(d=1), ZOmega(d=-1), ZOmega(b=1), k=1), 0.5),
-        (Sd, H, Z): (DyadicMatrix(ZOmega(c=-1), ZOmega(c=1), ZOmega(a=1), ZOmega(a=1), k=1), 0.75),
-        (Sd, H, Sd): (DyadicMatrix(ZOmega(d=1), ZOmega(b=-1), ZOmega(b=-1), ZOmega(d=1), k=1), 0.0),
+        (I,): DyadicMatrix(ZOmega(d=1), ZOmega(), ZOmega(), ZOmega(d=1)),
+        (H,): -DyadicMatrix(ZOmega(b=1), ZOmega(b=1), ZOmega(b=1), ZOmega(b=-1), k=1),
+        (S,): DyadicMatrix(ZOmega(a=-1), ZOmega(), ZOmega(), ZOmega(c=1)),
+        (X,): DyadicMatrix(ZOmega(), ZOmega(b=-1), ZOmega(b=-1), ZOmega()),
+        (Y,): DyadicMatrix(ZOmega(), ZOmega(d=-1), ZOmega(d=1), ZOmega()),
+        (Z,): DyadicMatrix(ZOmega(b=-1), ZOmega(), ZOmega(), ZOmega(b=1)),
+        (Sd,): DyadicMatrix(ZOmega(c=-1), ZOmega(), ZOmega(), ZOmega(a=1)),
+        (H, S): DyadicMatrix(ZOmega(c=-1), ZOmega(a=-1), ZOmega(c=-1), ZOmega(a=1), k=1),
+        (H, Z): DyadicMatrix(ZOmega(d=1), ZOmega(d=-1), ZOmega(d=1), ZOmega(d=1), k=1),
+        (H, Sd): DyadicMatrix(ZOmega(a=-1), ZOmega(c=-1), ZOmega(a=-1), ZOmega(c=1), k=1),
+        (S, H): DyadicMatrix(ZOmega(c=-1), ZOmega(c=-1), ZOmega(a=-1), ZOmega(a=1), k=1),
+        (S, X): DyadicMatrix(ZOmega(), ZOmega(c=-1), ZOmega(a=-1), ZOmega()),
+        (S, Y): DyadicMatrix(ZOmega(), ZOmega(a=1), ZOmega(c=1), ZOmega()),
+        (Z, H): DyadicMatrix(ZOmega(d=1), ZOmega(d=1), ZOmega(d=-1), ZOmega(d=1), k=1),
+        (Sd, H): DyadicMatrix(ZOmega(a=-1), ZOmega(a=-1), ZOmega(c=-1), ZOmega(c=1), k=1),
+        (S, H, S): DyadicMatrix(ZOmega(d=1), ZOmega(b=1), ZOmega(b=1), ZOmega(d=1), k=1),
+        (S, H, Z): DyadicMatrix(ZOmega(a=-1), ZOmega(a=1), ZOmega(c=1), ZOmega(c=1), k=1),
+        (S, H, Sd): DyadicMatrix(ZOmega(b=-1), ZOmega(d=-1), ZOmega(d=1), ZOmega(b=1), k=1),
+        (Z, H, S): DyadicMatrix(ZOmega(a=-1), ZOmega(c=1), ZOmega(a=1), ZOmega(c=1), k=1),
+        (Z, H, Z): DyadicMatrix(ZOmega(b=-1), ZOmega(b=1), ZOmega(b=1), ZOmega(b=1), k=1),
+        (Z, H, Sd): DyadicMatrix(ZOmega(c=-1), ZOmega(a=1), ZOmega(c=1), ZOmega(a=1), k=1),
+        (Sd, H, S): DyadicMatrix(ZOmega(b=-1), ZOmega(d=1), ZOmega(d=-1), ZOmega(b=1), k=1),
+        (Sd, H, Z): DyadicMatrix(ZOmega(c=-1), ZOmega(c=1), ZOmega(a=1), ZOmega(a=1), k=1),
+        (Sd, H, Sd): DyadicMatrix(ZOmega(d=1), ZOmega(b=-1), ZOmega(b=-1), ZOmega(d=1), k=1),
+    }
+    return {gate: SO3Matrix(su2) for gate, su2 in clifford_elems.items()}
+
+
+def _clifford_gates_to_SU2() -> dict:
+    r"""Returns a dictionary mapping single-qubit Clifford group elements to their corresponding SU(2) matrices
+    and global phase scaled by :math:`\pi^{-1}`."""
+    I, X, Y, Z = qml.I(0), qml.X(0), qml.Y(0), qml.Z(0)
+    H, S, Sd = qml.H(0), qml.S(0), qml.adjoint(qml.S(0))
+    return {
+        I: (DyadicMatrix(ZOmega(d=1), ZOmega(), ZOmega(), ZOmega(d=1)), 0.0),
+        H: (-DyadicMatrix(ZOmega(b=1), ZOmega(b=1), ZOmega(b=1), ZOmega(b=-1), k=1), 0.5),
+        S: (DyadicMatrix(ZOmega(a=-1), ZOmega(), ZOmega(), ZOmega(c=1)), 0.25),
+        X: (DyadicMatrix(ZOmega(), ZOmega(b=-1), ZOmega(b=-1), ZOmega()), 0.5),
+        Y: (DyadicMatrix(ZOmega(), ZOmega(d=-1), ZOmega(d=1), ZOmega()), 0.5),
+        Z: (DyadicMatrix(ZOmega(b=-1), ZOmega(), ZOmega(), ZOmega(b=1)), 0.5),
+        Sd: (DyadicMatrix(ZOmega(c=-1), ZOmega(), ZOmega(), ZOmega(a=1)), 0.75),
     }
-    return {gate: (SO3Matrix(su2), phase) for gate, (su2, phase) in clifford_elems.items()}
 
 
 @lru_cache
@@ -65,11 +81,21 @@ def _parity_transforms() -> dict:
             - The SO(3) matrix representation of its inverse.
             - The PennyLane gate representation of the transformation.
             - The global phase scaled by :math:`\pi^{-1}` of the transformation.
+
+    Developer Notes:
+
+        - The following dictionary maps the keys in the Matsumoto-Amano normal form to their:
+            1. SO(3) matrix representation (useful for obtaining the parity vector)
+            2. SO(3) matrix representation of their inverse (useful of reversing the operation)
+            3. PennyLane gate representation (useful for capturing normal form output)
+            4. Global phase difference from PL operators scaled by :math:`\pi^{-1}` (read below).
+
+        - We use the exact matrix representation of the T-gate instead of its SU(2) representation,
+          as representing :math:`\exp(\pm i\pi/8)` in the `ZOmega` ring seems non trivial. We take
+          this into account when computing the global phase in (4), where we compute the global
+          phase by checking the difference between the PL-gates implementation of the keys from
+          the `DyadicMatrix` implementation.
     """
-    # The following dictionary maps the keys in the Matsumoto-Amano normal form to their:
-    # 1. SO(3) matrix representation (useful for obtaining the parity vector)
-    # 2. SO(3) matrix representation of their inverse (useful of reversing the operation)
-    # 3. PennyLane gate representation (useful for capturing normal form output)
     transform_ops = {
         "C": (
             SO3Matrix(DyadicMatrix(ZOmega(d=1), ZOmega(), ZOmega(), ZOmega(d=1))),
@@ -81,19 +107,19 @@ def _parity_transforms() -> dict:
             SO3Matrix(DyadicMatrix(ZOmega(d=1), ZOmega(), ZOmega(), ZOmega(c=1))),
             SO3Matrix(DyadicMatrix(ZOmega(d=1), ZOmega(), ZOmega(), ZOmega(a=-1))),
             (qml.T(0),),
-            0.125,
+            0.0,
         ),
         "HT": (
             SO3Matrix(DyadicMatrix(ZOmega(d=1), ZOmega(c=1), ZOmega(d=1), ZOmega(c=-1), k=1)),
             SO3Matrix(DyadicMatrix(ZOmega(d=1), ZOmega(d=1), ZOmega(a=-1), ZOmega(a=1), k=1)),
             (qml.H(0), qml.T(0)),
-            0.625,
+            1.5,
         ),
         "SHT": (
             SO3Matrix(DyadicMatrix(ZOmega(d=1), ZOmega(c=1), ZOmega(b=1), ZOmega(a=-1), k=1)),
             SO3Matrix(DyadicMatrix(ZOmega(d=1), ZOmega(b=-1), ZOmega(a=-1), ZOmega(c=1), k=1)),
             (qml.S(0), qml.H(0), qml.T(0)),
-            0.875,
+            1.25,
         ),
     }
     return {
@@ -121,27 +147,48 @@ def _ma_normal_form(
         Tuple[qml.operation.Operator, float] | tuple[int, tuple[int, ...], int, float]: The decomposition of the SO(3) matrix into Matsumoto-Amano normal forms and acquired global phase.
     """
     parity_transforms = _parity_transforms()
-    clifford_elements = _clifford_group_to_SO3()
+    clifford_so3s = _clifford_group_to_SO3()
+    clifford_su2s = _clifford_gates_to_SU2()
 
     so3_op = deepcopy(op)
 
     # The following use lemmas from arXiv:1312.6584.
+    a, c, k = ZOmega(d=1), ZOmega(), 0  # Useful for global phase tracking.
     decomposition, rep_bits, g_phase = [], [], 0.0
     while (parity_vec := tuple(so3_op.parity_vec)) != (1, 1, 1):  # Fig. 1 and Lemma 6.3
         so3_val, op_gate, op_phase = parity_transforms[parity_vec]  # Lemma 4.10
         so3_op = so3_val @ so3_op  # Lemma 6.4
-        rep_bits.append(1 if parity_vec == (2, 0, 2) else 0)
         decomposition.extend(op_gate)
         g_phase += op_phase
 
+        if parity_vec == (2, 2, 0):  # T
+            c = c * ZOmega(c=1)
+            rep_bits.append(0)
+        elif parity_vec == (0, 2, 2):  # HT
+            a, c, k = ZOmega(b=-1) * (a + c), ZOmega(a=-1) * (a - c), k + 1
+            rep_bits.append(0)
+        else:  # SHT
+            ic = ZOmega(b=1) * c
+            a, c, k = ZOmega(c=-1) * (a + ic), ZOmega(b=-1) * (a - ic), k + 1
+            rep_bits.append(1)
+
     cl_index = -1
-    for clifford_index, (clifford_gate, clifford_data) in enumerate(clifford_elements.items()):
-        if clifford_data[0] == so3_op:
-            decomposition.extend(clifford_gate)
+    for clifford_index, (clifford_ops, clifford_so3) in enumerate(clifford_so3s.items()):
+        if clifford_so3 == so3_op:
             cl_index = clifford_index
-            g_phase += clifford_data[1]
+            for clf_op in clifford_ops:
+                decomposition.append(clf_op)
+                su2, gp = clifford_su2s[clf_op]
+                a, c, k = su2.a * a + su2.b * c, su2.c * a + su2.d * c, k + su2.k
+                g_phase -= gp
             break
 
+    # Extract the global phase from the decomposition from the
+    # tracked elements (`a` and `c`) of the Dyadic matrix.
+    su2mat = op.matrix
+    g_angle = -qml.math.angle(complex(su2mat.a) / complex(a) * _SQRT2 ** (k - su2mat.k))
+    g_phase = g_angle / PI - g_phase
+
     if not compressed:
         return decomposition, g_phase
 
diff --git a/pennylane/ops/op_math/decompositions/rings.py b/pennylane/ops/op_math/decompositions/rings.py
index 58677c6e0..595dd1d53 100644
--- a/pennylane/ops/op_math/decompositions/rings.py
+++ b/pennylane/ops/op_math/decompositions/rings.py
@@ -74,8 +74,10 @@ class ZSqrtTwo:
             return ZSqrtTwo(self.a * int(other), self.b * int(other))
         raise TypeError(f"Unsupported type {type(other)} for multiplication with ZSqrtTwo")
 
+    def __rsub__(self, other: ZSqrtTwo | int | float) -> ZSqrtTwo:
+        return other + (-self)
+
     __radd__ = __add__
-    __rsub__ = __sub__
     __rmul__ = __mul__
 
     def __pow__(self, power: int) -> ZSqrtTwo:
@@ -238,10 +240,10 @@ class ZOmega:
             return ZOmega(self.a, self.b, self.c, self.d + int(other))
         raise TypeError(f"Unsupported type {type(other)} for addition with ZOmega")
 
-    def __sub__(self, other):
+    def __sub__(self, other: ZOmega | int | float) -> ZOmega:
         return self + (-other)
 
-    def __rsub__(self, other):
+    def __rsub__(self, other: ZOmega | int | float) -> ZOmega:
         return other + (-self)
 
     __radd__ = __add__
diff --git a/pennylane/ops/op_math/decompositions/ross_selinger.py b/pennylane/ops/op_math/decompositions/ross_selinger.py
index ecd2af0b4..99e24b4d5 100644
--- a/pennylane/ops/op_math/decompositions/ross_selinger.py
+++ b/pennylane/ops/op_math/decompositions/ross_selinger.py
@@ -40,6 +40,30 @@ def _domain_correction(theta: float) -> tuple[float, ZOmega]:
         theta -= math.pi * 4
     abs_theta = abs(theta)
 
+    PI_8 = math.pi / 8
+    div_ = round(abs_theta / PI_8, 12)  # Find and round the quotient
+    ivl_ = int(div_)  # Integer part of the quotient
+    mod_ = round(abs(ivl_ * PI_8 - abs_theta), 12)  # Find and round the remainder
+
+    # Check if abs_theta is an odd multiple of PI_8
+    if div_ == ivl_ and mod_ == 0.0 and ivl_ % 2 == 1.0:
+        scale_vals = (
+            [
+                (ZOmega(d=1), ZOmega(d=1)),
+                (ZOmega(b=-1), ZOmega(c=1)),
+                (ZOmega(d=-1), ZOmega(b=1)),
+                (ZOmega(b=1), ZOmega(a=1)),
+            ]
+            if sign == -1
+            else [
+                (ZOmega(b=1), ZOmega(b=1)),
+                (ZOmega(d=-1), ZOmega(d=1)),
+                (ZOmega(b=-1), ZOmega(b=-1)),
+                (ZOmega(d=1), ZOmega(d=-1)),
+            ]
+        )
+        return (sign, ivl_), scale_vals[ivl_ // 2]
+
     if pi_vals[0] <= abs_theta < pi_vals[1]:  # pi/4 <= |theta| < 3pi/4
         return -sign * math.pi / 2, ZOmega(b=sign)
 
@@ -63,10 +87,6 @@ def rs_decomposition(op, epsilon, *, max_search_trials=20, max_factoring_trials=
 
     This algorithm produces a decomposition with :math:`O(3\text{log}_2(1/\epsilon)) + O(\text{log}_2(\text{log}_2(1/\epsilon)))` operations.
 
-    .. note::
-        Currently, the global phase :math:`\theta` returned by the decomposition might differ from the
-        true global phase :math:`\theta^{*}` by an additive factor of :math:`\pi`.
-
     Args:
         op (~pennylane.RZ | ~pennylane.PhaseShift): A :class:`~.RZ` or :class:`~.PhaseShift` gate operation.
         epsilon (float): The maximum permissible error.
@@ -102,7 +122,7 @@ def rs_decomposition(op, epsilon, *, max_search_trials=20, max_factoring_trials=
         matrix_rs = qml.prod(*reversed(ops)).matrix()
 
     When the function is run for a sufficient ``max_search_trials``, the output gate sequence
-    should implement the same operation approximately, up to a global phase.
+    should implement the same operation approximately, up to an :math:`\epsilon`-error.
 
     >>> qml.math.allclose(op.matrix(), matrix_rs, atol=1e-3)
     True
@@ -118,19 +138,26 @@ def rs_decomposition(op, epsilon, *, max_search_trials=20, max_factoring_trials=
         # Get the implemented angle with the domain correction and scaling factor for it.
         angle = -op.data[0] / 2
         shift, scale = _domain_correction(angle)
-
-        # Get the grid problem for the angle.
-        u_solutions = GridIterator(angle + shift, epsilon, max_trials=max_search_trials)
+        interface = qml.math.get_interface(angle)
+        phase = 0.0 if isinstance(op, qml.RZ) else angle
 
         u, t, k = ZOmega(d=1), ZOmega(), 0
-        for u_sol, k_val in u_solutions:
-            xi = ZSqrtTwo(2**k_val) - u_sol.norm().to_sqrt_two()
-            if (t_sol := _solve_diophantine(xi, max_trials=max_factoring_trials)) is not None:
-                u, t, k = u_sol * scale, t_sol * scale, k_val
-                break
+        if not isinstance(scale, ZOmega):  # Get solution for the ± (2k + 1) . PI / 4 case.
+            t_mat = DyadicMatrix(u, t, t, ZOmega(c=1)) @ DyadicMatrix(scale[0], t, t, u)
+            dyd_mat = t_mat * scale[1]
+            phase += math.pi / 8 if shift[0] == -1 else (8 - shift[1]) * math.pi / 8
+
+        else:  # Get the solution from the grid search solver.
+            u_solutions = GridIterator(angle + shift, epsilon, max_trials=max_search_trials)
+            for u_sol, k_val in u_solutions:
+                xi = ZSqrtTwo(2**k_val) - u_sol.norm().to_sqrt_two()
+                if (t_sol := _solve_diophantine(xi, max_trials=max_factoring_trials)) is not None:
+                    u, t, k = u_sol * scale, t_sol * scale, k_val
+                    break
+
+            dyd_mat = DyadicMatrix(u, -t.conj(), t, u.conj(), k=k)
 
         # Get the normal form of the decomposition.
-        dyd_mat = DyadicMatrix(u, -t.conj(), t, u.conj(), k=k)
         so3_mat = SO3Matrix(dyd_mat)
         decomposition, g_phase = _ma_normal_form(so3_mat)
 
@@ -148,9 +175,6 @@ def rs_decomposition(op, epsilon, *, max_search_trials=20, max_factoring_trials=
             _ = [qml.apply(op) for op in new_tape.operations]
 
     # TODO: Improve the global phase information to the decomposition.
-    interface = qml.math.get_interface(angle)
-    phase = 0.0 if isinstance(op, qml.RZ) else angle
-
     phase += qml.math.mod(g_phase, 2) * math.pi
     global_phase = qml.GlobalPhase(qml.math.array(phase, like=interface))
 


Test patch:
diff --git a/tests/ops/op_math/decompositions/test_normal_forms.py b/tests/ops/op_math/decompositions/test_normal_forms.py
index d4e4e94d0..d8796d954 100644
--- a/tests/ops/op_math/decompositions/test_normal_forms.py
+++ b/tests/ops/op_math/decompositions/test_normal_forms.py
@@ -19,6 +19,7 @@ import pytest
 
 import pennylane as qml
 from pennylane.ops.op_math.decompositions.normal_forms import (
+    _clifford_gates_to_SU2,
     _clifford_group_to_SO3,
     _ma_normal_form,
     _parity_transforms,
@@ -39,11 +40,10 @@ class TestNormalForms:
         assert isinstance(clifford_matrices, dict)
         assert len(clifford_matrices) == 24
 
-        for gates, (so3mat, phase) in clifford_matrices.items():
+        for gates, so3mat in clifford_matrices.items():
             # Check that the SO(3) matrix is orthogonal and has determinant 1
             assert isinstance(so3mat, SO3Matrix), "All gates should be SO3Matrix instances"
             assert so3mat.k == 0, f"Gate {gates} should have k=0"
-            assert isinstance(phase, float), "Phase should be a float."
 
             su2mat = _SU2_transform(qml.matrix(qml.prod(*gates)))[0]
             w, x, y, z = _quaternion_transform(su2mat)
@@ -58,6 +58,20 @@ class TestNormalForms:
                 ),
             ), "SO(3) matrix does not match expected form"
 
+    def test_clifford_gates_to_SU2(self):
+        """Test that the Clifford gates are correctly mapped to SU(2) matrices."""
+        clifford_matrices = _clifford_gates_to_SU2()
+        assert isinstance(clifford_matrices, dict)
+        assert len(clifford_matrices) == 7
+
+        for gate, (su2mat, phase) in clifford_matrices.items():
+            assert isinstance(su2mat, DyadicMatrix), "All gates should be DyadicMatrix instances"
+            assert isinstance(phase, float), "Phase should be a float."
+
+            assert qml.math.allclose(
+                su2mat.ndarray, _SU2_transform(qml.matrix(gate))[0]
+            ), "SU(2) matrix does not match expected form"
+
     def test_parity_transforms(self):
         """Test that the parity transforms are correctly defined."""
         parity_transforms = _parity_transforms()
@@ -72,7 +86,10 @@ class TestNormalForms:
                 gate, (qml.ops.Operation, qml.ops.op_math.Prod)
             ), "Each transform should have a gate"
             assert isinstance(phase, float)
-            assert qml.math.allclose(phase * math.pi, _SU2_transform(qml.matrix(gate))[1])
+            if parity_vec != (1, 1, 1):  # 0.125 here comes from the matrix form of the T-gate.
+                assert qml.math.allclose(
+                    phase, (0.125 - _SU2_transform(qml.matrix(gate))[1] / math.pi) % 2
+                )
 
             su2mat = _SU2_transform(qml.matrix(qml.adjoint(gate)))[0]
             w, x, y, z = _quaternion_transform(su2mat)
@@ -118,7 +135,7 @@ class TestNormalForms:
                     so3rep @= op_
 
         cl_list = list(clifford_elements.keys())
-        so3mat @= clifford_elements[cl_list[c]][0]
+        so3mat @= clifford_elements[cl_list[c]]
         for op_ in cl_list[c]:
             so3rep @= op_
 
diff --git a/tests/ops/op_math/decompositions/test_rings.py b/tests/ops/op_math/decompositions/test_rings.py
index 5e0f990a2..e22ba4190 100644
--- a/tests/ops/op_math/decompositions/test_rings.py
+++ b/tests/ops/op_math/decompositions/test_rings.py
@@ -68,6 +68,7 @@ class TestZSqrtTwo:
         assert z2 % 2 == ZSqrtTwo(1, 0)
         assert z1.to_omega() == ZOmega(a=-2, b=0, c=2, d=1)
         assert z2.to_omega() == ZOmega(a=-4, b=0, c=4, d=3)
+        assert 1 - ZSqrtTwo(1, 2) == ZSqrtTwo(0, -2)
 
     def test_arithmetic_errors(self):
         """Test that arithmetic operations raise errors for invalid types."""
diff --git a/tests/ops/op_math/decompositions/test_ross_selinger.py b/tests/ops/op_math/decompositions/test_ross_selinger.py
index 6fc30363f..82e365167 100644
--- a/tests/ops/op_math/decompositions/test_ross_selinger.py
+++ b/tests/ops/op_math/decompositions/test_ross_selinger.py
@@ -26,7 +26,8 @@ from pennylane.ops.op_math.decompositions.ross_selinger import _domain_correctio
     ("angle", "result"),
     [
         (math.pi / 3, (-math.pi / 2, ZOmega(b=1))),
-        (math.pi / 8, (0.0, ZOmega(d=1))),
+        (math.pi / 8, ((1.0, 1), (ZOmega(b=1), ZOmega(b=1)))),
+        (-3 * math.pi / 8, ((-1.0, 3), (ZOmega(b=-1), ZOmega(c=1)))),
         (math.pi / 2, (-math.pi / 2, ZOmega(b=1))),
         (5 * math.pi / 3, (-3 * math.pi / 2, ZOmega(b=-1))),
         (3 * math.pi / 2, (-3 * math.pi / 2, ZOmega(b=-1))),
@@ -36,7 +37,7 @@ from pennylane.ops.op_math.decompositions.ross_selinger import _domain_correctio
 def test_domain_correction(angle, result):
     """Test the functionality to create domain correction"""
     shift, scale = _domain_correction(angle)
-    assert qml.math.isclose(shift, result[0])
+    assert qml.math.allclose(shift, result[0])
     assert scale == result[1]
 
 
@@ -44,10 +45,14 @@ def test_domain_correction(angle, result):
     ("op", "epsilon"),
     [
         (qml.RZ(math.pi / 42, wires=[1]), 1e-4),
+        (qml.RZ(5 * math.pi / 4, wires=[1]), 1e-5),
         (qml.PhaseShift(math.pi / 7, wires=["a"]), 1e-3),
         (qml.RZ(math.pi / 3, wires=[1]), 1e-5),
+        (qml.RZ(-math.pi / 3, wires=[1]), 1e-3),
+        (qml.RZ(-math.pi / 4, wires=[1]), 1e-4),
         (qml.PhaseShift(-math.pi / 6, wires=[0]), 1e-3),
         (qml.RZ(-math.pi / 5, wires=[0]), 1e-4),
+        (qml.RZ(-math.pi / 7, wires=[1]), 1e-4),
         (qml.RZ(-math.pi / 8, wires=[2]), 1e-2),
     ],
 )
@@ -59,14 +64,7 @@ def test_ross_selinger(op, epsilon):
 
     matrix_rs = qml.matrix(qml.tape.QuantumScript(gates))
 
-    mat1 = qml.matrix(op)
-    mat2 = matrix_rs
-    phase = qml.math.divide(
-        mat1, mat2, out=qml.math.zeros_like(mat1, dtype=complex), where=mat1 != 0
-    )[qml.math.nonzero(qml.math.round(mat1, 10))]
-
-    assert qml.math.allclose(phase / phase[0], qml.math.ones(len(phase)), atol=1e-3)
-    assert qml.math.allclose(qml.matrix(op), phase[0] * matrix_rs, atol=epsilon)
+    assert qml.math.allclose(qml.matrix(op), matrix_rs, atol=epsilon)
     assert qml.prod(*gates, lazy=False).wires == op.wires