A faster implementation for float result.

[0382]

It seems that you just want float result, then https://github.com/0382/WignerSymbol is a faster implementation.

A simple benchmark code

#include "WignerSymbol.hpp"
#include <chrono>

using namespace util;

double compute_all_wigner_3j(int lmax)
{
    WignerSymbols wigner;
    wigner.reserve(lmax, "Jmax", 3);
    double sum = 0;
    for(int l1 = 0; l1 <= lmax; ++l1)
    {
        for(int l2 = 0; l2 <= lmax; ++l2)
        {
            for(int l3 = 0; l3 <= lmax; ++l3)
            {
                for(int m1 = -l1; m1 <= l1; ++m1)
                {
                    for(int m2 = -l2; m2 <= l2; ++m2)
                    {
                        for(int m3 = -l3; m3 <= l3; ++m3)
                        {
                            sum += wigner.f3j(2*l1, 2*l2, 2*l3, 2*m1, 2*m2, 2*m3);
                        }
                    }
                }
            }
        }
    }
    return sum;
}

int main()
{
    for(int lmax : {4, 8, 12, 16, 20})
    {
        auto now = std::chrono::system_clock::now();
        double sum = 0;
        for(int rep = 0; rep < 10; ++rep)
        {
            sum += compute_all_wigner_3j(lmax);
        }
        auto end = std::chrono::system_clock::now();
        auto dura = std::chrono::duration_cast<std::chrono::milliseconds>(end - now).count();
        std::cout << "lmax = " << lmax << ", time = " << dura / 10. << "ms, sum = " << sum << '\n';
    }
    return 0;
}

[Luthaf]

Thanks a lot for sharing this! It does look quite a bit faster indeed, I'll dig deeper to see how I could use the same algorithm here!

How far have you validated the coefficients though? I tried your code on the example in issue #7 (i.e. wigner.f3j(2*100, 2*300, 2*285, 2*2, 2*-2, 0)), and the result is -0.0 instead of the expected 0.001979165708981953.

For smaller j1/j2/j3 values, everything seems fine.

[0382]

It seems that my code makes something wrong. Thank you for the test, I will check it.

[0382]

I have debugged my code, and find the problem comes from vary large float numbers multiply overflow to Inf. I will not fix this, because my code is just designed for fast numeric calculation, and the large J is not a common angular momentum in real world physical system. I also test the GSL library, which also does also not give the right result. The algorithm in your code is better for large J. Thanks again for your test.

[Luthaf]

Sounds good! I'll add your code to the set of benchmarks in this repo as well.