  Starting from recent advances in the first-principles modeling of polarons,
variational polaron equations in the strong-coupling adiabatic approximation
are formulated in Bloch space. In this framework, polaron formation energy as
well as individual electron, phonon and electron-phonon contributions are
obtained. We suggest an efficient gradient-based optimization algorithm and
apply these equations to the generalized Fr\"ohlich model with anisotropic
non-degenerate electronic bands, both in two- and three-dimensional cases. The
effect of the divergence of Fr\"ohlich electron-phonon matrix elements at
$\Gamma$-point is treated analytically, improving the convergence with respect
to the sampling in reciprocal space. The whole methodology is validated by
obtaining the known asymptotic solution of the standard Fr\"ohlich model in
isotropic scenario and also by comparing our results with the Gaussian ansatz
approach, showing the difference between the numerically exact and Gaussian
trial wavefunctions. Additionally, decomposition of the energy into individual
terms allows one to recover the Pekar's 1:2:3:4 theorem, which is shown to be
valid even in the anisotropic case. We expect that the improvements in the
formalism and numerical implementation will be applicable beyond the large
polaron hypothesis inherent to Fr\"ohlich model.
