A commonly used prior distribution for the infinite-dimensional Bayesian formulation is the Whittle-Mat’ern prior, which involves a covariance operator based on a fractional elliptic operator. Due to computational considerations, the approach is computationally challenging for non-integer exponents. A similar computational bottleneck arises in generating samples using the SPDE approach to Gaussian random fields. We present two different solvers for efficiently applying the discretized covariance operator (and its square root) for all admissible values of the exponent: the first uses a multipreconditioned Krylov subspace method, and the second exploits a certain low-rank structure in the solution. We also discuss how the multipreconditioned solver can be used to accelerate the solution of linear Bayesian inverse problems to obtain the maximum a posteriori estimate and the approximate posterior variance. Numerical experiments demonstrate the performance and scalability of the solvers and their applicability to model and real-data inverse problems in tomography and a time-dependent heat equation. This is joint work with Harbir Antil (George Mason) and Hussam Al Daas (Rutherford Appleton Laboratory).