code: def zyz_decomposition(gate: Gate) -> Circuit:
      
      if gate.qubit_nb != 1:
          raise ValueError('Expected 1-qubit gate')
  
      q, = gate.qubits
  
      U = asarray(gate.asoperator())
      U /= np.linalg.det(U) ** (1/2)    # SU(2)
  
      if abs(U[0, 0]) > abs(U[1, 0]):
          theta1 = 2 * np.arccos(min(abs(U[0, 0]), 1))
      else:
          theta1 = 2 * np.arcsin(min(abs(U[1, 0]), 1))
  
      cos_halftheta1 = np.cos(theta1/2)
      if not np.isclose(cos_halftheta1, 0.0):
          phase = U[1, 1] / cos_halftheta1
          theta0_plus_theta2 = 2 * np.arctan2(np.imag(phase), np.real(phase))
      else:
          theta0_plus_theta2 = 0.0
  
      sin_halftheta1 = np.sin(theta1/2)
      if not np.isclose(sin_halftheta1, 0.0):
          phase = U[1, 0] / sin_halftheta1
          theta0_sub_theta2 = 2 * np.arctan2(np.imag(phase), np.real(phase))
      else:
          theta0_sub_theta2 = 0.0
  
      theta0 = (theta0_plus_theta2 + theta0_sub_theta2) / 2
      theta2 = (theta0_plus_theta2 - theta0_sub_theta2) / 2
  
      t0 = theta0/np.pi
      t1 = theta1/np.pi
      t2 = theta2/np.pi
  
      circ1 = Circuit()
      circ1 += TZ(t2, q)
      circ1 += TY(t1, q)
      circ1 += TZ(t0, q)
  
      return circ1
docs: Returns the Euler Z-Y-Z decomposition of a local 1-qubit gate.