  We calculate the non-singlet, the pure singlet contribution, and their
interference term, at $O(\alpha^2)$ due to electron-pair initial state
radiation to $e^+ e^-$ annihilation into a neutral vector boson in a direct
analytic computation without any approximation. The correction is represented
in terms of iterated incomplete elliptic integrals. Performing the limit $s \gg
m_e^2$ we find discrepancies with the earlier results of
Ref.~\cite{Berends:1987ab} and confirm results obtained in
Ref.~\cite{Blumlein:2011mi} where the effective method of massive operator
matrix elements has been used, which works for all but the power corrections in
$m^2/s$. In this way, we also confirm the validity of the factorization of
massive partons in the Drell-Yan process. We also add non-logarithmic terms at
$O(\alpha^2)$ which have not been considered in \cite{Berends:1987ab}. The
corrections are of central importance for precision analyzes in $e^+e^-$
annihilation into $\gamma^*/Z^*$ at high luminosity.
