  The construction of constraint-satisfying initial data is an essential
element for the numerical exploration of the dynamics of compact-object
binaries. While several codes have been developed over the years to compute
generic quasi-equilibrium configurations of binaries comprising either two
black holes, or two neutron stars, or a black hole and a neutron star, these
codes are often not publicly available or they provide only a limited
capability in terms of mass ratios and spins of the components in the binary.
We here present a new open-source collection of spectral elliptic solvers that
are capable of exploring the major parameter space of binary black holes
(BBHs), binary neutron stars (BNSs), and mixed binaries of black holes and
neutron stars (BHNSs). Particularly important is the ability of the
spectral-solver library to handle neutron stars that are either irrotational or
with an intrinsic spin angular momentum that is parallel to the orbital one. By
supporting both analytic and tabulated equations of state at zero or finite
temperature, the new infrastructure is particularly geared towards allowing for
the construction of BHNS and BNS binaries. For the latter, we show that the new
solvers are able to reach the most extreme corners in the physically plausible
space of parameters, including extreme mass ratios and spin asymmetries, thus
representing the most extreme BNS computed to date. Through a systematic series
of examples, we demonstrate that the solvers are able to construct
quasi-equilibrium and eccentricity-reduced initial data for BBHs, BNSs, and
BHNSs, achieving spectral convergence in all cases. Furthermore, using such
initial data, we have carried out evolutions of these systems from the inspiral
to after the merger, obtaining evolutions with eccentricities $\lesssim
10^{-4}-10^{-3}$, and accurate gravitational waveforms.
