legendre_T(int n, int l) coeffecients not working

[max91219]

Hi

I have tried to use the following function to transform some Legendre representation greens functions to imfreq however the coefficients appear to be wrong. I think this has something to do with using the following to get the values for negative n

T{-n,l} = T{n,l}^*

as I believe this should be:

T{-n,l} = T{-(n+1),l}^*

I have tested this by testing that the transformation is unitary as per equation E3 of the following paper Orthogonal polynomial representation of imaginary-time Green’s functions however when I do this using the current function in TRIQS I get the following result

__ CURRENT _ 0 0 0.594695 0 1 2.04668e-33 0 2 0.195572 1 0 2.04668e-33 1 1 0.507233 1 2 -2.06978e-33 2 0 0.195572 2 1 -2.06978e-33 2 2 0.905484 __ FIX _ 0 0 0.99998 0 1 0 0 2 -4.42049e-05 1 0 0 1 1 1 1 2 0 2 0 -4.42049e-05 2 1 0 2 2 0.999901

where we would expect a kronecker delta on the values of l and l' (the first two values printed on each line)

this is the code that was used in order to produce the above results

#include <triqs/utility/legendre.hpp>
#include <boost/math/special_functions/bessel.hpp>

inline std::complex<double> legendre_T_nl(int n, int l);

int main(int argc, char **argv) {

    std::complex<double> res;
    std::cout << "__________ CURRENT _________" << std::endl;

    for (int l = 0; l <= 2; l++) {
        for (int l_p = 0; l_p <= 2; l_p++) {
            for (int n = -10250; n <= 10250; n++) {
                res += std::conj(triqs::utility::legendre_T(n, l)) * triqs::utility::legendre_T(n, l_p);
            }
            std::cout << l << "    " << l_p << "    " <<  res.real() << std::endl;
            res = 0.0;
        }
    }

    std::cout << "__________ FIX _________" << std::endl;

    for (int l = 0; l <= 2; l++) {
        for (int l_p = 0; l_p <= 2; l_p++) {
            for (int n = -10250; n <= 10250; n++) {
                res += std::conj(legendre_T_nl(n, l)) * legendre_T_nl(n, l_p);
            }
            std::cout << l << "    " << l_p << "    " <<  res.real() << std::endl;
            res = 0.0;
        }
    }

    return 0;
}

inline std::complex<double> legendre_T_nl(int n, int l) {

    const std::complex<double> i_c(0.0, 1.0);
    const double pi = boost::math::constants::pi<double>();

    bool neg_n = false;
    if(n < 0) {
        n = std::abs(n+1);
        neg_n = true;
    }

    double minus_fac = n%2 == 0 ? 1 : -1;
    double arg = ((2.0 * n + 1.0) * pi) / 2.0;

    std::complex<double> res = minus_fac * pow(i_c, l + 1) * sqrt(2 * l + 1.0) * boost::math::sph_bessel(static_cast<unsigned int>(l), arg);

    return neg_n ? std::conj(res) : res;
}```

[krivenko]

Thank you for the fix!