  A scheme for rapidly and accurately computing solutions to boundary integral
equations (BIEs) on rotationally symmetric surfaces in three dimensions is
presented. The scheme uses the Fourier transform to reduce the original BIE
defined on a surface to a sequence of BIEs defined on a generating curve for
the surface. It can handle loads that are not necessarily rotationally
symmetric. Nystrom discretization is used to discretize the BIEs on the
generating curve. The quadrature used is a high-order Gaussian rule that is
modified near the diagonal to retain high-order accuracy for singular kernels.
The reduction in dimensionality, along with the use of high-order accurate
quadratures, leads to small linear systems that can be inverted directly via,
e.g., Gaussian elimination. This makes the scheme particularly fast in
environments involving multiple right hand sides. It is demonstrated that for
BIEs associated with Laplace's equation, the kernel in the reduced equations
can be evaluated very rapidly by exploiting recursion relations for Legendre
functions. Numerical examples illustrate the performance of the scheme; in
particular, it is demonstrated that for a BIE associated with Laplace's
equation on a surface discretized using 320 000 points, the set-up phase of the
algorithm takes 2 minutes on a standard desktop, and then solves can be
executed in 0.5 seconds.
