diff --git a/src/pyqasm/algorithms/__init__.py b/src/pyqasm/algorithms/__init__.py
new file mode 100644
index 00000000..91f8d6fc
--- /dev/null
+++ b/src/pyqasm/algorithms/__init__.py
@@ -0,0 +1,14 @@
+"""
+Sub module for quantum algorithms.
+
+Functions:
+----------
+    solovay_kitaev: Solovay-Kitaev algorithm for approximating unitary gates.
+
+"""
+
+from .solovay_kitaev import solovay_kitaev
+
+__all__ = [
+    "solovay_kitaev",
+]
diff --git a/src/pyqasm/algorithms/solovay_kitaev/__init__.py b/src/pyqasm/algorithms/solovay_kitaev/__init__.py
new file mode 100644
index 00000000..91f8d6fc
--- /dev/null
+++ b/src/pyqasm/algorithms/solovay_kitaev/__init__.py
@@ -0,0 +1,14 @@
+"""
+Sub module for quantum algorithms.
+
+Functions:
+----------
+    solovay_kitaev: Solovay-Kitaev algorithm for approximating unitary gates.
+
+"""
+
+from .solovay_kitaev import solovay_kitaev
+
+__all__ = [
+    "solovay_kitaev",
+]
diff --git a/src/pyqasm/algorithms/solovay_kitaev/basic_approximation.py b/src/pyqasm/algorithms/solovay_kitaev/basic_approximation.py
new file mode 100644
index 00000000..5c77177e
--- /dev/null
+++ b/src/pyqasm/algorithms/solovay_kitaev/basic_approximation.py
@@ -0,0 +1,98 @@
+"""
+Definition of the basic approximation algorithm for the Solovay-Kitaev theorem.
+"""
+
+import numpy as np
+
+from pyqasm.algorithms.solovay_kitaev.utils import TU2Matrix, get_tu2matrix_for_basic_approximation
+from pyqasm.elements import BasisSet
+
+
+def rescursive_traversal(
+    target_matrix, approximated_matrix, target_gate_set_list, current_depth, params
+):
+    """Recursively traverse the tree to find the best approximation of the target matrix.
+    Args:
+        target_matrix (np.ndarray): The target matrix to approximate.
+        approximated_matrix (TU2Matrix): The approximated matrix.
+        target_gate_set_list (list): The list of target gates to approximate.
+        current_depth (int): The current depth of the tree.
+        params (tuple): The parameters for the approximation.
+
+    Returns:
+        float: The closest difference between the target and approximated matrix.
+        TU2Matrix: The closest approximated matrix.
+        TU2Matrix: The best approx
+
+    """
+    accuracy, max_tree_depth, best_gate = params
+
+    if current_depth >= max_tree_depth:
+        return best_gate
+
+    for gate in target_gate_set_list:
+        if not approximated_matrix.can_multiple(gate):
+            continue
+        approximated_matrix_copy = approximated_matrix.copy()
+        approximated_matrix = approximated_matrix * gate
+
+        diff = approximated_matrix.distance(target_matrix)
+        if diff < accuracy:
+            best_gate = approximated_matrix.copy()
+            return best_gate
+
+        # Update the closest gate if the current one is closer
+        if diff < best_gate.distance(target_matrix):
+            best_gate = approximated_matrix.copy()
+
+        best_gate = rescursive_traversal(
+            target_matrix,
+            approximated_matrix.copy(),
+            target_gate_set_list,
+            current_depth + 1,
+            (accuracy, max_tree_depth, best_gate),
+        )
+        approximated_matrix = approximated_matrix_copy.copy()
+
+    return best_gate
+
+
+def basic_approximation(target_matrix, target_gate_set, accuracy=0.001, max_tree_depth=3):
+    """Approximate the target matrix using the basic approximation algorithm.
+
+    Args:
+        target_matrix (np.ndarray): The target matrix to approximate.
+        target_gate_set (BasisSet): The target gate set to approximate.
+        accuracy (float): The accuracy of the approximation.
+        max_tree_depth (int): The maximum depth of the tree.
+
+    Returns:
+        TU2Matrix: The approximated matrix.
+    """
+    approximated_matrix = TU2Matrix(np.identity(2), [], None, None)
+    target_gate_set_list = get_tu2matrix_for_basic_approximation(target_gate_set)
+    current_depth = 0
+    best_gate = TU2Matrix(np.identity(2), [], None, None)
+
+    params = (accuracy, max_tree_depth, best_gate)
+
+    best_gate = rescursive_traversal(
+        target_matrix, approximated_matrix.copy(), target_gate_set_list, current_depth, params
+    )
+
+    return best_gate
+
+    # result = None
+
+    # if best_gate:
+    #     result = best_gate.copy()
+    # else:
+    #     result = closest_gate.copy()
+
+    # return result
+
+
+if __name__ == "__main__":
+    target_matrix_U = np.array([[0.70711, 0.70711j], [0.70711j, 0.70711]])
+
+    print(basic_approximation(target_matrix_U, BasisSet.CLIFFORD_T, 0.0001, 3))
diff --git a/src/pyqasm/algorithms/solovay_kitaev/gc_decompose.py b/src/pyqasm/algorithms/solovay_kitaev/gc_decompose.py
new file mode 100644
index 00000000..99cddeb4
--- /dev/null
+++ b/src/pyqasm/algorithms/solovay_kitaev/gc_decompose.py
@@ -0,0 +1,59 @@
+import math
+from typing import List
+import numpy as np
+from pyqasm.algorithms.solovay_kitaev.utils import SU2Matrix
+
+def u_to_bloch(u):
+  """Compute angle and axis for a unitary."""
+
+  angle = np.real(np.arccos((u[0, 0] + u[1, 1]) / 2))
+  sin = np.sin(angle)
+  if sin < 1e-10:
+    axis = [0, 0, 1]
+  else:
+    nx = (u[0, 1] + u[1, 0]) / (2j * sin)
+    ny = (u[0, 1] - u[1, 0]) / (2 * sin)
+    nz = (u[0, 0] - u[1, 1]) / (2j * sin)
+    axis = [nx, ny, nz]
+  return axis, 2 * angle
+
+def Rotation(vparm: List[float], theta: float, name: str):
+  """Produce the single-qubit rotation operator."""
+
+  v = np.asarray(vparm)
+  if v.shape != (3,) or not math.isclose(v @ v, 1) or not np.all(np.isreal(v)):
+    raise ValueError('Rotation vector v must be a 3D real unit vector.')
+
+  return np.cos(theta / 2) * np.identity(2) - 1j * np.sin(theta / 2) * (v[0] * np.array([[0.0, 1.0], [1.0, 0.0]]) + v[1] * np.array([[0.0, -1.0j], [1.0j, 0.0]]) + v[2] * np.array([[1.0, 0.0], [0.0, -1.0]]))
+
+
+def RotationX(theta: float):
+  return Rotation([1.0, 0.0, 0.0], theta, 'Rx')
+
+
+def RotationY(theta: float):
+  return Rotation([0.0, 1.0, 0.0], theta, 'Ry')
+
+
+def RotationZ(theta: float):
+  return Rotation([0.0, 0.0, 1.0], theta, 'Rz')
+
+def gc_decomp(u):
+    axis, theta = u_to_bloch(u)
+    
+    phi = 2.0 * np.arcsin(np.sqrt(np.sqrt((0.5 - 0.5 * np.cos(theta / 2)))))
+    v = RotationX(phi)
+    if axis[2] > 0:
+        w = RotationY(2 * np.pi - phi)
+    else:
+        w = RotationY(phi)
+    
+    _, ud = np.linalg.eig(u)
+    
+    vwvdwd = v @ w @ v.conj().T @ w.conj().T
+    
+    s = ud @ vwvdwd.conj().T
+
+    v_hat = s @ v @ s.conj().T
+    w_hat = s @ w @ s.conj().T
+    return SU2Matrix(v_hat, []), SU2Matrix(w_hat, [])
\ No newline at end of file
diff --git a/src/pyqasm/algorithms/solovay_kitaev/optimizer.py b/src/pyqasm/algorithms/solovay_kitaev/optimizer.py
new file mode 100644
index 00000000..55ab294a
--- /dev/null
+++ b/src/pyqasm/algorithms/solovay_kitaev/optimizer.py
@@ -0,0 +1,94 @@
+"""
+Definition of the optimizer for the Solovay-Kitaev theorem.
+"""
+
+from pyqasm.algorithms.solovay_kitaev.utils import get_identity_weight_group_for_optimizer
+from pyqasm.elements import BasisSet
+
+
+def optimize_gate_sequence(seq: list[str], target_basis_set):
+    """Optimize a gate sequence based on the identity weight group.
+    Args:
+        seq (list[str]): The gate sequence to optimize.
+        target_basis_set (BasisSet): The target basis set.
+
+    Returns:
+        list[str]: The optimized gate sequence.
+    """
+    target_identity_weight_group = get_identity_weight_group_for_optimizer(target_basis_set)
+    for _ in range(int(1e6)):
+        current_group = None
+        current_weight = 0
+        start_index = 0
+        changed = False
+
+        for i, gate_name in enumerate(seq):
+            if gate_name not in target_identity_weight_group:
+                continue
+
+            gate = target_identity_weight_group[gate_name]
+            new_group = gate["group"]
+            new_weight = gate["weight"]
+
+            if current_group is None or new_group != current_group:
+                current_group = new_group
+                current_weight = new_weight
+                start_index = i
+            else:
+                current_weight += new_weight
+
+            if current_weight == 1:
+                seq = seq[:start_index] + seq[i + 1 :]
+                changed = True
+                break
+            if current_weight > 1:
+                remaining_weight = current_weight - 1
+                for key, value in target_identity_weight_group.items():
+                    if value["group"] == current_group and value["weight"] == remaining_weight:
+                        seq = seq[:start_index] + [key] + seq[i + 1 :]
+                        changed = True
+                        break
+                break
+
+        if not changed:
+            return seq
+
+    return seq
+
+
+if __name__ == "__main__":
+    s1 = ["s", "s", "s", "t", "t", "tdg", "sdg", "sdg", "sdg", "tdg", "s", "h", "s"]
+    s2 = [
+        "t",
+        "s",
+        "s",
+        "s",
+        "t",
+        "tdg",
+        "tdg",
+        "sdg",
+        "sdg",
+        "sdg",
+        "t",
+        "s",
+        "s",
+        "s",
+        "t",
+        "tdg",
+        "tdg",
+        "sdg",
+        "sdg",
+        "sdg",
+        "s",
+        "h",
+        "s",
+    ]
+    s3 = ["h", "s", "s", "t", "t", "s", "t"]  # ['h', 't']
+    s4 = ["h", "s", "s", "t", "t", "s", "h"]  # []
+    s5 = ["h", "s", "s", "t", "h", "h", "t", "s", "h", "t"]  # ['t']
+
+    print(optimize_gate_sequence(s1, BasisSet.CLIFFORD_T) == ["s", "h", "s"])
+    print(optimize_gate_sequence(s2, BasisSet.CLIFFORD_T) == ["s", "h", "s"])
+    print(optimize_gate_sequence(s3, BasisSet.CLIFFORD_T) == ["h", "t"])
+    print(optimize_gate_sequence(s4, BasisSet.CLIFFORD_T) == [])
+    print(optimize_gate_sequence(s5, BasisSet.CLIFFORD_T) == ["t"])
diff --git a/src/pyqasm/algorithms/solovay_kitaev/solovay_kitaev.py b/src/pyqasm/algorithms/solovay_kitaev/solovay_kitaev.py
new file mode 100644
index 00000000..349cc03d
--- /dev/null
+++ b/src/pyqasm/algorithms/solovay_kitaev/solovay_kitaev.py
@@ -0,0 +1,182 @@
+"""
+Definition of the Solovay-Kitaev algorithm.
+"""
+
+from typing import List, Tuple
+
+import numpy as np
+
+from pyqasm.algorithms.solovay_kitaev.basic_approximation import basic_approximation
+from pyqasm.algorithms.solovay_kitaev.gc_decompose import gc_decomp
+from pyqasm.algorithms.solovay_kitaev.optimizer import optimize_gate_sequence
+from pyqasm.algorithms.solovay_kitaev.utils import (
+    SU2Matrix,
+    get_su2matrix_for_solovay_kitaev_algorithm,
+)
+from pyqasm.elements import BasisSet
+count = 1093
+
+def group_commutator(a: SU2Matrix, b: SU2Matrix) -> SU2Matrix:
+    """Compute the group commutator [a,b] = aba^{-1}b^{-1}."""
+    return a * b * a.dagger() * b.dagger()
+
+
+def find_basic_approximation(
+    target_matrix: SU2Matrix, target_basis_set, use_optimization, accuracy=1e-6, max_tree_depth=10
+) -> SU2Matrix:
+    """Find the basic approximation of a target matrix.
+
+    Args:
+        target_matrix (SU2Matrix): The target matrix to approximate.
+        target_basis_set (BasisSet): The basis set to use for the approximation.
+        use_optimization (bool): Whether to use optimization to reduce the number of gates.
+        accuracy (float): The accuracy of the approximation.
+        max_tree_depth (int): The maximum depth of the tree.
+
+    Returns:
+        SU2Matrix: The basic approximation of the target matrix.
+    """
+    gates = basic_approximation(target_matrix, target_basis_set, accuracy, max_tree_depth)
+    if use_optimization:
+        gates.name = optimize_gate_sequence(gates.name, target_basis_set)
+    return SU2Matrix(gates.matrix, gates.name)
+
+
+def decompose_group_element(
+    target: SU2Matrix,
+    target_gate_set,
+    basic_gates: List[SU2Matrix],
+    depth,
+    accuracy,
+    use_optimization,
+) -> Tuple[List[SU2Matrix], float]:
+    """Decompose a group element into a sequence of basic gates.
+
+    Args:
+        target (SU2Matrix): The target group element.
+        target_gate_set (BasisSet): The target gate set.
+        basic_gates (List[SU2Matrix]): The basic gates to use for the approximation.
+        depth (int): The depth of the approximation.
+        accuracy (float): The accuracy of the approximation.
+        use_optimization (bool): Whether to use optimization to reduce the number of gates.
+
+    Returns:
+        Tuple[List[SU2Matrix], float]: The sequence of basic gates and the error.
+    """
+    # global count
+    # count -= 1
+    # print(count)
+    
+    if depth == 0:
+        best_approx = find_basic_approximation(
+            target.matrix, target_gate_set, use_optimization=use_optimization
+        )
+        return best_approx, target.distance(best_approx)
+
+    # Recursive approximation
+    prev_sequence, prev_error = decompose_group_element(
+        target, target_gate_set, basic_gates, depth - 1, accuracy, use_optimization
+    )
+
+    # If previous approximation is good enough, return it
+    # ERROR IS HARD CODED RIGHT NOW -> CHANGE THIS TO FIT USER-INPUT
+    if prev_error < accuracy:
+        return prev_sequence, prev_error
+
+    # error = target * prev_sequence.dagger()
+
+    # # Find Va and Vb such that their group commutator approximates the error
+    # best_v = None
+    # best_w = None
+    # best_error = float("inf")
+
+    # for v in basic_gates:
+    #     for w in basic_gates:
+    #         comm = group_commutator(v, w)
+    #         curr_error = error.distance(comm)
+    #         if curr_error < best_error:
+    #             best_error = curr_error
+    #             best_v = v
+    #             best_w = w
+
+    result = prev_sequence
+    
+    best_v, best_w = gc_decomp((target * prev_sequence.dagger()).matrix)
+
+    # Add correction terms
+    if best_v is not None and best_w is not None:
+        v_sequence, error = decompose_group_element(
+            best_v, target_gate_set, basic_gates, depth - 1, accuracy, use_optimization
+        )
+        w_sequence, error = decompose_group_element(
+            best_w, target_gate_set, basic_gates, depth - 1, accuracy, use_optimization
+        )
+
+        result = group_commutator(v_sequence, w_sequence) * prev_sequence
+
+    final_error = target.distance(result)
+    return result, final_error
+
+
+def solovay_kitaev(
+    target: np.ndarray, target_basis_set, depth: int = 3, accuracy=1e-6, use_optimization=True
+) -> List[np.ndarray]:
+    """
+    Main function to run the Solovay-Kitaev algorithm.
+
+    Args:
+        target: The target unitary matrix to approximate
+        target_basis_set: The basis set to use for the approximation
+        depth: The depth of the approximation
+        accuracy: The accuracy of the approximation
+        use_optimization: Whether to use optimization to reduce the number
+
+    Returns:
+        A list of gates that approximate the target unitary matrix
+    """
+    # Convert inputs to SU2Matrix objects
+    target_su2 = SU2Matrix(target, [])
+
+    basic_gates_su2 = get_su2matrix_for_solovay_kitaev_algorithm(target_basis_set)
+
+    # Run the decomposition
+    sequence, _ = decompose_group_element(
+        target_su2, target_basis_set, basic_gates_su2, depth, accuracy, use_optimization
+    )
+
+    if use_optimization:
+        sequence.name = optimize_gate_sequence(sequence.name, target_basis_set)
+        return sequence
+
+    return sequence
+
+
+# if __name__ == "__main__":
+#     target_matrix_U = np.array([[0.70711, 0.70711j], [0.70711j, 0.70711]])
+
+#     r0 = solovay_kitaev(target_matrix_U, BasisSet.CLIFFORD_T, depth=0)
+#     print(r0.name)  # Output: ['s', 'h', 's']
+
+#     r1 = solovay_kitaev(target_matrix_U, BasisSet.CLIFFORD_T, depth=1)
+#     print(
+#         r1.name
+#     )  # Output: ['s', 's', 's', 't', 't', 'tdg', 'sdg', 'sdg', 'sdg', 'tdg', 's', 'h', 's']
+
+#     r2 = solovay_kitaev(target_matrix_U, BasisSet.CLIFFORD_T, depth=2)
+#     print(r2.name)  # Output: ['t', 's', 's', 's', 't',
+#     #             'tdg', 'tdg', 'sdg', 'sdg', 'sdg',
+#     #             't', 's', 's', 's', 't',
+#     #             'tdg', 'tdg', 'sdg', 'sdg', 'sdg',
+#     #             's', 'h', 's']
+
+#     print(np.allclose(r0.matrix, r1.matrix))  # Output: True
+#     print(np.allclose(r1.matrix, r2.matrix))  # Output: True
+#     print(np.allclose(r2.matrix, r0.matrix))  # Output: True
+
+if __name__ == '__main__':
+    U = np.array([[-0.8987577 +0.j, -0.31607208+0.30386352j],[-0.14188982-0.41485163j,  0.79903828+0.41146473j]])
+    r0 = solovay_kitaev(U, BasisSet.CLIFFORD_T, depth=2)
+    print("----------------")
+    print(r0.distance(SU2Matrix(U, [])))
+    print(r0.matrix)
+    # print(r0.name)
diff --git a/src/pyqasm/algorithms/solovay_kitaev/testing.ipynb b/src/pyqasm/algorithms/solovay_kitaev/testing.ipynb
new file mode 100644
index 00000000..8f553fae
--- /dev/null
+++ b/src/pyqasm/algorithms/solovay_kitaev/testing.ipynb
@@ -0,0 +1,348 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import qiskit\n",
+    "from qiskit.circuit.random import random_circuit\n",
+    "from qiskit.quantum_info import Operator\n",
+    "import random"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Generating random circuits"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class CircuitInfo:\n",
+    "    def __init__(self, qc, resultant_unitary, gate_names):\n",
+    "        self.qc = qc\n",
+    "        self.resultant_unitary = resultant_unitary\n",
+    "        self.gate_names = gate_names"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "circuit_list = []\n",
+    "\n",
+    "for _ in range(10):\n",
+    "    qc = random_circuit(1, random.randint(1, 20), measure=False)\n",
+    "    resultant_unitary = Operator(qc).data\n",
+    "    gate_names = [instruction.operation.name for instruction in qc.data]\n",
+    "    \n",
+    "    circuit_list.append(CircuitInfo(qc, resultant_unitary, gate_names))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASMAAABuCAYAAABskXUrAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAACBdJREFUeJzt3QtMlecdx/EfBxQOIqbqVCheEARvIBtWJZ2dOs2mibaLWWtmrYvWJdu8dGXi2rlpdB3zss44s8Ul1ksbDWbZ6maadFl0yIwXHF5o1eJQjCJonfSCIAqc5X1NujTi6uUV/uc9309iMOd9fc855PD1eS+8T1QoFAoJADpYoKNfAAA4iBEAE4gRABOIEQATiBEAE4gRABOIEQATiBEAE4gRABOIEQATiBEAE4gRABOIEQATiBEAE4gRABOIEQATiBEAE4gRABOIEQATiBEAE4gRABOIEQATiBEAE4gRABOIEQATiBEAE4gRABOIEQATiBEAE4gRABOIEQATiBEAE4gRABOIEQATiBEAE4gRABOIEQATYjr6BfhZKCTdaFFYiYuWoqK8214oFFJDY7PCSXwwRlEefROc96+mJoWV2FjP3v/9IEaPkBOise8orJRMkYIefiqcECWM2aZwUn/wBXWJ7+TNxpqa1PzsbIWTmJ1bpbi4dn9edtMAmECMAJhAjACYQIwAmECMAJhAjACYQIwAmECMAJhAjACYQIwAmECMAJhAjACYQIwAmOC7GF29elUFBQVKT09XXFyc+vbtq0WLFun69euaO3eue2uEDRs2dPTLBODnW4gcO3ZMkydPVm1trbp06aKhQ4fq0qVLWr9+vSorK3Xt2jV3vZycHIWbfz19b/eXyfjFXnXNGie/eW1hrl59MUdzfr5Pm98+c8fyvZumKG9EL+XO2KX3/10nvym+ekWTDvxDvxqarZfTBre5Tue/7tSUXkl6e/RYhaMYP42Ipk6d6oYoPz9fy5YtU9euXd1lq1ev1pIlSxQTc/umWdnZ2Qo3A3705l2XNdWeVc2OZYpJ7Km4xzPlR8t/d1RTv9ZPr/94tP52oFrVlxs+W/bS88M07okk/WRdqS9DFCl8E6OFCxfq4sWLmj9/vtauXfu5Zc5u2/bt23X8+HGlpqYqMTFR4abHuOfbfLy1qUGnC/KkQLRSFxepU/ck+dGt5lbNXrpPh96apk3Lx+qb33/XfTxjQDe9tmCkDp64ojVbyjv6ZSLSjxmdOnVKRUVF6tmzpwoLC9tcJzc31/06YsQI+UnV+jlqrDqhlNmrlJg9QX529NR/VLjpuL7xZIrmTc9UIBClba895d4m1wlVa2uoo18iIn1ktGPHDrW2tmrmzJlKSEhoc51gMOi7GNX+aY3q/lmkx776nHo/k69IsPIPRzVtXD+tzR+lnME9NDqrl15ec0gVVR8rEjS0tOhquN1TO5JitGfPHvfr+PHj77qOswvnpxh9cuzvqn7zFQX7Z6n/gk2KFM3NIXcUVLpjmn7w3BCVlNVq3VvvKVKs+OB9948f+SJG58+fd7/279+/zeXNzc3av3//Q8do5MiR7gHyexXVOaje6+488/Owmi5X6ezaGYoOdlXaK39WdFwXz7adkTFIoZuNnm2vVZ2k7kvlpY/rb6rpZos6d4rWOyUX3FlYvDQoI0MB3fJkW8FAQCdz8uSVF/sN1PTkvm0um3yw2JPnyMjIUGNr6wP92z59+ujIkSORGyPnGiJHY2PbP0TO8STnbJtzds05gP2gnBBVV1ff8/qB2Hj1lrecA9aVhd9Sy/U6pS/drdikNE+371wK4TyHZ6I6S93lqc0rxrohOllZp6Xfy9HOd8/p7MVPPdt+zaVLUuimJ9uKj46WPLySJD0hQV//ktefqjs/A87uYHvzRYycGtfV1amsrEx5eZ//X6impkaLFy92/+6c0n+Y+aCc57kfzsjIa+c3zFPjuWNKnrlS3XIne7795ORkz0dGNZ5tTVrwnaEaPypZr64/ol17z6us6Bm9sWKsxs3xbk6opORkT0dG4SY5OfmhRkYRHaOJEye6Z9RWrVqlSZMmucNMR2lpqWbNmuWOiry42PF+h5/O3IVezpt2edfrurZvu7qNflp9vv1TPQoVFWc8nTftesMtz+ZNS++XqMJFI3W4/EOteuOEe/Zs+e/LVLjoCTdSv91+0pPnOVNR4dm8aaEbN8Ju3rSKigpFMW/ag3GuI+rRo4cuXLigYcOGKSsrS4MGDdKoUaM0cOBATZgwIewPXn96Yq8ubilQ7OOZSn1pW4fM+NmRnLe7ZeVTig5EafbS4s9O46/eXK7S9z50IzUw5fZFrghPvhgZpaSkqKSkxN0dKy4uVlVVlfurIBs3btS8efOUlpYW1jG6da1GZ9c8K7W26LG86fro8F/uum5wQLbiB4TfFeZfJH92lp78cm8V/OawTp/732l8J0rf/dm+R7K7hvblixg5hgwZot27d9/xeH19vRunQCCg4cOHKxzdqP5AzZ/c3tWs/eMv/++6STOW+S5Gg1O7aeUPv6IDx6/o11vvPI1/svKjR7K7hvYVFQp5fWLUlkOHDmnMmDHKzMzU6dOn2/W5vT5m1B5KpsjsMaP2Un/whYg+ZhSzcyvHjB6F8vLysN5FAyIFMQJgAjECYIJvDmB/0e+tAbDN9yMjAOGBGAEwgRgBMIEYATCBGAEwgRgBMIEYATCBGAEwgRgBMIEYATCBGAEwwff3M+pIznf2RvtPsvBQ4qJv3+LVK87Hq8G5sVMYiQ/GeHZbX/fHK9wmXYyN7ZDbGhMjACawmwbABGIEwARiBMAEYgTABGIEwARiBMAEYgTABGIEwARiBMAEYgTABGIEwARiBMAEYgTABGIEwARiBMAEYgTABGIEwARiBMAEYgTABGIEwARiBMAEYgTABGIEwARiBMAEYgRAFvwXHZ9MV2Hpx2oAAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 352.675x117.056 with 1 Axes>"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "circuit_list[0].qc.draw(output='mpl')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[ 0.70710678+0.j, -0.70710678+0.j],\n",
+       "       [-0.70710678+0.j, -0.70710678+0.j]])"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "circuit_list[0].resultant_unitary"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "['z', 'x', 'h']"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "circuit_list[0].gate_names"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Performing Testing"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import timeit\n",
+    "from functools import partial\n",
+    "\n",
+    "from pyqasm.algorithms import solovay_kitaev\n",
+    "from pyqasm.maps.gates import BasisSet"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 34,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def measure_result_and_time(target_matrix_U, target_gate_set):\n",
+    "    result = solovay_kitaev(target_matrix_U, target_gate_set, depth=20)\n",
+    "    \n",
+    "    partial_func = partial(solovay_kitaev, target_matrix_U, target_gate_set, 20)\n",
+    "    avg_exc_time = timeit.timeit(partial_func, number=1)/1\n",
+    "    return result, avg_exc_time"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 35,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def test_and_analyze_sk_algo(circuit_list, target_gate_set):\n",
+    "    for circuit in circuit_list:\n",
+    "        result, avg_exc_time = measure_result_and_time(circuit.resultant_unitary, target_gate_set)\n",
+    "        result_distance = np.linalg.norm(result.matrix - circuit.resultant_unitary, 2)\n",
+    "        \n",
+    "        print(\"------------------------------------------------\")\n",
+    "        print(f\"Execution time:     {avg_exc_time}\")\n",
+    "        print(f\"Circuit gates:      {circuit.gate_names}\")\n",
+    "        print(f\"Decomposed gates:   {result.name}\")\n",
+    "        print(f\"Distance:           {result_distance}\")\n",
+    "        print(f\"Is within accuracy: {result_distance < 1e-6}\")\n",
+    "        print(\"Orignal resultant matrix:\")\n",
+    "        print(circuit.resultant_unitary)\n",
+    "        print(\"Decomposed resultant unitary:\")\n",
+    "        print(result.matrix)\n",
+    "        print(\"------------------------------------------------\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 36,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "------------------------------------------------\n",
+      "Execution time:     2.4317583000083687\n",
+      "Circuit gates:      ['z', 'x', 'h']\n",
+      "Decomposed gates:   ['t', 't', 't', 't', 'h', 't', 't', 's']\n",
+      "Distance:           2.7488925102312213e-16\n",
+      "Is within accuracy: True\n",
+      "Orignal resultant matrix:\n",
+      "[[ 0.70710678+0.j -0.70710678+0.j]\n",
+      " [-0.70710678+0.j -0.70710678+0.j]]\n",
+      "Decomposed resultant unitary:\n",
+      "[[ 0.70710678+0.00000000e+00j -0.70710678+0.00000000e+00j]\n",
+      " [-0.70710678-2.53429949e-17j -0.70710678-5.55111512e-17j]]\n",
+      "------------------------------------------------\n",
+      "------------------------------------------------\n",
+      "Execution time:     4.228998799997498\n",
+      "Circuit gates:      ['r', 'u1']\n",
+      "Decomposed gates:   ['h', 't', 't', 't', 't', 't', 's', 'h', 'sdg', 'tdg', 'tdg', 'tdg', 'tdg', 'tdg', 'h', 's', 's', 'h', 's', 's', 'h', 's', 's', 'h']\n",
+      "Distance:           0.4046729693889623\n",
+      "Is within accuracy: False\n",
+      "Orignal resultant matrix:\n",
+      "[[-0.8987577 +0.j         -0.31607208+0.30386352j]\n",
+      " [-0.14188982-0.41485163j  0.79903828+0.41146473j]]\n",
+      "Decomposed resultant unitary:\n",
+      "[[-0.85355339+0.35355339j -0.14644661+0.35355339j]\n",
+      " [-0.14644661-0.35355339j  0.85355339+0.35355339j]]\n",
+      "------------------------------------------------\n",
+      "------------------------------------------------\n",
+      "Execution time:     2.119061999997939\n",
+      "Circuit gates:      ['z', 'sdg', 'sdg', 't', 's', 'sx', 'y']\n",
+      "Decomposed gates:   ['s', 'h', 't', 's', 's', 'h', 's', 's', 'h', 's']\n",
+      "Distance:           2.027528834450315e-16\n",
+      "Is within accuracy: True\n",
+      "Orignal resultant matrix:\n",
+      "[[-0.5       -0.5j         0.        +0.70710678j]\n",
+      " [-0.5       +0.5j        -0.70710678+0.j        ]]\n",
+      "Decomposed resultant unitary:\n",
+      "[[-5.00000000e-01-5.00000000e-01j -1.11855716e-17+7.07106781e-01j]\n",
+      " [-5.00000000e-01+5.00000000e-01j -7.07106781e-01+1.11855716e-17j]]\n",
+      "------------------------------------------------\n",
+      "------------------------------------------------\n",
+      "Execution time:     2.345279599991045\n",
+      "Circuit gates:      ['r', 's', 'rz', 'sdg']\n",
+      "Decomposed gates:   ['t', 'h', 's', 's', 'h', 's', 's', 's', 'h', 't']\n",
+      "Distance:           0.27439900654353205\n",
+      "Is within accuracy: False\n",
+      "Orignal resultant matrix:\n",
+      "[[ 0.20618486-0.55826326j -0.59488154+0.54032018j]\n",
+      " [ 0.59488154+0.54032018j  0.20618486+0.55826326j]]\n",
+      "Decomposed resultant unitary:\n",
+      "[[-9.59550172e-18-0.70710678j -5.00000000e-01+0.5j       ]\n",
+      " [ 5.00000000e-01+0.5j         5.25894755e-16+0.70710678j]]\n",
+      "------------------------------------------------\n",
+      "------------------------------------------------\n",
+      "Execution time:     19.87103229999775\n",
+      "Circuit gates:      ['u3', 'ry', 's', 'ry', 'sx', 'ry', 'ry', 'sdg', 'u', 'x', 'rx', 'sxdg', 'id', 'ry', 'u', 'u3']\n",
+      "Decomposed gates:   ['s', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 't', 'h', 't', 'h', 't', 's', 'h']\n",
+      "Distance:           0.2770351557535887\n",
+      "Is within accuracy: False\n",
+      "Orignal resultant matrix:\n",
+      "[[ 0.62808441+0.45021745j  0.59411116+0.22326252j]\n",
+      " [ 0.29575769-0.56155286j -0.17038302+0.75376084j]]\n",
+      "Decomposed resultant unitary:\n",
+      "[[ 0.70710678+5.00000000e-01j  0.5       -9.06080907e-18j]\n",
+      " [ 0.35355339-3.53553391e-01j -0.14644661+8.53553391e-01j]]\n",
+      "------------------------------------------------\n",
+      "------------------------------------------------\n",
+      "Execution time:     2.138298299993039\n",
+      "Circuit gates:      ['u2', 'ry', 'h']\n",
+      "Decomposed gates:   ['h', 's', 's', 'h', 't', 'h', 't', 's', 'h', 's']\n",
+      "Distance:           0.3227556363764241\n",
+      "Is within accuracy: False\n",
+      "Orignal resultant matrix:\n",
+      "[[ 0.92162512+0.21404063j -0.32184159+0.0348101j ]\n",
+      " [ 0.09123887+0.31059494j  0.55102236+0.76914277j]]\n",
+      "Decomposed resultant unitary:\n",
+      "[[ 0.85355339+0.35355339j -0.35355339-0.14644661j]\n",
+      " [ 0.14644661+0.35355339j  0.35355339+0.85355339j]]\n",
+      "------------------------------------------------\n",
+      "------------------------------------------------\n",
+      "Execution time:     8.068621699989308\n",
+      "Circuit gates:      ['u']\n",
+      "Decomposed gates:   ['s', 'h', 't', 'h', 't', 'h', 's', 's', 't']\n",
+      "Distance:           0.3480364791371961\n",
+      "Is within accuracy: False\n",
+      "Orignal resultant matrix:\n",
+      "[[ 0.81342817+0.j         -0.05970368-0.5785932j ]\n",
+      " [-0.24774313+0.52626795j -0.76763181-0.26908508j]]\n",
+      "Decomposed resultant unitary:\n",
+      "[[ 8.53553391e-01+0.14644661j  3.64046306e-17-0.5j       ]\n",
+      " [-3.53553391e-01+0.35355339j -7.07106781e-01-0.5j       ]]\n",
+      "------------------------------------------------\n",
+      "------------------------------------------------\n",
+      "Execution time:     18.017034499993315\n",
+      "Circuit gates:      ['y', 'u1', 'r', 'h', 'z', 'u1', 'p', 'u', 'u2', 's', 'u', 'sx', 't', 'x', 'u2', 'y']\n",
+      "Decomposed gates:   ['s', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'tdg', 'tdg', 'tdg', 'tdg', 's', 't', 'sdg', 'tdg', 'h', 't', 't', 'h', 's', 's', 's', 'h', 't', 't']\n",
+      "Distance:           0.39865290368990064\n",
+      "Is within accuracy: False\n",
+      "Orignal resultant matrix:\n",
+      "[[-0.06007544-0.01807924j -0.70686923+0.70455658j]\n",
+      " [ 0.92164167-0.38293696j -0.0063209 +0.06241765j]]\n",
+      "Decomposed resultant unitary:\n",
+      "[[-5.19907294e-17+7.49842011e-17j -7.07106781e-01+7.07106781e-01j]\n",
+      " [ 7.07106781e-01-7.07106781e-01j  3.89002693e-16-3.66009221e-16j]]\n",
+      "------------------------------------------------\n",
+      "------------------------------------------------\n",
+      "Execution time:     23.155897600008757\n",
+      "Circuit gates:      ['z', 'tdg', 'z', 'u2', 't', 'id', 'sxdg', 'sx', 'id', 's', 'y', 'id', 'sx', 't', 'y', 'ry', 'rz']\n",
+      "Decomposed gates:   ['s', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 't', 'tdg', 'h', 's', 's', 'h', 's', 's', 's', 'h', 's', 't']\n",
+      "Distance:           0.45184583063972\n",
+      "Is within accuracy: False\n",
+      "Orignal resultant matrix:\n",
+      "[[-0.25037996-0.72070377j -0.17476521-0.62237696j]\n",
+      " [ 0.54982081-0.33998975j -0.67202971+0.36122031j]]\n",
+      "Decomposed resultant unitary:\n",
+      "[[-1.36185234e-17-7.07106781e-01j -5.00000000e-01-5.00000000e-01j]\n",
+      " [ 7.07106781e-01+2.32873552e-16j -5.00000000e-01+5.00000000e-01j]]\n",
+      "------------------------------------------------\n",
+      "------------------------------------------------\n",
+      "Execution time:     30.734558600001037\n",
+      "Circuit gates:      ['u', 'rz', 'p', 's', 'z', 'z', 'p', 'sx', 'tdg', 't', 'rx', 'rx', 'rx']\n",
+      "Decomposed gates:   ['s', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'sdg', 'sdg', 's', 's', 'tdg', 'tdg', 'tdg', 'tdg', 's', 't', 'tdg', 'tdg', 'tdg', 's', 't', 'sdg', 'tdg', 'h', 't', 't', 'h', 's', 'h', 's', 't', 'h', 's']\n",
+      "Distance:           0.23391882033249825\n",
+      "Is within accuracy: False\n",
+      "Orignal resultant matrix:\n",
+      "[[ 0.53716092+0.49705851j -0.4588246 +0.50385611j]\n",
+      " [-0.15799653+0.66289372j  0.71324392+0.16398819j]]\n",
+      "Decomposed resultant unitary:\n",
+      "[[ 5.00000000e-01+5.00000000e-01j -5.00000000e-01+5.00000000e-01j]\n",
+      " [-5.98455758e-16+7.07106781e-01j  7.07106781e-01+1.35729377e-16j]]\n",
+      "------------------------------------------------\n"
+     ]
+    }
+   ],
+   "source": [
+    "test_and_analyze_sk_algo(circuit_list, BasisSet.CLIFFORD_T)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "qbraid-venv",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/src/pyqasm/algorithms/solovay_kitaev/utils.py b/src/pyqasm/algorithms/solovay_kitaev/utils.py
new file mode 100644
index 00000000..540eebba
--- /dev/null
+++ b/src/pyqasm/algorithms/solovay_kitaev/utils.py
@@ -0,0 +1,165 @@
+"""
+Definition of common utility functions for the Solovay-Kitaev algorithm.
+"""
+
+from typing import List
+
+import numpy as np
+
+from pyqasm.elements import BasisSet
+from pyqasm.maps.gates import GATE_OPT_DATA, SELF_INVERTING_ONE_QUBIT_OP_SET, ST_GATE_INV_MAP
+
+
+class SU2Matrix:
+    """Class representing a 2x2 Special Unitary matrix.
+
+    Used for the Solovay-Kitaev algorithm.
+    """
+
+    def __init__(self, matrix: np.ndarray, name: List[str]):
+        self.matrix = matrix
+        self.name = name
+
+    def __mul__(self, other: "SU2Matrix") -> "SU2Matrix":
+        """Calculates the dot product of two matrices, and stores them with the combined name.
+
+        Args:
+            other (SU2Matrix): The other matrix.
+
+        Returns:
+            SU2Matrix: The product of the two matrices.
+        """
+        matrix = np.dot(self.matrix, other.matrix)
+        name = self.name.copy()
+        name.extend(other.name)
+        return SU2Matrix(matrix, name)
+
+    def dagger(self) -> "SU2Matrix":
+        """Returns the conjugate transpose."""
+        matrix = self.matrix.conj().T
+        name = []
+        for n in self.name[::-1]:
+            name.append(self._get_dagger_gate_name(n))
+
+        return SU2Matrix(matrix, name)
+
+    def distance(self, other: "SU2Matrix") -> float:
+        """Calculates the operator norm distance between two matrices."""
+        return np.linalg.norm(self.matrix - other.matrix, 2)
+
+    def _get_dagger_gate_name(self, name: str):
+        """Returns the name of the dagger gate."""
+        if name in SELF_INVERTING_ONE_QUBIT_OP_SET:
+            return name
+
+        return ST_GATE_INV_MAP[name]
+
+    def __str__(self):
+        return f"name: {self.name}, matrix: {self.matrix}"
+
+
+class TU2Matrix:
+    """Class representing a 2x2 Traversal Unitary matrix.
+
+    Used for basic approximation in the Solovay-Kitaev algorithm.
+    """
+
+    def __init__(
+        self, matrix: np.ndarray, name: List[str], identity_group: str, identity_weight: float
+    ):
+        self.matrix = matrix
+        self.name = name
+        self.identity_group = identity_group
+        self.identity_weight = identity_weight
+
+    def __mul__(self, other: "TU2Matrix") -> "TU2Matrix":
+        """Calculates the dot product of two matrices, and stores them with the combined name.
+        Adds the identity weight if the identity group is the same.
+        Updates the identity group and weight if the identity group is different.
+        Args:
+            other (TU2Matrix): The other matrix.
+
+        Returns:
+            TU2Matrix: The product of the two matrices.
+        """
+        matrix = np.dot(self.matrix, other.matrix)
+        name = self.name.copy()
+        name.extend(other.name)
+        identity_weight = 0
+        identity_group = ""
+        if self.identity_group == other.identity_group:
+            identity_weight = self.identity_weight + other.identity_weight
+        else:
+            identity_group = other.identity_group
+            identity_weight = other.identity_weight
+
+        return TU2Matrix(matrix, name, identity_group, identity_weight)
+
+    def can_multiple(self, other: "TU2Matrix"):
+        """Checks if the two matrices can be multiplied.
+        If the identity group is the same, the identity weight should be less than 1.
+        """
+        if self.identity_group != other.identity_group:
+            return True
+
+        return self.identity_weight + other.identity_weight < 1
+
+    def distance(self, other):
+        """Calculates the operator norm distance between two matrices."""
+        return np.linalg.norm(self.matrix - other, 2)
+
+    def __str__(self):
+        return f"""name: {self.name},
+            matrix: {self.matrix},
+            group: {self.identity_group},
+            weight: {self.identity_weight}"""
+
+    def copy(self) -> "TU2Matrix":
+        """Returns a copy of the current matrix."""
+        return TU2Matrix(
+            self.matrix.copy(), self.name.copy(), self.identity_group, self.identity_weight
+        )
+
+
+def get_su2matrix_for_solovay_kitaev_algorithm(target_basis_set) -> List[SU2Matrix]:
+    """Returns a list of SU2Matrix objects for the given basis set.
+    This list is used for the Solovay-Kitaev algorithm.
+    """
+    gate_list = GATE_OPT_DATA[target_basis_set]
+    return [SU2Matrix(gate["matrix"], [gate["name"]]) for gate in gate_list]
+
+
+def get_tu2matrix_for_basic_approximation(target_basis_set) -> List[TU2Matrix]:
+    """Returns a list of TU2Matrix objects for the given basis set.
+    This list is used for the basic approximation algorithm.
+    """
+    whole_gate_list = GATE_OPT_DATA[target_basis_set]
+    required_gate_list = [gate for gate in whole_gate_list if gate["used_for_basic_approximation"]]
+
+    return [
+        TU2Matrix(
+            gate["matrix"], [gate["name"]], gate["identity"]["group"], gate["identity"]["weight"]
+        )
+        for gate in required_gate_list
+    ]
+
+
+def get_identity_weight_group_for_optimizer(target_basis_set):
+    """Returns the identity weight group for the given basis set.
+    This is used for the optimization of the gate sequence.
+    """
+    gate_list = GATE_OPT_DATA[target_basis_set]
+    identity_weight_group = {}
+
+    for gate in gate_list:
+        identity_weight_group[gate["name"]] = {
+            "group": gate["identity"]["group"],
+            "weight": gate["identity"]["weight"],
+        }
+
+    return identity_weight_group
+
+
+if __name__ == "__main__":
+    result = get_identity_weight_group_for_optimizer(BasisSet.CLIFFORD_T)
+    print(result)
diff --git a/src/pyqasm/decomposer.py b/src/pyqasm/decomposer.py
index fd9b337d..8c612e41 100644
--- a/src/pyqasm/decomposer.py
+++ b/src/pyqasm/decomposer.py
@@ -15,9 +15,15 @@
 import openqasm3.ast as qasm3_ast
 from openqasm3.ast import BranchingStatement, QuantumGate
 
+from pyqasm.algorithms import solovay_kitaev
 from pyqasm.exceptions import RebaseError
-from pyqasm.maps.decomposition_rules import DECOMPOSITION_RULES, AppliedQubit
-from pyqasm.maps.gates import BASIS_GATE_MAP
+from pyqasm.maps.decomposition_rules import (
+    DECOMPOSITION_RULES,
+    ROTATIONAL_LOOKUP_RULES,
+    AppliedQubit,
+)
+from pyqasm.maps.expressions import CONSTANTS_MAP
+from pyqasm.maps.gates import BASIS_GATE_MAP, get_target_matrix_for_rotational_gates
 
 
 class Decomposer:
@@ -26,13 +32,17 @@ class Decomposer:
     """
 
     @classmethod
-    def process_gate_statement(cls, gate_name, statement, target_basis_set):
+    def process_gate_statement(
+        cls, gate_name, statement, target_basis_set, depth=10, accuracy=1e-6
+    ):
         """Process the gate statement based on the target basis set.
 
         Args:
             gate_name: The name of the gate to process.
             statement: The statement to process.
             target_basis_set: The target basis set to rebase the module to.
+            depth: The depth of the approximation.
+            accuracy: The accuracy of the approximation.
 
         Returns:
             list: The processed gates based on the target basis set.
@@ -47,17 +57,12 @@ def process_gate_statement(cls, gate_name, statement, target_basis_set):
             processed_gates_list = [statement]
         elif gate_name in decomposition_rules:
             # Decompose the gates
-            processed_gates_list = cls._get_decomposed_gates(
-                decomposition_rules, statement, gate_name
-            )
+            rule_list = decomposition_rules[gate_name]
+            processed_gates_list = cls._get_decomposed_gates(rule_list, statement)
         elif gate_name in {"rx", "ry", "rz"}:
-            # Approximate parameterized gates using Solovay-Kitaev
-            # Example -
-            # approx_gates = solovay_kitaev_algo(
-            #     gate_name, statement.arguments[0].value, accuracy=0.01
-            # )
-            # return approx_gates
-            pass
+            processed_gates_list = cls._process_rotational_gate(
+                gate_name, statement, target_basis_set, depth, accuracy
+            )
         else:
             # Raise an error if the gate is not supported in the target basis set
             error = f"Gate '{gate_name}' is not supported in the '{target_basis_set} set'."
@@ -66,12 +71,14 @@ def process_gate_statement(cls, gate_name, statement, target_basis_set):
         return processed_gates_list
 
     @classmethod
-    def process_branching_statement(cls, branching_statement, target_basis_set):
+    def process_branching_statement(cls, branching_statement, target_basis_set, depth, accuracy):
         """Process the branching statement based on the target basis set.
 
         Args:
             branching_statement: The branching statement to process.
             target_basis_set: The target basis set to rebase the module to.
+            depth: The depth of the approximation.
+            accuracy: The accuracy of the approximation.
 
         Returns:
             BranchingStatement: The processed branching statement based on the target basis set.
@@ -87,7 +94,9 @@ def process_branching_statement(cls, branching_statement, target_basis_set):
                 )
                 if_block.extend(processed_gates_list)
             elif isinstance(statement, BranchingStatement):
-                if_block.append(cls.process_branching_statement(statement, target_basis_set))
+                if_block.append(
+                    cls.process_branching_statement(statement, target_basis_set, depth, accuracy)
+                )
             else:
                 if_block.append(statement)
 
@@ -99,7 +108,9 @@ def process_branching_statement(cls, branching_statement, target_basis_set):
                 )
                 else_block.extend(processed_gates_list)
             elif isinstance(statement, BranchingStatement):
-                else_block.append(cls.process_branching_statement(statement, target_basis_set))
+                else_block.append(
+                    cls.process_branching_statement(statement, target_basis_set, depth, accuracy)
+                )
             else:
                 else_block.append(statement)
 
@@ -108,26 +119,25 @@ def process_branching_statement(cls, branching_statement, target_basis_set):
         )
 
     @classmethod
-    def _get_decomposed_gates(cls, decomposition_rules, statement, gate):
+    def _get_decomposed_gates(cls, rule_list, statement):
         """Apply the decomposed gates based on the decomposition rules.
 
         Args:
-            decomposition_rules: The decomposition rules to apply.
+            rule_list: The decomposition rules to apply.
             statement: The statement to apply the decomposition rules to.
-            gate: The name of the gate to apply the decomposition rules to.
 
         Returns:
             list: The decomposed gates to be applied.
         """
         decomposed_gates = []
 
-        for rule in decomposition_rules[gate]:
+        for rule in rule_list:
             qubits = cls._get_qubits_for_gate(statement.qubits, rule)
             arguments = [qasm3_ast.FloatLiteral(value=rule["param"])] if "param" in rule else []
 
             new_gate = qasm3_ast.QuantumGate(
                 modifiers=[],
-                name=qasm3_ast.Identifier(name=rule["gate"]),
+                name=qasm3_ast.Identifier(name=rule if isinstance(rule, str) else rule["gate"]),
                 arguments=arguments,
                 qubits=qubits,
             )
@@ -159,3 +169,42 @@ def _get_qubits_for_gate(cls, qubits, rule):
             else:
                 qubits = [qubits[1]]
         return qubits
+
+    @classmethod
+    def _process_rotational_gate(cls, gate_name, statement, target_basis_set, depth, accuracy):
+        """Process the rotational gates based on the target basis set.
+
+        Args:
+            gate_name: The name of the gate.
+            statement: The statement to process.
+            target_basis_set: The target basis set to rebase the module to.
+            depth: The depth of the approximation.
+            accuracy: The accuracy of the approximation.
+
+        Returns:
+            list: The processed gates based on the target basis set.
+        """
+        theta = statement.arguments[0].value
+        processed_gates_list = []
+
+        # Use lookup table if ∅ is pi, pi/2 or pi/4
+        if theta in [CONSTANTS_MAP["pi"], CONSTANTS_MAP["pi"] / 2, CONSTANTS_MAP["pi"] / 4]:
+            rotational_lookup_rules = ROTATIONAL_LOOKUP_RULES[target_basis_set]
+            rule_list = rotational_lookup_rules[gate_name][theta]
+            processed_gates_list = cls._get_decomposed_gates(rule_list, statement)
+
+        # Use Solovay-Kitaev's Algorithm for gate approximation
+        else:
+            target_matrix = get_target_matrix_for_rotational_gates(gate_name, theta)
+            approximated_gates = solovay_kitaev(target_matrix, target_basis_set, depth, accuracy)
+
+            for approximated_gate_name in approximated_gates.name:
+                new_gate = qasm3_ast.QuantumGate(
+                    modifiers=[],
+                    name=qasm3_ast.Identifier(name=approximated_gate_name),
+                    arguments=[],
+                    qubits=statement.qubits,
+                )
+
+                processed_gates_list.append(new_gate)
+        return processed_gates_list
diff --git a/src/pyqasm/maps/decomposition_rules.py b/src/pyqasm/maps/decomposition_rules.py
index b2421eb1..e9e427e5 100644
--- a/src/pyqasm/maps/decomposition_rules.py
+++ b/src/pyqasm/maps/decomposition_rules.py
@@ -96,6 +96,28 @@ class AppliedQubit(Enum):
     },
 }
 
+ROTATIONAL_LOOKUP_RULES = {
+    BasisSet.CLIFFORD_T: {
+        "rz": {
+            CONSTANTS_MAP["pi"]: ["s", "s"],
+            CONSTANTS_MAP["pi"] / 2: ["s"],
+            CONSTANTS_MAP["pi"] / 4: ["t"],
+        },
+        # Rx(∅) = H.Rz(∅).H
+        "rx": {
+            CONSTANTS_MAP["pi"]: ["h", "s", "s", "h"],
+            CONSTANTS_MAP["pi"] / 2: ["h", "s", "h"],
+            CONSTANTS_MAP["pi"] / 4: ["h", "t", "h"],
+        },
+        # Ry(∅) = S†.H.Rz(∅).H.S
+        "ry": {
+            CONSTANTS_MAP["pi"]: ["sdg", "h", "s", "s", "h", "s"],
+            CONSTANTS_MAP["pi"] / 2: ["sdg", "h", "s", "h", "s"],
+            CONSTANTS_MAP["pi"] / 4: ["sdg", "h", "t", "h", "s"],
+        },
+    }
+}
+
 # """TODO: Implement the Solovay-Kitaev algorithm"""
 #
 # def solovay_kitaev_algo(gate_name, param, accuracy):
diff --git a/src/pyqasm/maps/gates.py b/src/pyqasm/maps/gates.py
index 8228916a..b4f92b04 100644
--- a/src/pyqasm/maps/gates.py
+++ b/src/pyqasm/maps/gates.py
@@ -1111,6 +1111,56 @@ def two_qubit_gate_op(
     BasisSet.CLIFFORD_T: {"h", "t", "s", "cx", "tdg", "sdg"},
 }
 
+# Gate Optimization Data is used by Solovay-Kitaev algorithm to optimize the gate set
+# The data is a dictionary with the basis set as the key and a list of dictionaries as the value
+# Each dictionary in the list contains the following keys:
+# - name: the name of the gate
+# - identity:
+#       The dot product of sequence of some gates becomes Identity matrix.
+#       The `identity` key is used to determine such sequences of gates.
+#       It stores the following keys:
+#           - group: the group to which the gate belongs.
+#           - weight: the weightage given to that gate.
+#       If gates belong to smae group, their weights can be added.
+#       Once the total weight reaches 1, the resultant dot product of sequnce
+#       becomes Identity Matrix
+# - matrix: the matrix representation of the gate
+# - used_for_basic_approximation: whether the gate is used for basic approximation
+GATE_OPT_DATA = {
+    BasisSet.CLIFFORD_T: [
+        {
+            "name": "h",
+            "identity": {"group": "h", "weight": 0.5},
+            "matrix": (1 / np.sqrt(2)) * np.array([[1, 1], [1, -1]]),
+            "used_for_basic_approximation": True,
+        },
+        {
+            "name": "s",
+            "identity": {"group": "s-t", "weight": 0.25},
+            "matrix": np.array([[1, 0], [0, 1j]]),
+            "used_for_basic_approximation": True,
+        },
+        {
+            "name": "t",
+            "identity": {"group": "s-t", "weight": 0.125},
+            "matrix": np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]]),
+            "used_for_basic_approximation": True,
+        },
+        {
+            "name": "sdg",
+            "identity": {"group": "sdg-tdg", "weight": 0.25},
+            "matrix": np.array([[1, 0], [0, 1j]]).conj().T,
+            "used_for_basic_approximation": False,
+        },
+        {
+            "name": "tdg",
+            "identity": {"group": "sdg-tdg", "weight": 0.125},
+            "matrix": np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]]).conj().T,
+            "used_for_basic_approximation": False,
+        },
+    ]
+}
+
 
 def map_qasm_op_to_callable(op_name: str) -> tuple[Callable, int]:
     """
@@ -1236,3 +1286,31 @@ def map_qasm_ctrl_op_to_callable(op_name: str, ctrl_count: int):
     raise ValidationError(
         f"Unsupported controlled QASM operation: {op_name} with {ctrl_count} controls"
     )
+
+
+def get_target_matrix_for_rotational_gates(gate_name, theta):
+    """
+    Get the target matrix for the rotational gates based on the gate name and theta.
+
+    Args:
+        gate_name: The name of the gate.
+        theta: The angle of rotation.
+
+    Returns:
+        np.ndarray: The target matrix for the rotational gates.
+    """
+    if gate_name == "rx":
+        target_matrix = np.array(
+            [
+                [np.cos(theta / 2), -1j * np.sin(theta / 2)],
+                [-1j * np.sin(theta / 2), np.cos(theta / 2)],
+            ]
+        )
+    elif gate_name == "ry":
+        target_matrix = np.array(
+            [[np.cos(theta / 2), -np.sin(theta / 2)], [np.sin(theta / 2), np.cos(theta / 2)]]
+        )
+    else:
+        target_matrix = np.array([[np.exp(-1j * theta / 2), 0], [0, np.exp(1j * theta / 2)]])
+
+    return target_matrix
diff --git a/src/pyqasm/modules/base.py b/src/pyqasm/modules/base.py
index c0a796d0..1cb2dce7 100644
--- a/src/pyqasm/modules/base.py
+++ b/src/pyqasm/modules/base.py
@@ -527,7 +527,7 @@ def unroll(self, **kwargs):
             self._unrolled_ast = Program(statements=[], version=self.original_program.version)
             raise err
 
-    def rebase(self, target_basis_set, in_place=True):
+    def rebase(self, target_basis_set, in_place=True, depth=10, accuracy=1e-6):
         """Rebase the AST to use a specified target basis set.
 
         Will unroll the module if not already done.
@@ -535,6 +535,8 @@ def rebase(self, target_basis_set, in_place=True):
         Args:
             target_basis_set: The target basis set to rebase the module to.
             in_place (bool): Flag to indicate if the rebase operation should be done in place.
+            depth: The depth of the approximation.
+            accuracy: The accuracy of the approximation.
 
         Returns:
             QasmModule: The module with the gates rebased to the target basis set.
@@ -555,7 +557,7 @@ def rebase(self, target_basis_set, in_place=True):
 
                 # Decompose the gate
                 processed_gates_list = Decomposer.process_gate_statement(
-                    gate_name, statement, target_basis_set
+                    gate_name, statement, target_basis_set, depth, accuracy
                 )
                 for processed_gate in processed_gates_list:
                     rebased_statements.append(processed_gate)
@@ -564,7 +566,9 @@ def rebase(self, target_basis_set, in_place=True):
                 # Recursively process the if_block and else_block
 
                 rebased_statements.append(
-                    Decomposer.process_branching_statement(statement, target_basis_set)
+                    Decomposer.process_branching_statement(
+                        statement, target_basis_set, depth, accuracy
+                    )
                 )
 
             else:
diff --git a/tests/qasm3/test_rebase.py b/tests/qasm3/test_rebase.py
index 7407ea9e..a4f1538b 100644
--- a/tests/qasm3/test_rebase.py
+++ b/tests/qasm3/test_rebase.py
@@ -139,6 +139,81 @@ def test_rebase_rotational_cx(input_gates, decomposed_gates):
             cx q[0], q[1];
             """,
         ),
