Consider the Poisson problem on a d-dimensional cube. It is well-known that, if the problem is discretized with linear finite elements on a uniform tensor product mesh, the resulting stiffness matrix can be diagonalized using the Fast Fourier Transform. This fact can be exploited to solve the linear system yielding O(N log N) complexity, where N represents the number of degrees of freedom. Such approach is referred to as a fast Poisson solver. In this talk, we show how to generalize this idea to the case of B-splines of arbitrary degree p. The resulting algorithm solves the linear system with O((N+p) \log N) complexity. This is achieved by first splitting the spline space into an outlier-free subspace and a subspace with low dimension. On the latter subspace, the eigenvectors of the problem are computed numerically. On the former subspace, on the other hand, the eigenvectors are approximated using interpolated sinusoidal functions. The resulting approximated eigendecomposition can be used as a preconditioner for the linear system, yielding extremely fast convergence independently of N and p.