Precise statistical analyses of astronomical data are the key to validating models of complex phenomena, such as dark matter and dark energy. In particular, spatial statistics are needed for large-scale sky catalogs. The n-point correlation functions provide a complete description of any point process and are widely used to understand astronomical data. However, the computational cost of estimating these functions scales as N^n for N data points. Furthermore, these expensive computations must be repeated many times at many different scales in order to gain a detailed picture of the correlation function and to estimate its variance. Since astronomy surveys contain hundreds of millions or billions of points (and are growing rapidly), these computations are infeasible. We present a new approach based on multidimensional trees to overcome these computational obstacles. We build on the previously most efficient algorithm (Gray and Moore, 2001, Moore, et al., 2001) which improved over the N^n scaling of a direct computation. In this work, we incorporate the computations at different scales along with the variance estimation directly. We can therefore achieve an order of magnitude speedup over the current state-of-the-art method. We show preliminary scaling results on a mock galaxy catalog.