  As demonstrated by Slepian et. al. in a sequence of classical papers, prolate
spheroidal wave functions (PSWFs) provide a natural and efficient tool for
computing with bandlimited functions defined on an interval. As a result, PSWFs
are becoming increasing popular in various areas in which such function occur -
this includes physics (e.g. wave phenomena, fluid dynamics), engineering (e.g.
signal processing, filter design), etc.
  To use PSWFs as a computational tool, one needs fast and accurate numerical
algorithms for the evaluation of PSWFs and related quantities, as well as for
the construction of quadratures, interpolation formulas, etc. Even though, for
the last half a century, substantial progress has been made in design of such
algorithms, the complexity of many of the existing algorithms, however, is at
least quadratic in the band limit $c$. For example, the evaluation of the $n$th
eigenvalue of the prolate integral operator requires at least $O(c^2)$
operations. Therefore, while the existing algorithms are quite satisfactory for
moderate values of $c$ (e.g. $c \leq 10^3$), they tend to be relatively slow
when $c$ is large (e.g. $c \geq 10^4$).
  In this paper, we describe several numerical algorithms for the evaluation of
PSWFs and related quantities, and design a class of PSWF-based quadratures for
the integration of bandlimited functions. Also, we perform detailed analysis of
the related properties of PSWFs. While the analysis is somewhat involved, the
resulting numerical algorithms are quite simple and efficient in practice. For
example, the evaluation of the $n$th eigenvalue of the prolate integral
operator requires $O(n+c)$ operations; also, the construction of related
accurate quadrature rules requires $O(c)$ operations.
  Our results are illustrated via several numerical experiments.
