Fun with Friedrich Bessel
in [Fortran]
|
From: steve on 27 Feb 2010 12:45 On Feb 27, 6:28 am, user1 <us...(a)example.net> wrote: > steve wrote: > > On Feb 26, 2:09 am, Arjan<arjan.van.d...(a)rivm.nl> wrote: > >>> If your compiler returns a value other than 0.19899990535769083, > >>> you may want to contact your OS/compiler vendor. > > >> I added a loop of 1001 iterations to let z run from 0 to twice the > >> critical value shown (for readers unfamiliar with Bessel-functions, > >> 2.4048255576957729_dp is the value where J_0(z) is 0). On MinGW/MSYS, > >> g95 and gfortran, this resulted in nice, smooth series for J_3(z), > >> with the exception of a spurious value 0 for z = > >> 2.4048255576957729_dp. On linux, g95 gave the same result as on MinGW/ > >> MSYS, but gfortran changed the 0 at z = 2.4048255576957729_dp into - > >> Infinity... On linux, ifort (11.0) and Portland (old version 6.1) > >> produced complete bogus series for all values of z, with a value 0 for > >> the critical z! Interestingly enough, these two compilers agreed on > >> their bogus... > > >> It looks as if it is even worse than suggested in the first message! > > >> A. > > > It is indeed much worse. If you have this problem, then it may > > occur at every zero of J0(x). A longer explanation is > > >http://gcc.gnu.org/ml/fortran/2010-02/msg00213.html > > > In addition, all odd order Bessel functions will by -Inf. > > > troutmask:kargl[217] cat testjn.c > > #include<stdio.h> > > #include<math.h> > > > int > > main(void) > > { > > double z; > > int i, n; > > > z = 2.4048255576957729; > > for (n = 2; n< 10; n++) > > printf("%d %e\n", n, jn(n,z)); > > return (0); > > } > > troutmask:kargl[218] Solaris-o z testjn.c -lm&& ./z > > 2 4.317548e-01 > > 3 -inf > > 4 4.069027e-02 > > 5 -inf > > 6 3.247664e-03 > > 7 -inf > > 8 7.495602e-05 > > 9 -inf > > > With the patch I posted to the gfortran list, one gets > > > troutmask:kargl[215] ./z > > 2 4.317548e-01 > > 3 1.989999e-01 > > 4 6.474667e-02 > > 5 1.638924e-02 > > 6 3.404818e-03 > > 7 6.006884e-04 > > 8 9.216579e-05 > > 9 1.251727e-05 > > > -- > > steve > > Is it a gfortran problem or a system supplied library ? I got correct > values (match above) using the Sun Studio compiler on OpenSolaris, and > strange results using Sun Studio compiler on Linux. You answered your own question as neither 'Sun Studio compiler on OpenSolaris' nor 'Sun Studio compiler on Linux' are gfortran. Others have reported that ifort and g95 have this problem. The short story is that jn(n,x) [aka besjn(n,x)] in libm is derived from code from Netlib's fdlibm e_jn.c. FreeBSD's svn repository reveals that this bug is at least 15.5 years old. Oh, and for the record, I discovered this issue while using good old Fortran. :-) -- steve
From: glen herrmannsfeldt on 27 Feb 2010 14:56 steve <kargls(a)comcast.net> wrote: > On Feb 27, 6:28?am, user1 <us...(a)example.net> wrote: (snip) >> Is it a gfortran problem or a system supplied library ? I got correct >> values (match above) using the Sun Studio compiler on OpenSolaris, and >> strange results using Sun Studio compiler on Linux. (snip) > The short story is that jn(n,x) [aka besjn(n,x)] in libm is > derived from code from Netlib's fdlibm e_jn.c. FreeBSD's svn > repository reveals that this bug is at least 15.5 years old. There have been some bessel function routines in many C libraries over the years. As far as I know, not part of the standard. My (somewhat old) FreeBSD system has j0(), j1(), jn(), and the man page says that they appeared in Version 7 AT&T UNIX. > Oh, and for the record, I discovered this issue while using > good old Fortran. :-) It is interesting that they have been associated with Unix/C longer than with Fortran... -- glen
From: steve on 27 Feb 2010 15:16 On Feb 27, 11:56 am, glen herrmannsfeldt <g...(a)ugcs.caltech.edu> wrote: > steve <kar...(a)comcast.net> wrote: > > On Feb 27, 6:28?am, user1 <us...(a)example.net> wrote: > > (snip) > > >> Is it a gfortran problem or a system supplied library ? I got correct > >> values (match above) using the Sun Studio compiler on OpenSolaris, and > >> strange results using Sun Studio compiler on Linux. > > (snip) > > > The short story is that jn(n,x) [aka besjn(n,x)] in libm is > > derived from code from Netlib's fdlibm e_jn.c. FreeBSD's svn > > repository reveals that this bug is at least 15.5 years old. > > There have been some bessel function routines in many C libraries > over the years. As far as I know, not part of the standard. > > My (somewhat old) FreeBSD system has j0(), j1(), jn(), and > the man page says that they appeared in Version 7 AT&T UNIX. > > > Oh, and for the record, I discovered this issue while using > > good old Fortran. :-) > > It is interesting that they have been associated with Unix/C > longer than with Fortran... > Well, it was impossible to implement general purpose special functions in Fortran until the IEEE 754 intrinsic module was standardized. There was no way to treat exceptional values and subnormal numbers. It was also impossible to use the bit twiddling techniques (used in at least fdlibm) until TRANSFER was standardized. At this point in time, I see no reason to re-invent the wheel, and re-implement libm in Fortran. -- steve
From: glen herrmannsfeldt on 27 Feb 2010 15:51 steve <kargls(a)comcast.net> wrote: (snip) > Well, it was impossible to implement general purpose special > functions in Fortran until the IEEE 754 intrinsic module was > standardized. There was no way to treat exceptional values > and subnormal numbers. It was also impossible to use the > bit twiddling techniques (used in at least fdlibm) until > TRANSFER was standardized. > At this point in time, I see no reason to re-invent the > wheel, and re-implement libm in Fortran. I have a book called "Computation of Special Functions" by Zhang with Fortran 77 routines for many functions, including Bessel functions. It was written in 1996. You might find it useful. I think I found it on eBay when I was looking for something else. The price was good enough, and I figured I would someday find use for it. So far I only used it for erf() when working with a defective implementation of erf() in VBA. -- glen
From: Richard Maine on 27 Feb 2010 15:54
steve <kargls(a)comcast.net> wrote: > Well, it was impossible to implement general purpose special > functions in Fortran until the IEEE 754 intrinsic module was > standardized. "Impossible" seems like a bit of an overstatement, particularly insomuch as I'm sure counterexamples exist. Not that I'm going to bother to search for suitable citations of such counterexamples. If you want to claim things like awkward to do efficiently, I wouldn't much argue, but "impossible" is implausible, even if you are imagining an awfully restrictive definition of "general purpose special function". One also suspects you have a particular set of special functions in mind, as there are "general purpose special functions" that are trivial. -- Richard Maine | Good judgment comes from experience; email: last name at domain . net | experience comes from bad judgment. domain: summertriangle | -- Mark Twain |