Fourier transforms are an often necessary component in many computational tasks, and can be computed efficiently through the fast Fourier transform (FFT) algorithm. However, many applications involve an underlying continuous signal, and a more natural choice would be to work with e.g. the Fourier series (FS) coefficients in order to avoid the additional overhead of translating between the analog and discrete domains. Unfortunately, there exists very little literature and tools for the manipulation of FS coefficients from discrete samples. This paper introduces a Python library called pyFFS for efficient FS coefficient computation, convolution, and interpolation. While the libraries SciPy and NumPy provide efficient functionality for discrete Fourier transform coefficients via the FFT algorithm, pyFFS addresses the computation of FS coefficients through what we call the fast Fourier series (FFS). Moreover, pyFFS includes an FS interpolation method based on the chirp Z-transform that can make it more than an order of magnitude faster than the SciPy equivalent when one wishes to perform interpolation. GPU support through the CuPy library allows for further acceleration, e.g. an order of magnitude faster for computing the 2-D FS coefficients of 1000 x 1000 samples and nearly two orders of magnitude faster for 2-D interpolation. As an application, we discuss the use of pyFFS in Fourier optics. pyFFS is available as an open source package at https://github.com/imagingofthings/pyFFS, with documentation at https://pyffs.readthedocs.io.