Wigner 3j returns 0s instead of actual input, but yields good inputs on permutation of the l's

[joeydumont]

Thanks for Jeff Zheng for pointing this out to me via email (@jeffzhen). This comes mostly verbatim for his e-mail to me.

I found some quite confusing bug where it returns 0 (where Mathematica returns non-zero). In all cases, the whole list of values for all l1 are 0 or -0. However, if I permute the 1,2,3 indices to say 2,3,1, I always get the correct answer. Here's one example:
[l1, l2, l3, m1, m2, m3] = [ 325 999 1221 280 899 -1179] Mathematica answer: 0.000163542255419
call (l1, l2, l3, m1, m2, m3): 0
call (l2, l3, m1, m2, m3): all 0's returned
call (l2, l3, l1, m2, m3, m1): 0.000163542255418

I appended more examples at the end. Would you mind take a quick look and see what is causing this? Since your code can clearly handle l numbers much larger than these, and it can return correct answer by just permuting the indices, my guess is that it is some benign bug?

More examples:
[ 693 896 1371 513 -838 325] 3.82641344761e-140 -0.0 0.0175918809377 3.82641344761e-140
[ 772 874 1231 442 756 -1198] 0.00171146251356 0.0 0.0 0.00171146251356
[ 842 996 1714 376 979 -1355] -9.12730017838e-14 -0.0 0.0 -9.12730017838e-14
[ 101 987 1084 -80 -935 1015] -0.00274486614091 -0.0 0.0 -0.00274486614091
[ 51 1003 978 -32 993 -961] -0.00561491486489 -0.0 0.0 -0.00561491486489
[ 217 1008 1107 -76 -987 1063] -0.00129163704248 -0.0 0.0 -0.00129163704248
[ 408 894 954 -114 -840 954] 0.000141186572123 0.0 0.0 0.000141186572123
[ 808 980 1734 341 954 -1295] -1.13581014918e-35 -0.0 0.0 -1.13581014918e-35
[ 827 1020 1570 590 902 -1492] -0.00091099434249 -0.0 0.0 -0.00091099434249

[joeydumont]

The issue seems to be related to the algorithm itself, not this particular implementation. For such high quantum numbers, the classical region (where the solution to the recurrence relations are oscillatory) are quite large. There doesn't seem to be an easy around this. For the time being, I have implemented a temporary hack that simply determines which quantum number is the largest and uses that one as l1. Note that will not solve the problem when using std::vector<double> wigner3j(...), as it is precisely this sequence that is unstable. The other indices are permuted accordingly (this uses the permutation symmetries of the Wigner 3j symbols).

There exists a similar algorithm in the literature that claims to remove the overflow problems associated with the Schulten+Gordon algorithm. I will try implementing it and see if this issue still appears.

• J. Luscombe and M. Luban, "Simplified recursive algorithm for Wigner 3j and 6j symbols," Phys. Rev. E 57, 7274–7277 (1998). DOI: 10.1103/PhysRevE.57.7274.

[jeffzhen]

Hi Joey,
Thanks a lot for spending time looking into this! If I'm understanding you
correctly, your current temporary fix only applies to the 3j function that
returns 1 number, but not to the function that return a whole vector for
all possible l1? We will be eagerly looking forward to your new algorithm,
since we actually use the vector version of the 3j function, rather than
the single value version.

Thank you so much for your great work again!
Jeff

On Wed, Jul 16, 2014 at 2:56 PM, Joey Dumont [email protected]
wrote:

The issue seems to be related to the algorithm itself, not this particular
implementation. For such high quantum numbers, the classical region (where
the solution to the recurrence relations are oscillatory) are quite large.
There doesn't seem to be an easy around this. For the time being, I have
implemented a temporary hack that simply determines which quantum number is
the largest and uses that one as l1. The other indices are permuted
accordingly (this uses the permutation symmetries of the Wigner 3j
symbols).

There exists a similar algorithm in the literature that claims to remove
the overflow problems associated with the Schulten+Gordon algorithm. I will
try implementing it and see if this issue still appears.

• J. Luscombe and M. Luban, "Simplified recursive algorithm for Wigner
  3j and 6j symbols," Phys. Rev. E 57, 7274–7277 (1998). DOI:
  10.1103/PhysRevE.57.7274 http:https://dx.doi.org/10.1103/PhysRevE.57.7274.

