Returns: A randomly sampled Density from the Hilbert–Schmidt
                ensemble of quantum states

    Ref: "Induced measures in the space of mixed quantum states" Karol
        Zyczkowski, Hans-Juergen Sommers, J. Phys. A34, 7111-7125 (2001)
        https://arxiv.org/abs/quant-ph/0012101

def random_density(qubits: Union[int, Qubits]) -> Density:
    """
    Returns: A randomly sampled Density from the Hilbert–Schmidt
                ensemble of quantum states

    Ref: "Induced measures in the space of mixed quantum states" Karol
        Zyczkowski, Hans-Juergen Sommers, J. Phys. A34, 7111-7125 (2001)
        https://arxiv.org/abs/quant-ph/0012101
    """
    N, qubits = qubits_count_tuple(qubits)
    size = (2**N, 2**N)
    ginibre_ensemble = (np.random.normal(size=size) +
                        1j * np.random.normal(size=size)) / np.sqrt(2.0)
    matrix = ginibre_ensemble @ np.transpose(np.conjugate(ginibre_ensemble))
    matrix /= np.trace(matrix)

    return Density(matrix, qubits=qubits)