  A precise calculation of the lepton anomalous magnetic moments (AMM) requires
an evaluation of the quantum electrodynamics (QED) Feynman diagrams up to five
independent loops. The complicated structure of ultraviolet (UV), infrared (IR)
and mixed divergences in the corresponding integrals makes it difficult to
calculate these high-order contributions in reasonable computer time frame.
  We demonstrate a method that eliminates all divergences point by point in
Feynman parametric space (before integration) and possesses a flexibility that
can be used for improving the precision of the numerical integration. This
flexibility is especially actual for calculating the contributions of the
Feynman diagrams with electron loops to the muon AMM. 3-loop and 4-loop
numerical test results are provided.
  The subtraction procedure is based on a forest formula with linear operators
applied to the Feynman amplitudes of UV divergent subdiagrams. It is similar to
BPHZ; the difference is in the linear operators used and in the way of
combining them. It is equivalent to the on-shell renormalization after
summation over diagrams: no residual renormalization is required.
