The computation of system properties which need information from all pairs of system elements (e.g. particles, mesh points) are examples for problems with quadratic computational complexity, i.e. they need O(N²) operations to compute the desired solution. In many scenarios it is not only the number of operations, but also the size in memory and communication costs, which increase in a non-linear way. Both characteristics are a strong limitation for the performance and scalability of the underlying method and therefore, numerical methods have been developed, which reduce the complexity to O(N log(N)) or O(N).  

The Fast Multipole Method, has been originally developed in the context of electrostatic interactions of complex particle systems. It combines optimal complexity O(N) and adjustable accuracy in one fast summation technique for a wide range of applications in modern numerical science. Meanwhile it has been identified as one of the Top10 algorithms of the last century [0] and a demonstrator for future Exascale systems [1] . It bridges theory and simulation for more than 30 years and has has entered many application domains due to its universal properties of reducing complexity. Applications range from energy computations in N-body systems, governed by Coulomb[2] or gravitational [3] interactions, electromagnetic scattering [4], boundary element methods [5] for Stokes, acoustic wave [6] and elastostatic problems [7] to preconditioners for dense linear algebra [8]. As a mesh-free method the FMM shows significant advantages over grid-based methods, especially for applications with spatial non-homogeneously distributed input. Flexibility of FMMs has been increased by the advent of kernel-independent implementations [9] and a priori error control schemes [10].