[Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic

[Steve Kargl]

On Fri, Aug 20, 2010 at 08:05:58AM +0200, Tobias Burnus wrote:
>  Steve Kargl wrote:
> >There may be a method to improve the values computed, but it
> >would substantially increase the computation.  Instead of
> >recursing on the actually estimates of the bessel function,
> >we could recurse on the significands by careful bookkeeping
> >of the exponents (which may require an additional array), e.g.,
> >
> [...]
> 
> I essentially do not use Bessel functions - and the few times I did only 
> for few N. I think gfortran's transformational function should use 
> something faster than doing N2-N1+1 calls to a full JN/YN function - but 
> if N2-N1+1 is a large number, currently the accuracy seems to be too low.
> 
> As you regularly seem to use Bessel functions (though mostly spherical 
> ones?), would you be satisfied by the current implementation in 
> gfortran?

Well, I may not be the best person to answer this question because
I tend to be suspicious of 'black boxes', so I write simply benchmarks
and test code.    

> Or would you use the book keeping version, which is in terms 
> of performance between the current recurrence version and the calling 
> the library function every time?

I think what you have now is good enough.  

> What are typical (larger) values for 
> N2-N1+1? Is this more likely 5 or 10 or 20 or 1000?

In classical field theory (both acoustic and EM), scattering
solutions are often written as infinite summations of partial
wave functions.  If one defines the dimensionless parameter
kd (k is the wave number and d is characteristic dimension),
then the seriesi usually can be truncated to N > int(kd) + m where 
m is usual around 20 or so.  For Mie scattering (eg, the cause
of rainbows), kd can be around 15000.  One needs to use
special techniques to do the summation, so one would most likely
use his own implementation.  For the acoustics research I do,
N is normally less than 500.

> I think the goal should be to have a version which is sufficiently 
> robust to be usable while at the same time being faster than calling the 
> elemental function N2-N1+1 times.
> 
> As you are a heavy user of Bessel functions, I would like if you could 
> decide on the algorithm - thus, what do you suggest?

I think what you have is sufficient.  Power users will use their
own specialized implementations.

-- 
Steve