The standard formula for squared Mahalanobis distance between two data points is
$$ D_{ij} = (x_1-x_2)^T \Sigma^{-1} (x_1-x_2) $$$$ D_{12} = (x_1-x_2)^T \Sigma^{-1} (x_1-x_2) $$
where $x_i$ is a $p \times 1$ vector corresponding to observation $i$. Typically, the covariance matrix is estimated from the observed data. Not counting matrix inversion, this operation requires $p^2+p$ multiplications and $p^2+2p$ additions, each repeated $n(n-1)/2$ times.
Consider the following derivation:
\begin{eqnarray*} D_{ij} &=& (x_1-x_2)^T \Sigma^{-1} (x_1-x_2) \\ &=& (x_1-x_2)^T \Sigma^{-\frac{1}{2}} \Sigma^{-\frac{1}{2}} (x_1-x_2) \\ &=& (x_1^T \Sigma^{-\frac{1}{2}} - x_2^T \Sigma^{-\frac{1}{2}}) (\Sigma^{-\frac{1}{2}}x_1 - \Sigma^{-\frac{1}{2}}x_2) \\ &=& (q_1^T - q_2^T)(q_1 - q_2) \end{eqnarray*}\begin{eqnarray*} D_{12} &=& (x_1-x_2)^T \Sigma^{-1} (x_1-x_2) \\ &=& (x_1-x_2)^T \Sigma^{-\frac{1}{2}} \Sigma^{-\frac{1}{2}} (x_1-x_2) \\ &=& (x_1^T \Sigma^{-\frac{1}{2}} - x_2^T \Sigma^{-\frac{1}{2}}) (\Sigma^{-\frac{1}{2}}x_1 - \Sigma^{-\frac{1}{2}}x_2) \\ &=& (q_1^T - q_2^T)(q_1 - q_2) \end{eqnarray*}
where $q_i = \Sigma^{-\frac{1}{2}}x_i$. Note that $x_i^T \Sigma^{-\frac{1}{2}} = (\Sigma^{-\frac{1}{2}} x_i)^T = q_i^T$. This relies on the fact that $\Sigma^{-\frac{1}{2}}$ is symmetric, which holds due to the fact that for any symmetric diagonalizable matrix $A = PEP^T$,
\begin{eqnarray*} A^{\frac{1}{2}^T} &=& (PE^{\frac{1}{2}}P^T)^T \\ &=& P^{T^T} E^{\frac{1}{2}^T} P^T \\ &=& PE^{\frac{1}{2}}P^T \\ &=& A^{\frac{1}{2}} \end{eqnarray*}
If we let $A=\Sigma^{-1}$, and note that $\Sigma^{-1}$ is symmetric, we see that that $\Sigma^{-\frac{1}{2}}$ must also be symmetric. If $X$ is the $n \times p$ matrix of observations and $Q$ is the $n \times p$ matrix such that the $i^{th}$ row of $Q$ is $q_i$, then $Q$ can be succinctly expressed as $X\Sigma^{-\frac{1}{2}}$. This and the previous results imply that
$$D_{k\ell} = \sum_{i=1}^p (Q_{ki}-Q_{\ell i})^2.$$ the only operations that are computed $n(n-1)/2$ times are $p$ multiplications and $2p$ additions (as opposed to the $p^2+p$ multiplications and $p^2+2p$ additions in the above method), resulting in an algorithm that is of computational complexity order $O(pn^2 + p^2n)$ instead of the original $O(p^2n^2)$.
require(ICSNP) # for pair.diff(), C implementation
fastPwMahal = function(data) {
# Calculate inverse square root matrix
invCov = solve(cov(data))
svds = svd(invCov)
invCovSqr = svds$u %*% diag(sqrt(svds$d)) %*% t(svds$u)
Q = data %*% invCovSqr
# Calculate distances
# pair.diff() calculates the n(n-1)/2 element-by-element
# pairwise differences between each row of the input matrix
sqrDiffs = pair.diff(Q)^2
distVec = rowSums(sqrDiffs)
# Create dist object without creating a n x n matrix
attr(distVec, "Size") = nrow(data)
attr(distVec, "Diag") = F
attr(distVec, "Upper") = F
class(distVec) = "dist"
return(distVec)
}