  Within the nonrelativistic QCD (NRQCD) factorization framework, we compute
the order-$\alpha_s^3$ perturbative corrections to the leptonic decay of the
$\Upsilon$ and $J/\psi$ with high numerical accuracy, at the lowest order in
velocity expansion. We confirm the existing three-loop results in literature.
Furthermore, we explicitly consider the complex-valued light-by-light (singlet)
contributions. We for the first time also consider the finite charm quark mass
effect for $\Upsilon$ leptonic decay, and obtain a new piece of 3-loop
contribution to the anomalous dimension related to the composite NRQCD bilinear
of vector current that arises from the nonzero charm quark mass. Based on the
complete three-loop NRQCD short-distance coefficients, we also present a
comprehensive phenomenological analysis for $\Upsilon(J/\psi)$ leptonic width.
