New extension for Associated Legendre Functions (ALFs)

New extension for Associated Legendre Functions (ALFs)

[Patrick Alken]

Hello all,

  Most current implementations of associated Legendre functions
(ALFs), including the current GSL version fail due to underflow
above degree and order (L,M) 1900.  A 2002 paper by Holmes and
Featherstone shows that with a very simple scaling modification
this limit can be extended to 2700 in double precision. I have
recently needed to compute high degree ALFs so I have implemented
their algorithm and made a new GSL extension if others are interested.

There are several other advantages in this extension over the
current GSL implementation:

1) I have included support for several normalization conventions:
   - Schmidt semi-normalized ALFs
   - Spherical harmonic normalized ALFs (GSL already has this)
   - Fully normalized ALFs
   - Unnormalized ALFs (GSL has this)

2) The user can choose whether to compute the Condon-Shortley
   phase of (-1)^m

3) The GSL array versions currently only compute ALFs for a fixed
   order m and multiple degrees l. This extension now efficiently
   computes all ALFs for a given maximum degree lmax for all
   (l,m) up to lmax.

4) The routines in this extension use a workspace which precomputes
   various factors in the recurrence relations which should lead
   to a (minor) speedup in the calculations.

This extension includes a complete test suite and documentation.
The docs contain a plot which shows how the current GSL
implementation fails for L,M = 2700.

I have put a link on the main GSL webpage for the extension.

Patrick Alken

Re: New extension for Associated Legendre Functions (ALFs)

[Brian Gough]

At Mon, 16 Mar 2009 15:31:37 -0600,
Patrick Alken wrote:
>   Most current implementations of associated Legendre functions
> (ALFs), including the current GSL version fail due to underflow
> above degree and order (L,M) 1900.  A 2002 paper by Holmes and
> Featherstone shows that with a very simple scaling modification
> this limit can be extended to 2700 in double precision. I have
> recently needed to compute high degree ALFs so I have implemented
> their algorithm and made a new GSL extension if others are interested.

Great work!

Incidentally, I found a copy of their paper online at
http://www.cage.curtin.edu.au/~geogrp/research.html
http://www.cage.curtin.edu.au/~will/20760279.pdf

-- 
Brian Gough

Re: New extension for Associated Legendre Functions (ALFs)

[Patrick Alken]

I should also mention that I played around with computing
ALFs using quad precision in the recurrence relations. In
doing this you can avoid underflow problems to degrees up
in the 10s of thousands (I only tried up to L,M = 20000).

But unfortunately I don't think sizeof(long double) is standard,
though I do think its guaranteed to be at least 10 bytes
by the latest ansi standard. So for now I left everything
in double precision.

Patrick

On Tue, Mar 17, 2009 at 11:28:52AM +0000, Brian Gough wrote:
> At Mon, 16 Mar 2009 15:31:37 -0600,
> Patrick Alken wrote:
> >   Most current implementations of associated Legendre functions
> > (ALFs), including the current GSL version fail due to underflow
> > above degree and order (L,M) 1900.  A 2002 paper by Holmes and
> > Featherstone shows that with a very simple scaling modification
> > this limit can be extended to 2700 in double precision. I have
> > recently needed to compute high degree ALFs so I have implemented
> > their algorithm and made a new GSL extension if others are interested.
> 
> Great work!
> 
> Incidentally, I found a copy of their paper online at
> http://www.cage.curtin.edu.au/~geogrp/research.html
> http://www.cage.curtin.edu.au/~will/20760279.pdf
> 
> -- 
> Brian Gough