moble/spherical_functions — pull request #15 by moble (2014-12-04 05:43:43+00:00)
Title: New sum for better stability

As mentioned in issue #12, the calculations suddenly get very bad for ell>29. I also noticed that the errors are growing exponentially from low ell, for certain rotations. I've replaced the key (alternating) sum with what should be a much stabler algorithm. The basic idea is that you can rewrite any rotation as the product of two rotations for which that alternating sum has a special form. In this case, the magnitudes of the R.a and R.b components are equal for each of the two rotations. We can evaluate the sums exactly offline, and store them for retrieval when needed (as the new `Delta` function).

I can still use Horner form for fast and accurate evaluations. In fact, I even managed to fold that loop onto itself, saving ~1/2 the iterations, by using a symmetry of the Delta.

These changes make the evaluation substantially slower (~20% slower for ell=8), but far more accurate. In particular, I've been able to tighten the tolerances on the tests, so that errors of roughly 2*ell times machine precision are expected.