  Based on the Beylkin-Cramer summation rule, we introduce a new fast algorithm
that enable us to explore the high order statistics efficiently in large data
sets. Central to this technique is to make decomposition both of fields and
operators within the framework of multi-resolution analysis (MRA), and realize
theirs discrete representations. Accordingly, a homogenous point process could
be equivalently described by a operation of a Toeplitz matrix on a vector,
which is accomplished by making use of fast Fourier transformation. The
algorithm could be applied widely in the cosmic statistics to tackle large data
sets. Especially, we demonstrate this novel technique using the spherical,
cubic and cylinder counts in cells respectively. The numerical test shows that
the algorithm produces an excellent agreement with the expected results.
Moreover, the algorithm introduces naturally a sharp-filter, which is capable
of suppressing shot noise in weak signals. In the numerical procedures, the
algorithm is somewhat similar to particle-mesh (PM) methods in N-body
simulations. As scaled with $O(N\log N)$, it is significantly faster than the
current particle-based methods, and its computational cost does not relies on
shape or size of sampling cells. In addition, based on this technique, we
propose further a simple fast scheme to compute the second statistics for
cosmic density fields and justify it using simulation samples. Hopefully, the
technique developed here allows us to make a comprehensive study of
non-Guassianity of the cosmic fields in high precision cosmology. A specific
implementation of the algorithm is publicly available upon request to the
author.
