  We describe an algorithm to generate temperature and polarization maps of the
cosmic microwave background radiation containing non-Gaussianity of arbitrary
local type. We apply an optimized quadrature scheme that allows us to predict
and control integration accuracy, speed up the calculations, and reduce memory
consumption by an order of magnitude. We generate 1000 non-Gaussian CMB
temperature and polarization maps up to a multipole moment of l_max = 1024. We
validate the method and code using the power spectrum and the fast cubic
(bispectrum) estimator and find consistent results. The simulations are
provided to the community.
