Prev: Treacherous Behavior of "gfortran" Compiler?
Next: Build project in CVF6.6c by using hotkey
From: Janne Blomqvist on 6 May 2008 13:32
On 2008-05-06, James Van Buskirk <not_valid(a)comcast.net> wrote:
> "Janne Blomqvist" <foo(a)bar.invalid> wrote in message
> news:slrng204af.eif.foo(a)vipunen.hut.fi...
>
>> On 2008-05-05, James Van Buskirk <not_valid(a)comcast.net> wrote:
>
>>> While we are on the
>>> subject of Bessel functions, did you notice that gfortran's
>>> initialization expression evaluation mechanism produced weird
>>> low-precision approximations of Bessel functions in the examples
>>> where it worked?
>
>> I suspect this will probably be fixed as a side effect of using MPFR
>> to evaluate bessel functions in initialization expressions, see
>
>> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36117
>
> But I think MPFR is already being used and that's why the trailing
> digits of the decimal expansion are zeros.

No, it was evaluated at runtime. That's also the reason why the
testcase you gave that needed the result at compile time failed.

You'll notice that I used the past tense above. Tobias Burnus posted a
patch for PR 36117 a little less than an hour ago, and it was
simultanously approved within minutes by FX and yours truly; finally
it was committed 5 minutes ago (as of this writing).

I guess the reason for the trailing zeros are probably that your
platform libc (as I understand it mingw uses the Microsoft C library)
doesn't provide long double (in C-speak, in Fortran on x86/x86_64 that
typically corresponds to 10 byte real) bessel functions, and mingw
then substitutes the double (8 byte real) versions instead.

You can see what happens by compiling with -fdump-tree-original. That
creates a file with a pseudo-C-like language which is what the
gfortran frontend feeds the gcc middle-end. Also see the -S compiler
option, which outputs the assembler output. The -fdump-tree-original
file shows calls to functions named something like "__builtin_j0l" and
so forth. On Linux, in the .s file one can then see that
e.g. __builtin_j0l has been translated to j0l, which is the long
double version of the j0 Bessel function. I guess that in your case,
this instead gets translated to j0, which is the double version; hence
the trailing zeros.

Now, with MPFR support initialization expressions will be fixed, but
if you still lack the long double libc functions needed for proper
runtime support I suppose it won't make it that much more useful for
high precision work.

--
Janne Blomqvist
From: James Van Buskirk on 6 May 2008 16:37
"James Van Buskirk" <not_valid(a)comcast.net> wrote in message
news:ed2dncmj0oLsK73VnZ2dnUVZ_o3inZ2d(a)comcast.com...

> I didn't see a PR for the f.p. control word.

It's http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323 and the 100
idiots who want PR323 'fixed' seem to have won out over sanity on
Win64 unlike all other platforms. Is there some way to 'unfix'
PR323 so that Win64 works like gawd & Billy Kahan intended?

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end


From: James Van Buskirk on 5 May 2008 17:53
"FX" <coudert(a)alussinan.org> wrote in message
news:fvk46j$1fuu$1(a)nef.ens.fr...

>> I would assume that the absence of BESSEL_J0, BESSEL_J1, BESSEL_Y0, and
>> BESSEL_Y1 for REAL(10) inputs for non-initialization expressions for
>> both Win32 and Win64 is a different issue altogether.

> This isn't, for now, treated as a bug: Windows C library doesn't provide
> these C99 intrinsic functions ({j,y}{0,1}l). It's been our policy until
> now that while we provide fallback versions of some C99 intrinsics when
> they're missing (for example, float variants or double variants), support
> for real types larger than double (real kinds 10 and 16) entirely relies
> on the system libraries for mathematical intrinsics. Maybe this should be
> documented somewhere, but I don't know really where. Suggestions as to a
> point where to insert this in our current documentation are welcome.

I know know where to put the documentation, either. How do you
document:

