recurseLegendrePoly computes associated Legendre functions recursively
function (n, x, Pnm_minus1, Pnm_minus2) {
    
    var Pnm = new Array(n+1);
    switch(n) {
        case 1:
            var x2 = numeric.mul(x,x);
            var P10 = x;
            var P11 = numeric.sqrt(numeric.sub(1,x2));
            Pnm[0] = P10;
            Pnm[1] = P11;
            break;
        case 2:
            var x2 = numeric.mul(x,x);
            var P20 = numeric.mul(3,x2);
            P20 = numeric.sub(P20,1);
            P20 = numeric.div(P20,2);
            var P21 = numeric.sub(1,x2);
            P21 = numeric.sqrt(P21);
            P21 = numeric.mul(3,P21);
            P21 = numeric.mul(P21,x);
            var P22 = numeric.sub(1,x2);
            P22 = numeric.mul(3,P22);
            Pnm[0] = P20;
            Pnm[1] = P21;
            Pnm[2] = P22;
            break;
        default:
            var x2 = numeric.mul(x,x);
            var one_min_x2 = numeric.sub(1,x2);
            // last term m=n
            var k = 2*n-1;
            var dfact_k = 1;
            if ((k % 2) == 0) {
                for (var kk=1; kk<k/2+1; kk++) dfact_k = dfact_k*2*kk;
            }
            else {
                for (var kk=1; kk<(k+1)/2+1; kk++) dfact_k = dfact_k*(2*kk-1);
            }
            Pnm[n] = numeric.mul(dfact_k, numeric.pow(one_min_x2, n/2));
            // before last term
            Pnm[n-1] = numeric.mul(2*n-1, numeric.mul(x, Pnm_minus1[n-1])); // P_{n(n-1)} = (2*n-1)*x*P_{(n-1)(n-1)}
            // three term recursence for the rest
            for (var m=0; m<n-1; m++) {
                var temp1 = numeric.mul( 2*n-1, numeric.mul(x, Pnm_minus1[m]) );
                var temp2 = numeric.mul( n+m-1, Pnm_minus2[m] );
                Pnm[m] = numeric.div( numeric.sub(temp1, temp2), n-m); // P_l = ( (2l-1)xP_(l-1) - (l+m-1)P_(l-2) )/(l-m)
            }
    }
    return Pnm;
}