docs: def kak_canonicalize_vector(x: float, y: float, z: float) -> KakDecomposition:
      """Canonicalizes an XX/YY/ZZ interaction by swap/negate/shift-ing axes.
  
      Args:
          x: The strength of the XX interaction.
          y: The strength of the YY interaction.
          z: The strength of the ZZ interaction.
  
      Returns:
          The canonicalized decomposition, with vector coefficients (x2, y2, z2)
          satisfying:
  
              0 ≤ abs(z2) ≤ y2 ≤ x2 ≤ π/4
              z2 ≠ -π/4
  
          Guarantees that the implied output matrix:
  
              g · (a1 ⊗ a0) · exp(i·(x2·XX + y2·YY + z2·ZZ)) · (b1 ⊗ b0)
  
          is approximately equal to the implied input matrix:
  
              exp(i·(x·XX + y·YY + z·ZZ))
      """
  
      phase = [complex(1)]  # Accumulated global phase.
      left = [np.eye(2)] * 2  # Per-qubit left factors.
      right = [np.eye(2)] * 2  # Per-qubit right factors.
      v = [x, y, z]  # Remaining XX/YY/ZZ interaction vector.
  
      # These special-unitary matrices flip the X, Y, and Z axes respectively.
      flippers = [
          np.array([[0, 1], [1, 0]]) * 1j,
          np.array([[0, -1j], [1j, 0]]) * 1j,
          np.array([[1, 0], [0, -1]]) * 1j
      ]
  
      # Each of these special-unitary matrices swaps two the roles of two axes.
      # The matrix at index k swaps the *other two* axes (e.g. swappers[1] is a
      # Hadamard operation that swaps X and Z).
      swappers = [
          np.array([[1, -1j], [1j, -1]]) * 1j * np.sqrt(0.5),
          np.array([[1, 1], [1, -1]]) * 1j * np.sqrt(0.5),
          np.array([[0, 1 - 1j], [1 + 1j, 0]]) * 1j * np.sqrt(0.5)
      ]
  
      # Shifting strength by ½π is equivalent to local ops (e.g. exp(i½π XX)∝XX).
      def shift(k, step):
          v[k] += step * np.pi / 2
          phase[0] *= 1j ** step
          right[0] = combinators.dot(flippers[k] ** (step % 4), right[0])
          right[1] = combinators.dot(flippers[k] ** (step % 4), right[1])
  
      # Two negations is equivalent to temporarily flipping along the other axis.
      def negate(k1, k2):
          v[k1] *= -1
          v[k2] *= -1
          phase[0] *= -1
          s = flippers[3 - k1 - k2]  # The other axis' flipper.
          left[1] = combinators.dot(left[1], s)
          right[1] = combinators.dot(s, right[1])
  
      # Swapping components is equivalent to temporarily swapping the two axes.
      def swap(k1, k2):
          v[k1], v[k2] = v[k2], v[k1]
          s = swappers[3 - k1 - k2]  # The other axis' swapper.
          left[0] = combinators.dot(left[0], s)
          left[1] = combinators.dot(left[1], s)
          right[0] = combinators.dot(s, right[0])
          right[1] = combinators.dot(s, right[1])
  
      # Shifts an axis strength into the range (-π/4, π/4].
      def canonical_shift(k):
          while v[k] <= -np.pi / 4:
              shift(k, +1)
          while v[k] > np.pi / 4:
              shift(k, -1)
  
      # Sorts axis strengths into descending order by absolute magnitude.
      def sort():
          if abs(v[0]) < abs(v[1]):
              swap(0, 1)
          if abs(v[1]) < abs(v[2]):
              swap(1, 2)
          if abs(v[0]) < abs(v[1]):
              swap(0, 1)
  
      # Get all strengths to (-¼π, ¼π] in descending order by absolute magnitude.
      canonical_shift(0)
      canonical_shift(1)
      canonical_shift(2)
      sort()
  
      # Move all negativity into z.
      if v[0] < 0:
          negate(0, 2)
      if v[1] < 0:
          negate(1, 2)
      canonical_shift(2)
  
      return KakDecomposition(
          global_phase=phase[0],
          single_qubit_operations_after=(left[1], left[0]),
          interaction_coefficients=(v[0], v[1], v[2]),
          single_qubit_operations_before=(right[1], right[0]))