Re: (very) slow recursive function


Ralf Schaa wrote:
Hi there,

for computing Bessel functions and Legendre polynomials there exists
recursive formulas (see for example Numerical Recipes in F77, p.172),
which can be implemented straightforward using only a few lines of
Fortran. Comparing with iterative (and much more sophisticated)
versions, the recursive versions are *really* slow.

The example below computes the Legendre polynomials recursively as
given in equations (6.8.7), (6.8.8) and (6.8.9) from Numerical Recipes
in F77, p.247

Only for 'l' up to 20 this version is suitable, for higher orders it is
really slow and takes ages. The iterative version from Numerical
Recipes however has no problems at all. Similar problems are
encountered with the computation of the modified Bessel functions with
fractional order ...

So my question: why is this so slow compared with the iterative
versions. Is this an 'inherently' problem of recursion? Is it my
program?

Cheers
-Ralf

Example:

pure recursive function legendrePoly(x,m,l) result(Plm)

! Computes the (associated) Legendre polynomial Plm(x)
! with (0 <= m <= l) and (-1 <= x <= 1) after the recursive
! formula in Numerical Recipes in Fortran 77 (Volume 1, p. 247).
! 'Ordinary' Legendre polynomals are given for m=0.
!
! Calls: pure recursive function ifact -> double factorial n!!

implicit none

real(dp) :: Plm
real(dp), intent(in) :: x
integer, intent(in) :: m
integer, intent(in) :: l

if(l == m) then
Plm = (-1)**m * ifact(2.0d0*m-1.0d0) * (1.0d0-x**2)**(m/2.0d0)
elseif(l == m-1) then
Plm = 0.0d0
else
Plm = (x*(2.0d0*l-1.0d0)*legendrePoly_(x,m,l-1) - &
(l+m-1.0d0)*legendrePoly_(x,m,l-2)) / (l-m)
endif
end function legendrePoly

First of all, an aside on the comment in your code: The condition m <=
l is neither necessary (this is a property of the Legendre functions)
nor honored during the course of the recursion.

If you wish to retain the benefits of recursion, but wish to improve
efficiency, a simple means is to maintain two associated global arrays,
one boolean (b) and another real (p). In the user code that calls Plm,
initialize the booleans to false before calling Plm. At the top of your
function Plm, check to see whether b(l) is true; if so, set Plm to p(l)
and return immediately. If not, use the recursion rule and update b(l)
and p(l) before returning.

If the user code makes more than one call to Plm, the initialization
needs to be done every time m or x is changed, or b and p will need to
become two or three dimensional arrays.

Alas, this bookkeeping improves the efficiency but makes the Plm
function impure.

.

Relevant Pages

  • (very) slow recursive function
    ... recursive formulas (see for example Numerical Recipes in F77, p.172), ... The example below computes the Legendre polynomials recursively as ... Is this an 'inherently' problem of recursion? ...
    (comp.lang.fortran)
  • Re: (very) slow recursive function
    ... for computing Bessel functions and Legendre polynomials there exists ... recursive formulas (see for example Numerical Recipes in F77, p.172), ... The example below computes the Legendre polynomials recursively as ... Is this an 'inherently' problem of recursion? ...
    (comp.lang.fortran)
  • Re: (very) slow recursive function
    ... recursive formulas (see for example Numerical Recipes in F77, p.172), ... The example below computes the Legendre polynomials recursively as ... Is this an 'inherently' problem of recursion? ... Do this 20 times and you have 2**20 executions of the ...
    (comp.lang.fortran)