+        (
+            "rx(pi) q[0];",
+            """
+            h q[0];
+            s q[0];
+            s q[0];
+            h q[0];
+            """,
+        ),
+        (
+            "rx(pi/2) q[0];",
+            """
+            h q[0];
+            s q[0];
+            h q[0];
+            """,
+        ),
+        (
+            "rx(pi/4) q[0];",
+            """
+            h q[0];
+            t q[0];
+            h q[0];
+            """,
+        ),
+        (
+            "ry(pi) q[0];",
+            """
+            sdg q[0];
+            h q[0];
+            s q[0];
+            s q[0];
+            h q[0];
+            s q[0];
+            """,
+        ),
+        (
+            "ry(pi/2) q[0];",
+            """
+            sdg q[0];
+            h q[0];
+            s q[0];
+            h q[0];
+            s q[0];
+            """,
+        ),
+        (
+            "ry(pi/4) q[0];",
+            """
+            sdg q[0];
+            h q[0];
+            t q[0];
+            h q[0];
+            s q[0];
+            """,
+        ),
+        (
+            "rz(pi) q[0];",
+            """
+            s q[0];
+            s q[0];
+            """,
+        ),
+        (
+            "rz(pi/2) q[0];",
+            """
+            s q[0];
+            """,
+        ),
+        (
+            "rz(pi/4) q[0];",
+            """
+            t q[0];
+            """,
+        ),
     ],
 )
 def test_rebase_clifford_t(input_gates, decomposed_gates):