code: def kronecker_decomposition(gate: Gate) -> Circuit:
      
      # An alternative approach would be to take partial traces, but
      # this approach appears to be more robust.
  
      if gate.qubit_nb != 2:
          raise ValueError('Expected 2-qubit gate')
  
      U = asarray(gate.asoperator())
      rank = 2**gate.qubit_nb
      U /= np.linalg.det(U) ** (1/rank)
  
      R = np.stack([U[0:2, 0:2].reshape(4),
                    U[0:2, 2:4].reshape(4),
                    U[2:4, 0:2].reshape(4),
                    U[2:4, 2:4].reshape(4)])
      u, s, vh = np.linalg.svd(R)
      v = vh.transpose()
      A = (np.sqrt(s[0]) * u[:, 0]).reshape(2, 2)
      B = (np.sqrt(s[0]) * v[:, 0]).reshape(2, 2)
  
      q0, q1 = gate.qubits
      g0 = Gate(A, qubits=[q0])
      g1 = Gate(B, qubits=[q1])
  
      if not gates_close(gate, Circuit([g0, g1]).asgate()):
          raise ValueError("Gate cannot be decomposed into two 1-qubit gates")
  
      circ = Circuit()
      circ += zyz_decomposition(g0)
      circ += zyz_decomposition(g1)
  
      assert gates_close(gate, circ.asgate())  # Sanity check
  
      return circ
docs: Decompose a 2-qubit unitary composed of two 1-qubit local gates.
  
      Uses the "Nearest Kronecker Product" algorithm. Will give erratic
      results if the gate is not the direct product of two 1-qubit gates.