Prev: interactive command editting
Next: Verifying function interface
From: mecej4 on 28 Feb 2010 01:07
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] cc -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

The problem is present with those compilers that use Bessel functions
from libm, such as GCC, GFortran, G77 on Linux-IA32, Linux-IA32, Linux
AMD-64, Linux-Itanium2 and Linux Z90. Oddly enough, GCC 4.3.2 on Cygwin
_does_not_ exhibit the problem.

On Windows, Microsoft VC has the problem. The Intel C compiler gives
correct results but warns (through macros in the Microsoft math.h header
file) that jn() is deprecated and should be replaced by _jn(). If one
does this, however, the Intel compiler, having only jn() instead of
_jn() in its libmmd.lib, uses the _jn() from Microsoft's libcxx library
and displays the erroneous results. Some reward for taking good advice!

As an independent check on your results, I used D.E. Amos's routine
ZBSUBS from TOMS 644 ( www.netlib.org/TOMS/644 ) with the following driver

program testjn
real*8 fnu,zr,zi
real*8 CYR(8),CYI(8)
zr = 2.4048255576957729d0 ! real part of arg
zi = 0d0 ! imag. part of arg
FNU = 2d0 ! starting order \nu = 2
KODE = 1
N = 8 ! number of integral values of \nu to compute for

call zbesj(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR)

write(*,*)' Ierr = ',ierr,' NZ = ',nz
write(*,'(1x,I2,1p,2D14.6)')(n+1,CYR(n),CYI(n),n=1,8)
end

and obtained results that agreed with your results to all digits. The
same perfect agreement is seen with results from Matlab and Maple.

-- mecej4
From: Richard Maine on 28 Feb 2010 01:30
mecej4 <mecej4_nspm(a)operamail_nspm_.com> wrote:

> As an independent check on your results, I used D.E. Amos's routine
> ZBSUBS from TOMS 644 ( www.netlib.org/TOMS/644 )

But that's written in Fortran, so, if I interpret Steve correctly, it is
"impossible" for it to be a "general purpose" "implementation"... or
something like that. Maybe someone should tell the ACM folk. :-)

P.S. The url appears to need the TOMS in lower case.

--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
From: FX on 28 Feb 2010 15:15
> My feeling is that the new functions only ended up in the Fortran
> standard after being standardized in the C and thus part of the systems
> C library (e.g. glibc).

It's not the case for all of them. At least ERFC_SCALED is a function
which has no equivalent in standard system libm's.

> Talking about math functions. Does anyone know what happened to the
> proposal to support the extra math functions to Fortran, which were
> specified in C's TR 24747 [1,2]? I remember seeing something about it
> on the J3 mailing list, though I have not heard anything lately. Does
> anyone know a (C) library / C compiler which implements those already?

The GNU libstdc++ has them (at least in the trunk version). They do not
rely on system implementations, and were written (mostly) by Edward
Smith-Rowland.

I'd they they could be used for gfortran (licence is identical to that of
libgfortran, i.e. GPL + exception).

> (The TR provides Laguerre/Legendre/Hermite polynominals, beta/zeta
> function, elliptical integrals, cylindrical/spherical Bessel/Neumann
> functions.)

It would really be great to have those in Fortran. I mean, really really.

--
FX
First  |  Prev  | 
Pages: 1 2 3 4 5 6
Prev: interactive command editting
Next: Verifying function interface