We study black hole linear perturbation theory in a four-dimensional Schwarzschild (anti) de Sitter background. When dealing with a positive cosmological constant, the corresponding spectral problem is solved systematically via the Nekrasov-Shatashvili functions or, equivalently, classical Virasoro conformal blocks. However, this approach can be more complicated to implement for certain perturbations if the cosmological constant is negative. For these cases, we propose an alternative method to set up perturbation theory for both small and large black holes in an analytical manner. Our analysis reveals a new underlying recursive structure that involves multiple polylogarithms. We focus on gravitational, electromagnetic, and conformally coupled scalar perturbations subject to Dirichlet and Robin boundary conditions. The low-lying modes of the scalar sector of gravitational perturbations and its hydrodynamic limit are studied in detail.