Re: [Patch,Fortran,committed] PR 45367/36158 - Increase numeric tolerance in bessel_{6,7}.f90

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

On Sun, Aug 22, 2010 at 11:17:19PM +0200, Tobias Burnus wrote:
>  H.J. Lu wrote:
> >gfortran.dg/bessel_7.f90 still fails with -m32 -mfpmath=sse:
> >
> >http://gcc.gnu.org/ml/gcc-testresults/2010-08/msg02282.html
> >http://gcc.gnu.org/ml/gcc-testresults/2010-08/msg02284.html
> 
> I think that's because we have different libm versions. And because the 
> recurrence algorithm is somehow not very stable for certain combinations 
> of N and X.
> 

The recurrence algorithm is stable.  What you are see
is a cancellation issue.  If j(n-1) is a small number
then the most significant bits in (2n/x)*j(n) and j(n+1)
are the same, so in j(n-1) = 2n/x * j(n) - j(n+1) one
loses precision.

When one refers to the stability of the recurrence
algorithm, one is normally interested in the propagation
of error.  For bessel functions, backwards recurrence
is stable for all x and n with the caveat that one needs
j(n) and j(n+1) to start the recurrence.  If n < x, then
forward recurrence is stable but once n exceeds x the
error starts to grow.
 
!
! Forward recurrence for bessel fcns.
!
program testing
  real j2, j0, j1
  j0 = besjn(0,2.)
  j1 = besjn(1,2.)
  do i = 1, 10
     j2 = (2 * i / 2.) * j1 - j0
     print *, j2, besjn(i+1,2.)
     j0 = j1
     j1 = j2
  end do
end program

laptop:kargl[214] gfc4x -o z test.f90 && ./z
  0.35283405      0.35283405    
  0.12894326      0.12894325    
  3.39957476E-02  3.39957215E-02
  7.03972578E-03  7.03962985E-03
  1.20288134E-03  1.20242906E-03
  1.77562237E-04  1.74944071E-04
  4.00543213E-05  2.21795508E-05
  1.42872334E-04  2.49234358E-06
  1.24579668E-03  2.51538637E-07
  1.23150945E-02  2.30428494E-08

There is a paper in an obscure applied math journal
that shows the error analysis and shows that the
stability of the recurrence.  I don't have that paper
handy, but I could probably dig it up if interested.

-- 
Steve