  We determine the infinite volume coefficients of the perturbative expansions
of the self-energies of static sources in the fundamental and adjoint
representations in SU(3) gluodynamics to order \alpha^{20} in the strong
coupling parameter \alpha. We use numerical stochastic perturbation theory,
where we employ a new second order integrator and twisted boundary conditions.
The expansions are obtained in lattice regularization with the Wilson action
and two different discretizations of the covariant time derivative within the
Polyakov loop. Overall, we obtain four different perturbative series. For all
of them the high order coefficients display the factorial growth predicted by
the conjectured renormalon picture, based on the operator product expansion.
This enables us to determine the normalization constants of the leading
infrared renormalons of heavy quark and heavy gluino pole masses and to
translate these into the modified minimal subtraction scheme (MS). We also
estimate the four-loop \beta-function coefficient of the lattice scheme.
