Re: (very) slow recursive function
- From: Louis Krupp <lkrupp@xxxxxxxxxxxxxxxxxxxxxxx>
- Date: Mon, 30 Oct 2006 01:34:55 -0700
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
Recursion is *generally* slower than iteration because function calls are *generally* slower than inline code.
Looking at your example in particular, I'm not sure what the underscore in legendrePoly_(x,...) does, but it looks like every non-trivial call to your function generates two more calls, so the number of calls is roughly 2**(l-m). Note also that for every l > m+1 (or something like that), you're computing legendrePoly_(x,m,l-2) twice: once for (x, m, l), and once for (x, m, l-1).
Then there's ifact; if you're computing the factorial recursively, that will slow you down, too.
I'm not sure what the compiler does with (-1) ** m, but if it generates code to do m multiplications of -1, it's going to be a bit slower than something like 1 - 2*mod(m, 2).
I believe that (1.0d0-x**2)**(m/2.0d0) is likely to be a bit slower than sqrt((1.0 - x**2)**m).
As a matter of style -- which you didn't ask about -- I would use 0.0 and 1.0 and 2.0 instead of 0.0d0 and so on in most places. The value is the same, and I think it's more readable without the d0.
Louis
.
- Follow-Ups:
- Re: (very) slow recursive function
- From: Arjen Markus
- Re: (very) slow recursive function
- References:
- (very) slow recursive function
- From: Ralf Schaa
- (very) slow recursive function
- Prev by Date: Re: READ-statement and a segmentation fault
- Next by Date: Re: Short records on unformatted reads
- Previous by thread: (very) slow recursive function
- Next by thread: Re: (very) slow recursive function
- Index(es):
Relevant Pages
|
