import torch


def quat_mul(q1: torch.Tensor, q2: torch.Tensor):
    """Computes qout = q1 * q2, where * is the Hamilton product between the two
    quaternions. Note that the Hamiltonian product is not commutative.

    Args:
        q1: Quaternions of shape ([*, ], 4)
        q2: Quaternions of shape ([*, ], 4)

    Returns:
        qout = q1*q2.
    """
    assert q1.shape[-1] == 4 and q2.shape[-1] == 4

    qout = torch.stack([
        q1[..., 0] * q2[..., 0] - q1[..., 1] * q2[..., 1] - q1[..., 2] * q2[..., 2] - q1[..., 3] * q2[..., 3],
        q1[..., 0] * q2[..., 1] + q1[..., 1] * q2[..., 0] + q1[..., 2] * q2[..., 3] - q1[..., 3] * q2[..., 2],
        q1[..., 0] * q2[..., 2] - q1[..., 1] * q2[..., 3] + q1[..., 2] * q2[..., 0] + q1[..., 3] * q2[..., 1],
        q1[..., 0] * q2[..., 3] + q1[..., 1] * q2[..., 2] - q1[..., 2] * q2[..., 1] + q1[..., 3] * q2[..., 0]
    ], dim=-1)
    return qout