That the multiwavelet representation of this integral operator is sparse can be readily seen from the multipole expansion of 1/r. Consider the interaction between two wavelets of order k in two boxes separated by r. Since the first k moments vanish, the interaction will decay as r−(2k+1) . It is also necessary to consider the interaction between wavelets and polynomials, which will decay at worst as r−(k+1) . We commonly employ wavelets of order 5-13, and, by increasing the order of the wavelets as we increase the required precision, it is never
necessary to include interactions beyond the first 26 nearest neighbors of a box. This immediately results in an O(N ) algorithm with guaranteed precision. Also, the iterative solution of equations has been replaced with a single, fast, sparse matrix-vector product.