Numerical evaluation of two-dimensional harmonic polylogarithms

The two-dimensional harmonic polylogarithms $\G(\vec{a}(z);y)$, a generalization of the harmonic polylogarithms, themselves a generalization of Nielsen's polylogarithms, appear in analytic calculations of multi-loop radiative corrections in quantum field theory. We present an algorithm for the numerical evaluation of two-dimensional harmonic polylogarithms, with the two arguments $y,z$ varying in the triangle $0\le y \le 1$, $ 0\le z \le 1$, $\ 0\le (y+z) \le 1$. This algorithm is implemented into a {\tt FORTRAN} subroutine {\tt tdhpl} to compute two-dimensional harmonic polylogarithms up to weight 4.