  We present a complete framework for numerical calculation of the power
spectrum and bispectrum in canonical inflation with an arbitrary number of
light or heavy fields. Our method includes all relevant effects at tree-level
in the loop expansion, including (i) interference between growing and decaying
modes near horizon exit; (ii) correlation and coupling between species near
horizon exit and on superhorizon scales; (iii) contributions from mass terms;
and (iv) all contributions from coupling to gravity. We track the evolution of
each correlation function from the vacuum state through horizon exit and the
superhorizon regime, with no need to match quantum and classical parts of the
calculation; when integrated, our approach corresponds exactly with the
tree-level Schwinger or 'in-in' formulation of quantum field theory. In this
paper we give the equations necessary to evolve all two- and three-point
correlation functions together with suitable initial conditions. The final
formalism is suitable to compute the amplitude, shape, and scale dependence of
the bispectrum in models with |fNL| of order unity or less, which are a target
for future galaxy surveys such as Euclid, DESI and LSST. As an illustration we
apply our framework to a number of examples, obtaining quantitatively accurate
predictions for their bispectra for the first time. Two accompanying reports
describe publicly-available software packages that implement the method.
