We develop a new type of orthogonal polynomial, the modified discrete Laguerre (MDL) polynomials, designed to accelerate the computation of bosonic Matsubara sums in statistical physics. The MDL polynomials lead to a rapidly convergent Gaussian quadrature scheme for Matsubara sums, and more generally for any sum $F(0)/2 + F(h) + F(2h) + cdots$ of exponentially decaying summands $F(nh) = f(nh)e^{-nhs}$ where $hs>0$. We demonstrate this technique for computation of finite-temperature Casimir forces arising from quantum field theory, where evaluation of the summand $F$ requires expensive electromagnetic simulations. A key advantage of our scheme, compared to previous methods, is that the convergence rate is nearly independent of the spacing $h$ (proportional to the thermodynamic temperature). We also prove convergence for any polynomially decaying $F$.