  We present an algorithm to compute correlation functions for systems with the
quantum numbers of many identical mesons from lattice quantum chromodynamics
(QCD). The algorithm is numerically stable and allows for the computation of
$n$-pion correlation functions for $n \in \{ 1, \dots, N\}$ using a single $N
\times N$ matrix decomposition, improving on previous algorithms. We apply the
algorithm to calculations of correlation functions with up to 6144 $\pi^+$s
using two ensembles of gauge field configurations generated with quark masses
corresponding to a pion mass $m_\pi = 170$ MeV and spacetime volumes of
$(4.4^3\times 8.8)\ {\rm fm}^4$ and $(5.8^3\times 11.6)\ {\rm fm}^4$. We also
discuss statistical techniques for the analysis of such systems, in which the
correlation functions vary over many orders of magnitude. In particular, we
observe that the many-pion correlation functions are well approximated by
log-normal distributions, allowing the extraction of the energies of these
systems. Using these energies, the large-isospin-density, zero-baryon-density
region of the QCD phase diagram is explored. A peak is observed in the energy
density at an isospin chemical potential $\mu_I\sim 1.5 m_\pi$, signalling the
transition into a Bose-Einstein condensed phase. The isentropic speed of sound
in the medium is seen to exceed the ideal-gas (conformal) limit ($c_s^2\leq
1/3$) over a wide range of chemical potential before falling towards the
asymptotic expectation at $\mu_I\sim 15 m_\pi$. These, and other thermodynamic
observables, indicate that the isospin chemical potential must be large for the
system to be well described by an ideal gas or perturbative QCD.
