[C] Legendre Polynomial (by recurrence) programming

A simple legendre polynomial computation in C/C++.

The code

/**
* Berenger (contact at berenger dot eu)
* This is the source code to construct the legendre polynome in C
* This is fast but you can improve the code by using pointer instead of
* accessing using index on the array and to compute (2*l-1) with a recurrence.
* Ref: Fast and accurate determination of the Wigner rotation matrices in FMM
*/

#include <cmath>
#include <cstdio>

// To use double or float
typedef double FReal;

// To access the P_l,m in a array
// {0,0}{1,0}{1,1}{2,0}{2,1}{2,2} ...
int atLm(const int l, const int m){
    // summation series over l + m => (l*(l+1))/2 + m
    return ((l*(l+1))>>1) + m;
}


/**
* Compute the legendre polynome by recurrence.
* If needed you can use:
* P_l,m (−x) = (−1)^m+l P_l,m(x)
* Also:
* P_l,-m (x) = (-1)^m * (l-m)!/(l+m)! * P_l,m(x)
*/
void computeLegendre(FReal legendre[], const FReal x, const int P){
    //This factor is reuse −sqrt(1 − x^2)
    const FReal factor = -sqrt(1.0-pow(x,2));

    // Init legendre
    legendre[atLm(0,0)] = 1.0;        // P_0,0(x) = 1
    // Easy values
    legendre[atLm(1,0)] = x;	  // P_1,0(x) = x
    legendre[atLm(1,1)] = factor;     // P_1,1(x) = −sqrt(1 − x^2)

    for(int l = 2; l <= P ; ++l ){
        for( int m = 0; m < l - 1 ; ++m ){
            // P_l,m = (2l-1)*x*P_l-1,m - (l+m-1)*x*P_l-2,m / (l-k)
            legendre[atLm(l,m)] = (FReal(2*l-1) * x * legendre[atLm(l-1,m)] - FReal( l + m - 1 ) * legendre[atLm(l-2,m)] )
                    / FReal( l - m );
        }
        // P_l,l-1 = (2l-1)*x*P_l-1,l-1
        legendre[atLm(l,l-1)] = FReal(2*l-1) * x * legendre[atLm(l-1,l-1)];
        // P_l,l = (2l-1)*factor*P_l-1,l-1
        legendre[atLm(l,l)] = FReal(2*l-1) * factor * legendre[atLm(l-1,l-1)];
    }
}

/**
* To test the legendre function
*/
int main(){
    static const int P = 2;
    FReal legendre[((P+2)*(P+1))/2];

    computeLegendre(legendre, 0.5, P);

    for(int l = 0 ; l <= P ; ++l){
        for(int m = 0 ; m <= l ; ++m){
            printf("P{%d,%d} = %lf t", l, m, legendre[atLm(l,m)]);
        }
        printf("n");
    }

    return 0;
}

Output

$ ./test
P{0,0} = 1.000000
P{1,0} = 0.500000 	P{1,1} = -0.866025
P{2,0} = -0.125000 	P{2,1} = -1.299038 	P{2,2} = 2.250000