diff --git a/src/qibo/gates/abstract.py b/src/qibo/gates/abstract.py
index e744165e97..92f74b2c8f 100644
--- a/src/qibo/gates/abstract.py
+++ b/src/qibo/gates/abstract.py
@@ -1,5 +1,6 @@
 import collections
 import json
+from math import pi
 from typing import List, Sequence, Tuple
 
 import sympy
@@ -359,18 +360,114 @@ def controlled_by(self, *qubits: int) -> "Gate":
             self.control_qubits = qubits
         return self
 
-    def decompose(self, *free) -> List["Gate"]:
+    def _base_decompose(self, *free, use_toffolis=True) -> List["Gate"]:
+        """Base decomposition for gates.
+
+        Returns a list containing the gate itself. Should be overridden by
+        subclasses that support decomposition to simpler gates.
+
+        Args:
+            free: Ids of free qubits to use for the gate decomposition.
+            use_toffolis: If ``True`` the decomposition contains only ``TOFFOLI`` gates.
+                If ``False`` a congruent representation is used for ``TOFFOLI`` gates.
+                See :class:`qibo.gates.TOFFOLI` for more details on this representation.
+
+        Returns:
+            list: Synthesis of the original gate in another gate set.
+        """
+        return [self.__class__(*self.init_args, **self.init_kwargs)]
+
+    @staticmethod
+    def _gates_cancel(g1, g2):
+        """Determines if two gates cancel each other.
+
+        Two gates are considered to cancel if:
+          - They are of the same type (class).
+          - They act on the same target and control qubits.
+          - For fixed gates (like H, CX, X, Y, Z, SWAP), they always cancel in pairs.
+          - For parametrized rotation gates (subclasses of _Rn_), their parameters sum to a multiple of 2π.
+
+        Note:
+            Multi-parameter gates are not currently supported by this check.
+
+        Args:
+            g1, g2: Gate instances to compare.
+
+        Returns:
+            bool: True if the gates cancel each other, False otherwise.
+        """
+        if g1.__class__ != g2.__class__:
+            return False
+
+        if g1.target_qubits != g2.target_qubits:
+            return False
+
+        if g1.control_qubits != g2.control_qubits:
+            return False
+
+        # Identity conditions for fixed gates
+        name = g1.name
+        if name in ("h", "cx", "x", "y", "z", "swap", "ecr", "ccx", "ccz"):
+            return True
+
+        # Check for parametrized rotation gates
+        if "_Rn_" in [base.__name__ for base in g1.__class__.__bases__]:
+            theta1 = g1.parameters[0]
+            theta2 = g2.parameters[0]
+            # Check if theta1 + theta2 is a multiple of 2π
+            return bool((theta1 + theta2) % (2 * pi) < 1e-8)
+
+        return False
+
+    def _control_mask_after_stripping(self, gates: List["Gate"]) -> List[bool]:
+        """Returns a mask indicating which gates should be controlled."""
+        left = 0
+        right = len(gates) - 1
+        mask = [True] * len(gates)
+        while left < right:
+            g1, g2 = gates[left], gates[right]
+            if self._gates_cancel(g1, g2):
+                mask[left] = False
+                mask[right] = False
+            left += 1
+            right -= 1
+        return mask
+
+    def decompose(self, *free, use_toffolis: bool = True) -> List["Gate"]:
         """Decomposes multi-control gates to gates supported by OpenQASM.
 
         Decompositions are based on `arXiv:9503016 <https://arxiv.org/abs/quant-ph/9503016>`_.
+        If the gate is already controlled, it recursively decomposes the base gate and updates
+        the control qubits accordingly.
 
         Args:
             free: Ids of free qubits to use for the gate decomposition.
+            use_toffolis(bool, optional): If ``True``, the decomposition contains only
+                :class:`qibo.gates.TOFFOLI` gates. If ``False``, a congruent
+                representation is used for :class:`qibo.gates.TOFFOLI` gates.
+                See :class:`qibo.gates.TOFFOLI` for more details on this representation.
 
         Returns:
             list: gates that have the same effect as applying the original gate.
         """
