Repository: PennyLaneAI/pennylane
Pull request #7691 (2025-06-16 16:14:33)

MCM-TT documentation update 
**Context:**

Now, Lightning devices have support for MCM with the Tree-Traversal algorithm. 

**Description of the Change:**
Update the documentation about the [MCM-TT](https://docs.pennylane.ai/en/stable/introduction/dynamic_quantum_circuits.html#tree-traversal-algorithm) 

**Benefits:**
Rendered version: https://xanaduai-pennylane--7691.com.readthedocs.build/en/7691/introduction/dynamic_quantum_circuits.html

**Possible Drawbacks:**

**Related GitHub Issues:**


[sc-93412]


Patch:
diff --git a/doc/introduction/dynamic_quantum_circuits.rst b/doc/introduction/dynamic_quantum_circuits.rst
index 12398eb2e..5f69c541b 100644
--- a/doc/introduction/dynamic_quantum_circuits.rst
+++ b/doc/introduction/dynamic_quantum_circuits.rst
@@ -266,7 +266,7 @@ scalings  with respect to the number of mid-circuit measurements (and shots) are
 +--------------------------+-------------------------------------------+-----------------------------------------------------------+-------------------------------------------+--------------+--------------+
 | Dynamic one-shot         | :gr:`\ ` :math:`\mathcal{O}(1)`           | :rd:`\ ` :math:`\mathcal{O}(n_{shots})`                   | :or:`\ ` finite differences\ :math:`{}^2` | :gr:`\ ` yes | :rd:`\ ` no  |
 +--------------------------+-------------------------------------------+-----------------------------------------------------------+-------------------------------------------+--------------+--------------+
-| Tree-traversal           | :or:`\ ` :math:`\mathcal{O}(n_{MCM}+1)`   | :or:`\ ` :math:`\mathcal{O}(min(n_{shots}, 2^{n_{MCM}}))` | :or:`\ ` finite differences\ :math:`{}^2` | :gr:`\ ` yes | :rd:`\ ` no  |
+| Tree-traversal           | :or:`\ ` :math:`\mathcal{O}(n_{MCM}+1)`   | :or:`\ ` :math:`\mathcal{O}(min(n_{shots}, 2^{n_{MCM}}))` | :or:`\ ` finite differences\ :math:`{}^2` | :gr:`\ ` yes | :gr:`\ ` yes |
 +--------------------------+-------------------------------------------+-----------------------------------------------------------+-------------------------------------------+--------------+--------------+
 
 
@@ -393,11 +393,19 @@ Since the counts of many nodes come out to be zero for shot-based simulations,
 it is often possible to ignore entire sub-trees, thereby reducing the computational
 cost.
 
-.. warning::
+.. note::
 
-    The tree-traversal algorithm is only supported by the
-    :class:`~.pennylane.devices.DefaultQubit` device, and currently does
-    not support just-in-time (JIT) compilation.
+    The tree-traversal algorithm is supported by the following devices:
+
+    * :class:`~.pennylane.devices.DefaultQubit`,
+    
+    * `lightning.qubit <https://docs.pennylane.ai/projects/lightning/en/stable/lightning_qubit/device.html>`_,
+    
+    * `lightning.gpu <https://docs.pennylane.ai/projects/lightning/en/stable/lightning_gpu/device.html>`_,
+    
+    * `lightning.kokkos <https://docs.pennylane.ai/projects/lightning/en/stable/lightning_kokkos/device.html>`_,
+    
+    and currently does not support just-in-time (JIT) compilation.
 
 .. _mcm_config:
 
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index fced6b309..7f16965ab 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -410,6 +410,45 @@
   See the documentation for more details.
   [(#7531)](https://github.com/PennyLaneAI/pennylane/pull/7531)
 
+* A new decomposition method for :func:`~.clifford_t_decomposition` is now available with `method="rs"`
+  (the [Ross-Selinger algorithm](https://arxiv.org/abs/1403.2975)) that produces orders of magnitude
+  less gates than `method="sk"` (the Solovay-Kitaev algorithm) in many cases. It is directly accessible
+  via :func:`~.ops.rs_decomposition` function.
+  [(#7588)](https://github.com/PennyLaneAI/pennylane/pull/7588)
+  [(#7641)](https://github.com/PennyLaneAI/pennylane/pull/7641)
+  [(#7611)](https://github.com/PennyLaneAI/pennylane/pull/7611)
+  [(#7711)](https://github.com/PennyLaneAI/pennylane/pull/7711)
+
+  The Ross-Selinger algorithm can drastically outperform the Solovay-Kitaev algorithm in many cases.
+  Consider this simple circuit:
+
+  ```python
+  @qml.qnode(qml.device("lightning.qubit", wires=2))
+  def circuit(x, y):
+
+      qml.RX(x, 0)
+      qml.CNOT([0, 1])
+      qml.RY(y, 0)
+
+      return qml.expval(qml.Z(0))
+
+  rs_circuit = qml.clifford_t_decomposition(circuit, method="rs")
+  sk_circuit = qml.clifford_t_decomposition(circuit, method="sk")
+
+  rs_specs = qml.specs(rs_circuit)(x, y)["resources"]
+  sk_specs = qml.specs(sk_circuit)(x, y)["resources"]
+  ```
+
+  Decomposing with `method="rs"` instead of `method="sk"` gives a significant reduction in overall 
+  gate counts, specifically the `qml.T` count:
+
+  ```pycon
+  >>> print(rs_specs.num_gates, sk_specs.num_gates)
+  267 48637
+  >>> print(rs_specs.gate_types['T'], sk_specs.gate_types['T'])
+  104 8507
+  ```
+
 <h3>Improvements 🛠</h3>
 
 * Adds a new `allocation` module containing `allocate` and `deallocate` instructions for requesting dynamic wires. This is currently
@@ -952,6 +991,11 @@ Here's a list of deprecations made this release. For a more detailed breakdown o
   off from the original tape by a global phase.
   [(#7619)](https://github.com/PennyLaneAI/pennylane/pull/7619)
 
+* Updated documentation for mid-circuit measurements using the Tree Traversal algorithm
+  to reflect supported devices and usage in analytic simulations,
+  in the :doc:`/introduction/dynamic_quantum_circuits` page.
+  [(#7691)](https://github.com/PennyLaneAI/pennylane/pull/7691)
+
 <h3>Contributors ✍️</h3>
 
 This release contains contributions from (in alphabetical order):
@@ -967,6 +1011,7 @@ Simone Gasperini,
 Korbinian Kottmann,
 Christina Lee,
 Anton Naim Ibrahim,
+Luis Alfredo Nuñez Meneses
 Oumarou Oumarou,
 Lee J. O'Riordan,
 Mudit Pandey,
diff --git a/pennylane/ops/op_math/__init__.py b/pennylane/ops/op_math/__init__.py
index d11a9c113..98c1d4b9d 100644
--- a/pennylane/ops/op_math/__init__.py
+++ b/pennylane/ops/op_math/__init__.py
@@ -89,6 +89,7 @@ Decompositions
     ~two_qubit_decomposition
     ~multi_qubit_decomposition
     ~sk_decomposition
+    ~rs_decomposition
 
 Control Decompositions
 ~~~~~~~~~~~~~~~~~~~~~~
@@ -127,6 +128,7 @@ from .controlled_ops import (
 from .decompositions import (
     one_qubit_decomposition,
     sk_decomposition,
+    rs_decomposition,
     two_qubit_decomposition,
     multi_qubit_decomposition,
     ctrl_decomp_bisect,
diff --git a/pennylane/ops/op_math/decompositions/__init__.py b/pennylane/ops/op_math/decompositions/__init__.py
index 4bd2689b8..2f3167af7 100644
--- a/pennylane/ops/op_math/decompositions/__init__.py
+++ b/pennylane/ops/op_math/decompositions/__init__.py
@@ -21,4 +21,5 @@ from .unitary_decompositions import (
     multi_qubit_decomposition,
 )
 from .solovay_kitaev import sk_decomposition
+from .ross_selinger import rs_decomposition
 from .controlled_decompositions import ctrl_decomp_bisect, ctrl_decomp_zyz
diff --git a/pennylane/ops/op_math/decompositions/grid_problems.py b/pennylane/ops/op_math/decompositions/grid_problems.py
new file mode 100644
index 000000000..614bf9942
--- /dev/null
+++ b/pennylane/ops/op_math/decompositions/grid_problems.py
@@ -0,0 +1,881 @@
+# 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.
+"""Contains 1-D and 2-D grid problem solving utilities."""
+
+from __future__ import annotations
+
+import math
+from copy import copy
+from functools import lru_cache
+from itertools import chain
+from typing import Iterable
+
+from pennylane.ops.op_math.decompositions.rings import _SQRT2, ZOmega, ZSqrtTwo
+
+_LAMBDA = 1 + math.sqrt(2)
+
+
+class Ellipse:
+    r"""A class representing an ellipse as a positive definite matrix.
+
+    The ellipse is defined by the positive definite matrix :math:`D`:
+
+    .. math::
+        D = \begin{bmatrix}
+            a & b \\
+            b & d
+        \end{bmatrix}
+
+    The matrix is related to the equation of the ellipse :math:`ax^2 + 2bxy + dy^2 = 1` by:
+
+    .. math::
+        a = \frac{\sin^2(\theta)}{m^2} + \frac{\cos^2(\theta)}{n^2}
+        b = \cos(\theta) \sin(\theta) \left(\frac{1}{m^2} - \frac{1}{n^2}\right)
+        d = \frac{\sin^2(\theta)}{m^2} + \frac{\cos^2(\theta)}{n^2}
+
+    where :math:`m` and :math:`n` are the lengths of the semi-major and semi-minor axes,
+    and :math:`\theta` is the angle of the ellipse.
+
+    Args:
+        D (list[float, float, float]): The elements of the positive definite matrix.
+        p (tuple[float, float]): The center of the ellipse.
+        axes (tuple[float, float]): The lengths of the semi-major and semi-minor axes.
+    """
+
+    def __init__(
+        self,
+        D: tuple[float, float, float],
+        p: tuple[float, float] = (0, 0),
+        axes: tuple[float, float] = (0, 0),
+    ):
+        self.a, self.b, self.d = D
+        self.p = p
+        self.axes = axes
+        # a = eλ^{-z}, d = eλ^z, b^2 = e^2 - 1; (Eq. 31, arXiv:1403.2975)
+        self.z = 0.5 * math.log2(self.d / self.a) / math.log2(_LAMBDA)
+        self.e = math.sqrt(self.a * self.d)
+
+    @classmethod
+    def from_axes(cls, p: tuple[float, float], theta: float, axes: tuple[float, float]):
+        """Create an ellipse from its axes and center point."""
+        axes = tuple(ax + 1e-20 for ax in axes)  # Avoid division by zero.
+        a = (math.sin(theta) / axes[1]) ** 2 + (math.cos(theta) / axes[0]) ** 2
+        b = math.cos(theta) * math.sin(theta) * (1 / (axes[0]) ** 2 - 1 / (axes[1]) ** 2)
+        c = (math.sin(theta) / axes[0]) ** 2 + (math.cos(theta) / axes[1]) ** 2
+        return cls([a, b, c], p, axes)
+
+    @classmethod
+    def from_region(cls, theta: float, epsilon: float, k: int = 0):
+        r"""Create an ellipse that bounds the region :math:`u | u \bullet z \geq 1 - \epsilon^2 / 2`,
+        with :math:`u \in \frac{1}{\sqrt{2}^k} \mathbb{Z}[\omega]` and :math`z = \exp{-i\theta / 2}`.
+        """
+        const = 1 - epsilon**2 / 2
+        shift = (1 - const) / 3
+        scale = (2 ** (k // 2)) * (_SQRT2 ** (k % 2))
+        semi_major = 2 * scale * shift
+        semi_minor = 2 * scale * math.sqrt((const + 0.5 * epsilon**2) * epsilon**2 / 3)
+        x, y = scale * (const + shift) * math.cos(theta), scale * (const + shift) * math.sin(theta)
+        return cls.from_axes((x, y), theta, (semi_major, semi_minor))
+
+    def __repr__(self) -> str:
+        """Return a string representation of the ellipse."""
+        return f"Ellipse(a={self.a}, b={self.b}, d={self.d}, p={self.p}, axes={self.axes})"
+
+    def __eq__(self, other: "Ellipse") -> bool:
+        """Check if the ellipses are equal."""
+        return (
+            self.a == other.a
+            and self.b == other.b
+            and self.d == other.d
+            and self.p == other.p
+            and self.axes == other.axes
+        )
+
+    @property
+    def discriminant(self) -> float:
+        """Calculate the discriminant of the characteristic polynomial associated with the ellipse."""
+        return (self.a + self.d) ** 2 - 4 * (self.a * self.d - self.b**2)
+
+    @property
+    def determinant(self) -> float:
+        """Calculate the determinant of the ellipse."""
+        return self.a * self.d - self.b**2
+
+    @property
+    def positive_semi_definite(self) -> bool:
+        """Check if the ellipse is positive semi-definite."""
+        return (self.a + self.d) + math.sqrt(self.discriminant) >= 0
+
+    @property
+    def uprightness(self) -> float:
+        """Calculate the uprightness of the ellipse (Eq. 32, arXiv:1403.2975)."""
+        return math.pi / (2 * self.e) ** 2
+
+    @staticmethod
+    def b_from_uprightness(uprightness: float) -> float:
+        """Calculate the b value of the ellipse from its uprightness (Eq. 33, arXiv:1403.2975)."""
+        return math.sqrt((math.pi / (4 * uprightness)) ** 2 - 1)
+
+    def contains(self, x: float, y: float) -> bool:
+        """Check if the point (x, y) is inside the ellipse."""
+        x_, y_ = x - self.p[0], y - self.p[1]
+        return (self.a * x_**2 + 2 * self.b * x_ * y_ + self.d * y_**2) <= 1
+
+    def normalize(self) -> tuple["Ellipse", float]:
+        """Normalize the ellipse to have a determinant of 1."""
+        s_val = 1 / math.sqrt(self.determinant)
+        return self.scale(scale=s_val), s_val
+
+    def scale(self, scale: float) -> "Ellipse":
+        """Scale the ellipse by a factor of scale."""
+        D = (self.a * scale, self.b * scale, self.d * scale)
+        axes = (self.axes[0] * math.sqrt(scale), self.axes[1] * math.sqrt(scale))
+        return Ellipse(D, self.p, axes)
+
+    def x_points(self, y: float) -> tuple[float, float]:
+        """Compute the x-points of the ellipse for a given y-value."""
+        y = float(y) - self.p[1]  # shift y to the origin
+        discriminant = y**2 * (self.b**2 - self.a * self.d) + self.a
+        if discriminant < 0:
+            raise ValueError(f"Point y={y} is outside the ellipse")
+
+        x0, d0 = self.p[0], math.sqrt(discriminant)
+        x1 = (-self.b * y - d0) / self.a
+        x2 = (-self.b * y + d0) / self.a
+        return x0 + x1, x0 + x2
+
+    def y_points(self, x: float) -> tuple[float, float] | None:
+        """Compute the y-points of the ellipse for a given x-value."""
+        x = float(x) - self.p[0]  # shift x to the origin
+        discriminant = (self.b * x) ** 2 - self.d * (self.a * x**2 - 1)
+        if discriminant < 0:
+            raise ValueError(f"Point x={x} is outside the ellipse")
+
+        y0, d0 = self.p[1], math.sqrt(discriminant)
+        y1 = (-self.b * x - d0) / self.d
+        y2 = (-self.b * x + d0) / self.d
+        return y0 + y1, y0 + y2
+
+    def bounding_box(self) -> tuple[float, float, float, float]:
+        """Return the bounding box of the ellipse in the form [[x0, x1], [y0, y1]]."""
+        denom = self.determinant
+        x, y = math.sqrt((self.d / denom)), math.sqrt((self.a / denom))
+        return (-x, x, -y, y)
+
+    def offset(self, offset: float) -> "Ellipse":
+        """Return the ellipse shifted by the offset."""
+        p_offset = (self.p[0] + offset, self.p[1] + offset)
+        return Ellipse((self.a, self.b, self.d), p_offset, self.axes)
+
+    def apply_grid_op(self, grid_op: GridOp) -> "Ellipse":
+        """Apply a grid operation :math:`G` to the ellipse :math:`E` as :math:`G^T E G`."""
+        ga, gb, gc, gd = grid_op.flatten
+
+        D = (
+            ga**2 * self.a + 2 * ga * gc * self.b + self.d * gc**2,
+            ga * gb * self.a + (ga * gd + gb * gc) * self.b + gc * gd * self.d,
+            gb**2 * self.a + 2 * gb * gd * self.b + self.d * gd**2,
+        )
+
+        p1, p2 = self.p
+        gda, gdb, gdc, gdd = grid_op.inverse().flatten
+
+        return Ellipse(D, (gda * p1 + gdb * p2, gdc * p1 + gdd * p2), self.axes)
+
+
+class EllipseState:
+    """A class representing a state as a pair of normalized ellipses.
+
+    This is based on the Definition A.1 of arXiv:1403.2975,
+    where the pair of ellipses are represented by real symmetric
+    positive semi-definite matrices of determinant :math:`1`.
+
+    Args:
+        e1 (Ellipse): The first ellipse.
+        e2 (Ellipse): The second ellipse.
+    """
+
+    def __init__(self, e1: Ellipse, e2: Ellipse):
+        self.e1 = e1
+        self.e2 = e2
+
+    def __repr__(self) -> str:
+        """Return a string representation of the state."""
+        return f"EllipseState(e1={self.e1}, e2={self.e2})"
+
+    @property
+    def skew(self) -> float:
+        """Calculate the skew of the state."""
+        # Uses Definition A.1 of arXiv:1403.2975
+        return self.e1.b**2 + self.e2.b**2
+
+    @property
+    def bias(self) -> float:
+        """Calculate the bias of the state."""
+        # Uses Definition A.1 of arXiv:1403.2975
+        return self.e2.z - self.e1.z
+
+    def skew_grid_op(self) -> GridOp:
+        """Calculate the special grid operation for the state for reducing the skew."""
+        # Uses Lemma A.5 (Step Lemma) of arXiv:1403.2975 for obtaining special grid op.
+        grid_op = GridOp.from_string("I")
+        state = EllipseState(self.e1, self.e2)
+        while (skew := state.skew) >= 15:
+            new_grid_op, state = state.reduce_skew()
+            grid_op = grid_op * new_grid_op
+            if state.skew > 0.9 * skew:  # pragma: no cover
+                raise ValueError(f"Skew was not decreased for state {state}")
+
+        return grid_op
+
+    def apply_grid_op(self, grid_op: GridOp) -> EllipseState:
+        """Apply a grid operation :math:`G` to the state."""
+        # Uses Definition A.3 of arXiv:1403.2975
+        return EllipseState(self.e1.apply_grid_op(grid_op), self.e2.apply_grid_op(grid_op.adj2()))
+
+    def apply_shift_op(self, k: int) -> tuple[EllipseState, int]:
+        """Apply a shift operator to the state."""
+        # Uses Definition A.6 and Lemma A.8 of arXiv:1403.2975
+        k = int(math.floor((1 - self.bias) / 2))
+        pk_pow, nk_pow = _LAMBDA**k, _LAMBDA**-k
+        e1, e2 = copy(self.e1), copy(self.e2)
+        e1.a, e1.d, e1.z = e1.a * pk_pow, e1.d * nk_pow, e1.z - k
+        e2.a, e2.d, e2.z = e2.a * nk_pow, e2.d * pk_pow, e2.z + k
+        e2.b *= (-1) ** k
+        e1.e, e2.e = math.sqrt(e1.a * e1.d), math.sqrt(e2.a * e2.d)
+        return EllipseState(e1, e2), k
+
+    # pylint: disable=too-many-branches
+    def reduce_skew(self) -> tuple[GridOp, EllipseState]:
+        """Reduce the skew of the state.
+
+        This uses Step Lemma described in Appendix A.6 of arXiv:1403.2975.
+
+        Returns:
+            tuple[GridOp, EllipseState]: A tuple containing the grid operation
+                and the state with reduced skew.
+        """
+        if any(not e.positive_semi_definite for e in (self.e1, self.e2)):  # pragma: no cover
+            raise ValueError("Ellipse is not positive semi-definite")
+
+        state, k = copy(self), 0
+        grid_op = GridOp.from_string("I")
+
+        if abs(state.bias) > 1:
+            state, k = state.apply_shift_op(k)
+
+        if state.e2.b < 0:
+            grid_op = grid_op * GridOp.from_string("Z")
+
+        if (state.e1.z + state.e2.z) < 0:
+            grid_op = grid_op * GridOp.from_string("X")
+
+        new_state = state.apply_grid_op(grid_op)
+        e1, e2 = new_state.e1, new_state.e2
+
+        if e1.b >= 0:
+            if -0.8 <= e1.z <= 0.8 and -0.8 <= e2.z <= 0.8:
+                grid_op = grid_op * GridOp.from_string("R")
+            elif e1.z <= 0.3 and e2.z >= 0.8:
+                grid_op = grid_op * GridOp.from_string("K")
+            elif e1.z >= 0.3 and e2.z >= 0.3:
+                c = min(e1.z, e2.z)
+                n = int(max(1, math.floor((_LAMBDA**c) / 2)))
+                grid_op = grid_op * (GridOp.from_string("A") ** n)
+            elif e1.z >= 0.8 and e2.z <= 0.3:
+                grid_op = grid_op * GridOp.from_string("K").adj2()
+            else:  # pragma: no cover
+                raise ValueError(f"Skew couldn't be reduced for the state {state}")
+        else:
+            if e1.z >= -0.8 and e1.z <= 0.8 and e2.z >= -0.8 and e2.z <= 0.8:
+                grid_op = grid_op * GridOp.from_string("R")
+            elif e1.z >= -0.2 and e2.z >= -0.2:
+                c = min(e1.z, e2.z)
+                n = int(max(1, math.floor((_LAMBDA**c) / 2)))
+                grid_op = grid_op * (GridOp.from_string("B") ** n)
+            else:  # pragma: no cover
+                raise ValueError(f"Skew couldn't be reduced for the state {state}")
+
+        if k != 0:
+            grid_op = grid_op.apply_shift_op(k)
+
+        return grid_op, self.apply_grid_op(grid_op)
+
+
+class GridOp:
+    r"""A class representing a grid operation on a 2D grid.
+
+    This follows Definition 5.10 and Lemma 5.11 of arXiv:1403.2975,
+    where a grid operator :math:`G` is described as a linear map,
+    which can be identifed with a :math:`2 \times 2` matrix as follows:
+
+    .. math::
+        G = \begin{pmatrix} a & b \\ c & d \end{pmatrix}.
+
+    Each entry :math:`m` of the matrix is of the form :math:`m = m_0 + m_1 / \sqrt{2}`,
+    where :math:`m_0, m_1 \in \mathbb{Z}`. They satisfy :math:`a_0+b_0+c_0+d_0 \equiv 0 (\mod 2)`
+    and :math:`a_1 \equiv b_1 \equiv c_1 \equiv d_1 (\mod 2)`.
+
+    Args:
+        a (tuple[int, int]): The a-coefficient of the grid operation.
+        b (tuple[int, int]): The b-coefficient of the grid operation.
+        c (tuple[int, int]): The c-coefficient of the grid operation.
+        d (tuple[int, int]): The d-coefficient of the grid operation.
+        check_valid (bool): If ``True``, the grid operation will be checked to be a valid
+            grid operation. Default is ``True``.
+
+    .. note::
+        The coefficients are given as tuples of the form (a_0, a_1), which corresponds
+        to an element being :math:`a = a_0 + a_1 / \sqrt{2}`.
+    """
+
+    # pylint: disable=too-many-arguments
+    def __init__(
+        self,
+        a: tuple[int, int],
+        b: tuple[int, int],
+        c: tuple[int, int],
+        d: tuple[int, int],
+        check_valid: bool = True,
+    ) -> None:
+        """Initialize the grid operation."""
+        self.a = a
+        self.b = b
+        self.c = c
+        self.d = d
+        if check_valid:
+            assert (a[0] + b[0] + c[0] + d[0]) % 2 == 0, "sum of a_0, b_0, c_0, d_0 must be even"
+            assert (
+                a[1] % 2 == b[1] % 2 == c[1] % 2 == d[1] % 2
+            ), "a_1, b_1, c_1, d_1 must have same parity"
+
+    def __repr__(self) -> str:
+        """Return a string representation of the grid operation."""
+        return f"GridOp(a={self.a}, b={self.b}, c={self.c}, d={self.d})"
+
+    def __eq__(self, other: "GridOp") -> bool:
+        """Check if the grid operations are equal."""
+        return self.a == other.a and self.b == other.b and self.c == other.c and self.d == other.d
+
+    def __mul__(self, other: "GridOp" | "ZOmega") -> "GridOp" | "ZOmega":
+        if isinstance(other, GridOp):
+            return self._mul_grid_op(other)
+        if isinstance(other, ZOmega):
+            return self._mul_z_omega(other)
+        raise TypeError(f"Cannot multiply GridOp with {type(other)}")
+
+    def __pow__(self, n: int) -> "GridOp":
+        """Raise the grid operator to a power."""
+        if self == self.from_string("I"):  # Identity:
+            return self
+        if self == self.from_string("A"):  # A
+            return GridOp((1, 0), (-2 * n, 0), (0, 0), (1, 0))
+        if self == self.from_string("B"):  # B
+            return GridOp((1, 0), (0, 2 * n), (0, 0), (1, 0))
+
+        x, res = self, GridOp((1, 0), (0, 0), (0, 0), (1, 0))
+        while n > 0:
+            if n % 2 == 1:
+                res = res * x
+            x = x * x
+            n //= 2
+        return res
+
+    def _mul_grid_op(self, other: "GridOp") -> "GridOp":
+        """Multiply the grid operator with another grid operator."""
+        a = (
+            self.a[0] * other.a[0]
+            + self.b[0] * other.c[0]
+            + (self.a[1] * other.a[1] + self.b[1] * other.c[1]) // 2,
+            self.a[0] * other.a[1]
+            + self.a[1] * other.a[0]
+            + self.b[0] * other.c[1]
+            + self.b[1] * other.c[0],
+        )
+        b = (
+            self.a[0] * other.b[0]
+            + self.b[0] * other.d[0]
+            + (self.a[1] * other.b[1] + self.b[1] * other.d[1]) // 2,
+            self.a[0] * other.b[1]
+            + self.a[1] * other.b[0]
+            + self.b[0] * other.d[1]
+            + self.b[1] * other.d[0],
+        )
+        c = (
+            self.c[0] * other.a[0]
+            + self.d[0] * other.c[0]
+            + (self.c[1] * other.a[1] + self.d[1] * other.c[1]) // 2,
+            self.c[0] * other.a[1]
+            + self.c[1] * other.a[0]
+            + self.d[0] * other.c[1]
+            + self.d[1] * other.c[0],
+        )
+        d = (
+            self.c[0] * other.b[0]
+            + self.d[0] * other.d[0]
+            + (self.c[1] * other.b[1] + self.d[1] * other.d[1]) // 2,
+            self.c[0] * other.b[1]
+            + self.c[1] * other.b[0]
+            + self.d[0] * other.d[1]
+            + self.d[1] * other.d[0],
+        )
+        return GridOp(a, b, c, d)
+
+    def _mul_z_omega(self, other: "ZOmega") -> "ZOmega":
+        """Multiply the grid operator with a ZOmega element."""
+        x1, y1 = other.d, other.b
+        x2, y2 = other.c - other.a, other.c + other.a
+        a_ = self.a[0] * x2 + self.a[1] * x1 + self.b[0] * y2 + self.b[1] * y1
+        c_ = self.c[0] * x2 + self.c[1] * x1 + self.d[0] * y2 + self.d[1] * y1
+        d_ = self.a[0] * x1 + self.b[0] * y1 + (self.a[1] * x2 + self.b[1] * y2) // 2
+        b_ = self.c[0] * x1 + self.d[0] * y1 + (self.c[1] * x2 + self.d[1] * y2) // 2
+        return ZOmega((c_ - a_) // 2, b_, (c_ + a_) // 2, d_)
+
+    @classmethod
+    def from_string(cls, string: str) -> GridOp:
+        """Return the grid operation from a string.
+
+        Args:
+            string: Supported string are: "I", "R", "A", "B", "K", "X", "Z".
+                (Fig 6, arXiv:1403.2975).
+
+        Returns:
+            GridOp: The grid operation corresponding to the string.
+        """
+        return _useful_grid_ops().get(string, "I")
+
+    @property
+    def determinant(self) -> tuple[float, float]:
+        r"""Calculate the determinant of the grid operation.
+
+        The determinant will be of form :math:`a + b / \sqrt{2}` and is given as tuple(a, b).
+        """
+        return (
+            self.a[0] * self.d[0]
+            - self.b[0] * self.c[0]
+            + (self.a[1] * self.d[1] - self.b[1] * self.c[1]) // 2,
+            self.a[0] * self.d[1]
+            - self.b[0] * self.c[1]
+            + self.a[1] * self.d[0]
+            - self.b[1] * self.c[0],
+        )
+
+    @property
+    def is_special(self) -> bool:
+        """Check if the grid operation is special based on Proposition 5.13 of arXiv:1403.2975."""
+        det1, det2 = self.determinant
+        return det2 == 0 and det1 in {1, -1}
+
+    @property
+    def flatten(self) -> list[float]:
+        """Flatten the grid operation based on description in Lemma 5.11 of arXiv:1403.2975."""
+        return [
+            self.a[0] + self.a[1] / _SQRT2,
+            self.b[0] + self.b[1] / _SQRT2,
+            self.c[0] + self.c[1] / _SQRT2,
+            self.d[0] + self.d[1] / _SQRT2,
+        ]
+
+    def inverse(self) -> "GridOp":
+        """Compute the inverse of the grid operation."""
+        det1, det2 = self.determinant
+        if det2 == 0 and det1 in {1, -1}:
+            return GridOp(
+                tuple(d // det1 for d in self.d),
+                tuple(-b // det1 for b in self.b),
+                tuple(-c // det1 for c in self.c),
+                tuple(a // det1 for a in self.a),
+            )
+        raise ValueError("Grid operator needs to be special to have an inverse.")
+
+    def transpose(self) -> "GridOp":
+        """Transpose the grid operation."""
+        return GridOp(
+            self.d,
+            self.b,
+            self.c,
+            self.a,
+        )
+
+    def adj2(self) -> "GridOp":
+        """Compute the sqrt(2)-conjugate of the grid operation."""
+        return GridOp(
+            (self.a[0], -self.a[1]),
+            (self.b[0], -self.b[1]),
+            (self.c[0], -self.c[1]),
+            (self.d[0], -self.d[1]),
+        )
+
+    def apply_to_ellipse(self, ellipse: Ellipse) -> Ellipse:
+        """Apply the grid operator to an ellipse based on Lemma A.4's proof in arXiv:1403.2975."""
+        return ellipse.apply_grid_op(self)
+
+    def apply_to_state(self, state: EllipseState) -> tuple[Ellipse, Ellipse]:
+        """Apply the grid operator to a state based on Definition A.3 of arXiv:1403.2975."""
+        return state.apply_grid_op(self)
+
+    def apply_shift_op(self, k: int) -> "GridOp":
+        """Apply a shift operator to the grid operation based on Lemma A.9 of arXiv:1403.2975."""
+        k, sign = abs(k), (-1) ** (k < 0)  # +1 if k > 0, -1 if k < 0
+        grid_op = self
+        for _ in range(k):
+            grid_op = GridOp(
+                (grid_op.a[1] + sign * grid_op.a[0], 2 * grid_op.a[0] + sign * grid_op.a[1]),
+                grid_op.b,
+                grid_op.c,
+                (grid_op.d[1] - sign * grid_op.d[0], 2 * grid_op.d[0] - sign * grid_op.d[1]),
+            )
+        return grid_op
+
+
+@lru_cache(maxsize=1)
+def _useful_grid_ops() -> dict["str", "GridOp"]:
+    """Generate a list of useful grid operations based on Fig 6 of arXiv:1403.2975."""
+    return {
+        "I": GridOp((1, 0), (0, 0), (0, 0), (1, 0)),
+        "R": GridOp((0, 1), (0, -1), (0, 1), (0, 1)),
+        "A": GridOp((1, 0), (-2, 0), (0, 0), (1, 0)),
+        "B": GridOp((1, 0), (0, 2), (0, 0), (1, 0)),
+        "K": GridOp((-1, 1), (0, -1), (1, 1), (0, 1)),
+        "X": GridOp((0, 0), (1, 0), (1, 0), (0, 0)),
+        "Z": GridOp((1, 0), (0, 0), (0, 0), (-1, 0)),
+    }
+
+
+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,
+    and a unit disk.
+
+    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.
+    """
+
+    def __init__(self, theta: float = 0.0, epsilon: float = 1e-3, max_trials: int = 20):
+        self.theta = theta
+        self.epsilon = epsilon
+        self.zval = math.cos(theta), math.sin(theta)
+        self.kmin = int(3 * math.log2(1 / epsilon) // 2)
+        self.max_trials = max_trials
+        self.target = 1 - self.epsilon**2 / 2
+
+    def __repr__(self) -> str:
+        """Return a string representation of the grid iterator."""
+        return f"GridIterator(theta={self.theta}, epsilon={self.epsilon}, max_trials={self.max_trials})"
+
+    def __iter__(self) -> Iterable[tuple[ZOmega, int]]:
+        """Iterate over the solutions to the scaled grid problem."""
+        # Warm start for an initial guess, where 14 is kmin for 1e-3.
+        k, i_ = min(self.kmin, 14), 6  # Give 6 trials for warm start.
+        e_, t_ = max(self.epsilon, 1e-3), min(self.target, 0.9999995)
+
+        e1 = Ellipse.from_region(self.theta, e_, k)  # Ellipse for the epsilon-region.
+        e2 = Ellipse.from_axes(p=(0, 0), theta=0, axes=(1, 1))  # Ellipse for the unit disk.
+        en, _ = e1.normalize()  # Normalize the epsilon-region.
+        grid_op = EllipseState(en, e2).skew_grid_op()  # Skew grid operation for the epsilon-region.
+
+        if k == self.kmin:
+            g_ = grid_op
+        else:
+            en_, _ = Ellipse.from_region(self.theta, self.epsilon, k).normalize()
+            g_ = EllipseState(en_, e2).skew_grid_op()
+
+        int_s, init_k = [ZOmega(d=1)], [k]  # Fallback solution.
+        guess_solutions = (  # Solutions for the trivial cases.
+            ZOmega(a=-1),
+            ZOmega(b=-1),
+            ZOmega(b=1),
+            ZOmega(c=1),
+            ZOmega(d=1),
+            ZOmega(d=-1),
+        )
+
+        for ix in range(self.max_trials):
+            # Update the radius of the unit disk.
+            radius = 2 ** (k // 2) * (math.sqrt(2) ** (k % 2))
+            e2 = Ellipse.from_axes(p=(0, 0), theta=0, axes=(radius, radius))
+            # Apply the grid operation to the state and solve the two-dimensional grid problem.
+            state = EllipseState(e1, e2).apply_grid_op(grid_op)
+            potential_solutions = self.solve_two_dim_problem(state)
+            try:
+                for solution in chain(guess_solutions, potential_solutions):
+                    # Normalize the solution and obtain the scaling exponent of sqrt(2).
+                    scaled_sol, kf = (grid_op * solution).normalize()
+
+                    complx_sol = complex(scaled_sol)
+                    sol_real, sol_imag = complx_sol.real, complx_sol.imag
+                    norm_zsqrt_two = float(scaled_sol.norm().to_sqrt_two())
+
+                    k_ = k - kf  # Update the scaling exponent of sqrt(2).
+                    dot_prod = (self.zval[0] * sol_real + self.zval[1] * sol_imag) / (
+                        2 ** (k_ // 2) * (math.sqrt(2) ** (k_ % 2))
+                    )
+                    # Check if the solution is follows the constraints of the target-region.
+                    if abs(norm_zsqrt_two) <= 2**k_:
+                        if dot_prod >= self.target:
+                            yield scaled_sol, k_
+                        elif dot_prod >= t_:
+                            int_s.append(scaled_sol)
+                            init_k.append(k_)
+
+                if ix == i_:
+                    k, e_, grid_op, t_ = max(self.kmin, k + 1), self.epsilon, g_, t_ / 10
+                else:
+                    k = k + 1
+
+                e1 = Ellipse.from_region(self.theta, e_, k)
+
+            except ValueError:  # pragma: no cover
+                break
+
+        for s, k in zip(int_s, init_k):
+            yield s, k
+
+    def solve_two_dim_problem(
+        self, state: EllipseState, num_points: int = 1000
+    ) -> Iterable[ZOmega]:
+        r"""Solve the grid problem for the state(E1, E2).
+
+        The solutions :math:`u \in Z[\omega]` are such that :math:`u \in E1` and
+        :math:`u.adj2() \in E2`, where ``adj2`` is :math:`\sqrt(2)` conjugation.
+
+        This is based on Proposition 5.21 and Theorem 5.18 of arXiv:1403.2975.
+
+        Args:
+            state: The state of the grid problem.
+            num_points: The number of points to use to determine if the rectangle is wider than the other.
+
+        Returns:
+            Iterable[ZOmega]: The list of solutions to the two dimensional grid problem.
+        """
+        e1, e2 = state.e1, state.e2
+        state2 = EllipseState(e1.offset(-1 / _SQRT2), e2.offset(1 / _SQRT2))
+
+        bbox1 = e1.bounding_box()
+        bbox11 = tuple(bb_ + e1.p[ix_ // 2] for ix_, bb_ in enumerate(bbox1))
+        bbox12 = tuple(bb_ - 1 / _SQRT2 for bb_ in bbox11)
+
+        bbox2 = e2.bounding_box()
+        bbox21 = tuple(bb_ + e2.p[ix_ // 2] for ix_, bb_ in enumerate(bbox2))
+        bbox22 = tuple(bb_ + 1 / _SQRT2 for bb_ in bbox21)
+
+        # Check if it is easier to solve problem for either of x-or-y-interval.
+        # Based on this, we can try to first solve for alpha-or-beta and then refine for other.
+        # If both of them are balanced, rely on doing naive search over both.
+        bbox_zip = list(zip((bbox11, bbox12), (bbox21, bbox22)))
+        num_x1, num_x2 = (self.bbox_grid_points(bb1[:2] + bb2[:2]) for bb1, bb2 in bbox_zip)
+        num_y1, num_y2 = (self.bbox_grid_points(bb1[2:] + bb2[2:]) for bb1, bb2 in bbox_zip)
+        num_b1 = [num_x1 > num_points * num_y1, num_y1 > num_points * num_x1]
+        num_b2 = [num_x2 > num_points * num_y2, num_y2 > num_points * num_x2]
+
+        # Solve the problem for the two cosets of ZOmega ring and add non-zero offset to odd one.
+        potential_sols1 = self.solve_upright_problem(state, bbox11, bbox21, num_b1, ZOmega())
+        potential_sols2 = self.solve_upright_problem(state2, bbox12, bbox22, num_b2, ZOmega(c=1))
+
+        for solution in chain(potential_sols1, potential_sols2):
+            sol1, sol2 = complex(solution), complex(solution.adj2())
+            x1, y1 = sol1.real, sol1.imag
+            x2, y2 = sol2.real, sol2.imag
+            if e1.contains(x1, y1) and e2.contains(x2, y2):
+                yield solution
+
+    # pylint:disable = too-many-arguments, too-many-branches
+    def solve_upright_problem(
+        self,
+        state: EllipseState,
+        bbox1: tuple[float],
+        bbox2: tuple[float],
+        num_b: list[bool],
+        shift: ZOmega,
+    ) -> Iterable[ZOmega]:
+        r"""Iterates over the solutions to the grid problem for two upright rectangles.
+
+        The solutions :math:`u \in Z[\omega]` are such that :math:`u \in A` and
+        :math:`u.adj2() \in B`, where ``adj2`` is :math:`\sqrt(2)` conjugation
+        and two rectangles :math:`A` and :math:`B`, form the subregions of
+        :math:`\mathbb{R}^2` of the form :math:`[x0, x1] \times [y0, y1]`.
+
+        Args:
+            state (State): The state of the grid problem.
+            bbox1 (tuple[float]): The bounding box of the first rectangle.
+            bbox2 (tuple[float]): The bounding box of the second rectangle.
+            num_b (list[bool]): Whether the second rectangle is wider than the first.
+            shift (ZOmega): The shift operator.
+
+        Returns:
+            Iterable[ZOmega]: The list of solutions to the upright grid problem for two rectangles.
+        """
+        e1, e2 = state.e1, state.e2
+        Ax0, Ax1, Ay0, Ay1 = bbox1
+        Bx0, Bx1, By0, By1 = bbox2
+
+        if num_b[0]:  # pragma: no cover # If it is easier to solve for beta first.
+            beta_solutions1 = self.solve_one_dim_problem(Ay0, Ay1, By0, By1)
+            for beta in beta_solutions1:
+                Ax0_tmp, Ax1_tmp = e1.x_points(beta)
+                Bx0_tmp, Bx1_tmp = e2.x_points(beta.adj2())
+                if Ax1_tmp - Ax0_tmp > 0 and Bx1_tmp - Bx0_tmp > 0:
+                    new_alpha_solutions = self.solve_one_dim_problem(
+                        Ax0_tmp, Ax1_tmp, Bx0_tmp, Bx1_tmp
+                    )
+                    for alpha in new_alpha_solutions:
+                        yield ZOmega.from_sqrt_pair(alpha, beta, shift)
+        elif num_b[1]:  # If it is easier to solve for alpha first.
+            alpha_solutions1 = self.solve_one_dim_problem(Ax0, Ax1, Bx0, Bx1)
+            for alpha in alpha_solutions1:
+                Ay0_tmp, Ay1_tmp = e1.y_points(alpha)
+                By0_tmp, By1_tmp = e2.y_points(alpha.adj2())
+                if Ay1_tmp - Ay0_tmp > 0 and By1_tmp - By0_tmp > 0:
+                    new_beta_solutions = self.solve_one_dim_problem(
+                        Ay0_tmp, Ay1_tmp, By0_tmp, By1_tmp
+                    )
+                    for beta in new_beta_solutions:
+                        yield ZOmega.from_sqrt_pair(alpha, beta, shift)
+        else:  # If both of them are balanced, solve for both and refine.
+            alpha_solutions1 = self.solve_one_dim_problem(Ax0, Ax1, Bx0, Bx1)
+            beta_solutions1 = self.solve_one_dim_problem(Ay0, Ay1, By0, By1)
+            found_beta1_solutions = []
+            for alpha in alpha_solutions1:
+                if len(found_beta1_solutions) == 0:
+                    for beta in beta_solutions1:
+                        found_beta1_solutions.append(beta)
+                        yield ZOmega.from_sqrt_pair(alpha, beta, shift)
+                else:
+                    for beta in found_beta1_solutions:
+                        yield ZOmega.from_sqrt_pair(alpha, beta, shift)
+
+    @staticmethod
+    def bbox_grid_points(bbox: tuple[float, float, float, float]) -> int:
+        """Count the number of grid points in a bounding box.
+
+        This gives an estimation on the expected number of solution within
+        the bounding box. This is based on the Lemma 16 of arXiv:1212.6253.
+
+        Args:
+            bbox (tuple[float, float, float, float]): The bounding box.
+
+        Returns:
+            int: The number of grid points in the bounding box.
+        """
+        d_ = math.log2(ZSqrtTwo(1, 1))
+        l1, l2 = ZSqrtTwo(1, 1), ZSqrtTwo(-1, 1)
+        d1, d2 = (bbox[1] - bbox[0], bbox[3] - bbox[2])
+
+        # Find the integer scaling factor for the x and y intervals, such that
+        # \delta_{1/2} \cdot (\lambda - 1)^{k_{1/2}} < 1, where \lambda= 1/√2.
+        k1, k2 = (int(math.floor(math.log2(d + 1e-20) / d_ + 1)) for d in (d1, d2))
+        if abs(k1) > abs(k2):  # If y-interval is wider than x-interval, swap.
+            bbox, k1, k2 = (bbox[2], bbox[3], bbox[0], bbox[1]), k2, k1
+
+        # Scale the x and y intervals to enter the intended interval.
+        # Look at `solve_one_dim_problem` for more details.
+        x_scale = float((l1 if k1 < 0 else l2) ** abs(k1))
+        y_scale = float((-1) ** k1 * (l2 if k1 < 0 else l1) ** abs(k1))
+
+        x0_scaled, x1_scaled = x_scale * bbox[0], x_scale * bbox[1]
+        y0_scaled, y1_scaled = sorted((y_scale * bbox[2], y_scale * bbox[3]))
+
+        # Check if we are indeed within the intended interval.
+        if x1_scaled - x0_scaled < 1 - _SQRT2:  # pragma: no cover
+            raise ValueError(f"Value should be larger than 1 - sqrt(2) for bbox {bbox}")
+
+        # Use the constraints x0 <= a + b * sqrt(2) <= x1 and y0 <= a - b * sqrt(2) <= y1
+        # to obtain the bounds on b and eliminate the variable a to obtain the bounds on b.
+        lower_bound_b = (x0_scaled - y1_scaled) / (2 * _SQRT2)
+        upper_bound_b = (x1_scaled - y0_scaled) / (2 * _SQRT2)
+        return 1 + int(upper_bound_b - lower_bound_b)
+
+    @staticmethod
+    def solve_one_dim_problem(
+        x0: float,
+        x1: float,
+        y0: float,
+        y1: float,
+    ) -> Iterable[ZSqrtTwo]:
+        r"""Iterates the solutions to the one dimensional grid problem given intervals :math:`[x0, x1]` and :math:`[y0, y1]`.
+
+        Given two real intervals :math:`[x0, x1]` and :math:`[y0, y1]`
+        such that :math:`\sqrt{(x1 - x0)*(y1 - y0)} >= (1 + \sqrt(2))`,
+        iterates over all solutions of the form :math:`a + b\sqrt(2)` such that
+        :math:`a + b\sqrt(2) \in [x0, x1]` and :math:`a - b\sqrt(2) \in [y0, y1]`.
+
+        This is based on the Lemmas 16 and 17 of arXiv:1212.6253.
+
+        Args:
+            x0 (float): The lower bound of the x-interval.
+            x1 (float): The upper bound of the x-interval.
+            y0 (float): The lower bound of the y-interval.
+            y1 (float): The upper bound of the y-interval.
+
+        Returns:
+            Iterable[ZSqrtTwo]: The list of solutions to the one dimensional grid problem.
+        """
+        d_ = math.log2(ZSqrtTwo(1, 1))
+        l1, l2 = ZSqrtTwo(1, 1), ZSqrtTwo(-1, 1)
+        d1, d2 = (x1 - x0, y1 - y0)
+
+        f_adj2 = False  # Check if we need to apply the sqrt(2) conjugation.
+        # Find the integer scaling factor for the x and y intervals, such that
+        # \delta_{1/2} \cdot (\lambda - 1)^{k_{1/2}} < 1, where \lambda= 1/√2.
+        k1, k2 = (int(math.floor(math.log2(d + 1e-20) / d_ + 1)) for d in (d1, d2))
+        if abs(k1) > abs(k2):  # If y-interval is wider than x-interval, swap.
+            f_adj2, k1, k2 = True, k2, k1
+            x0, x1, y0, y1 = y0, y1, x0, x1
+
+        # Turn the problem into a scaled grid problem,
+        # such that we get to solve the specific case of:
+        # -1 + √2 <= x1 - x0 < 1 and (x1 - x0)(y1 - y0) >= (1 + √2)^2.
+        s_scale = ZSqrtTwo(1, 1) ** abs(k1)
+        x_scale = float((l1 if k1 < 0 else l2) ** abs(k1))
+        y_scale = float((-1) ** k1 * (l2 if k1 < 0 else l1) ** abs(k1))
+
+        x0_scaled, x1_scaled = x_scale * x0, x_scale * x1
+        y0_scaled, y1_scaled = sorted((y_scale * y0, y_scale * y1))
+
+        # Check if we are solving the problem for the intended interval.
+        if x1_scaled - x0_scaled < 1 - _SQRT2:  # pragma: no cover
+            bbox = (x0_scaled, x1_scaled, y0_scaled, y1_scaled)
+            raise ValueError(f"Value should be larger than 1 - sqrt(2) for bbox {bbox}")
+
+        # Use the constraints y0 <= a - b * sqrt(2) <= y1 and x0 <= a + b * sqrt(2) <= x1
+        # to obtain the bounds on b and eliminate the variable a to obtain the bounds on b.
+        lower_bound_b = (x0_scaled - y1_scaled) / (2 * _SQRT2)
+        upper_bound_b = (x1_scaled - y0_scaled) / (2 * _SQRT2)
+
+        for b in range(int(math.floor(upper_bound_b)), int(math.ceil(lower_bound_b)) - 1, -1):
+            # Use the constraints x0 <= a + b * sqrt(2) <= x1 to obtain the bounds on a.
+            lower_bound_a = x0_scaled - b * _SQRT2
+            upper_bound_a = x1_scaled - b * _SQRT2
+            if upper_bound_a - lower_bound_a >= 1:  # pragma: no cover
+                raise ValueError(f"Value should be less than one for {(x0, x1, y0, y1)}")
+
+            # Check if the bounds on the interval contains an integer.
+            if math.ceil(lower_bound_a) == math.floor(upper_bound_a):
+                a = int(math.ceil(lower_bound_a))
+                # Check if the solution satisfies both bounds on x and y.
+                if (x0_scaled + y0_scaled <= 2 * a) and (2 * a <= x1_scaled + y1_scaled):
+                    alpha, beta = a + b * _SQRT2, a - b * _SQRT2
+                    # Check if the consecutive solutions are within the desired bounds.
+                    if x0_scaled <= alpha <= x1_scaled and y0_scaled <= beta <= y1_scaled:
+                        # Undo the scaling to obtain the solution.
+                        sol = ZSqrtTwo(a, b) / s_scale if k1 < 0 else ZSqrtTwo(a, b) * s_scale
+                        yield sol if not f_adj2 else sol.adj2()
diff --git a/pennylane/ops/op_math/decompositions/norm_solver.py b/pennylane/ops/op_math/decompositions/norm_solver.py
new file mode 100644
index 000000000..fdca4a5de
--- /dev/null
+++ b/pennylane/ops/op_math/decompositions/norm_solver.py
@@ -0,0 +1,416 @@
+# 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 module contains mathematical operations for solving diophantine equations."""
+
+import math
+from functools import lru_cache, singledispatch
+from random import randrange
+
+from pennylane.ops.op_math.decompositions.rings import ZOmega, ZSqrtTwo
+
+
+@singledispatch
+def _gcd(a, b):
+    """Compute the greatest common divisor of two integers."""
+    raise NotImplementedError(
+        f"GCD is not implemented for types {type(a)} and {type(b)}. "
+        "Please implement a custom GCD function for these types."
+    )  # pragma: no cover
+
+
+@_gcd.register(int)
+def _(elem1, elem2):
+    return math.gcd(elem1, elem2)
+
+
+@_gcd.register(ZSqrtTwo)
+@_gcd.register(ZOmega)
+def _(elem1, elem2):
+    while elem2 != 0:
+        elem1, elem2 = elem2, elem2 % elem1
+    return elem1
+
+
+@lru_cache(maxsize=100)
+def _prime_factorize(n: int, max_trials=1000, z_sqrt_two: bool = True) -> list[int] | None:
+    r"""Computes the prime factorization of a number :math:`n`.
+
+    This function uses a combination of the Brent's variant of Pollard's rho algorithm for integer factorization
+    and Miller-Rabin primality test to find the prime factors of :math:`n`. It handles trivial cases and uses a
+    stack-based approach to iteratively find factors until all prime factors are identified.
+
+    .. note::
+        The function returns ``None`` cases where prime factor might be 7,
+        as it cannot be expressed in the ring :math:`\mathbb{Z}[\sqrt{2}]`.
+
+    Args:
+        n (int): The number to factor.
+        max_tries (int): The maximum number of attempts to find a factor.
+        z_sqrt_two (bool): When ``True``, the function will only return the factorization
+            iff it can be expressed in :math:`\mathbb{Z}[\sqrt{2}]`. Therefore, if a prime
+            factor :math:`p ≡ 7 mod 8` is encountered, the function will return ``None``.
+            When ``False``, the function will return all factors. Default is ``True``.
+
+    Returns:
+        list[int] | None: A list of prime factors of :math:`n`, or ``None`` if no valid factors are found.
+    """
+    factors, stack = [], [n]
+    while len(stack) > 0:
+        p = stack.pop()
+        # Trivial case
+        if p <= 1:
+            continue
+        # Cannot split in Z[√2], if p ≡ 7 mod 8.
+        if _primality_test(p):
+            if z_sqrt_two and p % 8 == 7:
+                return None
+            factors.append(p)
+            continue
+        if (factor := _integer_factorize(p, max_trials)) is None:  # pragma: no cover
+            return None
+        if z_sqrt_two and factor % 7 == 0:  # pragma: no cover
+            return None
+        # If we have found an integer factor,
+        # push it and its complement onto the stack
+        stack.extend([factor, p // factor])
+
+    return sorted(factors)
+
+
+@lru_cache(maxsize=100)
+def _integer_factorize(n: int, max_tries=1000) -> int | None:
+    r"""Computes an integer factor of a number :math:`n`.
+
+    This function implements the `Brent's variant <https://doi.org/10.1007/BF01933190>`_
+    of the Pollard's rho algorithm for integer factorization.
+
+    Args:
+        n (int): The number to factor.
+        max_tries (int): The maximum number of attempts to find a factor.
+
+    Returns:
+        int | None: An integer factor of :math:`n`, or ``None`` if no factors are found.
+    """
+    if n <= 2:
+        return None
+
+    if n % 2 == 0:
+        return 2
+
+    # Main loop: retry with different parameters on failure
+    while max_tries > 0:
+        x, c, m = randrange(1, n), randrange(1, n), randrange(1, n)
+        g, r, q = 1, 1, 1
+
+        while g == 1:
+            y = x
+            # Process next `r` steps
+            for _ in range(r):
+                x = (x * x % n + c) % n
+
+            k = 0
+            while k < r and g == 1:
+                xs, rs = x, min(m, r - k)
+                # Process next `min(m, r-k)`` steps
+                for _ in range(rs):
+                    x = (x * x % n + c) % n
+                    q = (q * abs(y - x)) % n
+                g = _gcd(q, n)
+                k += m
+            r <<= 1
+
+        # Linear search for a factor if the above doesn't yield a result.
+        if g == n:
+            g, y = 1, xs
+            while g == 1:
+                y = (y * y % n + c) % n
+                g = _gcd(abs(x - y), n)
+
+        # If we found a valid integer factor, such that it is neither 1 nor n.
+        if g not in (1, n):
+            return g
+
+        max_tries -= 1
+
+    return None
+
+
+def _factorize_prime_zsqrt_two(p: int) -> list[ZSqrtTwo]:
+    r"""Find the factorization of a prime number :math:`p` in ring :math:`\mathbb{Z}[\sqrt{2}]`.
+
+    This uses theory from Appendix C.2 of `arXiv:1403.2975 <https://arxiv.org/abs/1403.2975>`_,
+    which includes the following results:
+
+        - Lemma C.7: Prime factorization of an integer prime :math:`p` in
+          :math:`\mathbb{Z}[\sqrt{2}]` consists of one or two factors.
+        - Lemma C.8: Prime factorization of :math:`\pm 2 = (0 + 1\sqrt{2})(0 \pm 1\sqrt{2})`
+          in :math:`\mathbb{Z}[\sqrt{2}]`.
+        - Lemma C.9: Prime factorization of :math:`p \equiv 3 \mod 8` and :math:`p \equiv 5 \mod 8`
+          will have one factor in :math:`\mathbb{Z}[\sqrt{2}]`.
+        - Lemma C.11: Prime factorization of :math:`p \equiv 1 \mod 8` and :math:`p \equiv 7 \mod 8`
+          will have two factors in :math:`\mathbb{Z}[\sqrt{2}]`.
+
+    Args:
+        n (int): The prime number for which to find the factorization.
+
+    Returns:
+        list[ZSqrtTwo]: A list of factors in the ring :math:`\mathbb{Z}[\sqrt{2}]`,
+            or `None` if no factorization exists.
+    """
+    if abs(p) == 2:
+        # Lemma C.8: ±2 = (0 + 1√2)(0 ± 1√2)
+        return [ZSqrtTwo(0, 1), ZSqrtTwo(0, (-1) ** (p < 0))]
+
+    if p % 8 in (3, 5):
+        # Lemma C.9: p = (p + 0√2) is prime in Z[√2]
+        return [ZSqrtTwo(p, 0)]
+
+    # Default case (Lemma C.11): p ≡ ±1 mod 8 |  p % 8 in (1, 7)
+    # Solve t^2 ≡ 2 (mod p) to split p
+    if (t := _sqrt_modulo_p(2, p)) is None:
+        return None
+    # Perform ring GCD to get (a + b√2)(a - b√2)
+    res = _gcd(ZSqrtTwo(p, 0), ZSqrtTwo(min(t, p - t), 1))
+    return [res, res.adj2()]
+
+
+# pylint: disable=too-many-return-statements
+def _factorize_prime_zomega(x: ZSqrtTwo, p: int) -> ZOmega | None:
+    r"""Find a prime factor of an element :math:`x` in the ring :math:`\mathbb{Z}[\omega]`,
+    where :math:`x` divides a prime integer :math:`p`.
+
+    This function uses theory from Appendix C.3 of `arXiv:1403.2975 <https://arxiv.org/abs/1403.2975>`_,
+    which includes the following results:
+
+        - Lemma C.13: Prime factorization of a prime :math:`p` in :math:`\mathbb{Z}[\sqrt{2}]`
+          has one or two factors in :math:`\mathbb{Z}[\omega]` consists.
+        - Lemma C.20-21: Prime factorization of a prime :math:`p` in :math:`\mathbb{Z}[\sqrt{2}]`\
+          is can be used as a possible solution to the Diophantine equation :math:`t^* t = \xi`
+          iff :math:`p == 2` or :math:`p \equiv 1, 3, 5 \mod 8`.
+        - Lemma C.20: For the cases, :math:`p \equiv 3 \mod 8` and :math:`p \equiv 1, 5 \mod 8`,
+          quadratic reciprocity in the proof of Lemma C.20 gives :math:`-2` and :math:`-1` as
+          square modulo :math:`p`, which are used for the factorization.
+
+    Args:
+        x (ZSqrtTwo): The element in the ring :math:`\mathbb{Z}[\sqrt{2}]`, such that x | p.
+        p (int): The prime number for which to find the factorization.
+
+    Returns:
+        ZOmega: Factor of :math:`p` in the ring :math:`\mathbb{Z}[\omega]`,
+            or ``None`` if no factorization exists.
+    """
+    # Basic cases
+    if p == 2:
+        return ZOmega(0, 0, 1, 1)
+
+    # p = 2k or p ≡ 7 mod 8, no factorization in Z[ω]
+    if ((a := p % 8) % 2 == 0) or a == 7:
+        return None
+
+    # p = 1, 5 mod 8, use h = sqrt(-1) mod p
+    if a in (1, 5):
+        if (h := _sqrt_modulo_p(-1, p)) is None:  # pragma: no cover
+            return None
+        return _gcd(ZOmega(0, 1, 0, h), ZOmega(-x.b, 0, x.b, x.a))
+
+    # Default case: a == 3
+    # p = 3 mod 8, use h = = sqrt(-2) mod p
+    if (h := _sqrt_modulo_p(-2, p)) is None:
+        return None
+    return _gcd(ZOmega(1, 0, 1, h), ZOmega(-x.b, 0, x.b, x.a))
+
+
+@lru_cache(maxsize=100)
+def _primality_test(n: int) -> bool:
+    r"""Determines whether an integer is prime or not.
+
+    This function implements the deterministic variant of `Miller-Rabin primality test
+    <https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test>`_.
+
+    Args:
+        n (int): The number to test for primality.
+
+    Returns:
+        bool: ``True`` if :math:`n` is likely prime, ``False`` otherwise.
+    """
+    if n < 2 or n == 4:
+        return False
+    if n < 4:
+        return True
+
+    # Small primes quick check
+    ps = {5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97}
+    if n in ps:
+        return True
+
+    if any(n % p == 0 for p in ps):
+        return False
+
+    # Factor n−1 as d·2^s with d odd
+    d, s = n - 1, 0
+    while d & 1 == 0:
+        d >>= 1
+        s += 1
+
+    # Deterministic bases for testing 64-bit range [https://miller-rabin.appspot.com/]
+    bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
+    for base in bases:
+        base = base if base < n else base % n
+        if base == 0 or base < 2:  # pragma: no cover
+            continue
+        x = pow(base, d, n)
+        if x in (1, n - 1):
+            continue
+        for _ in range(s - 1):
+            x = pow(x, 2, n)
+            if x == 1:  # pragma: no cover
+                return False
+            if x == n - 1:
+                break
+        if x != n - 1:
+            return False
+
+    return True
+
+
+@lru_cache(maxsize=100)
+def _legendre_symbol(a: int, p: int) -> int:
+    r"""Computes the Legendre symbol :math:`\left(\frac{a}{p}\right)`.
+
+    This function uses the definition of the Legendre symbol:
+
+    .. math::
+        \left(\frac{a}{p}\right) = a^{(p-1)/2} \mod p
+
+    Args:
+        a (int): The number to compute the Legendre symbol of.
+        p (int): The prime number.
+
+    Returns:
+        int: The Legendre symbol of :math:`a` modulo :math:`p`.
+    """
+    return pow(a, (p - 1) // 2, p)
+
+
+def _sqrt_modulo_p(n: int, p: int) -> int | None:
+    r"""Computes square root of :math:`n` under modulo :math:`p` if it exists.
+
+    This uses `Tonelli-Shanks algorithm <https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm>`_
+    to compute :math:`x`, such that :math:`x^2 \equiv n (mod p)`, where :math:`p` is an odd prime.
+
+    Args:
+        n (int): The number to compute the square root of.
+        p (int): The odd prime modulus.
+
+    Returns:
+        int or None: The square root of :math:`n` under modulo :math:`p`, or None if it does not exist.
+    """
+    # Trivial cases
+    if (a := n % p) == 0:
+        return 0
+    if p == 2:
+        return a
+    if _legendre_symbol(a, p) != 1 or p % 2 == 0:
+        return None
+
+    # Factor p−1 as q·2^s with q odd
+    q, s = p - 1, 0
+    while q & 1 == 0:
+        q >>= 1
+        s += 1
+
+    # (p - 1) = 2.q ==> 2 * k
+    if s == 1:
+        return pow(a, (p + 1) // 4, p)
+
+    # Find a quadratic non‐residue z,
+    # such that z^((p-1)/2) ≡ −1 mod p
+    for z in range(2, p):
+        if _legendre_symbol(z, p) == (p - 1):
+            break
+
+    r = pow(a, (q + 1) // 2, p)
+    c = pow(z, q, p)
+    t = pow(a, q, p)
+
+    m = s
+    while t % p != 1:
+        # Find least i in [1, m),
+        # such that t^(2^i) ≡ 1 mod p
+        ix, t2 = 1, pow(t, 2, p)
+        while t2 % p != 1 and ix < m:
+            t2, ix = pow(t2, 2, p), ix + 1
+
+        # Compute b = c^(2^(m-ix-1)) mod p
+        # and update the intial elements
+        if m - ix - 1 < 0:  # pragma: no cover
+            return None
+        b = pow(c, 1 << (m - ix - 1), p)
+        r = (r * b) % p
+        c = pow(b, 2, p)
+        t = (t * c) % p
+        m = ix
+
+    return r
+
+
+def _solve_diophantine(xi: ZSqrtTwo, max_trials: int = 1000) -> ZOmega | None:
+    r"""Solve the Diophantine equation :math:`t^* t = \xi` for :math:`t \in \mathbb{Z}[\omega]`
+    and :math:`\xi \in \mathbb{Z}[\sqrt{2}]`.
+
+    This function uses the theory from Appendix C and D (Proof of Lemma 8.4) of
+    `arXiv:1403.2975 <https://arxiv.org/abs/1403.2975>`_.
+
+    Args:
+        xi (ZSqrtTwo): An element of the ring :math:`\mathbb{Z}[\sqrt{2}]`.
+
+    Returns:
+        ZOmega | None: An element of the ring :math:`\mathbb{Z}[\omega]` that satisfies the equation, or ``None`` if no solution exists.
+    """
+    if xi.a == 0 and xi.b == 0:
+        return ZOmega(0, 0, 0, 0)
+
+    if (p := abs(xi)) < 2:
+        return None
+
+    if not (factors := _prime_factorize(p, max_trials)):
+        return None
+
+    scale, next_xi = ZOmega(d=1), xi
+    for factor in factors:
+        if (primes_zsqrt_two := _factorize_prime_zsqrt_two(factor)) is None:
+            return None  # pragma: no cover
+
+        for eta in primes_zsqrt_two:
+            # Scale the next_xi by the factor in Z[√2].
+            next_xi, next_ab = next_xi * eta.adj2(), abs(eta)
+            # Check if the next_xi is divisible by the element eta in Z[√2].
+            if next_xi.a % next_ab == 0 and next_xi.b % next_ab == 0:
+                next_xi = ZSqrtTwo(next_xi.a // next_ab, next_xi.b // next_ab)
+                if (t := _factorize_prime_zomega(eta, factor)) is None:
+                    return None  # pragma: no cover
+                scale *= t
+
+    # the remaining quotient should be divisible
+    s_val = (scale.conj() * scale).to_sqrt_two()
+    s_new, s_abs = (xi * s_val.adj2()), abs(s_val)
+    if any(x_ % s_abs != 0 for x_ in s_new.flatten):
+        return None
+
+    # the remaining quotient should be a unit in Z[√2]
+    t2 = xi / s_val
+    if abs(t2) ** 2 != 1:
+        return None  # pragma: no cover
+
+    return scale * t2.sqrt().to_omega()
diff --git a/pennylane/ops/op_math/decompositions/normal_forms.py b/pennylane/ops/op_math/decompositions/normal_forms.py
new file mode 100644
index 000000000..0e0f7b137
--- /dev/null
+++ b/pennylane/ops/op_math/decompositions/normal_forms.py
@@ -0,0 +1,151 @@
+# 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.
+"""Decompositions for Matsumoto-Amano normal forms in SO(3) matrices."""
+
+from copy import deepcopy
+from functools import lru_cache
+
+import pennylane as qml
+from pennylane.ops.op_math.decompositions.rings import 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}`."""
+    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),
+    }
+    return {gate: (SO3Matrix(su2), phase) for gate, (su2, phase) in clifford_elems.items()}
+
+
+@lru_cache
+def _parity_transforms() -> dict:
+    r"""Returns information required to perform the parity transforms used in the _ma_normal_form.
+
+    Returns:
+        dict: A dictionary mapping parity vectors to tuples containing:
+            - 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.
+    """
+    # 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))),
+            SO3Matrix(DyadicMatrix(ZOmega(d=1), ZOmega(), ZOmega(), ZOmega(d=1))),
+            (qml.I(0),),
+            0.0,
+        ),  # Identity is used as a placeholder to represent arbitrary Clifford group element.
+        "T": (
+            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,
+        ),
+        "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,
+        ),
+        "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,
+        ),
+    }
+    return {
+        tuple(so3_1.parity_vec): (so3_2, gate, phase)
+        for (so3_1, so3_2, gate, phase) in transform_ops.values()
+    }
+
+
+def _ma_normal_form(
+    op: SO3Matrix, compressed=False
+) -> tuple[qml.operation.Operator, float] | tuple[int, tuple[int, ...], int, float]:
+    r"""Decompose an SO(3) matrix into Matsumoto-Amano normal form.
+
+    A Matsumoto-Amano normal form - :math:`(T | \epsilon) (HT | SHT)^* \mathcal{C}`, consists of a rightmost
+    Clifford operator, followed by any number of syllables of the form ``HT`` or ``SHT``, followed by an
+    optional syllable ``T`` [`arXiv:1312.6584 <https://arxiv.org/abs/1312.6584>`_\ ].
+
+    Args:
+        op (SO3Matrix): The SO(3) matrix to decompose.
+        compressed (bool): If ``True``, the output will be a tuple containing information about the decomposition
+            in terms of bits and indices. If ``False``, the output will be a list of all gates in the decomposition.
+            Default is ``False``.
+
+    Returns:
+        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()
+
+    so3_op = deepcopy(op)
+
+    # The following use lemmas from arXiv:1312.6584.
+    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
+
+    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)
+            cl_index = clifford_index
+            g_phase += clifford_data[1]
+            break
+
+    if not compressed:
+        return decomposition, g_phase
+
+    t_bit = int(decomposition[0] == qml.T(0))
+    c_bit = max(0, cl_index)
+
+    return (t_bit, tuple(rep_bits[t_bit:]), c_bit, g_phase)
diff --git a/pennylane/ops/op_math/decompositions/rings.py b/pennylane/ops/op_math/decompositions/rings.py
new file mode 100644
index 000000000..a1767685f
--- /dev/null
+++ b/pennylane/ops/op_math/decompositions/rings.py
@@ -0,0 +1,689 @@
+# 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 module provides classes for algebraic prerequisites for Clifford+T decompositions"""
+
+
+from __future__ import annotations
+
+import math
+from copy import deepcopy
+from typing import List
+
+import numpy as np
+
+_SQRT2 = math.sqrt(2)
+_OMEGA = (1 + 1j) / _SQRT2
+
+
+class ZSqrtTwo:
+    r"""Represents the elements of the ring of integers adjoined with the square root of 2.
+
+    ..math::
+        \mathbb{Z}[\sqrt{2}] = \{ a + b\sqrt{2} \mid a, b \in \mathbb{Z} \}
+
+    Args:
+        a (int): The integer part of the element.
+        b (int): The coefficient of sqrt(2) in the element.
+    """
+
+    def __init__(self, a: int = 0, b: int = 0) -> None:
+        self.a = int(a)
+        self.b = int(b)
+
+    def __str__(self: ZSqrtTwo) -> str:
+        if self.b == 0:
+            return str(self.a)
+        if self.a == 0:
+            return f"{self.b}√2"
+        return f"{self.a} + {self.b}√2"
+
+    def __repr__(self: ZSqrtTwo) -> str:
+        return f"ZSqrtTwo(a={self.a}, b={self.b})"
+
+    def __float__(self: ZSqrtTwo) -> float:
+        return float(self.a) + float(self.b) * _SQRT2
+
+    def __add__(self, other: ZSqrtTwo | int | float) -> ZSqrtTwo:
+        if isinstance(other, ZSqrtTwo):
+            return ZSqrtTwo(self.a + other.a, self.b + other.b)
+        if isinstance(other, int) or (isinstance(other, float) and other.is_integer()):
+            return ZSqrtTwo(self.a + int(other), self.b)
+        raise TypeError(f"Unsupported type {type(other)} for addition with ZSqrtTwo")
+
+    def __sub__(self, other: ZSqrtTwo | int | float) -> ZSqrtTwo:
+        return self + (-other)
+
+    def __mul__(self, other: ZSqrtTwo | int | float) -> ZSqrtTwo:
+        if isinstance(other, ZSqrtTwo):
+            return ZSqrtTwo(
+                self.a * other.a + 2 * self.b * other.b,
+                self.a * other.b + self.b * other.a,
+            )
+        if isinstance(other, int) or (isinstance(other, float) and other.is_integer()):
+            return ZSqrtTwo(self.a * int(other), self.b * int(other))
+        raise TypeError(f"Unsupported type {type(other)} for multiplication with ZSqrtTwo")
+
+    __radd__ = __add__
+    __rsub__ = __sub__
+    __rmul__ = __mul__
+
+    def __pow__(self, power: int) -> ZSqrtTwo:
+        if power == 0:
+            return ZSqrtTwo(1, 0)
+        if power < 0:
+            raise ValueError(f"Negative powers {power} are unsupported for ZSqrtTwo.")
+        if isinstance(power, float) and not power.is_integer():
+            raise ValueError(f"Non-integer powers {power} are unsupported for ZSqrtTwo.")
+        result = self
+        while power > 1:
+            result *= self
+            power -= 1
+        return result
+
+    def __truediv__(self, other: ZSqrtTwo) -> ZSqrtTwo:
+        if isinstance(other, ZSqrtTwo):
+            return (self * other.adj2()) / abs(other)
+        if isinstance(other, int) or (isinstance(other, float) and other.is_integer()):
+            other = int(other)
+            if self.a % other == 0 and self.b % other == 0:
+                return ZSqrtTwo(self.a // other, self.b // other)
+        raise TypeError(f"Unsupported type {type(other)} for dividing ZSqrtTwo")
+
+    def __floordiv__(self, other: int) -> ZSqrtTwo:
+        if isinstance(other, int):
+            return ZSqrtTwo(self.a // other, self.b // other)
+        raise TypeError(f"Unsupported type {type(other)} for floor dividing ZSqrtTwo")
+
+    def __eq__(self, other) -> bool:
+        if isinstance(other, ZSqrtTwo):
+            return self.a == other.a and self.b == other.b
+        if isinstance(other, int) or (isinstance(other, float) and other.is_integer()):
+            return self.a == int(other) and self.b == 0
+        return np.isclose(float(self), float(other))
+
+    def __abs__(self) -> float:
+        return self.a**2 - 2 * self.b**2
+
+    def __neg__(self) -> ZSqrtTwo:
+        return ZSqrtTwo(-self.a, -self.b)
+
+    def __mod__(self, other: ZSqrtTwo | int | float) -> ZSqrtTwo:
+        if isinstance(other, int) or (isinstance(other, float) and other.is_integer()):
+            return ZSqrtTwo(self.a % int(other), self.b % int(other))
+
+        d = abs(other)
+        n1, n2 = (self.a * other.a - 2 * self.b * other.b), (self.b * other.a - self.a * other.b)
+        return self - ZSqrtTwo(round(n1 / d), round(n2 / d)) * other
+
+    @property
+    def flatten(self: ZSqrtTwo) -> List[int]:
+        """Flatten to a list."""
+        return [self.a, self.b]
+
+    def conj(self) -> ZSqrtTwo:
+        r"""Return the complex conjugate.
+
+        .. math::
+            (a + b\sqrt{2})^{\dagger} = a + b\sqrt{2}
+        """
+        return ZSqrtTwo(self.a, self.b)
+
+    def adj2(self) -> ZSqrtTwo:
+        r"""Return the adjoint, i.e., the root-2 conjugate.
+
+        .. math::
+            (a + b\sqrt{2})^{\bullet} = a - b\sqrt{2}
+        """
+        return ZSqrtTwo(self.a, -self.b)
+
+    def sqrt(self) -> ZSqrtTwo | None:
+        """Return the square root."""
+        r = math.isqrt(max((d := abs(self)), 0))
+        if r * r != d:
+            return None
+
+        xs = (math.isqrt((self.a + r) // 2), math.isqrt((self.a - r) // 2))
+        ys = (math.isqrt((self.a - r) // 4), math.isqrt((self.a + r) // 4))
+
+        res = None
+        for x, y in zip(xs, ys):
+            zrt = ZSqrtTwo(x, y)
+            if zrt * zrt == self:
+                res = zrt
+                break
+            art = zrt.adj2()
+            if art * art == self:  # pragma: no cover
+                res = art
+                break
+
+        return res
+
+    def to_omega(self) -> ZOmega:
+        """Convert to the an ring of integers adjoined with omega."""
+        return ZOmega(-self.b, 0, self.b, self.a)
+
+
+class ZOmega:
+    r"""Represents the elements of the ring of integers adjoined with :math:`\omega`, the eight root of unity.
+
+    ..math::
+        \mathbb{Z}[\omega] = \{ a\omega^3 + b\omega^2 + c\omega + d \mid a, b, c, d \in \mathbb{Z} \}
+
+    where :math:`\omega = (1 + i) / \sqrt{2}` is the eighth root of unity.
+
+    Args:
+        a (int): Coefficient of :math:`\omega^3`.
+        b (int): Coefficient of :math:`\omega^2`.
+        c (int): Coefficient of :math:`\omega`.
+        d (int): Constant term.
+    """
+
+    def __init__(self, a: int = 0, b: int = 0, c: int = 0, d: int = 0) -> None:
+        self.a = int(a)
+        self.b = int(b)
+        self.c = int(c)
+        self.d = int(d)
+
+    def __str__(self: ZOmega) -> str:
+        terms = []
+        if self.a:
+            terms.append(f"{self.a} ω^3")
+        if self.b:
+            terms.append(f"{self.b} ω^2")
+        if self.c:
+            terms.append(f"{self.c} ω")
+        if self.d:
+            terms.append(f"{self.d}")
+        return " + ".join(terms) if terms else "0"
+
+    def __repr__(self: ZOmega) -> str:
+        return f"ZOmega(a={self.a}, b={self.b}, c={self.c}, d={self.d})"
+
+    def __complex__(self) -> complex:
+        return complex(self.a * _OMEGA**3 + self.b * _OMEGA**2 + self.c * _OMEGA + self.d)
+
+    def __abs__(self: ZOmega) -> float:
+        a, b, c, d = (self.a, self.b, self.c, self.d)
+        return (a**2 + b**2 + c**2 + d**2) ** 2 - 2 * (a * b + b * c + c * d - d * a) ** 2
+
+    def __mul__(self, other: ZOmega | int | float) -> ZOmega:
+        if isinstance(other, ZOmega):
+            a, b, c, d = (self.a, self.b, self.c, self.d)
+            _a, _b, _c, _d = (other.a, other.b, other.c, other.d)
+            new_a = a * _d + b * _c + c * _b + d * _a
+            new_b = b * _d + c * _c + d * _b - a * _a
+            new_c = c * _d + d * _c - a * _b - b * _a
+            new_d = d * _d - a * _c - b * _b - c * _a
+            return ZOmega(new_a, new_b, new_c, new_d)
+        if isinstance(other, int) or (isinstance(other, float) and other.is_integer()):
+            other = int(other)
+            return ZOmega(self.a * other, self.b * other, self.c * other, self.d * other)
+        raise TypeError(f"Unsupported type {type(other)} for multiplication with ZOmega")
+
+    def __add__(self, other: ZOmega | int | float) -> ZOmega:
+        if isinstance(other, ZOmega):
+            return ZOmega(self.a + other.a, self.b + other.b, self.c + other.c, self.d + other.d)
+        if isinstance(other, int) or (isinstance(other, float) and other.is_integer()):
+            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):
+        return self + (-other)
+
+    __radd__ = __add__
+    __rsub__ = __sub__
+    __rmul__ = __mul__
+
+    def __neg__(self: ZOmega) -> ZOmega:
+        return ZOmega(-self.a, -self.b, -self.c, -self.d)
+
+    def __eq__(self, other: ZOmega) -> bool:
+        if isinstance(other, ZOmega):
+            return (
+                self.a == other.a and self.b == other.b and self.c == other.c and self.d == other.d
+            )
+        if isinstance(other, int):
+            return self.a == 0 and self.b == 0 and self.c == 0 and self.d == other
+        return np.isclose(complex(self), complex(other))
+
+    def __pow__(self, power: int) -> ZOmega:
+        if power < 0:
+            raise ValueError(f"Negative powers {power} are unsupported for ZOmega.")
+        if isinstance(power, float) and not power.is_integer():
+            raise ValueError(f"Non-integer powers {power} are unsupported for ZOmega.")
+        if power == 0:
+            return ZOmega(0, 0, 0, 1)
+        result = self
+        while power > 1:
+            result *= self
+            power -= 1
+        return result
+
+    def __truediv__(self, other: int | float) -> ZOmega:
+        if isinstance(other, int) or (isinstance(other, float) and other.is_integer()):
+            other = int(other)
+            if all(a % other == 0 for a in self.flatten):
+                return ZOmega(self.a / other, self.b / other, self.c / other, self.d / other)
+        raise TypeError(f"Unsupported type {type(other)} for division with ZOmega")
+
+    def __floordiv__(self, other: int) -> ZOmega:
+        if isinstance(other, int) or (isinstance(other, float) and other.is_integer()):
+            other = int(other)
+            return ZOmega(self.a // other, self.b // other, self.c // other, self.d // other)
+        raise TypeError(f"Unsupported type {type(other)} for floor division with ZOmega")
+
+    def __mod__(self, other: ZOmega) -> ZOmega:
+        d = abs(other)
+        n = self * other.conj() * ((other * other.conj()).adj2())
+        return ZOmega(*[(s + d // 2) // d for s in n.flatten]) * other - self
+
+    @classmethod
+    def from_sqrt_pair(cls, alpha: ZSqrtTwo, beta: ZSqrtTwo, shift: ZOmega) -> ZOmega:
+        """Return ``ZOmega`` element as :math:`A + 1j * B + shift`, where :math:`A` and
+        :math:`B` are ``ZSqrtTwo`` elements and ``shift`` is ``ZOmega`` element."""
+        return cls(beta.b - alpha.b, beta.a, beta.b + alpha.b, alpha.a) + shift
+
+    @property
+    def flatten(self: ZOmega) -> List[int]:
+        """Flatten to a list."""
+        return [self.a, self.b, self.c, self.d]
+
+    def conj(self: ZOmega) -> ZOmega:
+        r"""Return the complex conjugate.
+
+        .. math::
+            (a\omega^3 + b\omega^2 + c\omega + d)^{\dagger} = -c\omega^3 + b\omega^2 - a\omega + d
+        """
+        return ZOmega(-self.c, -self.b, -self.a, self.d)
+
+    def adj2(self: ZOmega) -> ZOmega:
+        r"""Return the adjoint, i.e., the root-2 conjugate.
+
+        .. math::
+            (a\omega^3 + b\omega^2 + c\omega + d)^{\bullet} = -a\omega^3 + b\omega^2 - c\omega + d
+        """
+        return ZOmega(-self.a, self.b, -self.c, self.d)
+
+    def norm(self: ZOmega) -> float:
+        """Return the norm squared."""
+        return self * self.conj()
+
+    def parity(self: ZOmega) -> int:
+        """Return the parity indicating structure of real and imaginary parts as a DyadicMatrix element."""
+        return (self.a + self.c) % 2
+
+    def to_sqrt_two(self: ZOmega) -> ZSqrtTwo:
+        """Convert to the ring of integers adjoined with the square root of 2.
+
+        Returns:
+            ZSqrtTwo: The corresponding element in the ZSqrtTwo ring.
+        """
+        if (self.c + self.a) == 0 and self.b == 0:
+            return ZSqrtTwo(self.d, (self.c - self.a) // 2)
+        raise ValueError("Cannot convert ZOmega to ZSqrtTwo.")
+
+    def normalize(self: ZOmega) -> tuple[ZOmega, int]:
+        """Normalize the ZOmega element and return the number of times 2 was factored out."""
+        res, ix = self, 0
+        while (res.a + res.c) % 2 == 0 and (res.b + res.d) % 2 == 0:
+            a = (res.b - res.d) // 2
+            b = (res.a + res.c) // 2
+            c = (res.b + res.d) // 2
+            d = (res.c - res.a) // 2
+            res = ZOmega(a, b, c, d)
+            ix += 1
+        return res, ix
+
+
+class DyadicMatrix:
+    r"""Represents the matrices over the ring :math:`\mathbb{D}[\omega]`,
+    the ring of dyadic fractions adjoined with :math:`\omega`.
+
+    The dyadic fractions :math:`\mathbb{D} = \mathbb{Z}[\frac{1}{2}]` are defined as
+    :math:`\mathbb{D} = \{ a / 2^k \mid a \in \mathbb{Z}, k \in  \{0\} \cup \mathbb{N}\}`. This gives:
+
+    .. math::
+        \mathbb{D}[omega] = \mathbb{Z}[\frac{1}{\sqrt{2}}, i] = \{ a\omega^3 + b\omega^2 + c\omega + d \mid a, b, c, d \in \mathbb{Z}[\frac{1}{\sqrt{2}}] \}
+
+    The `~pennylane.ZOmega` (or :math:`\mathbb{Z}[\omega]) represents a subset of :math:`\mathbb{D}[\omega]`,
+    and therefore can be used to construct the elements of a ``DyadicMatrix``, which is represented as:
+
+    .. math::
+        \frac{1}{\sqrt{2}^k}
+        \begin{pmatrix}
+        a_{00} & a_{01} \\
+        a_{10} & a_{11}
+        \end{pmatrix},
+
+    where :math:`a_{ij} \in \mathbb{D}[\omega]` and :math:`k \in \mathbb{Z}`.
+
+    Args:
+        a (ZOmega): Element at position (0, 0) of the matrix.
+        b (ZOmega): Element at position (0, 1) of the matrix.
+        c (ZOmega): Element at position (1, 0) of the matrix.
+        d (ZOmega): Element at position (1, 1) of the matrix.
+        k (int): Optional integer to scale the matrix by a factor of :math:`1 / \sqrt{2}^k`.
+    """
+
+    # pylint:disable = too-many-positional-arguments, too-many-arguments
+    def __init__(self, a: ZOmega, b: ZOmega, c: ZOmega, d: ZOmega, k: int = 0) -> None:
+        self.a = a
+        self.b = b
+        self.c = c
+        self.d = d
+        self.k = k
+        self.normalize()
+
+    def __str__(self: DyadicMatrix) -> str:
+        return f"[[{self.a}, {self.b}], [{self.c}, {self.d}]]"
+
+    def __repr__(self: DyadicMatrix) -> str:
+        return f"DyadicMatrix(a={self.a}, b={self.b}, c={self.c}, d={self.d}, k={self.k})"
+
+    def __neg__(self: DyadicMatrix) -> DyadicMatrix:
+        return DyadicMatrix(-self.a, -self.b, -self.c, -self.d, self.k)
+
+    def __eq__(self: DyadicMatrix, other: DyadicMatrix) -> bool:
+        """Check if two dyadic matrices are equal."""
+        return (
+            self.a == other.a
+            and self.b == other.b
+            and self.c == other.c
+            and self.d == other.d
+            and self.k == other.k
+        )
+
+    def __mul__(self: DyadicMatrix, other: int | ZOmega) -> DyadicMatrix:
+        """Multiply the matrix by an integer."""
+        if isinstance(other, float) and other.is_integer():
+            other = int(other)
+
+        return DyadicMatrix(
+            self.a * other,
+            self.b * other,
+            self.c * other,
+            self.d * other,
+            self.k,
+        )
+
+    def __add__(self: DyadicMatrix, other: int | float | complex | DyadicMatrix) -> DyadicMatrix:
+        """Add two dyadic matrices."""
+        if isinstance(other, int) or (isinstance(other, float) and other.is_integer()):
+            other = DyadicMatrix(ZOmega(d=int(other)), ZOmega(), ZOmega(), ZOmega(d=int(other)))
+        if isinstance(other, complex) and other.real.is_integer() and other.imag.is_integer():
+            z_omega = ZOmega(b=other.imag, d=other.real)
+            other = DyadicMatrix(z_omega, ZOmega(), ZOmega(), z_omega)
+        if not isinstance(other, DyadicMatrix):
+            raise TypeError(f"Unsupported type {type(other)} for addition with DyadicMatrix")
+        # Ensure k is the maximum of both matrices
+        A, B = (self, other) if self.k >= other.k else (other, self)
+        k_scale, k_parity = int(math.pow(2, (A.k - B.k) // 2)), (A.k - B.k) % 2
+        b_elems = []
+        for b_elem in B.flatten:
+            a, b, c, d = [s * k_scale for s in b_elem.flatten]
+            if k_parity != 0:  # sqrt(2) factor
+                a, b, c, d = [(b - d), (c + a), (b + d), (c - a)]
+            b_elems.append(ZOmega(a, b, c, d))
+
+        return DyadicMatrix(
+            a=A.a + b_elems[0], b=A.b + b_elems[1], c=A.c + b_elems[2], d=A.d + b_elems[3], k=A.k
+        )
+
+    def __matmul__(self: DyadicMatrix, other: DyadicMatrix) -> DyadicMatrix:
+        """Multiply two dyadic matrices."""
+        if not isinstance(other, DyadicMatrix):
+            raise TypeError(
+                f"Unsupported type {type(other)} for matrix multiplication with DyadicMatrix"
+            )
+        return DyadicMatrix(
+            a=self.a * other.a + self.b * other.c,
+            b=self.a * other.b + self.b * other.d,
+            c=self.c * other.a + self.d * other.c,
+            d=self.c * other.b + self.d * other.d,
+            k=self.k + other.k,
+        )
+
+    @property
+    def ndarray(self: DyadicMatrix) -> np.ndarray:
+        """Convert the matrix to a NumPy array."""
+        return (_SQRT2**-self.k) * np.array(
+            [
+                [complex(self.a), complex(self.b)],
+                [complex(self.c), complex(self.d)],
+            ],
+            dtype=np.complex128,
+        )
+
+    @property
+    def flatten(self: DyadicMatrix) -> List[ZOmega]:
+        """Flatten the matrix elements to a list."""
+        return [self.a, self.b, self.c, self.d]
+
+    def conj(self: DyadicMatrix) -> DyadicMatrix:
+        """Return the conjugate of the matrix."""
+        return DyadicMatrix(
+            self.a.conj(),
+            self.b.conj(),
+            self.c.conj(),
+            self.d.conj(),
+            self.k,
+        )
+
+    def adj2(self: DyadicMatrix) -> DyadicMatrix:
+        """Return the root-2 adjoint of the matrix."""
+        return DyadicMatrix(
+            self.a.adj2(),
+            self.b.adj2(),
+            self.c.adj2(),
+            self.d.adj2(),
+            self.k,
+        )
+
+    def mult2k(self: DyadicMatrix, k: int) -> DyadicMatrix:
+        """Multiply the matrix by :math:`2^k`, i.e., an integer power of 2."""
+        if k == 0:
+            return self
+        k_val = min(0, self.k - 2 * k)
+        k_scale = abs(k_val % 2)
+        e_scale = int(math.pow(2, (k_scale - k_val) // 2))
+        return DyadicMatrix(
+            self.a * e_scale,
+            self.b * e_scale,
+            self.c * e_scale,
+            self.d * e_scale,
+            self.k + k_scale,
+        )
+
+    def normalize(self: DyadicMatrix) -> None:
+        """Reduce the k value of the dyadic matrix.
+
+        Example:
+            >>> A = DyadicMatrix(ZOmega(d=2), ZOmega(d=2), ZOmega(d=2), ZOmega(d=2), k = 4)
+            >>> A.normalize()
+            >>> A
+            DyadicMatrix(a=1, b=1, c=1, d=1, k=2)
+        """
+        # a, b, c, d = self.flatten
+        # Factoring 2: Derived using (1 - w^4) [a', b', c', d'] => [a, b, c, d]
+        if all(ZOmega() == s for s in self.flatten):
+            self.k = 0
+            return
+
+        while all(np.allclose([_s % 2 for _s in s.flatten], 0) for s in self.flatten):
+            self.k -= 2
+            self.a, self.b, self.c, self.d = [s // 2 for s in self.flatten]
+
+        # Factoring sqrt(2): Derived using (w - w^3) [a', b', c', d'] => [a, b, c, d]
+        sqrt2_flag = True
+        while sqrt2_flag and self.k > 0:
+            for s in self.flatten:
+                if (s.a + s.c) % 2 or (s.b + s.d) % 2:
+                    sqrt2_flag = False
+                    break
+            if sqrt2_flag:
+                self.k -= 1
+                elements = self.flatten
+                for i in range(4):
+                    a, b, c, d = elements[i].flatten
+                    elements[i] = ZOmega((b - d) // 2, (c + a) // 2, (b + d) // 2, (c - a) // 2)
+                self.a, self.b, self.c, self.d = elements
+
+
+class SO3Matrix:
+    r"""Represents the :math:`SO(3)` matrices over the ring :math:`\mathbb{D}[\sqrt{2}]`,
+    the ring of dyadic integers adjoined with :math:`\sqrt{2}`.
+
+    The `~pennylane.ZSqrtTwo` (or :math:`\mathbb{Z}[\sqrt{2}]) represents a subset of this ring,
+    and can be used to construct its elements. The matrix form is usually represented as:
+
+    .. math::
+
+        \frac{1}{\sqrt{2}^k}
+        \begin{pmatrix}
+        a_{00} & a_{01} & a_{11} \\
+        a_{10} & a_{11} & a_{12} \\
+        a_{20} & a_{21} & a_{22}
+        \end{pmatrix},
+
+    where :math:`a_{ij} \in \mathbb{Z}[\sqrt{2}]` and :math:`k \in \mathbb{Z}`.
+
+    Args:
+        matrix (DyadicMatrix): The :class:`~pennylane.DyadicMatrix` matrix from which the :math:`SO(3)` matrix is derived.
+    """
+
+    def __init__(self, matrix: DyadicMatrix) -> None:
+        """Initialize the SO(3) matrix with a dyadic matrix and an integer k."""
+        self.matrix = matrix
+        self.k = 0
+        self.so3mat = self.from_matrix(matrix)
+        self.normalize()
+
+    def __str__(self: SO3Matrix) -> str:
+        """Return a string representation of the SO(3) matrix."""
+        elements = self.flatten
+        str_repr = "["
+        for i in range(3):
+            str_repr += f"[{elements[i * 3]}, {elements[i * 3 + 1]}, {elements[i * 3 + 2]}], \n"
+        str_repr = str_repr.rstrip(", \n") + "]" + (f" * 1 / √2^{self.k}" if self.k else "")
+        return str_repr
+
+    def __repr__(self: SO3Matrix) -> str:
+        """Return a string representation of the SO(3) matrix."""
+        return f"SO3Matrix(matrix={self.matrix}, k={self.k})"
+
+    def __matmul__(self, other: SO3Matrix) -> SO3Matrix:
+        res = deepcopy(self)
+        res.matrix = self.matrix @ other.matrix
+        res.k = self.k + other.k
+
+        us_self, us_other = self.flatten, other.flatten
+        res.so3mat = [
+            [
+                us_self[0] * us_other[0] + us_self[1] * us_other[3] + us_self[2] * us_other[6],
+                us_self[0] * us_other[1] + us_self[1] * us_other[4] + us_self[2] * us_other[7],
+                us_self[0] * us_other[2] + us_self[1] * us_other[5] + us_self[2] * us_other[8],
+            ],
+            [
+                us_self[3] * us_other[0] + us_self[4] * us_other[3] + us_self[5] * us_other[6],
+                us_self[3] * us_other[1] + us_self[4] * us_other[4] + us_self[5] * us_other[7],
+                us_self[3] * us_other[2] + us_self[4] * us_other[5] + us_self[5] * us_other[8],
+            ],
+            [
+                us_self[6] * us_other[0] + us_self[7] * us_other[3] + us_self[8] * us_other[6],
+                us_self[6] * us_other[1] + us_self[7] * us_other[4] + us_self[8] * us_other[7],
+                us_self[6] * us_other[2] + us_self[7] * us_other[5] + us_self[8] * us_other[8],
+            ],
+        ]
+        res.normalize()
+        return res
+
+    def __eq__(self: SO3Matrix, other: SO3Matrix) -> bool:
+        return self.k == other.k and all(x == y for (x, y) in zip(self.flatten, other.flatten))
+
+    @property
+    def flatten(self: SO3Matrix) -> List[ZOmega]:
+        """Flatten the matrix to a 1D NumPy array."""
+        return [l for row in self.so3mat for l in row]
+
+    @property
+    def ndarray(self: SO3Matrix) -> np.ndarray:
+        """Convert the matrix to a NumPy array."""
+        matrix = np.array(list(map(float, self.flatten)))
+        return (_SQRT2**-self.k) * matrix.reshape(3, 3)
+
+    @property
+    def parity_mat(self: SO3Matrix) -> np.ndarray:
+        """Return the parity of the SO(3) matrix."""
+        return np.array([[x.a % 2 for x in row] for row in self.so3mat], dtype=np.int8)
+
+    @property
+    def parity_vec(self: SO3Matrix) -> np.ndarray:
+        """Return the permutation vector of the SO(3) matrix."""
+        return np.sum(self.parity_mat, axis=1)
+
+    def from_matrix(self, matrix: DyadicMatrix) -> List[List[ZSqrtTwo]]:
+        """Return the SO(3) matrix as a list of lists."""
+        su2_elems, k = matrix.flatten, 2 * matrix.k
+        if any(s.parity for s in su2_elems):
+            z_sqrt2 = [(ZSqrtTwo((s.c - s.a), s.d), ZSqrtTwo((s.c + s.a), s.b)) for s in su2_elems]
+            k += 2
+        else:  # pragma: no cover
+            z_sqrt2 = [
+                (ZSqrtTwo(s.d, (s.c - s.a) // 2), ZSqrtTwo(s.b, (s.c + s.a) // 2))
+                for s in su2_elems
+            ]
+
+        a_, b_, c_, d_ = z_sqrt2
+        so3_mat = [
+            [
+                a_[0] * d_[0] + a_[1] * d_[1] + b_[0] * c_[0] + b_[1] * c_[1],
+                a_[1] * d_[0] + b_[0] * c_[1] - b_[1] * c_[0] - a_[0] * d_[1],
+                a_[0] * c_[0] + a_[1] * c_[1] - b_[0] * d_[0] - b_[1] * d_[1],
+            ],
+            [
+                a_[0] * d_[1] - a_[1] * d_[0] + b_[0] * c_[1] - b_[1] * c_[0],
+                a_[0] * d_[0] + a_[1] * d_[1] - b_[0] * c_[0] - b_[1] * c_[1],
+                a_[0] * c_[1] - a_[1] * c_[0] - b_[0] * d_[1] + b_[1] * d_[0],
+            ],
+            [
+                2 * (a_[0] * b_[0] + a_[1] * b_[1]),
+                2 * (a_[1] * b_[0] - a_[0] * b_[1]),
+                a_[0] ** 2 + a_[1] ** 2 - b_[0] ** 2 - b_[1] ** 2,
+            ],
+        ]
+        self.k = k
+        return so3_mat
+
+    def normalize(self: SO3Matrix) -> None:
+        """Reduce the k value of the SO(3) matrix.
+
+        Example:
+            >>> A = DyadicMatrix(ZOmega(d=2), ZOmega(d=2), ZOmega(d=2), ZOmega(d=2), k = 4) * 2
+            >>> B = SO3Matrix(A @ A)
+            >>> B.normalize()
+            SO3Matrix(matrix=[[1, 1], [1, 1]], k=-6)
+        """
+        elements = self.flatten
+        if all(s.a == 0 and s.b == 0 for s in elements):
+            self.k = 0
+            return
+        # Factoring 2: Derived using (a + b . √2) => 2 . (a//2 + b//2 . √2)
+        while all(s.a % 2 == 0 and s.b % 2 == 0 for s in elements):
+            self.k -= 2
+            elements = [s // 2 for s in elements]
+        # Factoring sqrt(2): Derived using (a + b . √2) => √2 (b + a // 2 . √2)
+        while all(s.a % 2 == 0 for s in elements) and self.k > 0:
+            self.k -= 1
+            elements = [ZSqrtTwo(s.b, s.a // 2) for s in elements]
+        self.so3mat = [elements[i : i + 3] for i in range(0, len(elements), 3)]
diff --git a/pennylane/ops/op_math/decompositions/ross_selinger.py b/pennylane/ops/op_math/decompositions/ross_selinger.py
new file mode 100644
index 000000000..4410ae96f
--- /dev/null
+++ b/pennylane/ops/op_math/decompositions/ross_selinger.py
@@ -0,0 +1,152 @@
+# 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.
+"""Ross-Selinger (arXiv:1403.2975v3) implementation for approximate Pauli-Z rotation gate decomposition."""
+import math
+
+import pennylane as qml
+from pennylane.ops.op_math.decompositions.grid_problems import GridIterator
+from pennylane.ops.op_math.decompositions.norm_solver import _solve_diophantine
+from pennylane.ops.op_math.decompositions.normal_forms import _ma_normal_form
+from pennylane.ops.op_math.decompositions.rings import DyadicMatrix, SO3Matrix, ZOmega, ZSqrtTwo
+from pennylane.queuing import QueuingManager
+
+
+def _domain_correction(theta: float) -> tuple[float, ZOmega]:
+    r"""Return shifts for the angle :math:`\theta` for it to be in the interval :math:`[-\pi/4, \pi/4]` and the corresponding scaling factor for the matrix elements.
+
+    Args:
+        theta (float): The angle to shift.
+
+    Returns:
+        tuple[float, ZOmega]: The domain shift and the scaling factor.
+    """
+    pi_vals = [(2 * i + 1) * math.pi / 4 for i in range(4)]  # pi/4, 3pi/4, 5pi/4, 7pi/4
+
+    sign = 1
+    theta %= 4 * math.pi
+    if theta > 2 * math.pi:
+        sign = -1
+        theta -= math.pi * 4
+    abs_theta = abs(theta)
+
+    if pi_vals[0] <= abs_theta < pi_vals[1]:  # pi/4 <= |theta| < 3pi/4
+        return -sign * math.pi / 2, ZOmega(b=sign)
+
+    if pi_vals[1] <= abs_theta < pi_vals[2]:  # 3pi/4 <= |theta| < 5pi/4
+        return -sign * math.pi, ZOmega(d=-1)
+
+    if pi_vals[2] <= abs_theta < pi_vals[3]:  # 5pi/4 <= |theta| < 7pi/4
+        return -sign * 3 * math.pi / 2, ZOmega(b=-sign)
+
+    return 0.0, ZOmega(d=1)  # -pi/4 <= |theta| < pi/4 / 7pi/4 <= |theta| < 8pi/4
+
+
+def rs_decomposition(op, epsilon, *, max_trials=20):
+    r"""Approximate a phase shift rotation gate in the Clifford+T basis using the `Ross-Selinger algorithm <https://arxiv.org/abs/1403.2975>`_.
+
+    This method implements the Ross-Selinger decomposition algorithm that approximates any arbitrary
+    phase shift rotation gate with :math:`\epsilon > 0` error. The procedure exits when the approximation error
+    becomes less than :math:`\epsilon`, or when ``max_trials`` attempts have been made for solution search.
+    In the latter case, the approximation error could be :math:`\geq \epsilon`.
+
+    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.
+
+    Keyword Args:
+        max_trials (int): The maximum number of attempts to find a solution while performing the grid search according to the the Algorithm 7.6,
+            in the `arXiv:1403.2975v3 <https://arxiv.org/abs/1403.2975>`_. Default is ``20``.
+
+    Returns:
+        list[~pennylane.operation.Operation]: A list of gates in the Clifford+T basis set that approximates the given
+        operation along with a final global phase operation. The operations are in the circuit-order.
+
+    Raises:
+        ValueError: If the given operator is not a :class:`~.RZ` or :class:`~.PhaseShift` gate.
+
+    **Example**
+
+    Suppose one would like to decompose :class:`~.RZ` with rotation angle :math:`\phi = \pi/3`:
+
+    .. code-block:: python3
+
+        import numpy as np
+        import pennylane as qml
+
+        op  = qml.RZ(np.pi/3, wires=0)
+        ops = qml.ops.rs_decomposition(op, epsilon=1e-3)
+
+        # Get the approximate matrix from the ops
+        matrix_rs = qml.prod(*reversed(ops)).matrix()
+
+    When the function is run for a sufficient ``max_trials``, the output gate sequence
+    should implement the same operation approximately, up to a global phase.
+
+    >>> qml.math.allclose(op.matrix(), matrix_rs, atol=1e-3)
+    True
+
+    """
+
+    with QueuingManager.stop_recording():
+
+        # Check for length of wires in the operation
+        if not isinstance(op, (qml.RZ, qml.PhaseShift)):
+            raise ValueError(f"Operator must be a RZ or PhaseShift gate, got {op}")
+
+        # 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_trials)
+
+        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)) is not None:
+                u, t, k = u_sol * scale, t_sol * scale, k_val
+                break
+
+        # 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)
+
+        # Remove inverses if any in the decomposition and handle trivial case
+        new_tape = qml.tape.QuantumScript(decomposition)
+
+    # Map the wires to that of the operation and queue
+    if queuing := QueuingManager.recording():
+        QueuingManager.remove(op)
+
+    if (op_wire := op.wires[0]) != 0:
+        [new_tape], _ = qml.map_wires(new_tape, wire_map={0: op_wire}, queue=True)
+    else:
+        if queuing:
+            _ = [qml.apply(op) for op in new_tape.operations]
+
+    # TODO: Add 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))
+
+    # Return the gates from the mapped tape and global phase
+    return new_tape.operations + [global_phase]
diff --git a/pennylane/transforms/decompositions/clifford_t_transform.py b/pennylane/transforms/decompositions/clifford_t_transform.py
index 5b6ba4cfc..b3796e7eb 100644
--- a/pennylane/transforms/decompositions/clifford_t_transform.py
+++ b/pennylane/transforms/decompositions/clifford_t_transform.py
@@ -20,6 +20,7 @@ from itertools import product
 
 import pennylane as qml
 from pennylane.ops import Adjoint
+from pennylane.ops.op_math.decompositions.ross_selinger import rs_decomposition
 from pennylane.ops.op_math.decompositions.solovay_kitaev import sk_decomposition
 from pennylane.queuing import QueuingManager
 from pennylane.tape import QuantumScript, QuantumScriptBatch
@@ -360,9 +361,13 @@ class _CachedCallable:
                 self.decompose_fn = lru_cache(maxsize=cache_size)(
                     partial(sk_decomposition, epsilon=epsilon, **method_kwargs)
                 )
+            case "rs":
+                self.decompose_fn = lru_cache(maxsize=cache_size)(
+                    partial(rs_decomposition, epsilon=epsilon, **method_kwargs)
+                )
             case _:
                 raise NotImplementedError(
-                    f"Currently we only support Solovay-Kitaev ('sk') decomposition, got {method}"
+                    f"Currently we only support Solovay-Kitaev ('sk') and Ross-Selinger ('rs') decompositions, got {method}"
                 )
 
         self.method = method
@@ -414,11 +419,14 @@ def clifford_t_decomposition(
     Then, the leftover single qubit :class:`~.RZ` operations are approximated in the Clifford+T basis with
     :math:`\epsilon > 0` error. By default, we use the Solovay-Kitaev algorithm described in
     `Dawson and Nielsen (2005) <https://arxiv.org/abs/quant-ph/0505030>`_ for this.
+    Alternatively, the Ross-Selinger algorithm described in `Ross and Selinger (2016) <https://arxiv.org/abs/1403.2975v3>`_
+    can be used by setting the ``method`` to ``"rs"``.
 
     Args:
         tape (QNode or QuantumTape or Callable): The quantum circuit to be decomposed.
         epsilon (float): The maximum permissible operator norm error of the complete circuit decomposition. Defaults to ``0.0001``.
-        method (str): Method to be used for Clifford+T decomposition. Default value is ``"sk"`` for Solovay-Kitaev.
+        method (str): Method to be used for Clifford+T decomposition. Default value is ``"sk"`` for Solovay-Kitaev. Alternatively,
+            the Ross-Selinger algorithm can be used with ``"rs"``.
         cache_size (int): The size of the cache built for the decomposition function based on the angle. Defaults to ``1000``.
         **method_kwargs: Keyword argument to pass options for the ``method`` used for decompositions.
 
@@ -432,11 +440,15 @@ def clifford_t_decomposition(
         **max_depth** (int), **basis_set** (list[str]), **basis_length** (int) -- arguments for the ``"sk"`` method,
         where the decomposition is performed using the :func:`~.sk_decomposition` method.
 
+    - Ross-Selinger decomposition --
+        **max_trials** (int) -- argument for the ``"rs"`` method, where the decomposition is performed using the
+        :func:`~.rs_decomposition` method.
+
     Raises:
         ValueError: If a gate operation does not have a decomposition when required.
         NotImplementedError: If chosen decomposition ``method`` is not supported.
 
-    .. seealso:: :func:`~.sk_decomposition` for Solovay-Kitaev decomposition.
+    .. seealso:: :func:`~.rs_decomposition` and :func:`~.sk_decomposition` for Ross-Selinger and Solovay-Kitaev decomposition methods, respectively.
 
     **Example**
 


Test patch:
diff --git a/tests/ops/op_math/decompositions/test_grid_problems.py b/tests/ops/op_math/decompositions/test_grid_problems.py
new file mode 100644
index 000000000..1fcd6be7d
--- /dev/null
+++ b/tests/ops/op_math/decompositions/test_grid_problems.py
@@ -0,0 +1,202 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for the grid problems for Clifford+T decomposition."""
+
+# pylint: disable = too-few-public-methods
+
+import inspect
+import math
+
+import pytest
+
+from pennylane import math as pmath
+from pennylane.ops.op_math.decompositions.grid_problems import (
+    Ellipse,
+    EllipseState,
+    GridIterator,
+    GridOp,
+)
+from pennylane.ops.op_math.decompositions.rings import ZOmega
+
+
+class TestEllipse:
+    """Tests for the Ellipse class."""
+
+    def test_init_and_repr(self):
+        """Test that Ellipse initializes correctly and has correct representations."""
+        ellipse = Ellipse((1, 0, 1), (4, 5), (6, 7))
+        assert ellipse.a == 1 and ellipse.b == 0 and ellipse.d == 1
+        assert ellipse.p == (4, 5) and ellipse.axes == (6, 7)
+        assert repr(ellipse) == "Ellipse(a=1, b=0, d=1, p=(4, 5), axes=(6, 7))"
+        assert str(ellipse) == "Ellipse(a=1, b=0, d=1, p=(4, 5), axes=(6, 7))"
+        assert ellipse.determinant == 1 and ellipse.discriminant == 0
+        assert ellipse.bounding_box() == (-1, 1, -1, 1)
+        assert ellipse.offset(1) == Ellipse((1, 0, 1), p=(5, 6), axes=(6, 7))
+        assert ellipse.b_from_uprightness(0.2) == math.sqrt((math.pi / 0.8) ** 2 - 1)
+        assert ellipse.positive_semi_definite
+        assert ellipse.x_points(4) == (4.0, 4.0)
+        assert ellipse.y_points(5) == (5.0, 5.0)
+        assert pmath.allclose(ellipse.uprightness, 0.7853981633974483)
+
+    def test_raise_error(self):
+        """Test that Ellipse raises an error when the ellipse is not valid."""
+        ellipse = Ellipse((1, 0, 1), (4, 5), (6, 7))
+        with pytest.raises(ValueError, match="is outside the ellipse"):
+            ellipse.x_points(24)
+        with pytest.raises(ValueError, match="is outside the ellipse"):
+            ellipse.y_points(24)
+
+
+class TestEllipseState:
+    """Tests for the EllipseState class."""
+
+    def test_init_and_repr(self):
+        """Test that EllipseState initializes correctly and has correct representations."""
+        e1 = Ellipse((1, 0, 1), (4, 5), (6, 7))
+        e2 = Ellipse((2, 1, 2), (2, 3), (4, 5))
+        state = EllipseState(e1, e2)
+        assert state.e1 == e1
+        assert state.e2 == e2
+        assert repr(state) == f"EllipseState(e1={e1}, e2={e2})"
+        assert state.skew == 1.0 and state.bias == 0.0
+
+
+class TestGridOp:
+    """Tests for the GridOp class."""
+
+    def test_init_and_repr(self):
+        """Test that GridOp initializes correctly and has correct representations."""
+        grid_op = GridOp.from_string("I")
+        assert repr(grid_op) == "GridOp(a=(1, 0), b=(0, 0), c=(0, 0), d=(1, 0))"
+        assert grid_op**3 == grid_op
+        assert GridOp.from_string("R") ** 3 == GridOp(a=(0, -1), b=(0, -1), c=(0, 1), d=(0, -1))
+        assert GridOp.from_string("K").is_special
+        assert GridOp(a=(0, 1), b=(0, -1), c=(1, 1), d=(-1, 1)).transpose() == GridOp(
+            a=(-1, 1), b=(0, -1), c=(1, 1), d=(0, 1)
+        )
+
+        grid_op = GridOp.from_string("B")
+        e1 = Ellipse((1, 0, 1), (4, 5), (6, 7))
+        e2 = Ellipse((2, 1, 2), (2, 3), (4, 5))
+        state = EllipseState(e1, e2)
+        state1 = grid_op.apply_to_state(state)
+        state2 = grid_op.inverse().apply_to_state(state1)
+        assert state.e1 == state2.e1 and state.e2 == state2.e2
+        assert state1.e1 == grid_op.apply_to_ellipse(e1)
+
+    def test_raise_error(self):
+        """Test that GridOp raises an error when an operation is not permitted."""
+        with pytest.raises(TypeError, match="Cannot multiply GridOp with"):
+            _ = GridOp.from_string("B") * 3
+
+        with pytest.raises(ValueError, match="Grid operator needs to be special"):
+            _ = GridOp(a=(1, 1), b=(-2, -1), c=(1, 1), d=(-1, 4), check_valid=False).inverse()
+
+
+class TestGridIterator:
+    """Tests for the GridIterator class."""
+
+    # pylint: disable = too-many-arguments
+    @pytest.mark.parametrize(
+        "x0, x1, y0, y1, num",
+        [
+            (8.9, 9.5, -21, -18, 2),
+            (246.023423, 248.5823575862261, 778, 779.0106829464769, 3),
+            (13734300, 13734500, -13874089.232, -13874089.181, 6),
+        ],
+    )
+    def test_one_dim_problem(self, x0, x1, y0, y1, num):
+        """Test that the one dimensional grid problem is solved correctly."""
+        gitr = GridIterator()
+        sols = gitr.solve_one_dim_problem(x0, x1, y0, y1)
+        assert inspect.isgenerator(sols)
+        assert gitr.bbox_grid_points((x0, x1, y0, y1)) == num
+
+        ix = 0
+        for sol in sols:
+            ix, s1, s2 = ix + 1, float(sol), float(sol.adj2())
+            assert x0 <= s1 <= x1
+            assert y0 <= s2 <= y1
+        assert 0 < ix <= num
+
+    @pytest.mark.parametrize(
+        "bbox1, bbox2, res",
+        [
+            ((5, 6, 4, 5), (2, 3, -1, 0), (1, 2, 3, 4)),
+            ((-4, -3.8, 2.2, 2.4), (-9, -8, 2.5, 4.2), (-2, 3, 1, -6)),
+        ],
+    )
+    def test_upright_problem(self, bbox1, bbox2, res):
+        """Test that the upright grid problem is solved correctly."""
+        D, num_b, shifts = (1, 0, 1), [0, 0], [ZOmega(), ZOmega(c=1)]
+        bbox3 = tuple(bb_ - 1 / math.sqrt(2) for bb_ in bbox1)
+        bbox4 = tuple(bb_ + 1 / math.sqrt(2) for bb_ in bbox2)
+        state = EllipseState(Ellipse(D), Ellipse(D))
+
+        gitr = GridIterator()
+        sols1 = gitr.solve_upright_problem(state, bbox1, bbox2, num_b, shifts[0])
+        sols2 = gitr.solve_upright_problem(state, bbox3, bbox4, num_b, shifts[1])
+        assert inspect.isgenerator(sols1) and inspect.isgenerator(sols2)
+        assert ZOmega(*res) in list(sols1) + list(sols2)
+
+    @pytest.mark.parametrize(
+        "e1, e2, res",
+        [
+            (
+                Ellipse((3, 0.5, 1.0), (-1.7, 13.95)),
+                Ellipse((3, 0.3, 0.3), (-12.3, -7.9)),
+                (4, 3, 12, -7),
+            ),
+            (
+                Ellipse((2, 0.6, 0.9), (-1.8, 14.93)),
+                Ellipse((2, 0.4, 0.2), (-11.3, -6.9)),
+                (4, 4, 11, -6),
+            ),
+        ],
+    )
+    def test_two_dim_problem(self, e1, e2, res):
+        """Test that the two dimensional grid problem is solved correctly."""
+        state = EllipseState(e1, e2)
+        gitr = GridIterator()
+        sols = gitr.solve_two_dim_problem(state)
+        assert inspect.isgenerator(sols)
+        assert ZOmega(*res) in list(sols)
+
+    @pytest.mark.parametrize(
+        "theta, epsilon",
+        [
+            (0.0, 1e-3),
+            (math.pi / 4, 1e-4),
+            (math.pi / 8, 1e-4),
+            (math.pi / 3, 1e-3),
+            (math.pi / 5, 1e-5),
+            (0.392125483789636, 1e-6),
+            (0.6789684841313233, 1e-3),
+            (0.056202026824044335, 1e-5),
+            (0.21375826964510297, 1e-4),
+            (0.5549739238125396, 1e-6),
+            (-0.454645364564563, 1e-3),
+            (-0.5549739238125396, 1e-2),
+        ],
+    )
+    def test_grid_iterator(self, theta, epsilon):
+        """Test that the two dimensional grid problem is solved correctly."""
+        grid_sols = GridIterator(theta=theta, epsilon=epsilon)
+        assert repr(grid_sols) == f"GridIterator(theta={theta}, epsilon={epsilon}, max_trials=20)"
+        assert hasattr(grid_sols, "__iter__")
+        u_sol, k_sol = next(iter(grid_sols))
+        assert u_sol is not None
+        u_sol_complex = complex(u_sol) / math.sqrt(2) ** k_sol
+        assert pmath.allclose(u_sol_complex.real, pmath.cos(theta), atol=epsilon)
+        assert pmath.allclose(u_sol_complex.imag, pmath.sin(theta), atol=epsilon)
diff --git a/tests/ops/op_math/decompositions/test_norm_solvers.py b/tests/ops/op_math/decompositions/test_norm_solvers.py
new file mode 100644
index 000000000..0d950877e
--- /dev/null
+++ b/tests/ops/op_math/decompositions/test_norm_solvers.py
@@ -0,0 +1,160 @@
+# 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.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for the norm solver functions for Clifford+T decomposition."""
+
+from functools import reduce
+
+import pytest
+from flaky import flaky
+
+from pennylane.ops.op_math.decompositions.norm_solver import (
+    _factorize_prime_zomega,
+    _factorize_prime_zsqrt_two,
+    _integer_factorize,
+    _primality_test,
+    _prime_factorize,
+    _solve_diophantine,
+    _sqrt_modulo_p,
+)
+from pennylane.ops.op_math.decompositions.rings import ZOmega, ZSqrtTwo
+
+
+class TestFactorization:
+    """Tests for the factorization functions."""
+
+    @pytest.mark.parametrize(
+        "num, expected",
+        [
+            (28, None),  # 28 = 2^2 * 7 (7 is not included)
+            (1, []),
+            (0, []),
+            (100, [2, 2, 5, 5]),
+            (60, [2, 2, 3, 5]),
+            (53, [53]),  # 53 is prime
+        ],
+    )
+    def test_prime_factorize(self, num, expected):
+        """Test prime factorization."""
+        assert _prime_factorize(num) == expected
+
+    @flaky(max_runs=5)
+    @pytest.mark.parametrize(
+        "num, expected",
+        [
+            (28, [2, 4, 7]),
+            (1, [None]),
+            (0, [None]),
+            (100, [2, 4, 5, 10, 20, 50]),
+            (60, [2, 3, 4, 5, 15, 20, 30]),
+            (47, [None]),  # 47 is prime
+        ],
+    )
+    def test_integer_factorize(self, num, expected):
+        """Test integer factorization."""
+        assert _integer_factorize(num) in expected
+
+    @pytest.mark.parametrize(
+        "num, expected",
+        [
+            (2, [ZSqrtTwo(0, 1), ZSqrtTwo(0, 1)]),
+            (3, [ZSqrtTwo(3, 0)]),
+            (7, [ZSqrtTwo(3, 1), ZSqrtTwo(3, -1)]),
+            (29, [ZSqrtTwo(29, 0)]),
+            (47, [ZSqrtTwo(7, 1), ZSqrtTwo(7, -1)]),
+            (15, None),
+        ],
+    )
+    def test_factorize_prime_zsqrt_two(self, num, expected):
+        """Test factorization of prime Z-sqrt(2)."""
+        assert _factorize_prime_zsqrt_two(num) == expected
+        if expected is not None:
+            assert reduce(lambda x, y: x * y, expected, ZSqrtTwo(1, 0)) == num
+
+    @pytest.mark.parametrize(
+        "num, expected",
+        [
+            (3, ZOmega(d=3)),
+            (27, None),
+            (5, ZOmega(b=-2, d=-1)),
+            (7, None),
+            (11, ZOmega(d=11)),
+            (13, ZOmega(b=-2, d=-3)),
+        ],
+    )
+    def test_factorize_prime_zomega(self, num, expected):
+        """Test factorization of prime Z-omega."""
+        zsqrt_two = _factorize_prime_zsqrt_two(num)[-1]
+        assert _factorize_prime_zomega(zsqrt_two, num) == expected
+        if expected is not None:
+            assert num in (expected, expected * expected.conj())
+
+    @pytest.mark.parametrize(
+        "num, expected",
+        [
+            (2, True),
+            (4, False),
+            (5, True),
+            (29, True),
+            (561, False),  # 561 is not prime (3*11*17)
+            (1729, False),  # 1729 is not prime (1*7*13*19)
+            (7901, True),  # 7901 is prime
+            (41041, False),  # 41041 is not prime (7*11*13*41)
+            (101 * 431, False),  # 101 and 431 are prime, but product is not
+        ],
+    )
+    def test_primality_test(self, num, expected):
+        """Test primality test."""
+        assert _primality_test(num) == expected
+
+    @pytest.mark.parametrize(
+        "nums, expected",
+        [
+            ((3, 2), 1),
+            ((0, 1), 0),
+            ((4, 5), 3),
+            ((9, 11), 3),
+            ((16, 17), 4),
+            ((25, 29), 24),
+            ((-2, 7), None),
+            ((-1, 6), None),  # 6 is not prime
+            ((56, 101), 37),
+            ((3, 4), None),  # 4 is not prime
+        ],
+    )
+    def test_sqrt_modulo_p(self, nums, expected):
+        """Test square root modulo p."""
+        assert _sqrt_modulo_p(*nums) == expected
+        if expected is not None and nums[1] != 2:
+            assert expected**2 % nums[1] == nums[0]
+
+    @pytest.mark.parametrize(
+        "num, expected",
+        [
+            (ZSqrtTwo(0, 0), ZOmega()),
+            (ZSqrtTwo(0, 1), None),
+            (ZSqrtTwo(2, 1), ZOmega(a=0, b=0, c=1, d=1)),
+            (ZSqrtTwo(2, -1), ZOmega(a=1, b=-1, c=0, d=0)),
+            (ZSqrtTwo(7, 0), None),
+            (ZSqrtTwo(23, 0), None),
+            (ZSqrtTwo(7, 2), ZOmega(a=1, b=1, c=1, d=2)),
+            (ZSqrtTwo(17, 0), None),
+            (ZSqrtTwo(5, 2), ZOmega(a=-2, b=-1, c=0, d=0)),
+            (ZSqrtTwo(13, 6), ZOmega(a=0, b=2, c=3, d=0)),
+        ],
+    )
+    def test_solve_diophantine(self, num, expected):
+        """Test `solve_diophantine` solves diophantine equation."""
+        assert _solve_diophantine(num) == expected
+        if expected is not None:
+            assert (expected.conj() * expected).to_sqrt_two() == num
diff --git a/tests/ops/op_math/decompositions/test_normal_forms.py b/tests/ops/op_math/decompositions/test_normal_forms.py
new file mode 100644
index 000000000..d4e4e94d0
--- /dev/null
+++ b/tests/ops/op_math/decompositions/test_normal_forms.py
@@ -0,0 +1,135 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for the normal forms for Clifford+T decomposition."""
+
+import math
+
+import pytest
+
+import pennylane as qml
+from pennylane.ops.op_math.decompositions.normal_forms import (
+    _clifford_group_to_SO3,
+    _ma_normal_form,
+    _parity_transforms,
+)
+from pennylane.ops.op_math.decompositions.rings import DyadicMatrix, SO3Matrix, ZOmega
+from pennylane.ops.op_math.decompositions.solovay_kitaev import (
+    _quaternion_transform,
+    _SU2_transform,
+)
+
+
+class TestNormalForms:
+    """Tests for the Matsumoto-Amano normal forms."""
+
+    def test_clifford_group_to_SO3(self):
+        """Test that the Clifford group elements are correctly mapped to SO(3) matrices."""
+        clifford_matrices = _clifford_group_to_SO3()
+        assert isinstance(clifford_matrices, dict)
+        assert len(clifford_matrices) == 24
+
+        for gates, (so3mat, phase) 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)
+            assert qml.math.allclose(
+                so3mat.ndarray,
+                qml.math.array(
+                    [
+                        [1 - 2 * (y**2 + z**2), 2 * (x * y - w * z), 2 * (x * z + w * y)],
+                        [2 * (x * y + w * z), 1 - 2 * (x**2 + z**2), 2 * (y * z - w * x)],
+                        [2 * (x * z - w * y), 2 * (y * z + w * x), 1 - 2 * (x**2 + y**2)],
+                    ]
+                ),
+            ), "SO(3) matrix does not match expected form"
+
+    def test_parity_transforms(self):
+        """Test that the parity transforms are correctly defined."""
+        parity_transforms = _parity_transforms()
+        assert isinstance(parity_transforms, dict)
+        assert len(parity_transforms) == 4
+
+        parity_vecs = [(1, 1, 1), (2, 2, 0), (0, 2, 2), (2, 0, 2)]
+        for ix, (parity_vec, (so3mat, gates, phase)) in enumerate(parity_transforms.items()):
+            gate = qml.prod(*gates)
+            assert isinstance(so3mat, SO3Matrix), "All transform should have SO3Matrix instances"
+            assert isinstance(
+                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])
+
+            su2mat = _SU2_transform(qml.matrix(qml.adjoint(gate)))[0]
+            w, x, y, z = _quaternion_transform(su2mat)
+            so3mat2 = qml.math.array(
+                [
+                    [1 - 2 * (y**2 + z**2), 2 * (x * y - w * z), 2 * (x * z + w * y)],
+                    [2 * (x * y + w * z), 1 - 2 * (x**2 + z**2), 2 * (y * z - w * x)],
+                    [2 * (x * z - w * y), 2 * (y * z + w * x), 1 - 2 * (x**2 + y**2)],
+                ]
+            )
+            assert qml.math.allclose(
+                so3mat.ndarray, so3mat2
+            ), "SO(3) matrix does not match expected form"
+            assert qml.math.allclose(parity_vecs[ix], parity_vec)
+
+    @pytest.mark.parametrize("klen", [1, 2, 5, 6])
+    def test_ma_normal_form(self, klen):
+        """Test the Matsumoto-Amano normal form decomposition."""
+        clifford_elements = _clifford_group_to_SO3()
+        a = int(klen % 2)
+        b = qml.math.random.randint(2, size=klen)
+        c = qml.math.random.randint(len(clifford_elements))
+
+        if a:
+            so3mat = SO3Matrix(DyadicMatrix(ZOmega(d=1), ZOmega(), ZOmega(), ZOmega(c=1)))
+            so3rep = qml.T(0)
+        else:
+            so3mat = SO3Matrix(DyadicMatrix(ZOmega(d=1), ZOmega(), ZOmega(), ZOmega(d=1)))
+            so3rep = qml.I(0)
+
+        for b_ in b:
+            if b_:
+                so3mat @= SO3Matrix(
+                    DyadicMatrix(ZOmega(d=1), ZOmega(c=1), ZOmega(b=1), ZOmega(a=-1), k=1)
+                )
+                for op_ in (qml.S(0), qml.H(0), qml.T(0)):
+                    so3rep @= op_
+            else:
+                so3mat @= SO3Matrix(
+                    DyadicMatrix(ZOmega(d=1), ZOmega(c=1), ZOmega(d=1), ZOmega(c=-1), k=1)
+                )
+                for op_ in (qml.H(0), qml.T(0)):
+                    so3rep @= op_
+
+        cl_list = list(clifford_elements.keys())
+        so3mat @= clifford_elements[cl_list[c]][0]
+        for op_ in cl_list[c]:
+            so3rep @= op_
+
+        (t_bit, rep_bits, c_bit, phase) = _ma_normal_form(so3mat, compressed=True)
+
+        assert t_bit == a, "T bit does not match expected value"
+        assert (rep_bits == b).all(), "Representation bits do not match expected values"
+        assert c_bit == c, "Clifford bit does not match expected value"
+        assert isinstance(phase, float), "Phase should be a float"
+
+        decomposition, phase = _ma_normal_form(so3mat, compressed=False)
+        assert qml.equal(
+            qml.prod(*decomposition), qml.prod(*so3rep[not a :])
+        ), "Decomposition does not match expected operator"
diff --git a/tests/ops/op_math/decompositions/test_rings.py b/tests/ops/op_math/decompositions/test_rings.py
new file mode 100644
index 000000000..050ee5b67
--- /dev/null
+++ b/tests/ops/op_math/decompositions/test_rings.py
@@ -0,0 +1,228 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for the algebraic prerequisites such as rings and matrices for Clifford+T decomposition."""
+
+import numpy as np
+import pytest
+
+from pennylane.ops.op_math.decompositions.rings import DyadicMatrix, SO3Matrix, ZOmega, ZSqrtTwo
+
+# pylint: disable=too-few-public-methods
+
+
+class TestZSqrtTwo:
+    """Tests for the ZSqrtTwo class."""
+
+    @pytest.mark.parametrize(
+        "a, b",
+        [(0, 0), (2, 1), (3, -2), (-4, 0), (0, 5)],
+    )
+    def test_init_and_repr(self, a, b):
+        """Test that ZSqrtTwo initializes correctly and has correct representations."""
+        z_sqrt_two = ZSqrtTwo(a, b)
+        assert z_sqrt_two.a == a
+        assert z_sqrt_two.b == b
+
+        assert repr(z_sqrt_two) == f"ZSqrtTwo(a={a}, b={b})"
+        # pylint:disable=condition-evals-to-constant
+        assert (
+            str(z_sqrt_two)
+            == (f"{a}" if a else "" + " + " if a and b else "" + f"{b}√2" if b else "")
+            or "0"
+        )
+
+    def test_arithmetic_operations(self):
+        """Test arithmetic operations on ZSqrtTwo."""
+        z1 = ZSqrtTwo(1, 2)
+        z2 = ZSqrtTwo(3, 4)
+
+        assert z1 + z2 == ZSqrtTwo(4, 6) == float(z1) + float(z2)
+        assert z2 + 10 == ZSqrtTwo(13, 4) == float(z2) + 10
+        assert z2 + 2.0 == ZSqrtTwo(5, 4) == float(z2) + 2.0
+        assert z1 - z2 == ZSqrtTwo(-2, -2) == float(z1) - float(z2)
+        assert z1 * z2 == ZSqrtTwo(19, 10) == float(z1) * float(z2)
+        assert z2 * 10 == ZSqrtTwo(30, 40) == float(z2) * 10
+        assert z2 * 2.0 == ZSqrtTwo(6, 8) == float(z2) * 2.0
+        assert ZSqrtTwo(14, 7) / z1 == ZSqrtTwo(2, 3)  # Computed manually
+        assert z1**0 == ZSqrtTwo(1, 0) == float(z1) ** 0
+        assert z1**2 == ZSqrtTwo(9, 4) == float(z1) ** 2
+        assert ZSqrtTwo(5, 0) == 5
+        assert z1.sqrt() is None  # Not all ZSqrtTwo instances have a square root
+        assert z2.sqrt() is None  # Not all ZSqrtTwo instances have a square root
+        assert ZSqrtTwo(30, 6).sqrt() is None
+        assert (z1**2).sqrt() == z1
+        assert z1.conj() == ZSqrtTwo(1, 2)
+        assert z2.adj2() == ZSqrtTwo(3, -4) and z1.adj2() == ZSqrtTwo(1, -2)
+        assert z2 // 2 == ZSqrtTwo(1, 2)
+        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)
+
+    def test_arithmetic_errors(self):
+        """Test that arithmetic operations raise errors for invalid types."""
+        z1 = ZSqrtTwo(1, 2)
+        z2 = complex(3, 4)
+
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = z1 + z2
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = z1 * z2
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = z1 / z2
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = z1 // z2
+        with pytest.raises(ValueError, match="Negative powers"):
+            _ = z1**-2
+        with pytest.raises(ValueError, match="Non-integer powers"):
+            _ = z1**1.5
+
+
+class TestZOmega:
+    """Tests for the ZOmega class."""
+
+    @pytest.mark.parametrize(
+        "a, b, c, d",
+        [
+            (0, 0, 0, 0),
+            (1, 2, 3, 4),
+            (-1, -2, -3, -4),
+            (5, -6, 7, -8),
+        ],
+    )
+    def test_init_and_repr(self, a, b, c, d):
+        """Test that ZOmega initializes correctly and has correct representations."""
+        z_omega = ZOmega(a, b, c, d)
+        assert z_omega.a == a
+        assert z_omega.b == b
+        assert z_omega.c == c
+        assert z_omega.d == d
+        string_repr = [
+            f"{coeff} {var}".strip()
+            for coeff, var in zip([a, b, c, d], ["ω^3", "ω^2", "ω", ""])
+            if coeff
+        ]
+        assert str(z_omega) == " + ".join(string_repr) if string_repr else "0"
+        assert repr(z_omega) == f"ZOmega(a={a}, b={b}, c={c}, d={d})"
+
+    def test_arithmetic_operations(self):
+        """Test arithmetic operations on ZOmega."""
+        z1 = ZOmega(1, 2, 3, 4)
+        z2 = ZOmega(5, 6, 7, 8)
+
+        assert z1 + z2 == ZOmega(6, 8, 10, 12) == complex(z1) + complex(z2)
+        assert z2 + 10 == ZOmega(5, 6, 7, 18) == complex(z2) + 10
+        assert z2 + 2.0 == ZOmega(5, 6, 7, 10) == complex(z2) + 2.0
+        assert z1 - z2 == ZOmega(-4, -4, -4, -4) == complex(z1) - complex(z2)
+        assert z1 * z2 == ZOmega(60, 56, 36, -2) == complex(z1) * complex(z2)
+        assert z2 * 10 == ZOmega(50, 60, 70, 80) == complex(z2) * 10
+        assert z2 * 2.0 == ZOmega(10, 12, 14, 16) == complex(z2) * 2.0
+        assert z1**2 == ZOmega(20, 24, 20, 6) == complex(z1) ** 2
+        assert z2**0 == ZOmega(0, 0, 0, 1) == complex(z2) ** 0
+        assert z1.conj() == ZOmega(-3, -2, -1, 4) == complex(z1).conjugate()
+        assert z2.adj2() == ZOmega(-5, 6, -7, 8) and z1.adj2() == ZOmega(-1, 2, -3, 4)
+        assert (z2 * 4) / 2 == ZOmega(5 * 2, 6 * 2, 7 * 2, 8 * 2)
+        assert z2 // 2 == ZOmega(5 // 2, 6 // 2, 7 // 2, 8 // 2)
+        assert abs(z2) == 14788 and abs(z1) == 388  # Computed manually
+        assert ZOmega(0, 0, 0, 5) == 5
+        assert z1.parity() == 0 and z2.parity() == 0  # Both z1 and z2 have even parity
+        assert z2.norm() == abs(complex(z2)) * abs(complex(z2.conj()))
+        assert (z1 - ZOmega(a=2, b=2, c=2)).to_sqrt_two() == ZSqrtTwo(a=4, b=1)
+
+    def test_arithmetic_errors(self):
+        """Test that arithmetic operations raise errors for invalid types."""
+        z1 = ZOmega(1, 2, 3, 4)
+        z2 = "Dyadic(1, 2)"
+
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = z1 + z2
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = z1 * z2
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = z1 / z2
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = z1 // z2
+        with pytest.raises(ValueError, match="Negative powers"):
+            _ = z1**-2
+        with pytest.raises(ValueError, match="Cannot convert ZOmega to ZSqrtTwo"):
+            _ = z1.to_sqrt_two()
+        with pytest.raises(ValueError, match="Non-integer powers"):
+            _ = z1**1.5
+
+
+class TestDyadicMatrix:
+    """Tests for the DyadicMatrix class."""
+
+    def test_init_and_repr(self):
+        """Test that DyadicMatrix initializes correctly and has correct representations."""
+        z1 = ZOmega(1, 2, 3, 4)
+        z2 = ZOmega(5, 6, 7, 8)
+        dyadic_matrix = DyadicMatrix(z1, z2, z1, z2)
+
+        assert repr(dyadic_matrix) == f"DyadicMatrix(a={z1}, b={z2}, c={z1}, d={z2}, k=0)"
+        assert str(dyadic_matrix) == f"[[{z1}, {z2}], [{z1}, {z2}]]"
+        assert dyadic_matrix * 2 == DyadicMatrix(2 * z1, 2 * z2, 2 * z1, 2 * z2)
+        assert np.allclose(
+            dyadic_matrix.ndarray,
+            np.array([[complex(z1), complex(z2)], [complex(z1), complex(z2)]]),
+        )
+        assert dyadic_matrix.conj() == DyadicMatrix(z1.conj(), z2.conj(), z1.conj(), z2.conj())
+        assert dyadic_matrix.adj2() == DyadicMatrix(z1.adj2(), z2.adj2(), z1.adj2(), z2.adj2())
+        assert dyadic_matrix * 2 == dyadic_matrix.mult2k(k=1)
+        assert dyadic_matrix * 1.0 == dyadic_matrix.mult2k(k=0)
+        assert dyadic_matrix + 2.0 == DyadicMatrix(z1, z2, z1, z2) + (2 + 0j) == dyadic_matrix + 2
+
+        z3 = DyadicMatrix(z1, z2, z1, z2, k=2)
+        z4, z5 = ZOmega(a=-3, b=6, c=9, d=3), ZOmega(a=-3, b=18, c=21, d=3)
+        assert dyadic_matrix + z3 == DyadicMatrix(z4, z5, z4, z5, k=1)
+
+    def test_arithmetic_errors(self):
+        """Test that arithmetic operations raise errors for invalid types."""
+        z1 = ZOmega(1, 2, 3, 4)
+        z2 = "Dyadic(1, 2)"
+        dyadic_matrix = DyadicMatrix(z1, z2, z1, z2)
+
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = dyadic_matrix + z2
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = dyadic_matrix * z2
+        with pytest.raises(TypeError, match="Unsupported type"):
+            _ = dyadic_matrix @ z2
+
+
+class TestSO3Matrix:
+    """Tests for the SO3Matrix class."""
+
+    def test_init_and_repr(self):
+        """Test that SO3Matrix initializes correctly and has correct representations."""
+        z1 = ZOmega(1, 2, 3, 4)
+        z2 = ZOmega(5, 6, 7, 8)
+        dyadic_matrix = DyadicMatrix(z1, z2, z1, z2)
+        so3_matrix = SO3Matrix(dyadic_matrix)
+
+        assert repr(so3_matrix) == f"SO3Matrix(matrix={dyadic_matrix}, k={so3_matrix.k})"
+
+        so3mat, so3k = so3_matrix.so3mat, so3_matrix.k
+        str_repr = "["
+        for i in range(3):
+            str_repr += f"[{so3mat[i][0]}, {so3mat[i][1]}, {so3mat[i][2]}], \n"
+        str_repr = str_repr.rstrip(", \n") + "]" + (f" * 1 / √2^{so3k}" if so3k else "")
+        assert str(so3_matrix) == str_repr
+
+        assert np.allclose(
+            so3_matrix.parity_mat, np.array([[1, 0, 0], [0, 0, 0], [1, 0, 0]])
+        )  # Computed manually
+        assert np.allclose(so3_matrix.parity_vec, [1, 0, 1])  # Computed manually
+
+        z3 = ZOmega()
+        assert np.allclose(SO3Matrix(DyadicMatrix(z3, z3, z3, z3)).ndarray, 0.0)
diff --git a/tests/ops/op_math/decompositions/test_ross_selinger.py b/tests/ops/op_math/decompositions/test_ross_selinger.py
new file mode 100644
index 000000000..2703be661
--- /dev/null
+++ b/tests/ops/op_math/decompositions/test_ross_selinger.py
@@ -0,0 +1,97 @@
+# Copyright 2025 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for Ross-Selinger decomposition's implementation."""
+
+import math
+
+import pytest
+
+import pennylane as qml
+from pennylane.ops.op_math.decompositions.rings import ZOmega
+from pennylane.ops.op_math.decompositions.ross_selinger import _domain_correction, rs_decomposition
+
+
+@pytest.mark.parametrize(
+    ("angle", "result"),
+    [
+        (math.pi / 3, (-math.pi / 2, ZOmega(b=1))),
+        (math.pi / 8, (0.0, ZOmega(d=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))),
+        (math.pi, (-math.pi, ZOmega(d=-1))),
+    ],
+)
+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 scale == result[1]
+
+
+@pytest.mark.parametrize(
+    ("op", "epsilon"),
+    [
+        (qml.RZ(math.pi / 42, wires=[1]), 1e-4),
+        (qml.PhaseShift(math.pi / 7, wires=["a"]), 1e-3),
+        (qml.RZ(math.pi / 3, wires=[1]), 1e-5),
+        (qml.PhaseShift(-math.pi / 6, wires=[0]), 1e-3),
+        (qml.RZ(-math.pi / 5, wires=[0]), 1e-4),
+        (qml.RZ(-math.pi / 8, wires=[2]), 1e-2),
+    ],
+)
+def test_ross_selinger(op, epsilon):
+    """Test Ross-Selinger decomposition method with specified max-depth"""
+    with qml.queuing.AnnotatedQueue() as q:
+        gates = rs_decomposition(op, epsilon=epsilon)
+    assert q.queue == gates
+
+    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.prod(*gates, lazy=False).wires == op.wires
+
+
+def test_epsilon_value_effect():
+    """Test that different epsilon values create different decompositions."""
+    op = qml.RZ(math.pi / 5, 0)
+    decomp_with_error = rs_decomposition(op, 1e-4)
+    decomp_less_error = rs_decomposition(op, 1e-2)
+    assert len(decomp_with_error) > len(decomp_less_error)
+
+
+def test_warm_start():
+    """Test that warm start is working."""
+    op = qml.RZ(math.pi / 8, 0)
+    decomp_with_error = rs_decomposition(op, 1e-10, max_trials=100)
+    decomp_less_error = rs_decomposition(op, 1e-3, max_trials=100)
+    assert len(decomp_with_error) == len(decomp_less_error)
+
+
+def test_exception():
+    """Test operation wire exception in Ross-Selinger"""
+    op = qml.SingleExcitation(1.0, wires=[1, 2])
+
+    with pytest.raises(
+        ValueError,
+        match=r"Operator must be a RZ or PhaseShift gate",
+    ):
+        rs_decomposition(op, epsilon=1e-4, max_trials=1)
diff --git a/tests/transforms/test_cliffordt_transform.py b/tests/transforms/test_cliffordt_transform.py
index d75830921..12f28a8a7 100644
--- a/tests/transforms/test_cliffordt_transform.py
+++ b/tests/transforms/test_cliffordt_transform.py
@@ -155,6 +155,27 @@ class TestCliffordCompile:
         )
         qml.math.isclose(res1, tape_fn([res2]), atol=1e-2)
 
+    @pytest.mark.parametrize("circuit", [circuit_1, circuit_2, circuit_3])
+    def test_decomposition_with_rs(self, circuit):
+        """Test decomposition for the Clifford transform with Ross-Selinger method."""
+
+        old_tape = qml.tape.make_qscript(circuit)()
+
+        [new_tape], tape_fn = clifford_t_decomposition(old_tape, method="rs")
+
+        assert all(
+            isinstance(op, _CLIFFORD_PHASE_GATES)
+            or isinstance(getattr(op, "base", None), _CLIFFORD_PHASE_GATES)
+            for op in new_tape.operations
+        )
+
+        dev = qml.device("default.qubit")
+        transform_program = dev.preprocess_transforms()
+        res1, res2 = qml.execute(
+            [old_tape, new_tape], device=dev, transform_program=transform_program
+        )
+        qml.math.isclose(res1, tape_fn([res2]), atol=1e-2)
+
     def test_qnode_decomposition(self):
         """Test decomposition for the Clifford transform applied to a QNode."""
 
@@ -453,7 +474,7 @@ class TestCliffordCompile:
 
         with pytest.raises(
             NotImplementedError,
-            match=r"Currently we only support Solovay-Kitaev \('sk'\) decomposition",
+            match=r"Currently we only support Solovay-Kitaev \('sk'\) and Ross-Selinger \('rs'\) decompositions",
         ):
             decomposed_qfunc()