  We develop and demonstrate an acceleration of the Liu & Tegmark quadratic
estimator formalism for inverse variance foreground subtraction and power
spectrum estimation in 21 cm tomography from O(N^3) to O(N log N), where N is
the number of voxels of data. This technique makes feasible the megavoxel scale
analysis necessary for current and upcoming radio interferometers by making
only moderately restrictive assumptions about foreground models and survey
geometry. We exploit iterative and Monte Carlo techniques and the symmetries of
the foreground covariance matrices to quickly estimate the 21 cm brightness
temperature power spectrum, P(k_parallel, k_perpendicular), the Fisher
information matrix, the error bars, the window functions, and the bias. We also
extend the Liu & Tegmark foreground model to include bright point sources with
known positions in a way that scales as O[(N log N)(N point sources)] <
O(N^5/3). As a first application of our method, we forecast error bars and
window functions for the upcoming 128-tile deployment of the Murchinson
Widefield Array, showing that 1000 hours of observation should prove
sufficiently sensitive to detect the power spectrum signal from the Epoch of
Reionization.