-        return [self.__class__(*self.init_args, **self.init_kwargs)]
+        if self.is_controlled_by:
+            # Step 1: Error check with all controls/targets
+            if set(free) & set(self.qubits):
+                raise_error(
+                    ValueError,
+                    "Cannot decompose multi-controlled ``X`` gate if free "
+                    "qubits coincide with target or controls.",
+                )
+            # Step 2: Decompose base gate without controls
+            base_gate = self.__class__(*self.init_args, **self.init_kwargs)
+            decomposed = base_gate._base_decompose(*free, use_toffolis=use_toffolis)
+            mask = self._control_mask_after_stripping(decomposed)
+            for bool_value, gate in zip(mask, decomposed):
+                if bool_value:
+                    gate.is_controlled_by = True
+                    gate.control_qubits += self.control_qubits
+            return decomposed
+        return self._base_decompose(*free, use_toffolis=use_toffolis)
 
     def matrix(self, backend=None):
         """Returns the matrix representation of the gate.
diff --git a/src/qibo/gates/gates.py b/src/qibo/gates/gates.py
index 1ed4cdd6f5..853dfad54f 100644
--- a/src/qibo/gates/gates.py
+++ b/src/qibo/gates/gates.py
@@ -82,7 +82,7 @@ def controlled_by(self, *q):
             gate = super().controlled_by(*q)
         return gate
 
-    def decompose(self, *free, use_toffolis=True):
+    def _base_decompose(self, *free, use_toffolis=True):
         """Decomposes multi-control ``X`` gate to one-qubit, ``CNOT`` and ``TOFFOLI`` gates.
 
         Args:
@@ -133,12 +133,12 @@ def decompose(self, *free, use_toffolis=True):
             m1 = n // 2
             free1 = controls[m1:] + (target,) + tuple(free[1:])
             x1 = self.__class__(free[0]).controlled_by(*controls[:m1])
-            part1 = x1.decompose(*free1, use_toffolis=use_toffolis)
+            part1 = x1._base_decompose(*free1, use_toffolis=use_toffolis)
 
             free2 = controls[:m1] + tuple(free[1:])
             controls2 = controls[m1:] + (free[0],)
             x2 = self.__class__(target).controlled_by(*controls2)
-            part2 = x2.decompose(*free2, use_toffolis=use_toffolis)
+            part2 = x2._base_decompose(*free2, use_toffolis=use_toffolis)
 
             decomp_gates = [*part1, *part2]
 
@@ -286,7 +286,7 @@ def clifford(self):
     def qasm_label(self):
         return "sx"
 
-    def decompose(self):
+    def _base_decompose(self, *free, use_toffolis=True):
         """Decomposition of :math:`\\sqrt{X}` up to global phase.
 
         A global phase difference exists between the definitions of
@@ -336,7 +336,7 @@ def clifford(self):
     def qasm_label(self):
         return "sxdg"
 
-    def decompose(self):
+    def _base_decompose(self, *free, use_toffolis=True):
         """Decomposition of :math:`(\\sqrt{X})^{\\dagger}` up to global phase.
 
         A global phase difference exists between the definitions of
@@ -789,7 +789,7 @@ def _dagger(self) -> "Gate":
             self.target_qubits[0], theta, phi
         )  # pylint: disable=E1130
 
-    def decompose(self):
+    def _base_decompose(self, *free, use_toffolis=True):
         """Decomposition of Phase-:math:`RX` gate."""
         from qibo.transpiler.decompositions import (  # pylint: disable=C0415
             standard_decompositions,
@@ -1052,7 +1052,7 @@ def _dagger(self) -> "Gate":
         theta, lam, phi = tuple(-x for x in self.parameters)  # pylint: disable=E1130
         return self.__class__(self.target_qubits[0], theta, phi, lam)
 
-    def decompose(self) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis=True) -> List[Gate]:
         """Decomposition of :math:`U_{3}` up to global phase.
 
         A global phase difference exists between the definitions of
@@ -1154,7 +1154,7 @@ def clifford(self):
     def qasm_label(self):
         return "cx"
 
-    def decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
         q0, q1 = self.control_qubits[0], self.target_qubits[0]
         return [self.__class__(q0, q1)]
 
@@ -1194,7 +1194,7 @@ def clifford(self):
     def qasm_label(self):
         return "cy"
 
-    def decompose(self) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis=True) -> List[Gate]:
         """Decomposition of :math:`\\text{CY}` gate.
 
         Decompose :math:`\\text{CY}` gate into :class:`qibo.gates.SDG` in
