  We present a new algorithm to rapidly compute the two-point (2PCF),
three-point (3PCF) and n-point (n-PCF) correlation functions in roughly O(N log
N) time for N particles, instead of O(N^n) as required by brute force
approaches. The algorithm enables an estimate of the full 3PCF for as many as
10^6 galaxies. This technique exploits node-to-node correlations of a recursive
bisectional binary tree. A balanced tree construction minimizes the depth of
the tree and the worst case error at each node. The algorithm presented in this
paper can be applied to problems with arbitrary geometry.
  We describe the detailed implementation to compute the two point function and
all eight components of the 3PCF for a two-component field, with attention to
shear fields generated by gravitational lensing. We also generalize the
algorithm to compute the n-point correlation function for a scalar field in k
dimensions where n and k are arbitrary positive integers.
