  Three--loop contributions to massive QED vacuum polarization are evaluated by
a combination of analytical and numerical techniques. The first three Taylor
coefficients, at small $q^2$, are obtained analytically, using
$d$\/--dimensional recurrence relations. Combining these with analytical input
at threshold, and at large $q^2$, an accurate Pad\'e approximation is obtained,
for all $q^2$. Inserting this in the one--loop diagram for the muon anomalous
magnetic moment, we find reasonable agreement with four--loop,
single--electron--loop, muon--anomaly contributions, recently re--evaluated by
Kinoshita, using 8--dimensional Monte--Carlo integration. We believe that our
new method is at least two orders of magnitude more accurate than the
Monte--Carlo approach, whose uncertainties appear to have been underestimated,
by a factor of 6.
