  We present a new flexible, fast and accurate way to implement massive
neutrinos, warm dark matter and any other non-cold dark matter relics in
Boltzmann codes. For whatever analytical or numerical form of the phase-space
distribution function, the optimal sampling in momentum space compatible with a
given level of accuracy is automatically found by comparing quadrature methods.
The perturbation integration is made even faster by switching to an approximate
viscous fluid description inside the Hubble radius, which differs from previous
approximations discussed in the literature. When adding one massive neutrino to
the minimal cosmological model, CLASS becomes just 1.5 times slower, instead of
about 5 times in other codes (for fixed accuracy requirements). We illustrate
the flexibility of our approach by considering a few examples of standard or
non-standard neutrinos, as well as warm dark matter models.
