Legendre polynomials

The built-in legendre() calculates the Legendre polynomials calculated ALL the orders for a given degree. If you only need a given order, this is a waste of memory and computing time (especially for large blocks of data). The function legendreP(l,m,x) is a drop-in substitute for legendre(l,x), except that it only calculates the required order.

The polynomial coefficients are calculated analytically. The coefficients are computed recursively from earlier coefficients (to avoid computing multiple factorials). The polynomial is computed using exactly l/2 multiplications and additions for any structure of the data array x.

The polynomial coefficients can be very large for moderate/high degrees, leading to reduced precision. For these cases, it is better to use recursive formulas. A second optional output estimates the error, based on the largest polynomial coefficient.

The numerical returns were checked for every order for all degrees up to m using an array of size 128^1, filled with random numbers between -1 and 1; the difference compared with the built-in were within numerical errors (relative error ~1e-15). Some of the orders were checked for all degrees up to l=20 as well.

This function is much faster than the built in, especially for larger degrees; for 100^3 array, the built-in legendre(8,x) takes about 9.5 sec, while legendreP(8,6,x) takes only 0.15sec. Even calculating all 9 orders takes only about 1.5 sec.

Associated Legendre polynomials are defined for 0<=m<=l, and are 0 by definition if |m|>l, rather than returning an error. For |x|>1, returns NaN rather than error; comment this line if you don't need this check.

There is no additional scaling is done (e.g., Schmidt), but these should not affect significantly the run time. Uses only standard matlab features, should work with any version of matlab; tested on 2008b, 2009b, 2010a.