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

[Tobias Burnus]

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

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?

Tobias

PS: As Steve has already found out: PR 36158 contains a run-time 
implementation of the transformational BESSEL_JN/YN(N1, N2, X). Note: I 
would like to use the same algorithm for both the library and the 
simply.c version.