The implementation of a method to compute two-loop vertex diagrams with two massless on-shell legs in analytical form is discussed. The strategy is to map the integrals to massless three-loop two-point functions which can be computed with the help of the integration-by-parts algorithm. The known results for the two-loop quark form factor in QCD as well as the three-gluon vertex with two on-shell gluons could be reproduced. As an application we present the virtual corrections to Higgs production in gluon-gluon fusion up to order \alpha_s^4.