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

[Steve Kargl <sgk at troutmask dot apl dot washington dot edu>]

On Fri, Aug 20, 2010 at 10:32:21AM +0200, Tobias Burnus wrote:
>  On 08/20/2010 08:05 AM, Tobias Burnus wrote:
> > Steve Kargl wrote:
> >>There may be a method to improve the values computed, but it
> >>would substantially increase the computation.
> 
> I just have realized that I did a stupid comparison for the calculated 
> values! I should have compared
> abs ((YN-YN2)/YN) with epsilon - without the "/YN" the comparison is not 
> useful and will fail if abs(YN) is a large value - even though the 
> relative error is small. The following was tested against the run-time 
> version libgfortran (PR 36158) and the results seem to be sufficiently 
> accurate. Thus, I think one can live with the current algorithm.
> 
> The only think which we could handle explicitly is the -INF vs. NAN, cf. 
> below. Do you think that makes sense?
> 
> Test program (needs the patch from the PR for libgfortran as X is not 
> known at compile time):
> 
> implicit none
> real,parameter :: values(*) = [0.5, 1.0, 0.9, 
> 1.8,2.0,3.0,4.0,4.25,8.0,34.53, 475.78]
> integer, parameter :: Nmax = 100
> real :: rec(0:Nmax), lib(0:Nmax)
> integer :: i
> 
> do i = 1, ubound(values,dim=1)
>   call compare(values(i))
> end do
> contains
> subroutine compare(X)
> integer :: i
> real X
> rec = BESSEL_YN(0, Nmax, X)
> lib = [ (BESSEL_YN(i, X), i=0,Nmax) ]
> print *, 'YN for X = ', X, ' -- Epsilon = ',epsilon(x)
> do i = 0, Nmax
>   print '(i2,2e17.9,e12.2,f14.10,2l3)', i, rec(i), lib(i),&
>         rec(i)-lib(i), ((rec(i)-lib(i))/rec(i))/epsilon(x),&
>         abs((rec(i)-lib(i))/rec(i))<  epsilon(x),&
>         abs((rec(i)-lib(i))/rec(i))<  3*epsilon(x)
> end do
> end
> end
> 
> 
> Result as follows. (Details see attachment -- there also for JN).
> 
> a) As soon as the result is -Infinity at N, the N+2 step produces NaN
> instead of -INF - and thus a completely different result.
> One could add a check of the kind:
> 
> if (last2 == -INF)
>   setall remaining items to -INF
> 
> What do you think about adding such a check for YN?
> 

Returning an INF is probably preferable to an NaN,
but I'm not going to require you to implement this
before committing your patch.

> b) Until the result is -INF, the deviations are:
> x = 0.0: Exact (special case)
> x = 0.5:<  epsilon  (up to i = 25, above: -Inf)
> x = 1.0:<  epsilon  (up to i = 29 above: -Inf)
> x = 1.8:<  epsilon until  i=11, then<  3*epsilon (up to i = 34 above: -Inf)
> x = 2.0:<  epsilon  (up to i = 34 above: -Inf)
> x = 3.0:<  epsilon until  i=5,
>          <  3*epsilon until i=17
>          <  8*epsilon (up to i = 39 above: -Inf)
> x = 4.0:<  epsilon  (up to i = 43, above: -Inf)
> x = 4.25:<  epsilon until  i=34, then<  2*epsilon (up to i = 44 above: -Inf)
> x = 8.0:<  epsilon  (up to i = 55, above: -Inf)
> x = 34.53:<  epsilon until  i=5
>            <  6*epsilon until i=25
>            <  22*epsilon until i=100
> x = 475.78:<  epsilon until  i=7
>              <  2*epsilon until i=26
>              <  15*epsilon until i=62
>              <  28*epsilon until i=82
>              <  130*epsilon until i=100
 
The rather limited range of i before one hits an INF is
why I was surprised J3 included the transformational form
in F2008.  IIRC, you raised a concern with (or request to 
removei from?) F2008 over the transformational form and
J3 rejected it saying something about recursion is a
common method to compute besseli functions.  I think J3
underestimated the care that is needed and the fact that
the range is rather limited without specialized algorithms.

-- 
steve