Reconfigure calculation of Wigner D matrices

Based on my re-derivation of the D matrices, I think I understand
better how to deal with roundoff, overflow, and underflow.

Also, I realized that it's easy to pull some of those if statements
out of the loops, which should help improve the speed slightly.

Finally, I've tightened the tolerance, so that 1.e-15 is considered
epsilon, instead of 1.e-14.  This should improve the accuracy of the
results in certain cases with large ell values.
Comments
Add first nontrivial test of WignerD