We have developed a symbolic algebra approach to automatically produce, verify, and optimize computer code for the Fast Multipole Method (FMM) operators. This approach allows for flexibility in choosing a basis set and kernel, and can generate computer code for any expansion order in multiple languages. The procedure is implemented in the publicly available Python program Mosaic. Optimizations performed at the symbolic level through algebraic manipulations significantly reduce the number of mathematical operations compared with a straightforward implementation of the equations. We find that the optimizer is able to eliminate 20%–80% of the floating-point operations and for the expansion orders p≤ 10 it changes the observed scaling properties. We present our approach using three variants of the operators with the Cartesian basis set for the harmonic potential kernel 1∕ r, including the use of totally symmetric and traceless multipole tensors.