  We present an efficient separable approach to the estimation and
reconstruction of the bispectrum and the trispectrum from observational (or
simulated) large scale structure data. This is developed from general CMB
(poly-)spectra methods which exploit the fact that the bispectrum and
trispectrum in the literature can be represented by a separable mode expansion
which converges rapidly (with $n_\textrm{max}={\cal{O}}(30)$ terms). With an
effective grid resolution $l_\textrm{max}$ (number of particles/grid points
$N=l_\textrm{max}^3$), we present a bispectrum estimator which requires only
${\cal O}(n_\textrm{max} \times l_\textrm{max}^3)$ operations, along with a
corresponding method for direct bispectrum reconstruction. This method is
extended to the trispectrum revealing an estimator which requires only ${\cal
O}(n_\textrm{max}^{4/3} \times l_\textrm{max}^3)$ operations. The complexity in
calculating the trispectrum in this method is now involved in the original
decomposition and orthogonalisation process which need only be performed once
for each model. However, for non-diagonal trispectra these processes present
little extra difficulty and may be performed in ${\cal O}(l_\textrm{max}^4)$
operations. A discussion of how the methodology may be applied to the
quadspectrum is also given. An efficient algorithm for the generation of
arbitrary nonGaussian initial conditions for use in N-body codes using this
separable approach is described. This prescription allows for the production of
nonGaussian initial conditions for arbitrary bispectra and trispectra. A brief
outline of the key issues involved in parameter estimation, particularly in the
non-linear regime, is also given.