C:\gfortran\clf\bessel_test>type test1.f90
module mykinds
implicit none
integer, parameter :: ep = selected_real_kind(18,4931)
integer, parameter :: qp_preferred = selected_real_kind(33,4931)
integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1-sign(1,qp_preferred))/2*ep
end module mykinds

program test1
use mykinds
implicit none
real(qp), parameter :: x = 3
integer, parameter :: n = 1
real(kind(x)), parameter :: xj0 = BESSEL_J0(x)
real(kind(x)), parameter :: xj1 = BESSEL_J1(x)
real(kind(x)), parameter :: xjn = BESSEL_JN(n,x)
real(kind(x)), parameter :: xy0 = BESSEL_Y0(x)
real(kind(x)), parameter :: xy1 = BESSEL_Y1(x)
real(kind(x)), parameter :: xyn = BESSEL_YN(n,x)

write(*,*) 'xj0 = ', xj0
write(*,*) 'xj1 = ', xj1
write(*,*) 'xjn = ', xjn
write(*,*) 'xy0 = ', xy0
write(*,*) 'xy1 = ', xy1
write(*,*) 'xyn = ', xyn
end program test1

C:\gfortran\clf\bessel_test>gfortran test1.f90 -otest1

C:\gfortran\clf\bessel_test>test1
xj0 = -0.26005195490193345000
xj1 = 0.33905895852593645000
xjn = 0.33905895852593645000
xy0 = 0.37685001001279040000
xy1 = 0.32467442479180000000
xyn = 0.32467442479180000000

.... works, but not:

C:\gfortran\clf\bessel_test>type test2.f90
module mykinds
implicit none
integer, parameter :: ep = selected_real_kind(18,4931)
integer, parameter :: qp_preferred = selected_real_kind(33,4931)
integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1-sign(1,qp_preferred))/2*ep
end module mykinds

program test2
use mykinds
implicit none
real(qp), parameter :: x = 3
integer, parameter :: n = 1
real(kind(x)), parameter :: xj0 = BESSEL_J0(x)
real(kind(x)), parameter :: xj1 = BESSEL_J1(x)
real(kind(x)), parameter :: xjn = BESSEL_JN(n,x)
real(kind(x)), parameter :: xy0 = BESSEL_Y0(x)
real(kind(x)), parameter :: xy1 = BESSEL_Y1(x)
real(kind(x)), parameter :: xyn = BESSEL_YN(n,x)
real(int(xj0+kind(1.0)+1)) yj0
real(int(xj1+kind(1.0))) yj1
real(int(xjn+kind(1.0))) yjn
real(int(xy0+kind(1.0))) yy0
real(int(xy1+kind(1.0))) yy1
real(int(xyn+kind(1.0))) yyn

yj0 = 42
yj1 = 42
yjn = 42
yy0 = 42
yy1 = 42
yyn = 42
write(*,*) 'yj0 = ', yj0
write(*,*) 'yj1 = ', yj1
write(*,*) 'yjn = ', yjn
write(*,*) 'yy0 = ', yy0
write(*,*) 'yy1 = ', yy1
write(*,*) 'yyn = ', yyn
end program test2

C:\gfortran\clf\bessel_test>gfortran test2.f90 -otest2
f951.exe: internal compiler error: Segmentation fault
Please submit a full bug report,
with preprocessed source if appropriate.
See <http://gcc.gnu.org/bugs.html> for instructions.

Or even:

C:\gfortran\clf\bessel_test>type test3.f90
module mykinds
implicit none
integer, parameter :: ep = selected_real_kind(18,4931)
integer, parameter :: qp_preferred = selected_real_kind(33,4931)
integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1-sign(1,qp_preferred))/2*ep
end module mykinds