—
Reply to this email directly or view it on GitHub
#1 (comment)
.

[joeydumont]

That is correct. I will try to implement the new algorithm this week, but I am a little short on time.

If you need such high quantum numbers, however, it might be worth looking into some approximation schemes or semiclassical results that might be applicable in your specific case.

[jeffzhen]

Hi Joey,
Thanks for the suggestion! I'd be happy to implement some approximated
algorithms, if you could kindly point me to some publication that contain a
numerical recipe. I'm so poorly knowledged in this field that I don't know
what the right key words are to search for. I'd be happy to contribute such
code to your repository as an alternative option for large l.

Thanks a lot!
Jeff

On Wed, Jul 16, 2014 at 5:09 PM, Joey Dumont [email protected]
wrote:

That is correct. I will try to implement the new algorithm this week, but
I am a little short on time.

If you need such high quantum numbers, however, it might be worth looking
into some approximation schemes or semiclassical results that might be
applicable in your specific case.

—
Reply to this email directly or view it on GitHub
#1 (comment)
.

[jeffzhen]

Hi Joey,
I just realized that I might have misunderstood your suggestion. If so,
please ignore my message.

Thanks a lot!
Jeff

On Wed, Jul 16, 2014 at 5:16 PM, Jeff Zheng [email protected] wrote:

Hi Joey,
Thanks for the suggestion! I'd be happy to implement some approximated
algorithms, if you could kindly point me to some publication that contain a
numerical recipe. I'm so poorly knowledged in this field that I don't know
what the right key words are to search for. I'd be happy to contribute such
code to your repository as an alternative option for large l.

Thanks a lot!
Jeff

On Wed, Jul 16, 2014 at 5:09 PM, Joey Dumont [email protected]
wrote:

That is correct. I will try to implement the new algorithm this week, but
I am a little short on time.

If you need such high quantum numbers, however, it might be worth looking
into some approximation schemes or semiclassical results that might be
applicable in your specific case.

—
Reply to this email directly or view it on GitHub
#1 (comment)
.

[joeydumont]

Hi,

I don't know of any specific numerical recipe, but you might these (and the references therein) useful. The first one might have some useful series approximations and such, while the others are more particular to the current algorithm, but may contain useful references.

• K. Schulten, "Semiclassical approximations to 3j- and 6j-coefficients for quantum-mechanical coupling of angular momenta," J. Math. Phys. 16, 1971 (1975).
• K. Schulten, "Exact recursive evaluation of 3j- and 6j-coefficients for quantum-mechanical coupling of angular momenta," J. Math. Phys. 16, 1961 (1975).
• J. Luscombe and M. Luban, "Simplified recursive algorithm for Wigner 3j and 6j symbols," Phys. Rev. E 57, 7274–7277 (1998).

[joeydumont]

It turns out that the problems were caused by the sgn() function. The vector function should now return proper values.

[jeffzhen]

Hi Joey,
Exciting news! While it fixes a few cases in my tests, maybe I'm being
dumb, but it seem to be returning the same 0s as before for some of my test
cases. Would you mind testing out this specific combination, which I am
still getting zeros:
WignerSymbols::wigner3j(529, 992, 1243, 196, -901, 705)

which I suppose the vector version should be
WignerSymbols::wigner3j(992, 1243, 196, -901, 705)

Thanks a lot!
Jeff

On Fri, Jul 18, 2014 at 4:00 PM, Joey Dumont [email protected]
wrote:

It turns out that the problems were caused by the sgn() function. The
vector function should now return proper values.

—
Reply to this email directly or view it on GitHub
#1 (comment)
.

[joeydumont]

Since this is caused by another problem, I've opened up a new issue. See #2.

[jeffzhen]

Hi Joey,
I can't tell you how grateful I am for your constant effort in fixing the
bugs! I am very much looking forward to your upcoming fix!

Thanks a lot!
Jeff

On Fri, Jul 18, 2014 at 6:27 PM, Joey Dumont [email protected]
wrote:

Since this is caused by another problem, I've opened up a new issue. See
#2 #2.

—
Reply to this email directly or view it on GitHub
#1 (comment)
.