  We consider Quantum Chromodynamics with external vector, axial-vector, scalar
and pseudo-scalar currents and compute three-loop corrections to the
corresponding vertex function taking into account massive quarks. We consider
all non-singlet contributions as well as those singlet contributions where the
external current couples to a massive quark loop. We apply a semi-numerical
method which is based on expansions around singular and regular kinematical
points. They are matched at intermediate values of the squared partonic
center-of-mass energy $s$ which allows to cover the whole kinematic range for
negative and positive values of $s$. Our method permits a systematic increase
of the precision by varying the expansion depth and the choice of the
intermediate matching points. In our current set-up we have at least seven
significant digits for the finite contribution of all form factors. We present
our results as a combination of series expansions and interpolation functions
which allows for a straightforward use in practical applications.