program test3
use mykinds
implicit none
real, parameter :: x = 3
integer, parameter :: n = 1
real(kind(x)), parameter :: xj0 = BESSEL_J0(x)
real(kind(x)), parameter :: xj1 = BESSEL_J1(x)
real(kind(x)), parameter :: xjn = BESSEL_JN(n,x)
real(kind(x)), parameter :: xy0 = BESSEL_Y0(x)
real(kind(x)), parameter :: xy1 = BESSEL_Y1(x)
real(kind(x)), parameter :: xyn = BESSEL_YN(n,x)
real(int(xj0+kind(1.0)+1)) yj0
real(int(xj1+kind(1.0))) yj1
real(int(xjn+kind(1.0))) yjn
real(int(xy0+kind(1.0))) yy0
real(int(xy1+kind(1.0))) yy1
real(int(xyn+kind(1.0))) yyn

yj0 = 42
yj1 = 42
yjn = 42
yy0 = 42
yy1 = 42
yyn = 42
write(*,*) 'yj0 = ', yj0
write(*,*) 'yj1 = ', yj1
write(*,*) 'yjn = ', yjn
write(*,*) 'yy0 = ', yy0
write(*,*) 'yy1 = ', yy1
write(*,*) 'yyn = ', yyn
end program test3

C:\gfortran\clf\bessel_test>gfortran test3.f90 -otest3
test3.f90:20.28:

real(int(xj0+kind(1.0)+1)) yj0
1
Error: Constant expression required at (1)
test3.f90:21.26:

real(int(xj1+kind(1.0))) yj1
1
Error: Constant expression required at (1)
test3.f90:22.26:

real(int(xjn+kind(1.0))) yjn
1
Error: Constant expression required at (1)
test3.f90:23.26:

real(int(xy0+kind(1.0))) yy0
1
Error: Constant expression required at (1)
test3.f90:24.26:

real(int(xy1+kind(1.0))) yy1
1
Error: Constant expression required at (1)
test3.f90:25.26:

real(int(xyn+kind(1.0))) yyn
1
Error: Constant expression required at (1)
test3.f90:27.6:

yj0 = 42
1
Error: Symbol 'yj0' at (1) has no IMPLICIT type
test3.f90:28.6:

yj1 = 42
1
Error: Symbol 'yj1' at (1) has no IMPLICIT type
test3.f90:29.6:

yjn = 42
1
Error: Symbol 'yjn' at (1) has no IMPLICIT type
test3.f90:30.6:

yy0 = 42
1
Error: Symbol 'yy0' at (1) has no IMPLICIT type
test3.f90:31.6:

yy1 = 42
1
Error: Symbol 'yy1' at (1) has no IMPLICIT type
test3.f90:32.6:

yyn = 42
1
Error: Symbol 'yyn' at (1) has no IMPLICIT type

And then there is the circumstance that:

C:\gfortran\clf\bessel_test>type test4.f90
module mykinds
implicit none
integer, parameter :: ep = selected_real_kind(18,4931)
integer, parameter :: qp_preferred = selected_real_kind(33,4931)
integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1-sign(1,qp_preferred))/2*ep
end module mykinds

program test4
use mykinds
implicit none
real(qp), parameter :: x = 3
integer, parameter :: n = 1

write(*,*) 'BESSEL_J0(x) = ', BESSEL_J0(x)
write(*,*) 'BESSEL_J1(x) = ', BESSEL_J1(x)
write(*,*) 'BESSEL_JN(n,x) = ', BESSEL_JN(n,x)
write(*,*) 'BESSEL_Y0(x) = ', BESSEL_Y0(x)
write(*,*) 'BESSEL_Y1(x) = ', BESSEL_Y1(x)
write(*,*) 'BESSEL_YN(n,x) = ', BESSEL_YN(n,x)
end program test4

C:\gfortran\clf\bessel_test>gfortran test4.f90 -otest4

C:\gfortran\clf\bessel_test>test4
BESSEL_J0(x) = -0.26005195490193345000
BESSEL_J1(x) = 0.33905895852593645000
BESSEL_JN(n,x) = 0.33905895852593645000
BESSEL_Y0(x) = 0.37685001001279040000
BESSEL_Y1(x) = 0.32467442479180000000
BESSEL_YN(n,x) = 0.32467442479180000000

