This is the mail archive of the
`gcc-patches@gcc.gnu.org`
mailing list for the GCC project.

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |

Other format: | [Raw text] |

*From*: Steve Kargl <sgk at troutmask dot apl dot washington dot edu>*To*: Tobias Burnus <burnus at net-b dot de>*Cc*: Tobias Burnus <tobias dot burnus at physik dot fu-berlin dot de>, gcc-patches at gcc dot gnu dot org, fortran at gcc dot gnu dot org*Date*: Fri, 20 Aug 2010 13:13:43 -0700*Subject*: Re: [Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic*References*: <20100818103725.GA1023@physik.fu-berlin.de> <20100819033446.GA46540@troutmask.apl.washington.edu> <4C6CDD82.7090106@net-b.de> <20100820010452.GA52279@troutmask.apl.washington.edu> <4C6E1B46.2010004@net-b.de> <4C6E3D95.2080401@net-b.de>

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

**Follow-Ups**:

**References**:**Re: [Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic***From:*Tobias Burnus

**Re: [Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic***From:*Steve Kargl

**Re: [Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic***From:*Tobias Burnus

**Re: [Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic***From:*Steve Kargl

**Re: [Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic***From:*Tobias Burnus

**Re: [Patch,Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic***From:*Tobias Burnus

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |