Re: (very) slow recursive function
- From: "robin" <robin_v@xxxxxxxxxxx>
- Date: Tue, 26 Dec 2006 23:09:59 GMT
"Ralf Schaa" <schaa@xxxxxxxxxxxxxxxx> wrote in message
news:1162174976.068293.26740@xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
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?
No. But note that each time your procedure is executed
and the final assignment is executed, TWO recursive calls
are executed. Thus, on consecutive invocations,
you have 1, 2, 4, 8, 16, ..., invocations.
Do this 20 times and you have 2**20 executions of the
procedure.
Is it my program?
Partly.
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
(-1)**M is an inefficient way to change sign.
* ifact(2.0d0*m-1.0d0)
2.0D0*M is an inefficient way to double M (which is INTEGER).
(2*m-1) requires no conversions to float.
Similar remarks for expression for Plm below.
* (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
.
- Prev by Date: Re: Bit Concatenation
- Next by Date: Re: Popularity contests
- Previous by thread: Re: (very) slow recursive function
- Index(es):
Relevant Pages
|
|
