Change computation of core loop in Wigner D

By reorganizing the core loop calculating the sum for Wigner D matrices, we can
eliminate the need to index binomial coefficients twice per iteration.
Instead, we can factor out the initial values of those coefficients, and
convert the remaining factorials into multiplications at each iteration.  This
makes the overall loops a couple percent faster for complete D matrices, and
around 15% faster for SWSHs.  But this should make it significantly easier
to switch to gmpy2 to calculate the values at high precision for ell>29.

Note that this also required changing the overall coefficients to incorporate
those initial-value binomials.

The formulas behind these changes are explained in the gh-pages documentation:
<http://moble.github.io/spherical_functions/PolynomialsWithFactorials.html>.