code: def sk_decomposition(
      op,
      epsilon,
      *,
      max_depth=5,
      basis_set=("T", "T*", "H"),
      basis_length=10,
  ):
      r"""Approximate an arbitrary single-qubit gate in the Clifford+T basis using the `Solovay-Kitaev algorithm <https://arxiv.org/abs/quant-ph/0505030>`_.
  
      This method implements the Solovay-Kitaev decomposition algorithm that approximates any single-qubit
      operation with :math:`\epsilon > 0` error. The procedure exits when the approximation error
      becomes less than :math:`\epsilon`, or when ``max_depth`` approximation passes have been made. In the
      latter case, the approximation error could be :math:`\geq \epsilon`.
  
      This algorithm produces a decomposition with :math:`O(\text{log}^{3.97}(1/\epsilon))` operations.
  
      Args:
          op (~pennylane.operation.Operation): A single-qubit gate operation.
          epsilon (float): The maximum permissible error.
  
      Keyword Args:
          max_depth (int): The maximum number of approximation passes. A smaller :math:`\epsilon` would generally require
              a greater number of passes. Default is ``5``.
          basis_set (list[str]): Basis set to be used for the decomposition and building an approximate set internally.
              It accepts the following gate terms: ``['X', 'Y', 'Z', 'H', 'T', 'T*', 'S', 'S*']``, where ``*`` refers
              to the gate adjoint. Default value is ``['T', 'T*', 'H']``.
          basis_length (int): Maximum expansion length of Clifford+T sequences in the internally-built approximate set.
              Default is ``10``.
  
      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 acts on more than one wires.
  
      **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)
  
          # Get the gate decomposition in ['T', 'T*', 'H']
          ops = qml.ops.sk_decomposition(op, epsilon=1e-3)
  
          # Get the approximate matrix from the ops
          matrix_sk = qml.prod(*reversed(ops)).matrix()
  
      When the function is run for a sufficient ``depth`` with a good enough approximate set,
      the output gate sequence should implement the same operation approximately.
  
      >>> qml.math.allclose(op.matrix(), matrix_sk, atol=1e-3)
      True
  
      """
      # Check for length of wires in the operation
      if len(op.wires) != 1:
          raise ValueError(
              f"Operator must be a single qubit operation, got {op} acting on {op.wires} wires."
          )
  
      with QueuingManager.stop_recording():
          # Build the approximate set with caching
          (
              approx_set_ids,
              approx_set_mat,
              approx_set_gph,
              approx_set_qat,
          ) = _approximate_set(basis_set, max_length=basis_length)
  
          # Build the k-d tree with the current approximation set for querying in the base case
          kd_tree = KDTree(qml.math.array(approx_set_qat))
  
          # Obtain the SU(2) and quaternion for the operation
          op_matrix = op.matrix()
          interface = qml.math.get_deep_interface(op_matrix)
          gate_mat, gate_gph = _SU2_transform(
              qml.math.unwrap(op_matrix)
          )
          gate_qat = _quaternion_transform(gate_mat)
  
          def _solovay_kitaev(umat, n, u_n1_ids, u_n1_mat):
              """Recursive method as given in the Section 3 of arXiv:0505030"""
  
              if not n:
                  # Check the approximate gate in our approximate set
                  seq_node = qml.math.array(
                      [_quaternion_transform(umat)]
                  )
                  _, [index] = kd_tree.query(seq_node, workers=-1)
                  return approx_set_ids[index], approx_set_mat[index]
  
              # Get the decomposition for the remaining unitary: U @ U'.adjoint()
              v_n, w_n = _group_commutator_decompose(
                  umat @ qml.math.conj(qml.math.transpose(u_n1_mat))
              )
  
              # Get the approximation for the residual commutator unitaries: V and W
              c_ids_mats = []
              for c_n in [v_n, w_n]:
                  # Get the approximation for each commutator iteratively: C --> C'
                  c_n1_ids, c_n1_mat = None, None
                  for i in range(n):
                      c_n1_ids, c_n1_mat = _solovay_kitaev(
                          c_n, i, c_n1_ids, c_n1_mat
                      )
  
                  # Get the adjoints C' --> C'.adjoint()
                  c_n1_ids_adj = [
                      qml.adjoint(gate, lazy=False)
                      for gate in reversed(c_n1_ids)
                  ]
                  c_n1_mat_adj = qml.math.conj(
                      qml.math.transpose(c_n1_mat)
                  )
  
                  # Store the decompositions and matrices for C'
                  c_ids_mats.append(
                      [c_n1_ids, c_n1_mat, c_n1_ids_adj, c_n1_mat_adj]
                  )
  
              # Get the V' and W'
              (
                  v_n1_ids,
                  v_n1_mat,
                  v_n1_ids_adj,
                  v_n1_mat_adj,
              ) = c_ids_mats[0]
              (
                  w_n1_ids,
                  w_n1_mat,
                  w_n1_ids_adj,
                  w_n1_mat_adj,
              ) = c_ids_mats[1]
  
              # Build the operations and their SU(2) for return
              approx_ids = (
                  u_n1_ids
                  + w_n1_ids_adj
                  + v_n1_ids_adj
                  + w_n1_ids
                  + v_n1_ids
              )
              approx_mat = (
                  v_n1_mat
                  @ w_n1_mat
                  @ v_n1_mat_adj
                  @ w_n1_mat_adj
                  @ u_n1_mat
              )
  
              return approx_ids, approx_mat
  
          # If we have it already, use that, otherwise proceed for decomposition
          _, [index] = kd_tree.query(
              qml.math.array([gate_qat]), workers=-1
          )
          decomposition, u_prime = (
              approx_set_ids[index],
              approx_set_mat[index],
          )
  
          # Iterate until max_depth while doing an epsilon-error comparision
          for depth in range(max_depth):
              # For a SU(2) matrix Hilbert-Schmidt norm is √(|α|^2 + |β|^2),
              # which is simply the L2-norm for the first row of the matrix.
              if qml.math.norm(gate_mat[0] - u_prime[0]) <= epsilon:
                  break
              # Approximate the residual with the approximation from previous iteration
              decomposition, u_prime = _solovay_kitaev(
                  gate_mat, depth + 1, decomposition, u_prime
              )
  
          # Remove inverses if any in the decomposition and handle trivial case
          [new_tape], _ = qml.transforms.cancel_inverses(
              QuantumScript(decomposition or [qml.Identity(0)])
          )
  
      # Map the wires to that of the operation and queue
      [map_tape], _ = qml.map_wires(
          new_tape, wire_map={0: op.wires[0]}, queue=True
      )
  
      # Get phase information based on the decomposition effort
      phase = (
          approx_set_gph[index] - gate_gph
          if depth or qml.math.allclose(gate_gph, 0.0)
          else 0.0
      )
      global_phase = qml.GlobalPhase(
          qml.math.array(phase, like=interface)
      )
  
      # Return the gates from the mapped tape and global phase
      return map_tape.operations + [global_phase]
  
file_path: pennylane/pennylane/ops/op_math/decompositions/solovay_kitaev.py
imports: import math
  import warnings
  from functools import lru_cache
  from scipy.spatial import KDTree
  import pennylane
  from pennylane.queuing import QueuingManager
  from pennylane.tape import QuantumScript
label_desc: Use image node_id 6 for calling a global function with example usage: sk_decomposition(op, epsilon) and returns: unknown, approx_ids,  approx_mat, approx_set_ids,  approx_set_mat