....compiles, as does:

C:\gfortran\clf\bessel_test>type test5.f90
module mykinds
implicit none
integer, parameter :: ep = selected_real_kind(18,4931)
integer, parameter :: qp_preferred = selected_real_kind(33,4931)
integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1-sign(1,qp_preferred))/2*ep
end module mykinds

module funcs
use mykinds
implicit none
contains
subroutine sub(n,x)
integer n
real x

write(*,*) 'BESSEL_J0(x) = ', BESSEL_J0(x)
write(*,*) 'BESSEL_J1(x) = ', BESSEL_J1(x)
write(*,*) 'BESSEL_JN(n,x) = ', BESSEL_JN(n,x)
write(*,*) 'BESSEL_Y0(x) = ', BESSEL_Y0(x)
write(*,*) 'BESSEL_Y1(x) = ', BESSEL_Y1(x)
write(*,*) 'BESSEL_YN(n,x) = ', BESSEL_YN(n,x)
end subroutine sub
end module funcs

program test5
use mykinds
use funcs
implicit none
real, parameter :: x = 3
integer, parameter :: n = 1

call sub(n,x)
end program test5

C:\gfortran\clf\bessel_test>gfortran test5.f90 -otest5

C:\gfortran\clf\bessel_test>test5
BESSEL_J0(x) = -0.26005197
BESSEL_J1(x) = 0.33905897
BESSEL_JN(n,x) = 0.33905897
BESSEL_Y0(x) = 0.37685001
BESSEL_Y1(x) = 0.32467443
BESSEL_YN(n,x) = 0.32467443

But:

C:\gfortran\clf\bessel_test>type test6.f90
module mykinds
implicit none
integer, parameter :: ep = selected_real_kind(18,4931)
integer, parameter :: qp_preferred = selected_real_kind(33,4931)
integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
(1-sign(1,qp_preferred))/2*ep
end module mykinds

module funcs
use mykinds
implicit none
contains
subroutine sub(n,x)
integer n
real(qp) x

write(*,*) 'BESSEL_J0(x) = ', BESSEL_J0(x)
write(*,*) 'BESSEL_J1(x) = ', BESSEL_J1(x)
write(*,*) 'BESSEL_JN(n,x) = ', BESSEL_JN(n,x)
write(*,*) 'BESSEL_Y0(x) = ', BESSEL_Y0(x)
write(*,*) 'BESSEL_Y1(x) = ', BESSEL_Y1(x)
write(*,*) 'BESSEL_YN(n,x) = ', BESSEL_YN(n,x)
end subroutine sub
end module funcs

program test6
use mykinds
use funcs
implicit none
real(qp), parameter :: x = 3
integer, parameter :: n = 1

call sub(n,x)
end program test6

C:\gfortran\clf\bessel_test>gfortran test6.f90 -otest6
C:\DOCUME~1\ADMINI~1\LOCALS~1\Temp/ccUFaaaa.o:test6.f90:(.text+0x8a):
undefined
reference to `_j0l'
C:\DOCUME~1\ADMINI~1\LOCALS~1\Temp/ccUFaaaa.o:test6.f90:(.text+0x13b):
undefined
reference to `_j1l'
C:\DOCUME~1\ADMINI~1\LOCALS~1\Temp/ccUFaaaa.o:test6.f90:(.text+0x1fc):
undefined
reference to `_jnl'
C:\DOCUME~1\ADMINI~1\LOCALS~1\Temp/ccUFaaaa.o:test6.f90:(.text+0x2ad):
undefined
reference to `_y0l'
C:\DOCUME~1\ADMINI~1\LOCALS~1\Temp/ccUFaaaa.o:test6.f90:(.text+0x35e):
undefined
reference to `_y1l'
C:\DOCUME~1\ADMINI~1\LOCALS~1\Temp/ccUFaaaa.o:test6.f90:(.text+0x41f):
undefined
reference to `_ynl'
collect2: ld returned 1 exit status

