  The complexity of quantum many-body problems scales exponentially with the
size of the system, rendering any finite size scaling analysis a formidable
challenge. This is particularly true for methods based on the full
representation of the wave function, where one simply accepts the enormous
Hilbert space dimensions and performs linear algebra operations, e.g., for
finding the ground state of the Hamiltonian. If the system satisfies an
underlying symmetry where an operator with degenerate spectrum commutes with
the Hamiltonian, it can be block-diagonalized, thus reducing the complexity at
the expense of additional bookkeeping. At the most basic level, required for
Krylov space techniques (like the Lanczos algorithm) it is necessary to
implement a matrix-vector product of a block of the Hamiltonian with arbitrary
block-wavefunctions, potentially without holding the Hamiltonian block in
memory. An efficient implementation of this operation requires the calculation
of the position of an arbitrary basis vector in the canonical ordering of the
basis of the block. We present here an elegant and powerful, multi-dimensional
approach to this problem for the $U(1)$ symmetry appearing in problems with
particle number conservation. Our divide-and-conquer algorithm uses multiple
subsystems and hence generalizes previous approaches to make them scalable. In
addition to the theoretical presentation of our algorithm, we provide DanceQ, a
flexible and modern - header only - C++20 implementation to manipulate,
enumerate, and map to its index any basis state in a given particle number
sector as open source software under https://DanceQ.gitlab.io/danceq.
