  We present a new method to estimate three-point correlations in Cosmic
Microwave Background maps. Our Fast Fourier Transform based implementation
estimates three-point functions using all possible configurations (triangles)
at a controlled resolution. The speed of the technique depends both on the
resolution and the total number of pixels $N$. The resulting $N\log N$ scaling
is substantially faster than naive methods with prohibitive $N^3$ scaling. As
an initial application, we measure three-point correlation functions in the
First Year Data Release of the Wilkinson Anisotropy Probe. We estimate 336
cross-correlations of any triplet of maps from the 8 differential assemblies,
scanning altogether 2.6 million triangular configurations. Estimating
covariances from Gaussian signal plus realistic noise simulations, we perform a
null-hypothesis testing with regards to the Gaussianity of the Cosmic Microwave
Background. Our main result is that at the three-point level WMAP is fully
consistent with Gaussianity. To quantify the level of possible deviations, we
introduce false discovery rate analysis, a novel statistical technique to
analyze for three-point measurements. This confirms that the data are
consistent with Gaussianity at better than 1-$\sigma$ level when jointly
considering all configurations. We constrain a specific non-Gaussian model
using the quadratic approximation of weak non-Gaussianities in terms of the
$f_{NLT}$ parameter, for which we construct an estimator from the the
three-point function. We find that using the skewness alone is more
constraining than a heuristic suboptimal combination of all our results; our
best estimate is $f_{NLT} = -110\pm 150$ assuming a $\Lambda$CDM concordance
model.