....fails to link. I would have thought it simpler in some sense
to just make all this consistent on a given platform rather than
attempt to document all the exceptions. While we are on the
subject of Bessel functions, did you notice that gfortran's
initialization expression evaluation mechanism produced weird
low-precision approximations of Bessel functions in the examples
where it worked?

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end


From: FX on 6 May 2008 17:04
> Is there some way to 'unfix' PR323 so that Win64 works like gawd &
> Billy Kahan intended?

If you expect any kind of answer, you need to provide a clear,
self-contained question.

--
FX
From: James Van Buskirk on 5 May 2008 20:24
"James Van Buskirk" <not_valid(a)comcast.net> wrote in message
news:9tednQ43AIB-HILVnZ2dnUVZ_tSknZ2d(a)comcast.com...

[Snip examples of borderline syntax of Bessel functions]

There is another place that is even more fun: the standard itself is
broken regarding the SHAPE argument to RESHAPE. OK, I see that
N1723.pdf has fixed the error in N1601.pdf here.

C:\gfortran\clf\bessel_test>type test7.f90
module funcs
implicit none
private
public sub
interface sub
module procedure sub1, sub2, sub3, sub4, sub5, sub6, sub7
end interface sub
contains
subroutine sub1(x)
integer x(:)

write(*,'(a)') 'Welcome to subroutine sub1!'
write(*,'(a,7(i0:", "))') 'Shape(x) = ', shape(x)
end subroutine sub1

subroutine sub2(x)
integer x(:,:)

write(*,'(a)') 'Welcome to subroutine sub2!'
write(*,'(a,7(i0:", "))') 'Shape(x) = ', shape(x)
end subroutine sub2

subroutine sub3(x)
integer x(:,:,:)

write(*,'(a)') 'Welcome to subroutine sub3!'
write(*,'(a,7(i0:", "))') 'Shape(x) = ', shape(x)
end subroutine sub3

subroutine sub4(x)
integer x(:,:,:,:)

write(*,'(a)') 'Welcome to subroutine sub4!'
write(*,'(a,7(i0:", "))') 'Shape(x) = ', shape(x)
end subroutine sub4

subroutine sub5(x)
integer x(:,:,:,:,:)

write(*,'(a)') 'Welcome to subroutine sub5!'
write(*,'(a,7(i0:", "))') 'Shape(x) = ', shape(x)
end subroutine sub5

subroutine sub6(x)
integer x(:,:,:,:,:,:)

write(*,'(a)') 'Welcome to subroutine sub6!'
write(*,'(a,7(i0:", "))') 'Shape(x) = ', shape(x)
end subroutine sub6

subroutine sub7(x)
integer x(:,:,:,:,:,:,:)

write(*,'(a)') 'Welcome to subroutine sub7!'
write(*,'(a,7(i0:", "))') 'Shape(x) = ', shape(x)
end subroutine sub7
end module funcs

program test7
use funcs
implicit none
real(10), parameter :: x = 3
integer i

call sub(reshape([0],[(3,i=1,int(5+BESSEL_J0(x)))],pad=[0]))
end program test7

C:\gfortran\clf\bessel_test>gfortran test7.f90 -otest7
test7.f90:65.24:

call sub(reshape([0],[(3,i=1,int(5+BESSEL_J0(x)))],pad=[0]))
1
Error: 'shape' argument of 'reshape' intrinsic at (1) must be an array of
consta
nt size

So this is another limitation of the Bessel functions. I see that
gfortran's error message, in the first place incorrect for f08, also
is not in tune with the f08 update from an 'an array of constant
size' to the rather more complicated but accurate statement that
SIZE(x), where x is that SHAPE actual argument, must be an
initialization expression.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end


 |  Next  |  Last
Pages: 1 2
Prev: Treacherous Behavior of "gfortran" Compiler?
Next: Build project in CVF6.6c by using hotkey