We present a new formalism, together with efficient numerical methods, to
directly calculate the CMB bispectrum today from a given primordial bispectrum
using the full linear radiation transfer functions. Unlike previous analyses
which have assumed simple separable ansatze for the bispectrum, this work
applies to a primordial bispectrum of almost arbitrary functional form, for
which there may have been both horizon-crossing and superhorizon contributions.
We employ adaptive methods on a hierarchical triangular grid and we establish
their accuracy by direct comparison with an exact analytic solution, valid on
large angular scales. We demonstrate that we can calculate the full CMB
bispectrum to greater than 1% precision out to multipoles l<1800 on reasonable
computational timescales. We plot the bispectrum for both the superhorizon
('local') and horizon-crossing ('equilateral') asymptotic limits, illustrating
its oscillatory nature which is analogous to the CMB power spectrum