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

Tobias Burnus burnus@net-b.de
Sun Aug 22 09:55:00 GMT 2010


  As HJ and Dominique reported, for Linux/ia32 and 
powerpc-apple-darwin9, the tolerance values were too tight for the 
comparison of the recurrence algorithm with the library version.

I increased now the tolerance based on Dominique's values and cross 
checked using -m32 on x86-64-linux. I hope it will now work on most 
systems, but we still might need to tweak the values (or use xfail) for 
some architectures as the initial values are slightly different and as 
the architectures might do different rounding or take different code paths.

Thanks to Dominique for testing and adjusting the values - and to HJ for 
keeping an eye on regressions. Committed as Rev. 163454/163455.

  * * *

Note of warning to users: Be careful when applying the recurrence 
algorithm (i.e. when using BESSEL_JN/YN(N1, N2, X): The results can be 
quite different! For instance, using
   BESSEL_JN(0, 100, 475.78)
the deviation is < 6*epsilon(0.0) - except for the following N, where it 
is (much) larger:
N=63: < 9*epsilon(0.0)
N=91: < 13*epsilon(0.0)
N=49: < 326*epsilon(0.0)

Note in particular that checking the last item in the recurrence 
algorithm - for JN that's lowest, here: N=0, does not help: For the 
example above the difference to the value returned by the library for JN 
is just 0.5*epsilon while for N=49 the difference is much larger! (The 
reason is that for N=49 the result is -0.8E-04 while the other values 
are of the order 0.xE-01 and 0.xE-02.)
.
Tobias
-------------- next part --------------
A non-text attachment was scrubbed...
Name: commit.diff
Type: text/x-patch
Size: 2427 bytes
Desc: not available
URL: <http://gcc.gnu.org/pipermail/fortran/attachments/20100822/86123e50/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: commit2.diff
Type: text/x-patch
Size: 1633 bytes
Desc: not available
URL: <http://gcc.gnu.org/pipermail/fortran/attachments/20100822/86123e50/attachment-0001.bin>


More information about the Fortran mailing list