@@ -1247,7 +1247,7 @@ def hamming_weight(self):
     def qasm_label(self):
         return "cz"
 
-    def decompose(self) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis=True) -> List[Gate]:
         """Decomposition of :math:`\\text{CZ}` gate.
 
         Decompose :math:`\\text{CZ}` gate into :class:`qibo.gates.H` in
@@ -1292,7 +1292,7 @@ def __init__(self, q0, q1):
     def qasm_label(self):
         return "csx"
 
-    def decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
         """"""
         from qibo.transpiler.decompositions import (  # pylint: disable=C0415
             standard_decompositions,
@@ -1336,7 +1336,7 @@ def __init__(self, q0, q1):
     def qasm_label(self):
         return "csxdg"
 
-    def decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
         """"""
         from qibo.transpiler.decompositions import (  # pylint: disable=C0415
             standard_decompositions,
@@ -1457,7 +1457,7 @@ def hamming_weight(self):
     def qasm_label(self):
         return "cry"
 
-    def decompose(self) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis=True) -> List[Gate]:
         """Decomposition of :math:`\\text{CRY}` gate."""
         from qibo.transpiler.decompositions import (  # pylint: disable=C0415
             standard_decompositions,
@@ -1856,10 +1856,10 @@ def hamming_weight(self):
     def qasm_label(self):
         return "fswap"
 
-    def decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
         """"""
         q0, q1 = self.target_qubits
-        return [X(q1)] + GIVENS(q0, q1, np.pi / 2).decompose() + [X(q0)]
+        return [X(q1)] + GIVENS(q0, q1, np.pi / 2)._base_decompose() + [X(q0)]
 
 
 class fSim(ParametrizedGate):
@@ -2185,7 +2185,7 @@ def __init__(self, q0, q1, theta, trainable=True):
     def hamming_weight(self):
         return _is_hamming_weight_given_angle(self.parameters[0])
 
-    def decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
         """"""
         from qibo.transpiler.decompositions import (  # pylint: disable=C0415
             standard_decompositions,
@@ -2227,7 +2227,7 @@ def __init__(self, q0, q1, theta, trainable=True):
     def hamming_weight(self):
         return True
 
-    def decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
         """Decomposition of :math:`\\text{R_{XX-YY}}` up to global phase.
 
         This decomposition has a global phase difference with respect to
@@ -2354,7 +2354,7 @@ def _dagger(self) -> "Gate":
         """"""
         return self.__class__(*self.target_qubits, -self.parameters[0])
 
