The virtual two-loop corrections for Higgs production in
gluon fusion are calculated analytically in QCD for arbitrary Higgs and
quark masses. Both scalar and pseudo-scalar Higgs bosons are
considered. The results are obtained by expanding the known
one-dimensional integral representation in terms of m_H/m_q, and
matching it with a suitably chosen ansatz of Harmonic
Polylogarithms. This ansatz is motivated by the known analytic result
for the Higgs decay rate into two photons. The method also allows us to
check this result and to extend it to the pseudo-scalar decay rate.