Code for review: New wigner coupling coefficient code and rotation matrix code

[Jonathan Underwood <j dot underwood at open dot ac dot uk>]

Dear All,

I've finally put some code together to:

(a) improve the code for calculating wigner angular momentum coupling 
coefficients (3j, 6j and 9j), replacing the code in specfunc/coupling.c.

(b) added code for the calculation of the Wigner reduced rotation 
matrices.

This code can be downloaded from <http://physics.open.ac.uk/~ju83>. 
wigner.c contains the code, gsl_sf_wigner.h is the header file, and there 
is a modified test_sf.c including tests for this code. Should compile 
cleanly when added to gsl-1.5 released code, with additions to Makefile.am 
as appropriate. I haven't tested with CVS, but can't see why there would 
be a problem.

Justification:
(b) is easily justifiable - it's new :) The algorithm for the calculation 
of the reduced wigner rotation is very similar to that for the calculation 
of the 3j, 6j, and 9j symbols, hence it is bundled together.

(a) Why replce coupling.c ? coupling.c has a checkered history. The 
original algorithms involved multiplications of factorials. However it was 
noted that this became innaccurate for large angular momenta (see BUGS for 
gsl-1.3 and previous bug reports). This led to the contribution of new 
code for the calculation of 3j symbols in terms of binomial coefficients. 
However, the 6j routine will suffer from the same problem as the same 
strategy is used in coupling.c for their calculation; that hasn't yet been 
remedied. The 9j calculation is done as a sum over 6j's, and so the 9j 
code will also be affected. Calculation in terms of binomials, while 
correct, is a little more complicated and less efficient. What I've done 
with the wigner.c code is reverted to the original expressions for the nj 
symbols in terms of factorials (as found in eg. Angular Momentum by Zare) 
but carried out their calculation in terms of summing logarithms of 
factorials rather than multiplying factorials, thus avoiding overflow 
errors etc. I've used equivalent code for many years (my own code), and 
this is also the strategy outlined in the fortran code in Angular Momentum 
by Zare, and also Thompsons out-of-print book. This code is I believe 
accurate and efficient at high J (certainly the 3j calculation does not 
suffer the previously reported bug). Also in the process I noticed a bug 
in coupling.c in that no testing is done for (m,j) being both integer, or 
both half-integer - this should be undefined, and not zero. Furthermore, 
the coupling code returns 0.0 for the cases when |m| > j, whereas it 
should really return NaN - this is fixed in wigner.c. See the comments 
I've added to the test_sf.c code for both coupling and wigner for more info.

However, I'm not entirely sure wigner.c correctly analyses the errors 
currently, and this is where a second expert eye is needed to look at the 
code. I tried to work it out from coupling.c, but could make no sense of 
that in terms of error calculation (I'm not sure it rigoriously uses the 
prescription for errors for GSL as Brian previously listed on this list). 
I have propagated errors as best I could, but there are a couple of 
ambiguities:
(a) Quadrature - Brian's post previously said to ignore quadrature for 
addition and subtraction, so I've also done this for multplication - 
correct/wrong?

(b) There isn't a version of gsl_sf_pow_int_e that takes into account the 
errors in the input value, so I've used the fact that for a^b, the error 
in the result is b times the error in a - correct? I've bolted this on for 
now, but it might be an idea to add error propagation functions for 
integer powers too?

Haven't yet done any documentation, but happy to do that if this code is 
included in gsl (which i think it should be). All comments, critiques etc 
welcome.

Best wishes
Jonathan