-    def decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
         """Decomposition of RBS gate according to `ArXiv:2109.09685
         <https://arxiv.org/abs/2109.09685>`_."""
         from qibo.transpiler.decompositions import (  # pylint: disable=C0415
@@ -2412,7 +2412,7 @@ def _dagger(self) -> "Gate":
         """"""
         return self.__class__(*self.target_qubits, -self.parameters[0])
 
-    def decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
         """Decomposition of RBS gate according to `ArXiv:2109.09685
         <https://arxiv.org/abs/2109.09685>`_."""
         from qibo.transpiler.decompositions import (  # pylint: disable=C0415
@@ -2453,7 +2453,7 @@ def __init__(self, q0, q1):
     def clifford(self):
         return True
 
-    def decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
         """Decomposition of :math:`\\textup{ECR}` gate up to global phase.
 
         A global phase difference exists between the definitions of
@@ -2506,7 +2506,7 @@ def __init__(self, q0, q1, q2):
     def qasm_label(self):
         return "ccx"
 
-    def decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis: bool = True) -> List[Gate]:
         """Decomposition of :math:`\\text{TOFFOLI}` gate.
 
         Decompose :math:`\\text{TOFFOLI}` gate into :class:`qibo.gates.CNOT` gates,
@@ -2536,8 +2536,6 @@ def congruent(self, use_toffolis: bool = True) -> List[Gate]:
             List with ``RY`` and ``CNOT`` gates that have the same effect as
             applying the original ``TOFFOLI`` gate.
         """
-        if use_toffolis:
-            return self.decompose()
 
         control0, control1 = self.control_qubits
         target = self.target_qubits[0]
@@ -2592,7 +2590,7 @@ def hamming_weight(self):
     def qasm_label(self):
         return "ccz"
 
-    def decompose(self) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis=True) -> List[Gate]:
         """Decomposition of :math:`\\text{CCZ}` gate.
 
         Decompose :math:`\\text{CCZ}` gate into :class:`qibo.gates.H` in
@@ -2712,7 +2710,7 @@ def __init__(
     def hamming_weight(self):
         return len(self.init_args[0]) == len(self.init_args[1])
 
-    def decompose(self) -> List[Gate]:
+    def _base_decompose(self, *free, use_toffolis=True) -> List[Gate]:
         """Decomposition of :math:`\\text{gRBS}` gate.
 
         Decompose :math:`\\text{gRBS}` gate into :class:`qibo.gates.X`, :class:`qibo.gates.CNOT`,
diff --git a/tests/test_gates_abstract.py b/tests/test_gates_abstract.py
index 8dccfd0997..8bdca6920f 100644
--- a/tests/test_gates_abstract.py
+++ b/tests/test_gates_abstract.py
@@ -2,13 +2,15 @@
 `qibo/gates/gates.py`."""
 
 import json
+from math import pi
 from typing import Optional
 
 import pytest
 
 from qibo import gates, matrices
 from qibo.config import PRECISION_TOL
-from qibo.gates import abstract
+from qibo.gates import Gate, abstract
+from qibo.models import Circuit
 
 
 @pytest.mark.parametrize(
@@ -422,6 +424,76 @@ def test_decompose():
     assert isinstance(decomp_gates[0], gates.H)
 
 
+def test_decompose_controlled(backend):
+    target = gates.H(0).controlled_by(1)
+    decomp = target.decompose()
+
+    assert len(decomp) == 1
+    assert isinstance(decomp[0], gates.H)
+    assert decomp[0].control_qubits == (1,)
+
+    circuit_1 = Circuit(2)
+    circuit_2 = circuit_1.copy(deep=True)
+    circuit_1.add(target)
+    circuit_2.add(decomp)
+
+    backend.assert_circuitclose(circuit_1, circuit_2)
+
+
+def test_decompose_controlled_optimized(backend):
+    target = gates.RBS(1, 2, 0.1).controlled_by(0)
+    decomp = target.decompose()
+
+    assert len(decomp) == 6
+
+    controls_on_zero = sum([0 in g.control_qubits for g in decomp])
+    assert controls_on_zero == 2
+
+    circuit_1 = Circuit(3)
+    circuit_2 = circuit_1.copy(deep=True)
+    circuit_1.add(target)
+    circuit_2.add(decomp)
+
+    backend.assert_circuitclose(circuit_1, circuit_2)
+
+
+@pytest.mark.parametrize(
+    "g1, g2, expected",
+    [
+        (gates.X(0), gates.X(0), True),
+        (gates.X(0), gates.X(1), False),
+        (gates.X(0), gates.Y(0), False),
+        (gates.Y(1), gates.Y(1), True),
+        (gates.Z(2), gates.Z(2), True),
+        (gates.SX(0), gates.SX(0), False),
+        (gates.RX(0, pi), gates.RX(0, -pi), True),
+        (gates.RX(0, pi / 2), gates.RX(0, pi / 3), False),
+        (gates.RY(1, 0.5), gates.RY(1, -0.5), True),
+        (gates.RZ(2, 2 * pi), gates.RZ(2, -2 * pi), True),
+        (gates.CNOT(0, 1), gates.CNOT(0, 1), True),
+        (gates.CNOT(0, 2), gates.CNOT(1, 2), False),
+        (gates.TOFFOLI(0, 1, 2), gates.TOFFOLI(0, 1, 2), True),
+        (gates.CCZ(0, 1, 2), gates.CCZ(0, 1, 2), True),
+    ],
+)
+def test_gates_cancel(g1, g2, expected):
+    result = Gate._gates_cancel(g1, g2)
+    assert result is expected
+
+
+def test_toffoli_congruent(backend):
+    congruent = gates.TOFFOLI(0, 1, 2)
+
+    circuit = Circuit(3)
+    circuit.add(congruent.congruent())
+    congruent = circuit.unitary(backend)
+
+    target = backend.matrices.TOFFOLI
+    target[4, 4] = -1
+
+    assert backend.calculate_matrix_norm(congruent - target) < 1e-8
+
+
 def test_special_gate():
     from qibo.gates.abstract import SpecialGate