When dealing with angular-momentum coupling one often faces sums that involve (at least) one angular momentum appearing in a 6j symbol. If these sums have to be calculated numerically it is thus desirable to be able to accurately calculate sets of 6j symbols where one of the angular momenta varies from its minimum to its maximum allowed value.
To this end, Schulten and Gordon [1] derived a linear three-term recurrence that links the values of neigboring 6j symbols. Starting from an arbitrary value for the first and the last 6j symbol they perform forward and backward recursion, matching both sequences in a central region. Imposing a normalization condition on the sequence obtained this way they recover the sequence of 6j symbols. The rather complicated procedure is necessary beacuse the values of the 6j symbols vary over many orders of magnitude in the classically-forbidden regions, causing numerical instabilities in a simple forward or backward recursion.
A simpler procedure was proposed by Luscombe and Luban [2]. Instead of calculating the symbols directly they target the ratios of neighboring 6j symbols. This collapses the three-term to a two-term recurrence — albeit a nonlinear one — and one can recover the values of the 6j symbols by successively multiplying the ratios and then normalizing the sequence. An issue with this procedure is that the ratios can become infinite if there is a 6j equal to zero in the sequence. However, these zeros only occur in the classically-allowed region where the three-term recurrence is numerically well-behaved.
This program combines both approaches: throughout the classically-forbidden regions the sequence is built up from both ends using the recurrence for the ratios. Once the sequences stop being monotonic (i.e. once the absolute value of the ratios stop being smaller or larger than one) the classically-allowed region is reached and the program connects both endpoints by a simple three-term forward recursion.
The 6j symbols calculated with this program are fairly accurate. Comparing results to Mathematica (which can evaluate 6j's to arbitrary precision) for the three test cases given in Figure 4 of Ref. [1] yields a maximum absolute and relative deviation of 8.674 x 10-19 and 1.274 x 10-13, respectively. The maximum errors from a random sampling of thousands of 6j symbols with j < 85/2 were of the same order of magnitude.
Compared to calculating the 6j symbols separately with GSL performance is improved six- to sevenfold, taking about 0.2μs per symbol on a single core of a Core i5-3470 CPU (3.20GHz). This figure is obtained from three sequences of 6j symbols, each consisting of 49 elements, and dividing the time spent in the calculation routines by the number of symbols calculated.
The code can be obtained here: sixj.c. The file consists of the calculation routine itself, a few utility functions and a main function that acts as a test program. See the comments for details.
The code is released under MIT license:
Copyright (c) 2014 Roland Wirth
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
If you have comments or suggestions, feel free to write me (remove the ham from the address).