Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Are f.p. manipulation functions only used in initialization?

62 views
Skip to first unread message

James Van Buskirk

unread,
Apr 26, 2008, 4:43:52 PM4/26/08
to
I was writing some code with gfortran and found that spacing only
worked in initialization expressions. Doesn't anyone else use the
floating point manipulation functions in ordinary expressions and
specification expressions? Testing with gfortran I found that
RRSPACING, SCALE, SET_EXPONENT, and, SPACING were broken like this,
while EXPONENT, FRACTION, and NEAREST seemed to work:

C:\gfortran\james\archpi>c:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran -v
Using built-in specs.
Target: x86_64-pc-mingw32
Configured with:
../../trunk/configure --prefix=/home/FX/irun64 --build=i586-pc-
mingw32 --target=x86_64-pc-mingw32 --with-gmp=/home/FX/local --enable-languages=
c,fortran --disable-werror --disable-nls --enable-threads
Thread model: win32
gcc version 4.4.0 20080421 (experimental) [trunk revision 134506] (GCC)

C:\gfortran\james\archpi>type spacing_bug.f90
program spacing_bug
implicit none
real, parameter :: x1 = 3.14159265358979
integer, parameter :: ulps1 = 2
real, parameter :: y1 = ulps1*spacing(x1) ! Only f03
real, parameter :: y2 = rrspacing(x1) ! Only f03
real, parameter :: y3 = scale(x1,ulps1) ! Only f03
real, parameter :: y4 = set_exponent(x1,ulps1) ! Only f03
real x
integer ulps

write(*,*) 'ulps1*spacing(x1) = ', y1
write(*,*) 'rrspacing(x1) = ', y2
write(*,*) 'scale(x1,ulps1) = ', y3
write(*,*) 'set_exponent(x1,ulps1) = ', y4
x = x1
ulps = ulps1
call sub(x,ulps)
end program spacing_bug

subroutine sub(x,ulps)
implicit none
real x
integer ulps
integer j1(int(ulps*spacing(x)))
integer j2(int(rrspacing(x))/100000)
integer j3(int(scale(x,ulps)))
integer j4(int(set_exponent(x,ulps)))

write(*,*) 'ulps*spacing(x) = ', ulps*spacing(x)
write(*,*) 'rrspacing(x) = ', rrspacing(x)
write(*,*) 'scale(x,ulps) = ', scale(x,ulps)
write(*,*) 'set_exponent(x,ulps) = ', set_exponent(x,ulps)
write(*,*) 'size(j1) = ', size(j1)
write(*,*) 'size(j2) = ', size(j2)
write(*,*) 'size(j3) = ', size(j3)
write(*,*) 'size(j4) = ', size(j4)
end subroutine sub

C:\gfortran\james\archpi>c:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran
spacin
g_bug.f90 -ospacing_bug

C:\gfortran\james\archpi>spacing_bug
ulps1*spacing(x1) = 4.76837158E-07
rrspacing(x1) = 13176795.
scale(x1,ulps1) = 12.566371
set_exponent(x1,ulps1) = 3.1415927
ulps*spacing(x) = 2.0000000
rrspacing(x) = 3.1415927
scale(x,ulps) = 3.1415927
set_exponent(x,ulps) = 0.78539819
size(j1) = 2
size(j2) = 0
size(j3) = 3
size(j4) = 0

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


Steven G. Kargl

unread,
Apr 26, 2008, 5:52:30 PM4/26/08
to
In article <6_6dnatLa7meCY7V...@comcast.com>,

"James Van Buskirk" <not_...@comcast.net> writes:
> I was writing some code with gfortran and found that spacing only
> worked in initialization expressions. Doesn't anyone else use the
> floating point manipulation functions in ordinary expressions and
> specification expressions? Testing with gfortran I found that
> RRSPACING, SCALE, SET_EXPONENT, and, SPACING were broken like this,
> while EXPONENT, FRACTION, and NEAREST seemed to work:
>

(snip)



> C:\gfortran\james\archpi>c:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran
> spacin
> g_bug.f90 -ospacing_bug
>
> C:\gfortran\james\archpi>spacing_bug
> ulps1*spacing(x1) = 4.76837158E-07
> rrspacing(x1) = 13176795.
> scale(x1,ulps1) = 12.566371
> set_exponent(x1,ulps1) = 3.1415927
> ulps*spacing(x) = 2.0000000
> rrspacing(x) = 3.1415927
> scale(x,ulps) = 3.1415927
> set_exponent(x,ulps) = 0.78539819
> size(j1) = 2
> size(j2) = 0
> size(j3) = 3
> size(j4) = 0
>

mobile:kargl[203] ./z
ulps1*spacing(x1) = 4.7683716E-07
rrspacing(x1) = 1.3176795E+07
scale(x1,ulps1) = 12.56637
set_exponent(x1,ulps1) = 3.141593
ulps*spacing(x) = 4.7683716E-07
rrspacing(x) = 1.3176795E+07
scale(x,ulps) = 12.56637
set_exponent(x,ulps) = 3.141593
size(j1) = 0
size(j2) = 131
size(j3) = 12
size(j4) = 3

Appears to work for me. Perhaps, it's an OS problem.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

James Van Buskirk

unread,
Apr 26, 2008, 6:24:59 PM4/26/08
to
"Steven G. Kargl" <ka...@troutmask.apl.washington.edu> wrote in message
news:fv086u$roe$1...@gnus01.u.washington.edu...

> Appears to work for me. Perhaps, it's an OS problem.

Thank you for testing. Maybe an OS problem or a regression. What
version are you using? I'm on 20080421.

I ran into another problem that may be with the documentation. In:

http://gcc.gnu.org/onlinedocs/gcc/i386-and-x86_002d64-Options.html#i386-and-x86_002d64-Options

It says:

"-mpc32
-mpc64
-mpc80
Set 80387 floating-point precision to 32, 64 or 80 bits. When
-mpc32 is specified, the significands of results of floating-point
operations are rounded to 24 bits (single precision); -mpc64 rounds
the the significands of results of floating-point operations to 53
bits (double precision) and -mpc80 rounds the significands of results
of floating-point operations to 64 bits (extended double precision),
which is the default. When this option is used, floating-point
operations in higher precisions are not available to the programmer
without setting the FPU control word explicitly.

Setting the rounding of floating-point operations to less than the
default 80 bits can speed some programs by 2% or more. Note that some
mathematical libraries assume that extended precision (80 bit)
floating-point operations are enabled by default; routines in such
libraries could suffer significant loss of accuracy, typically through
so-called "catastrophic cancellation", when this option is used to set
the precision to less than extended precision."

My reading of this is that precision control should be set to
b'11' = Extended Precision, but I must be reading the wrong part
of the manual because I get:

C:\gfortran\james\archpi>c:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran -v
Using built-in specs.
Target: x86_64-pc-mingw32
Configured with:
../../trunk/configure --prefix=/home/FX/irun64 --build=i586-pc-
mingw32 --target=x86_64-pc-mingw32 --with-gmp=/home/FX/local --enable-languages=
c,fortran --disable-werror --disable-nls --enable-threads
Thread model: win32
gcc version 4.4.0 20080421 (experimental) [trunk revision 134506] (GCC)

C:\gfortran\james\archpi>type test_status.f90
program test_status
implicit none
interface
function status()
integer(selected_int_kind(4)) status
end function status
subroutine control(cw)
integer(selected_int_kind(4)) cw
end subroutine control
end interface
real(10) x
real(10) y
integer(selected_int_kind(2)) array(10)
integer(selected_int_kind(4)) old_control

x = 1
x = 4*atan(x)
y = x-nearest(x,-1.0_10)
y = x+y
array = transfer(x,array)
write(*,'(10z2.2)') array(size(array):1:-1)
array = transfer(y,array)
write(*,'(10z2.2)') array(size(array):1:-1)
old_control = status()
write(*,'(z4.4)') old_control
call control(ior(old_control,int(z'0300',kind(old_control))))
x = 1
x = 4*atan(x)
y = x-nearest(x,-1.0_10)
y = x+y
array = transfer(x,array)
write(*,'(10z2.2)') array(size(array):1:-1)
array = transfer(y,array)
write(*,'(10z2.2)') array(size(array):1:-1)
write(*,'(z4.4)') status()
call control(old_control)
end program test_status

C:\gfortran\james\archpi>type status.s
.text
.globl _status_
.def _status_; .scl 2; .type 32; .endef
_status_:
fstcw 8(%rsp)
movzwl 8(%rsp), %eax
ret
.globl _control_
.def _control_; .scl 2; .type 32; .endef
_control_:
fldcw (%rcx)
ret

C:\gfortran\james\archpi>c:\gfortran\win64\bin\x86_64-pc-mingw32-gfortran -mpc80
test_status.f90 status.s -otest_status

C:\gfortran\james\archpi>test_status
4000C90FDAA22168C000
4000C90FDAA22168C000
027F
4000C90FDAA22168C235
4000C90FDAA22168C236
037F

So you can see that gfortran isn't setting precision control to b'11'
even with the -mpc80 switch in effect so that real(10) arithmetic
isn't working until I set the x87 control word by hand. Probably I'm
just missing something in the documentation somewhere but it can sure
be hard to find it from a standing start.

Steven G. Kargl

unread,
Apr 26, 2008, 6:54:20 PM4/26/08
to
In article <QdCdnaPuP40qNo7V...@comcast.com>,

"James Van Buskirk" <not_...@comcast.net> writes:
> "Steven G. Kargl" <ka...@troutmask.apl.washington.edu> wrote in message
> news:fv086u$roe$1...@gnus01.u.washington.edu...
>
>> Appears to work for me. Perhaps, it's an OS problem.
>
> Thank you for testing. Maybe an OS problem or a regression. What
> version are you using? I'm on 20080421.

On i386-*-freebsd, I have
mobile:kargl[203] gfc --version
GNU Fortran (GCC) 4.2.3 20080130 (prerelease)

On amd64-*-freebsd, I have

troutmask:kargl[201] gfc --version
GNU Fortran (GCC) 4.2.3 20071024 (prerelease)

and

troutmask:kargl[202] ../sgk/work/4x/bin/gfortran --version
GNU Fortran (GCC) 4.4.0 20080415 (experimental)


All versions give the same results.


>
> I ran into another problem that may be with the documentation. In:
>
> http://gcc.gnu.org/onlinedocs/gcc/i386-and-x86_002d64-Options.html#i386-and-x86_002d64-Options
>
> It says:
>
> "-mpc32
> -mpc64
> -mpc80

First time I've seen these flags. These appear to be new in 4.4
(at least 4.2.x does not recognize the options).

> Set 80387 floating-point precision to 32, 64 or 80 bits. When
> -mpc32 is specified, the significands of results of floating-point
> operations are rounded to 24 bits (single precision); -mpc64 rounds
> the the significands of results of floating-point operations to 53
> bits (double precision) and -mpc80 rounds the significands of results
> of floating-point operations to 64 bits (extended double precision),
> which is the default. When this option is used, floating-point
> operations in higher precisions are not available to the programmer
> without setting the FPU control word explicitly.
>
> Setting the rounding of floating-point operations to less than the
> default 80 bits can speed some programs by 2% or more. Note that some
> mathematical libraries assume that extended precision (80 bit)
> floating-point operations are enabled by default; routines in such
> libraries could suffer significant loss of accuracy, typically through
> so-called "catastrophic cancellation", when this option is used to set
> the precision to less than extended precision."

I doubt gfortran honors these flags, particular with anything that
gfortran can constant fold or the results of a numeric inquiry function.

(program snipped)

> So you can see that gfortran isn't setting precision control to b'11'
> even with the -mpc80 switch in effect so that real(10) arithmetic
> isn't working until I set the x87 control word by hand. Probably I'm
> just missing something in the documentation somewhere but it can sure
> be hard to find it from a standing start.

The ChangeLog entry for these option suggest that they may be linux
(Unix-like OS?) specific.

2007-04-03 Uros Bizjak <ubi...@gmail.com>

* config.gcc (i[34567]86-*-linux*): Add i386/t-crtpc to tm-file.
(x86_64-*-linux*): Ditto.
* config/i386/i386.opt (mpc): New option.
* config/i386/i386.c (overrride_options): Handle
ix87_precision_string.
* config/i386/crtprec.c: New file.
* config/i386/t-crtpc: Ditto.
* config/i386/linux.h (ENDFILE_SPEC): Add handling of -mpc32, -mpc64
and -mpc80 options.
* config/i386/linux64.h (ENDFILE_SPEC): Ditto.
* config/i386/t-linux64 (EXTRA_MULTILIB_PARTS): Add
crtprec32.o, crtprec64.o and crtprec80.o.

* doc/invoke.texi (Machine Dependent Options): Add -mpc32, -mpc64
and -mpc80 options.
(i386 and x86-64 Options): Document -mpc32, -mpc64 and -mpc80 options.


--
Steve
http://troutmask.apl.washington.edu/~kargl/

Anony

unread,
Apr 26, 2008, 7:25:57 PM4/26/08
to

"James Van Buskirk" <not_...@comcast.net> wrote in message
news:6_6dnatLa7meCY7V...@comcast.com...

On Windows, my result (4.4.0 20080425) is as:

C:\TEMP\fortran>a


ulps1*spacing(x1) = 4.76837158E-07
rrspacing(x1) = 13176795.
scale(x1,ulps1) = 12.566371
set_exponent(x1,ulps1) = 3.1415927

ulps*spacing(x) = 4.76837158E-07
rrspacing(x) = 13176795.
scale(x,ulps) = 12.566371
set_exponent(x,ulps) = 3.1415927

James Van Buskirk

unread,
Apr 26, 2008, 10:10:58 PM4/26/08
to
"Anony" <invali...@equation.com> wrote in message
news:9SOQj.2034$e26.671@trnddc02...

> On Windows, my result (4.4.0 20080425) is as:

> C:\TEMP\fortran>a
> ulps1*spacing(x1) = 4.76837158E-07
> rrspacing(x1) = 13176795.
> scale(x1,ulps1) = 12.566371
> set_exponent(x1,ulps1) = 3.1415927
> ulps*spacing(x) = 4.76837158E-07
> rrspacing(x) = 13176795.
> scale(x,ulps) = 12.566371
> set_exponent(x,ulps) = 3.1415927
> size(j1) = 0
> size(j2) = 131
> size(j3) = 12
> size(j4) = 3

C:\gfortran\james\archpi>c:\gcc_equation\bin\x86_64-pc-mingw32-gfortran -v
Built by Equation Solution (http://www.Equation.com).


Using built-in specs.
Target: x86_64-pc-mingw32
Configured with:

../gcc-4.4-20080425-mingw/configure --host=x86_64-pc-mingw32 --
build=x86_64-unknown-linux-gnu --target=x86_64-pc-mingw32 --prefix=/home/gfortra
n/gcc-home/binary/mingw32/native/x86_64/gcc/4.4-20080425 --with-gmp=/home/gfortr
an/gcc-home/binary/mingw32/native/x86_64/gmp --with-mpfr=/home/gfortran/gcc-home
/binary/mingw32/native/x86_64/mpfr --with-sysroot=/home/gfortran/gcc-home/binary
/mingw32/cross/x86_64/gcc/4.4-20080425 --with-gcc --with-gnu-ld --with-gnu-as
--
disable-shared --disable-nls --disable-tls --enable-languages=c,fortran --enable
-threads=win32 --enable-libgomp --disable-win32-registry
Thread model: win32
gcc version 4.4.0 20080425 (experimental) (GCC)

C:\gfortran\james\archpi>c:\gcc_equation\bin\x86_64-pc-mingw32-gfortran
spacing_
bug.f90 -ospacing_bug

C:\gfortran\james\archpi>spacing_bug
ulps1*spacing(x1) = 4.76837158E-07
rrspacing(x1) = 13176795.
scale(x1,ulps1) = 12.566371
set_exponent(x1,ulps1) = 3.1415927
ulps*spacing(x) = 2.0000000
rrspacing(x) = 3.1415927
scale(x,ulps) = 3.1415927
set_exponent(x,ulps) = 0.78539819
size(j1) = 2
size(j2) = 0
size(j3) = 3
size(j4) = 0

So I presume your test was carried out on 32-bit Windows. As you
can see, same date build fails on 64-bit Windows. I recall that in
the earliest versions of gfortran for 64-bit Windows, all the single
precision intrinsics were missing. This may be a continuation of
that problem, but now for a subset of intrinsics rather than a
subset of KINDs. Steven Kargl's intuition about platform specificity
seems to have been accurate.

James Van Buskirk

unread,
Apr 26, 2008, 10:55:30 PM4/26/08
to
"James Van Buskirk" <not_...@comcast.net> wrote in message
news:QdCdnaPuP40qNo7V...@comcast.com...

> I ran into another problem that may be with the documentation.

Here's another item from the docs:

http://gcc.gnu.org/onlinedocs/gfortran/Code-Gen-Options.html#Code-Gen-Options

"-frepack-arrays
In some circumstances GNU Fortran may pass assumed shape array
sections via a descriptor describing a noncontiguous area of memory.
This option adds code to the function prologue to repack the data
into a contiguous block at runtime.

This should result in faster accesses to the array. However it can
introduce significant overhead to the function call, especially when
the passed data is noncontiguous."

It is not mentioned here that this is one of the options that makes
the program behave contrary to the standard:

C:\gfortran\james\archpi>type repack_test.f90
module test1
use ISO_C_BINDING, only: C_PTR, C_LOC
implicit none
contains
function point(x)
real, intent(in), target :: x(:)
type(C_PTR) point
real, pointer :: p

p => x(2)
point = C_LOC(p)
end function point
end module test1

program test2
use test1
use ISO_C_BINDING
implicit none
real, target :: x(7)
type(C_PTR) cp1, cp2
integer(C_INTPTR_T) ip

x = 42
cp1 = C_LOC(x(3))
cp2 = point(x(::2))
write(*,'(2(z16.16:1x))') transfer(cp1,ip), transfer(cp2,ip)
write(*,'(L1)') C_ASSOCIATED(cp1,cp2)
end program test2


C:\gfortran\james\archpi>c:\gcc_equation\bin\x86_64-pc-mingw32-gfortran -std=f20
03 repack_test.f90 -orepack_test

C:\gfortran\james\archpi>repack_test
000000000022FE48 000000000022FE48
T

C:\gfortran\james\archpi>c:\gcc_equation\bin\x86_64-pc-mingw32-gfortran -frepack
-arrays -std=f2003 repack_test.f90 -orepack_test

C:\gfortran\james\archpi>repack_test
000000000022FE48 00000000004A5A04
F

It would be nice to have a warning either in the documentation for
-frepack-arrays or at compile time when -std=f2003 is in force
that the result is generated code that is inconsistent with the
standard.

Tobias Burnus

unread,
Apr 27, 2008, 4:25:20 AM4/27/08
to
On Apr 27, 4:55 am, "James Van Buskirk" <not_va...@comcast.net> wrote:
> "James Van Buskirk" <not_va...@comcast.net> wrote in messagenews:QdCdnaPuP40qNo7V...@comcast.com...
> "-frepack-arrays

>
> It is not mentioned here that this is one of the options that makes
> the program behave contrary to the standard:

Confirmed. The problem is that in the function "point" a temporary
array is created to contain the repacked data of the dummy argument
"x". Thus, the pointer points to the temporary array and not to the
dummy array.

Solution: Do not repack dummy arguments with the TARGET attribute. I
filled now http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36059

Thanks for the bug report. (Maybe I should re-read at some point the
gfortran documentation, there are options I have never heard of.)

Tobias

James Van Buskirk

unread,
Apr 27, 2008, 5:32:59 AM4/27/08
to
"Tobias Burnus" <bur...@net-b.de> wrote in message
news:24b0ceab-d3e9-4b68...@c65g2000hsa.googlegroups.com...

> Confirmed. The problem is that in the function "point" a temporary
> array is created to contain the repacked data of the dummy argument
> "x". Thus, the pointer points to the temporary array and not to the
> dummy array.

> Solution: Do not repack dummy arguments with the TARGET attribute. I
> filled now http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36059

I think that both the actual argument and the dummy argument must
have the target attribute for a problem to arise. If the actual
argument doesn't have the target attribute any pointers that think
they are associated with it have in fact undefined association
status, and same for the dummy. Since the interface must be explicit
if the dummy argument has the target attribute (N1601.pdf, section
12.3.1.1 (2)(a)) the compiler knows when it is deciding whether to
repack the array (i.e. when compiling the caller) whether or not
both actual and dummy arguments have the target attribute.

Another possibility would be to tolerate the nonstandard behavior
but call attention to the possibility in the documentation. There
are examples of compilers that do this. Actually gfortran is more
benign in this regard because the nonstandard switch is not the
default.

> Thanks for the bug report. (Maybe I should re-read at some point the
> gfortran documentation, there are options I have never heard of.)

Thank you for your attention to this notification. Keep up the
good (volunteer) work.

Anony

unread,
Apr 27, 2008, 8:34:46 AM4/27/08
to

"James Van Buskirk" <not_...@comcast.net> wrote in message
news:eZudnecxz7swfY7V...@comcast.com...

>
>
> So I presume your test was carried out on 32-bit Windows. As you
> can see, same date build fails on 64-bit Windows. I recall that in
> the earliest versions of gfortran for 64-bit Windows, all the single
> precision intrinsics were missing. This may be a continuation of
> that problem, but now for a subset of intrinsics rather than a
> subset of KINDs. Steven Kargl's intuition about platform specificity
> seems to have been accurate.
>

Yes. I was on 32-bit windows when I copied and pasted your program. Sorry,
not to see you are on win64.

It is nice to know there are people working on win64. Win64 seems not ready
yet. There are some problems, even affecting binutuls. RANLIB and AR also
have problem (It cannot be reproduced everytime.). There was one user who
complined "hang forever" without providing a code example (I never
eperienced the problem, and I have no idea).

Charles Coldwell

unread,
Apr 27, 2008, 9:00:43 AM4/27/08
to
"James Van Buskirk" <not_...@comcast.net> writes:

> Setting the rounding of floating-point operations to less than the
> default 80 bits can speed some programs by 2% or more. Note that some
> mathematical libraries assume that extended precision (80 bit)
> floating-point operations are enabled by default; routines in such
> libraries could suffer significant loss of accuracy, typically through
> so-called "catastrophic cancellation", when this option is used to set
> the precision to less than extended precision."

Notably, on GNU/Linux (i.e. GNU glibc on Linux/x86)
/usr/include/fpu_control.h defines

#define _FPU_GETCW(cw) __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw))
#define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw))

and

/* precision control */
#define _FPU_EXTENDED 0x300 /*libm requires double extended precision. */
#define _FPU_DOUBLE 0x200
#define _FPU_SINGLE 0x0

The comment next to _FPU_EXTENDED should be written in a bold,
blinking, red font. You will screw up all of the libm routines (such
as cos, sin, exp, log) if you change the FPU precision. Since Fortran
intrinsics are usually (always?) implemented in terms of libm
routines, you will threfore screw up the corresponding intrinsics.
About the only reason to change the FPU precision is if you are
writing your own libm, or software floating point (such as some of the
arbitrary precision work done by David H. Bailey and Yozo Hida), or
something similar.

> C:\gfortran\james\archpi>test_status
> 4000C90FDAA22168C000
> 4000C90FDAA22168C000
> 027F
> 4000C90FDAA22168C235
> 4000C90FDAA22168C236
> 037F
>
> So you can see that gfortran isn't setting precision control to b'11'
> even with the -mpc80 switch in effect so that real(10) arithmetic
> isn't working until I set the x87 control word by hand.

I think what is happening here is an OS dependency. It must be that
gfortran/Win32 is linking to a Win32 math library that depends on
having the FPU in double precision, not double extended precision.
You might check to see what happens to some of the trigonometric
intrinsics after you diddle the FPU mode.

So perhaps the double extended precision can't work on Win32 because
either the FPU is in double extended precision and the math library is
borken, or the FPU is in single precision and real*10 is borken.

Chip

--
Charles M. "Chip" Coldwell
"Turn on, log in, tune out"
GPG Key ID: 852E052F
GPG Key Fingerprint: 77E5 2B51 4907 F08A 7E92 DE80 AFA9 9A8F 852E 052F

FX

unread,
Apr 27, 2008, 2:37:33 PM4/27/08
to
Hi Anony,

> It is nice to know there are people working on win64. Win64 seems not
> ready yet. There are some problems, even affecting binutuls. RANLIB and
> AR also have problem

Whatever problems you are experiencing, please report them to the
respective projects; this is how you can, as a user, contribute to the
progress of these projects.

--
FX

FX

unread,
Apr 27, 2008, 2:39:17 PM4/27/08
to
Hi James,

> Testing with gfortran I found that
> RRSPACING, SCALE, SET_EXPONENT, and, SPACING were broken like this,
> while EXPONENT, FRACTION, and NEAREST seemed to work

Thanks for your example, I will try to reduce it and understand what's
happening; if it happens to not be a Fortran front-end issue, I'll alert
our great Win64 developers to it.

--
FX

James Van Buskirk

unread,
Apr 27, 2008, 3:53:14 PM4/27/08
to
"Charles Coldwell" <cold...@gmail.com> wrote in message
news:m2wsmj1...@gmail.com...

> "James Van Buskirk" <not_...@comcast.net> writes:

> Notably, on GNU/Linux (i.e. GNU glibc on Linux/x86)
> /usr/include/fpu_control.h defines

> #define _FPU_GETCW(cw) __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw))
> #define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw))

You can see why I gave my floating control word code in GAS: at
least the reader had some kind of chance to puzzle out what I
was doing. Easier IMO to look at Intels documentation such as

http://download.intel.com/design/processor/manuals/253666.pdf

and write to the processor in assembly.

> and

> /* precision control */
> #define _FPU_EXTENDED 0x300 /*libm requires double extended precision. */
> #define _FPU_DOUBLE 0x200
> #define _FPU_SINGLE 0x0

Documented in

http://download.intel.com/design/processor/manuals/253665.pdf

> The comment next to _FPU_EXTENDED should be written in a bold,
> blinking, red font. You will screw up all of the libm routines (such
> as cos, sin, exp, log) if you change the FPU precision. Since Fortran
> intrinsics are usually (always?) implemented in terms of libm
> routines, you will threfore screw up the corresponding intrinsics.

Um, well, not really. x86_64 is a different animal from x86, and
at least SSE2 is guaranteed on the former platform, so in principle
one could get by without any use of x87. If you look at the gnu docs
I quoted,

http://gcc.gnu.org/onlinedocs/gcc/i386-and-x86_002d64-Options.html#i386-and-x86_002d64-Options

it says that -mpc80 is the default anyway, so at least the documentation
is in error here. In fact this whole page is garbage for Win64; I
don't know about other x86_64 platforms. For example the -mregparm=num
paragraph states that by default no parameters are passed in registers
whereas in fact four are in Win64 and lots more are in Linux x86_64.
See http://www.agner.org/optimize/calling_conventions.pdf for summary
information on default calling conventions on these platforms. This
page on gcc documentation needs a ?Redo from start.

> About the only reason to change the FPU precision is if you are
> writing your own libm, or software floating point (such as some of the
> arbitrary precision work done by David H. Bailey and Yozo Hida), or
> something similar.

You know that precision control only applies to addition, subtraction,
multiplication, division, and square root, right? Division and square
root can be accelerated by changing precision control, addition,
subtraction, and multiplication are pipelined operations that proceed
at the same pace (as measured by latency and throughput) regardless
of precision control. ISTR that LF95 messed with precision control
before doing a lot of square roots at single precision. It would
have been faster to do pipelined Goldschmidt iterations, but that
would have been too much fun :)

>> C:\gfortran\james\archpi>test_status
>> 4000C90FDAA22168C000
>> 4000C90FDAA22168C000
>> 027F
>> 4000C90FDAA22168C235
>> 4000C90FDAA22168C236
>> 037F

>> So you can see that gfortran isn't setting precision control to b'11'
>> even with the -mpc80 switch in effect so that real(10) arithmetic
>> isn't working until I set the x87 control word by hand.

> I think what is happening here is an OS dependency. It must be that
> gfortran/Win32 is linking to a Win32 math library that depends on
> having the FPU in double precision, not double extended precision.
> You might check to see what happens to some of the trigonometric
> intrinsics after you diddle the FPU mode.

The x87 transcendental functions are documented to run at extended
precision independent of precision control so unless better argument
reduction than the hardware does is requied, precision control doesn't
matter for them. If fancy argument reduction is going to happen, the
overhead of changing the precision control is a drop in the bucket.

> So perhaps the double extended precision can't work on Win32 because
> either the FPU is in double extended precision and the math library is
> borken, or the FPU is in single precision and real*10 is borken.

I think you mean double precision, not single precision above.
Real*8 would be broken if the FPU were in single precision mode
and the software actually used the FPU for double precision
operations, which for the most part is doesn't do in x86_84
mode. I'm not sure what's going on here, but as you have noted
the value of the precision control bits is contrary to documentation
both of us have been able to find and this is what is preventing
real*10 operations (which I can't find much documentation about
in the gfortran manual) from working in any kind of expected
fashion. I don't know if there is something down there in gcc
which requires setting pc = b'10' so that it is intentionally set
up that way (like MicroSoft libraries) or if it's just an accident
that can be fixed by setting pc = b'11' as I did in my example.

Steven G. Kargl

unread,
Apr 27, 2008, 4:30:26 PM4/27/08
to
In article <fv2h8l$2p4t$3...@nef.ens.fr>,

It's an OS issue. Look at -fdump-tree-original.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

FX

unread,
Apr 27, 2008, 6:12:09 PM4/27/08
to
> It's an OS issue. Look at -fdump-tree-original.

That's what I intend to do as soon as I get some time. There are two
things I suspect can happen, one of them is a GCC/gfortran issue (ie the
builtins themselves being wrongly defined, for example long vs. long
long) or an OS library issue. Wait and see!

--
FX

Charles Coldwell

unread,
Apr 28, 2008, 6:46:00 AM4/28/08
to
"James Van Buskirk" <not_...@comcast.net> writes:

>
> C:\gfortran\james\archpi>spacing_bug
> ulps1*spacing(x1) = 4.76837158E-07
> rrspacing(x1) = 13176795.
> scale(x1,ulps1) = 12.566371
> set_exponent(x1,ulps1) = 3.1415927
> ulps*spacing(x) = 2.0000000
> rrspacing(x) = 3.1415927
> scale(x,ulps) = 3.1415927
> set_exponent(x,ulps) = 0.78539819
> size(j1) = 2
> size(j2) = 0
> size(j3) = 3
> size(j4) = 0

$ gfortran --version
GNU Fortran (GCC) 4.1.2 20070925 (Red Hat 4.1.2-33)
Copyright (C) 2007 Free Software Foundation, Inc.

GNU Fortran comes with NO WARRANTY, to the extent permitted by law.
You may redistribute copies of GNU Fortran
under the terms of the GNU General Public License.
For more information about these matters, see the file named COPYING

$ uname -m
x86_64

$ gfortran spacing_bug.f90 && ./a.out


ulps1*spacing(x1) = 4.7683716E-07
rrspacing(x1) = 1.3176795E+07
scale(x1,ulps1) = 12.56637
set_exponent(x1,ulps1) = 3.141593
ulps*spacing(x) = 4.7683716E-07
rrspacing(x) = 1.3176795E+07
scale(x,ulps) = 12.56637
set_exponent(x,ulps) = 3.141593

size(j1) = 0
size(j2) = 131
size(j3) = 12
size(j4) = 3

So perhaps a regression, or perhaps it never worked on Win64?

Greg Lindahl

unread,
Apr 28, 2008, 8:25:49 PM4/28/08
to
In article <m2wsmj1...@gmail.com>,
Charles Coldwell <cold...@gmail.com> wrote:

>The comment next to _FPU_EXTENDED should be written in a bold,
>blinking, red font. You will screw up all of the libm routines (such
>as cos, sin, exp, log) if you change the FPU precision.

Since many fortran compilers allow you to manipulate the x87
precision, I don't think this is the case -- one commonly used
commercial compiler sets it to 64 bits *by default*.

If you want exact IEEE answers, this is only one of many things you
need to worry about.

Also, the x87 is used a lot less on x86-64 than x86.

I've also seen comments that glibc's libm only works in 1 rounding
mode. That is changed a lot less frequently than the x87 precision.

-- greg

James Van Buskirk

unread,
Apr 30, 2008, 1:17:01 PM4/30/08
to
"FX" <cou...@alussinan.org> wrote in message
news:fv2h8l$2p4t$3...@nef.ens.fr...

> Hi James,

I further tried all functions that can accept a single REAL
variable as input. I was surprised to see in the docs that BESSEL_J0,
BESSEL_J1, BESSEL_Y0, and BESSEL_Y1, aren't elemental. Is this what
the rule is for the new f08 functions? These Bessel functions failed
to link, requiring _j0l, _j1l, _y0l, and _y1l respectively.

The manual says DFLOAT needs an INTEGER actual argument but in fact
it will takes REAL*4, REAL*8, or REAL*10. LOG_GAMMA was out of
alphabetic order in the manual. Is it true the gfortran's
implementation of NINT doesn't take a KIND optional argument?
Amusing that there is an INT2 and INT8 but no INT1 nor INT4. I
suppose you don't want to add extensions like that if you don't
really have to. Similar with SNGL accepting REAL*8 but not
REAL*10 input.

I didn't see what the default KIND for the result of CEILING was
in the manual.

ERF, ERFC, GAMMA, and LOG_GAMMA returned NaN for REAL*10 inputs in
my test. EXPONENT returned garbage for REAL*10 and FRACTION
return its input for REAL*10. SPACING and RRSPACING were hosed in
my test so that I didn't capture any output.

Test is at http://home.comcast.net/~kmbtib/Fortran_stuff/funr1.ZIP
Oops, I forgot to include a README: you can run the tests by
deleting funr1_gfc.out and funr1_gfc.err and then running
funr1_gfc.bat. If not on Windows, you may have to edit those
*.bat files to make it all go. :)

Dan Nagle

unread,
Apr 30, 2008, 6:18:33 PM4/30/08
to
Hello,

On 2008-04-30 13:17:01 -0400, "James Van Buskirk" <not_...@comcast.net> said:

> I further tried all functions that can accept a single REAL
> variable as input. I was surprised to see in the docs that BESSEL_J0,
> BESSEL_J1, BESSEL_Y0, and BESSEL_Y1, aren't elemental.

The j0, j1, y0, y1 functions are elemental.
The confusion come from the jn and yn functions.

There are two versions of these. One is elemental.
That one is the straight y( :) = jn( n, x( :)) version.

The transformational version is j/yn( n1, n2, x) version.
This returns all j/y's between order n1 and order n2 inclusive.
This version was added because many libraries compute
the Bessel functions via recurrence relations, and you
may as well get a bunch at once. The idea is that if
you have an array of coefficients, coef( 1: m), then
you can use the jn or yn with n2 - n1 + 1 = m, and pass
the result and your coef array to a polynomial evaluation
routine with high precision and efficiency.

HTH

--
Cheers!

Dan Nagle

Dan Nagle

unread,
Apr 30, 2008, 6:20:37 PM4/30/08
to
Hello,

On 2008-04-30 13:17:01 -0400, "James Van Buskirk" <not_...@comcast.net> said:

> LOG_GAMMA was out of
> alphabetic order in the manual.

In the standard, IIRC, the intrinsics are listed in ASCII order.

YMMV

--
Cheers!

Dan Nagle

FX

unread,
Apr 30, 2008, 6:54:21 PM4/30/08
to
> The transformational version is j/yn( n1, n2, x) version. This returns
> all j/y's between order n1 and order n2 inclusive. This version was
> added because many libraries compute the Bessel functions via
> recurrence relations, and you may as well get a bunch at once.

I had been wondering why you added these, so I get the idea. Frankly, I
don't think it's that good an idea: this is a very specialized topic, so
people who require this will probably have their own libraries anyway,
and I really don't like that idea of adding special case after special
case. History already gave us a long and complicated list of intrinsics,
adding yet another one that has two different forms is, in my opinion, a
pain.

--
FX

Dan Nagle

unread,
Apr 30, 2008, 8:14:51 PM4/30/08
to
Hello,

On 2008-04-30 18:54:21 -0400, "FX" <cou...@alussinan.org> said:

> I had been wondering why you added these, so I get the idea. Frankly, I
> don't think it's that good an idea: this is a very specialized topic, so
> people who require this will probably have their own libraries anyway,
> and I really don't like that idea of adding special case after special
> case. History already gave us a long and complicated list of intrinsics,
> adding yet another one that has two different forms is, in my opinion, a
> pain.

Well, these functions are in standard C, or in most C implementations,
anyway. So why not have them in Fortran?

--
Cheers!

Dan Nagle

Greg Lindahl

unread,
Apr 30, 2008, 8:53:08 PM4/30/08
to
In article <fvatat$16j$1...@nef.ens.fr>, FX <cou...@alussinan.org> wrote:

> History already gave us a long and complicated list of intrinsics,
> adding yet another one that has two different forms is, in my opinion, a
> pain.

He hid them in a module, right? That's the right way to do it in the
modern era, not creating yet more namespace pollution.

-- greg


James Van Buskirk

unread,
Apr 30, 2008, 11:37:03 PM4/30/08
to
"Dan Nagle" <dann...@verizon.net> wrote in message
news:2008043018183343658-dannagle@verizonnet...

> On 2008-04-30 13:17:01 -0400, "James Van Buskirk" <not_...@comcast.net>
> said:

>> I further tried all functions that can accept a single REAL
>> variable as input. I was surprised to see in the docs that BESSEL_J0,
>> BESSEL_J1, BESSEL_Y0, and BESSEL_Y1, aren't elemental.

> The j0, j1, y0, y1 functions are elemental.
> The confusion come from the jn and yn functions.

Actually I was looking at the gfortran documentation:
http://gcc.gnu.org/onlinedocs/gfortran/BESSEL_005fJ0.html#BESSEL_005fJ0
and you are referring to N1723.pdf. You can see there that gfortran
says for BESSEL_J0:

"Arguments:
X The type shall be REAL(*) and it shall be scalar."

and N1723.pdf says:

"Case (i): The result of BESSEL_JN(N,X) is scalar."

But that clearly isn't the right thing to say because
BESSEL_JN(N,X) is elemental, so maybe the verbiage should
be borrowed from ATAN2 which doesn't say anything like that,
so it's clear that ATAN2([1.0,2.0,3.0],1.0) and
ATAN2(1.0,[1.0,2.0,3.0]) have shape [3].

> There are two versions of these. One is elemental.
> That one is the straight y( :) = jn( n, x( :)) version.

> The transformational version is j/yn( n1, n2, x) version.
> This returns all j/y's between order n1 and order n2 inclusive.
> This version was added because many libraries compute
> the Bessel functions via recurrence relations, and you
> may as well get a bunch at once. The idea is that if
> you have an array of coefficients, coef( 1: m), then
> you can use the jn or yn with n2 - n1 + 1 = m, and pass
> the result and your coef array to a polynomial evaluation
> routine with high precision and efficiency.

Now, does transformational mean that BESSEL_JN(N1,N2,X) acts like
SPREAD if X is an array, or does it mean that N1, N2, and X must
all really be scalars in this form of the function? Here is an
example of the kind of confusion caused by the verbiage in N1723.pdf,
section 13.7.22:

C:\gfortran\clf\bessel_test>type bessel_test.f90
program bessel_test
implicit none
integer N(3)
! real(10) x(3)
real x(3)

N = [1,2,3]
x = [1,2,3]
write(*,*) bessel_j0(x(1))
write(*,*) bessel_j0(x)
write(*,*) bessel_jn(n(1),x(1))
write(*,*) bessel_jn(n,x(1))
write(*,*) bessel_jn(n(1),x)
write(*,*) bessel_jn(n,x)
! write(*,*) bessel_jn(n(1),n(3),x(1))
end program bessel_test

C:\gfortran\clf\bessel_test>gfortran bessel_test.f90 -obessel_test

C:\gfortran\clf\bessel_test>bessel_test
0.76519769
0.76519769 0.22389078 -0.26005197
0.44005057
0.44005057 0.11490349 1.95633546E-02
0.44005057 0.57672483 0.33905897
0.44005057 0.35283402 0.30906272

So gfortran treats the elemental case of BESSEL_J0 elementally
even though its documentation says that X shall be scalar, and
it treats the elemental cases of BESSEL_JN elementally even though
its docs says N and X shall be scalar, and N1723.pdf says the result
of BESSEL_JN(N,X) is scalar.

gfortran doesn't document nor does it implement the transformational
flavor of BESSEL_JN(N1,N2,X), partly because the N1723 documentation
is vague (I think) and partly because transformational intrinsics
are a royal PITB. Look at how vendors are more than happy to
implement all sorts of elemental intrinsics of perhaps marginal
utility to keep perhaps isolated customers happy, but you don't
see any rush to create, say, a FIRST transformational intrinsic
no matter how useful it might be. In this case, N1723.pdf doesn't
restrict the KIND of N1 and N2 nor of X so that if by the time
N1723.pdf becomes law, gfortran may have 5 integer KINDs and 4
real KINDs, so they would have to have in place a testing regime
for 5X5X4 = 100 different flavors of this transformational
intrinsic, just assuming it doesn't start acting like SPREAD which
we can tell for sure by examining N1723.pdf.

As I was saying before I starting ranting, gfortran doesn't
implement the transformational form of BESSEL_JN, so that if
the line in the above program that invokes it like that is
uncommented, it doesn't compile. A real failure for gfortran
is the fact that the REAL(10) flavors of the Bessel functions
are implemented in the compiler but not in the library so that
if the real(10) declaration of x were in force, the program would
compile but fail at the linking stage -- something that should
be documented in the manual and with compiler error messages or
simply fixed.

Definitely some rough edges in the draft and gfortran's initial
implementation can be seen here. Par for the course at such an
early stage I should say.

Steven G. Kargl

unread,
May 1, 2008, 1:05:52 AM5/1/08
to
In article <0LqdnagZta99p4TV...@comcast.com>,

"James Van Buskirk" <not_...@comcast.net> writes:

> A real failure for gfortran
> is the fact that the REAL(10) flavors of the Bessel functions
> are implemented in the compiler but not in the library so that
> if the real(10) declaration of x were in force, the program would
> compile but fail at the linking stage -- something that should
> be documented in the manual and with compiler error messages or
> simply fixed.

Send a documentation patch to gfortran developers.

'or simply fixed'. I'm speechless. Implementing functions that
appear in libm is certainly not simple. Well, I guess it could
be simple if one does not care about the quality of implementation.

--
Steve
http://troutmask.apl.washington.edu/~kargl/

FX

unread,
May 1, 2008, 4:06:56 AM5/1/08
to
>> I had been wondering why you added these, so I get the idea. Frankly,
>> I don't think it's that good an idea: this is a very specialized
>> topic, so people who require this will probably have their own
>> libraries anyway, and I really don't like that idea of adding special
>> case after special case. History already gave us a long and
>> complicated list of intrinsics, adding yet another one that has two
>> different forms is, in my opinion, a pain.
>
> Well, these functions are in standard C, or in most C implementations,
> anyway. So why not have them in Fortran?

I'm not complaining about having one form of Bessel functions, I'm
complaining about the bessel_{j,y}n having two different forms: my
critics above apply to the second form (the transformational one) of
these. I agree having that the elemental form is something that exists in
C and as a libm function usually, so that's fine, but the
transformational one, well, that's just yet another special case of the
language: Fortran is already a difficult language to know well, because
there are so many special cases, let's not add to that.

--
FX

Dan Nagle

unread,
May 1, 2008, 6:19:46 AM5/1/08
to
Hello,

On 2008-04-30 23:37:03 -0400, "James Van Buskirk" <not_...@comcast.net> said:

> Look at how vendors are more than happy to
> implement all sorts of elemental intrinsics of perhaps marginal
> utility to keep perhaps isolated customers happy, but you don't
> see any rush to create, say, a FIRST transformational intrinsic
> no matter how useful it might be.

Well, we added a findloc() intrinsic, actually in response
to a public comment from f03. (We didn't have time
to add it to f03.)

--
Cheers!

Dan Nagle

FX

unread,
May 2, 2008, 4:38:14 AM5/2/08
to
>> Testing with gfortran I found that RRSPACING, SCALE, SET_EXPONENT,
>> and, SPACING were broken like this, while EXPONENT, FRACTION, and
>> NEAREST seemed to work

I isolated the failure, which turned out to be in the mingw-win64
libraries, reported it to them and they fixed it.

--
FX

James Van Buskirk

unread,
May 3, 2008, 3:16:15 PM5/3/08
to
"FX" <cou...@alussinan.org> wrote in message
news:fvejtl$bb3$1...@nef.ens.fr...

Thanks. The fix didn't make it into the 20080502 snapshot, though.
Maybe next week? Are the problems with ERF, ERFC, ERFC_SCALED,
EXPONENT, FRACTION, GAMMA, and LOG_GAMMA for REAL(10) inputs on
Win64 also addressed by this fix, or are those a separate issue?
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.

FX

unread,
May 4, 2008, 6:46:43 AM5/4/08
to
> 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.

As a different issue, you might wonder why we allow them to compile and
wait for them to fail to link: it's because there is no way to know at
compilation time which math library will be used later (this is true in
all cases, and especially important in the case of cross-compilers), so
we don't want to put too strong a restriction on power users who link
with a different library.

--
FX

James Van Buskirk

unread,
May 5, 2008, 5:53:39 PM5/5/08
to
"FX" <cou...@alussinan.org> wrote in message
news:fvk46j$1fuu$1...@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?

James Van Buskirk

unread,
May 5, 2008, 8:24:42 PM5/5/08
to
"James Van Buskirk" <not_...@comcast.net> wrote in message
news:9tednQ43AIB-HILV...@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.

Janne Blomqvist

unread,
May 6, 2008, 4:05:35 AM5/6/08
to
On 2008-05-05, James Van Buskirk <not_...@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


--
Janne Blomqvist

James Van Buskirk

unread,
May 6, 2008, 6:01:57 AM5/6/08
to
"Janne Blomqvist" <f...@bar.invalid> wrote in message
news:slrng204a...@vipunen.hut.fi...

> 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. A normal approximation
wouldn't have that property.

BTW, I did for single-real-argument initialization expressions
what http://home.comcast.net/~kmbtib/Fortran_stuff/funr1.ZIP did
for single-real-argument ordinary expressions. The Bessel functions
failed as noted upthread, but also:

C:\gfortran\james\num_initr1>type mykinds.f90
module mykinds
implicit none
integer, parameter :: ck1 = kind('x')
integer, parameter :: ik1 = selected_int_kind(2)
integer, parameter :: ik2 = selected_int_kind(4)
integer, parameter :: ik4 = selected_int_kind(9)
integer, parameter :: ik8 = selected_int_kind(18)
integer, parameter :: Lk1 = ik1
integer, parameter :: Lk2 = ik2
integer, parameter :: Lk4 = ik4
integer, parameter :: Lk8 = ik8
integer, parameter :: sp = selected_real_kind(6,37)
integer, parameter :: dp = selected_real_kind(15,307)


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

type mytype
integer(ik4) data
end type mytype
end module mykinds

C:\gfortran\james\num_initr1>gfortran -c mykinds.f90

C:\gfortran\james\num_initr1>type isnan.f90
program isnan_prog
use mykinds
implicit none
logical, parameter :: lsp = isnan(1.0_sp)
logical, parameter :: ldp = isnan(1.0_dp)
logical, parameter :: lqp = isnan(1.0_qp)
real(merge(sp,dp,lsp)) xsp
real(merge(sp,dp,ldp)) xdp
real(merge(sp,dp,lqp)) xqp

write(*,*) lsp, ldp, lqp
write(*,*) kind(xsp), kind(xdp), kind(xqp)
end program isnan_prog
C:\gfortran\james\num_initr1>gfortran isnan.f90 -oisnan
isnan.f90:7.24:

real(merge(sp,dp,lsp)) xsp


1
Error: Constant expression required at (1)

isnan.f90:8.24:

real(merge(sp,dp,ldp)) xdp


1
Error: Constant expression required at (1)

isnan.f90:9.24:

real(merge(sp,dp,lqp)) xqp


1
Error: Constant expression required at (1)

isnan.f90:12.33:

write(*,*) kind(xsp), kind(xdp), kind(xqp)
1
Error: Symbol 'xdp' at (1) has no IMPLICIT type
isnan.f90:12.44:

write(*,*) kind(xsp), kind(xdp), kind(xqp)
1
Error: Symbol 'xqp' at (1) has no IMPLICIT type
isnan.f90:12.22:

write(*,*) kind(xsp), kind(xdp), kind(xqp)
1
Error: Symbol 'xsp' at (1) has no IMPLICIT type

C:\gfortran\james\num_initr1>type erfc_scaled.f90
program erfc_scaled_prog
use mykinds
implicit none
real(sp), parameter :: xsp = 2.0_sp
real(dp), parameter :: xdp = 2.0_dp
real(qp), parameter :: xqp = 2.0_qp
integer, parameter :: ksp = kind(1.0_sp)+0*erfc_scaled(xsp)
integer, parameter :: kdp = kind(1.0_dp)+0*erfc_scaled(xdp)
integer, parameter :: kqp = kind(1.0_qp)+0*erfc_scaled(xqp)
real(ksp) ysp
real(kdp) ydp
real(kqp) yqp

write(*,*) kind(ysp), kind(ydp), kind(yqp)
end program erfc_scaled_prog

C:\gfortran\james\num_initr1>gfortran erfc_scaled.f90 -oerfc_scaled
erfc_scaled.f90:10.11:

real(ksp) ysp


1
Error: Constant expression required at (1)

erfc_scaled.f90:11.11:

real(kdp) ydp


1
Error: Constant expression required at (1)

erfc_scaled.f90:12.11:

real(kqp) yqp


1
Error: Constant expression required at (1)

erfc_scaled.f90:14.33:

write(*,*) kind(ysp), kind(ydp), kind(yqp)
1
Error: Symbol 'ydp' at (1) has no IMPLICIT type
erfc_scaled.f90:14.44:

write(*,*) kind(ysp), kind(ydp), kind(yqp)
1
Error: Symbol 'yqp' at (1) has no IMPLICIT type
erfc_scaled.f90:14.22:

write(*,*) kind(ysp), kind(ydp), kind(yqp)
1
Error: Symbol 'ysp' at (1) has no IMPLICIT type

C:\gfortran\james\num_initr1>type loc.f90
program loc_prog
use mykinds
use ISO_C_BINDING, only: C_INTPTR_T
implicit none
real(sp) r
real(sp), parameter :: pr = 1.0_sp
integer(C_INTPTR_T), parameter :: lpr = loc(pr)
integer(C_INTPTR_T), parameter :: lr = loc(r)

write(*,*) lpr, lr
write(*,*) loc(x=r)
end program loc_prog

C:\gfortran\james\num_initr1>gfortran loc.f90 -oloc
loc.f90:7.47:

integer(C_INTPTR_T), parameter :: lpr = loc(pr)
1
Error: 'array' argument of 'loc' intrinsic at (1) must be a variable
loc.f90:8.41:

integer(C_INTPTR_T), parameter :: lr = loc(r)
1
Error: Intrinsic function 'loc' at (1) is not permitted in an initialization
exp
ression
loc.f90:10.17:

write(*,*) lpr, lr
1
Error: Symbol 'lpr' at (1) has no IMPLICIT type
loc.f90:10.21:

write(*,*) lpr, lr
1
Error: Symbol 'lr' at (1) has no IMPLICIT type
loc.f90:11.13:

write(*,*) loc(x=r)
1
Error: Can't find keyword named 'x' in call to 'loc' at (1)

C:\gfortran\james\num_initr1>type sizeof.f90
program sizeof_prog
use mykinds
use ISO_C_BINDING, only: C_SIZE_T
implicit none
integer(C_SIZE_T), parameter :: isp = sizeof(1.0_sp)

write(*,*) isp
end program sizeof_prog

C:\gfortran\james\num_initr1>gfortran sizeof.f90 -osizeof
sizeof.f90:5.40:

integer(C_SIZE_T), parameter :: isp = sizeof(1.0_sp)
1
Error: Intrinsic function 'sizeof' at (1) is not permitted in an
initialization
expression
sizeof.f90:7.17:

write(*,*) isp
1
Error: Symbol 'isp' at (1) has no IMPLICIT type

Discussion:

ERFC_SCALED fails very much like the Bessel functions do for
initialization expressions, and ISNAN seems to as well, except
that it may be MERGE that is causing the problems:

C:\gfortran\james\num_initr1>type merge.f90
program merge_prog
use mykinds
implicit none
logical, parameter :: l = .FALSE.
real(merge(sp,dp,l)) r

r = 42
write(*,*) r
end program merge_prog

C:\gfortran\james\num_initr1>gfortran merge.f90 -omerge
merge.f90:5.22:

real(merge(sp,dp,l)) r


1
Error: Constant expression required at (1)

merge.f90:7.4:

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

So MERGE is definitely a problem for f03 because it's an elemental
function, not a transformational one. Further testing has shown
that ISNAN is also a problem:

C:\gfortran\james\num_initr1>type transfer1.f90
program isnan_prog
use mykinds
implicit none
logical, parameter :: lsp = .FALSE.
real(sp+0*transfer(lsp,0)) xsp

write(*,*) lsp
write(*,*) kind(xsp)
end program isnan_prog
C:\gfortran\james\num_initr1>gfortran transfer1.f90 -otransfer1

C:\gfortran\james\num_initr1>transfer1
F
4

C:\gfortran\james\num_initr1>type transfer.f90
program isnan_prog
use mykinds
implicit none
logical, parameter :: lsp = isnan(1.0_sp)
logical, parameter :: ldp = isnan(1.0_dp)
logical, parameter :: lqp = isnan(1.0_qp)
real(sp+0*transfer(lsp,0)) xsp
real(dp+0*transfer(ldp,0)) xdp
real(qp+0*transfer(lqp,0)) xqp

write(*,*) lsp, ldp, lqp
write(*,*) kind(xsp), kind(xdp), kind(xqp)
end program isnan_prog
C:\gfortran\james\num_initr1>gfortran transfer.f90 -otransfer
transfer.f90:7.28:

real(sp+0*transfer(lsp,0)) xsp


1
Error: Constant expression required at (1)

transfer.f90:8.28:

real(dp+0*transfer(ldp,0)) xdp


1
Error: Constant expression required at (1)

transfer.f90:9.28:

real(qp+0*transfer(lqp,0)) xqp


1
Error: Constant expression required at (1)

transfer.f90:12.33:

write(*,*) kind(xsp), kind(xdp), kind(xqp)
1
Error: Symbol 'xdp' at (1) has no IMPLICIT type
transfer.f90:12.44:

write(*,*) kind(xsp), kind(xdp), kind(xqp)
1
Error: Symbol 'xqp' at (1) has no IMPLICIT type
transfer.f90:12.22:

write(*,*) kind(xsp), kind(xdp), kind(xqp)
1
Error: Symbol 'xsp' at (1) has no IMPLICIT type

Now, LOC has somewhat different problems. The first error for LOC
is actually OK because I attempted to use it on a named constant,
contrary to the documentation. However, the documentation says that
the keyword for the argument is 'X' but the error message says it's
'ARRAY'. The last error message again shows that the documentation
is incorrect regarding the keyword for LOC's dummy argument. The
second error message says that LOC can't be used in initialization
expressions. Mmmm... I guess that's fair because the result can't
be known until link time but it's needed at compile time. So the
only complaint I have about LOC is the keyword: the compiler and
the documentation should agree.

SIZEOF is banned from initialization expression, but it shouldn't be
because it's the kind of thing you need if you're going to transfer
an instance of a derived type to an array that isn't an automatic
data object, SIZEOF would be the easiest thing to use to dimension
the destination array. Well, I always use SIZE and TRANSFER for this
task, but some people prefer SIZEOF.

In summary, ISNAN, ERFC_SCALED, and MERGE seem to have the same
sort of problems in initialization expressions as do the Bessel
functions: they are allowed in initialization expressions but not
evaluated in time to be used in KIND numbers.

LOC is OK except that the keyword doesn't agree with the
documentation.

SIZEOF is forbidden in initialization expressions, which is OK
because it's an extension (and what is its class: can it be
considered transformational?) but it seems illogical nonetheless.

I didn't notice any singularities regarding the other single-real-
argument intrinsics for initialization expressions.

Janne Blomqvist

unread,
May 6, 2008, 1:32:14 PM5/6/08
to
On 2008-05-06, James Van Buskirk <not_...@comcast.net> wrote:
> "Janne Blomqvist" <f...@bar.invalid> wrote in message
> news:slrng204a...@vipunen.hut.fi...
>
>> On 2008-05-05, James Van Buskirk <not_...@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

James Van Buskirk

unread,
May 6, 2008, 3:50:09 PM5/6/08
to
"Janne Blomqvist" <f...@bar.invalid> wrote in message
news:slrng215g...@vipunen.hut.fi...

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

I didn't think this was the case because when the REAL(10) versions
of the Bessel functions get invoked at runtime, you get a link error.

> 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).

Great! It's good to see that some of the problems I have been
finding have been fixed. Did ISNAN, ERFC_SCALED, and MERGE, which
seem to have the same issue regarding initialization expressions
as BESSEL_J0, BESSEL_J1, BESSEL_JN, BESSEL_Y0, BESSEL_Y1, and
BESSEL_YN get fixed automatically by this patch or do I have to
hunt down all of the instances of this behavior one by one?

> 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.

That's a step that mingw doesn't carry out for expressions that can't
be constant folded. In that case what you see is link failures or
stub functions that return incorrect values as seen in

http://home.comcast.net/~kmbtib/Fortran_stuff/funr1.ZIP

> 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.

I have violated my own rule about deducing what the compiler is doing
from HLL tests rather than examining its actual output. HLL tests
proved quite misleading in this instance, don't you think?

The trailing zeros seemed mysterious to me but they can be
explained by a combination of my forgetting to set the x87 control
word precision control to b'11' and a lucky input value. Setting
(actually resetting) the control word and rolling the dice again
we get:

C:\gfortran\clf\bessel_test>type finit.s
.text
.globl _xfinit_
.def _xfinit_; .scl 2; .type 32; .endef
_xfinit_:
finit
ret

C:\gfortran\clf\bessel_test>type test8.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 test8
use mykinds
implicit none
real(qp), parameter :: x = 3.5


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)

call xfinit


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

end program test8

C:\gfortran\clf\bessel_test>gfortran test8.f90 finit.s -otest8

C:\gfortran\clf\bessel_test>test8
xj0 = -0.38012773998726337739
xj1 = 0.13737752736232718572
xjn = 0.13737752736232718572
xy0 = 0.18902194392082650675
xy1 = 0.41018841788751188289
xyn = 0.41018841788751188289

So mystery solved! I didn't see a PR for the f.p. control word.
Precision control has to be b'11' for REAL(10) arithmetic to behave
normally or be output normally and it's annoying to have to do this
by hand every time. Microsoft definitely wants to set it at b'10' for
Win64, but I don't know whether this is only a prejudice against
accurate results or if it breaks some library that gfortran uses
if it's not set up this way. gfortran uses the 128-byte red zone
after all, and this is to my mind a much more dangerous violation
of the Win64 way of doing things, especially for multithreaded
programs.

> 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.

There are several things that have to happen for an intrinsic to be
fully supported in f03. It needs to be implemented for all
combinations of KINDs of inputs and for initialization, specification,
and ordinary expressions. Some accuracy and even speed is desirable.
I'm going through the list of intrinsics as time permits and for now
just checking existence for all KINDs and for initialization and
ordinary expressions. I haven't done specification expressions yet
because so far they seem to be going the same way as ordinary
expressions and gfortran has showstopper issues that have been
outstanding for quite a while that might make it more difficult to
implement a testing regime for specification expressions. I haven't
touched the issue of array arguments. It seems Dick is doing
something in that area.

James Van Buskirk

unread,
May 6, 2008, 4:37:04 PM5/6/08
to
"James Van Buskirk" <not_...@comcast.net> wrote in message
news:ed2dncmj0oLsK73V...@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?

FX

unread,
May 6, 2008, 5:04:28 PM5/6/08
to
> 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

James Van Buskirk

unread,
May 6, 2008, 5:46:35 PM5/6/08
to
"FX" <cou...@alussinan.org> wrote in message
news:fvqh4r$1n5d$1...@nef.ens.fr...

C:\gfortran\clf\kahan>type finit.s


.text
.globl _xfinit_
.def _xfinit_; .scl 2; .type 32; .endef
_xfinit_:

.globl XFINIT
.def XFINIT; .scl 2; .type 32; .endef
XFINIT:
finit
ret
.globl _getcw_
.def _getcw_; .scl 2; .type 32; .endef
_getcw_:
.globl GETCW
.def GETCW; .scl 2; .type 32; .endef
GETCW:
xorl %eax, %eax
push %rax
fstcw (%rsp)
pop %rax
ret

C:\gfortran\clf\kahan>type kahan.f90
program kahan
implicit none
integer, parameter :: C_INT16_T = selected_int_kind(4)
integer(C_INT16_T), external :: getcw

write(*,'(a,z4.4)') 'Initial control word = ', getcw()
call xfinit
write(*,'(a,z4.4)') 'After reset control word = ', getcw()
end program kahan

C:\gfortran\clf\kahan>as finit.s -ofinit.o

C:\gfortran\clf\kahan>gfortran kahan.f90 finit.o -okahan

C:\gfortran\clf\kahan>kahan
Initial control word = 027F
After reset control word = 037F

C:\gfortran\clf\kahan>ifort kahan.f90 finit.o
Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.
ifort: Command line warning: unrecognized source type 'finit.o'; object file
ass
umed

Microsoft (R) Incremental Linker Version 8.00.40310.39
Copyright (C) Microsoft Corporation. All rights reserved.

-out:kahan.exe
-subsystem:console
kahan.obj
finit.o

C:\gfortran\clf\kahan>kahan
Initial control word = 027F
After reset control word = 037F

The problem is that on Win64, the default value for the floating
point control word seems to be z'027F' which sets precision control
to 53 bits. This is not an issue for ifort because its extra
precision KIND is quadruple precision. For gfortran this is a big
problem because its extra precision KIND is extended precision and
so needs precision control set to 64 bits. As can be seen from the
above example, if gfortran were sporting enough to just issue an
FINIT at program start, precision control would be set to 64 bits
and life could go on from there. It's rather silly to set precision
control at 53 bits as Microsoft seems to want to do because on Win64
you are guaranteed SSE2 so that if you wanted to avoid 80-bit (i.e.
64-bit precision) temporaries you could just use SSE2 for your
calculations and your code would be more performant as well.

Setting the hardware default, z'037F', as by FINIT on program entry
seems to me to make sense on Win64, but I am not sure whether some
Microsoft library that gfortran needs will behave in an unexpected
manner given this setting. For sure REAL(10) arithmetic is busted
without this setting. Microsoft doesn't know if they are on foot or
horseback regarding x87 on Win64. For e.g. an early version of ml64
considered mmx operations an error because the programmers for that
assembler had heard a rumor that Win64 wouldn't preserve x87 state
across task switches. In fact it's only for ring-0 code that x87
state doesn't get preserved, so for ordinary ring-3 code x87 and so
mmx operations are OK.

The gcc docs state that 64-bit precision (as opposed to 53-bit
precision) is the default but these docs are hopelessly out of date
for Win64.

Ron Ford

unread,
May 6, 2008, 6:34:28 PM5/6/08
to

What James is doing looks so exciting. I wish I could follow it better. I
understand his kinds module:


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

This addresses how to test for what a machine gives you for quad as opposed
to double precision, and sets qp to the greater. Main doesn't look too
bad:


program test8
use mykinds
implicit none
real(qp), parameter :: x = 3.5
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)

call xfinit
write(*,*) 'xj0 = ', xj0
write(*,*) 'xj1 = ', xj1
write(*,*) 'xjn = ', xjn
write(*,*) 'xy0 = ', xy0
write(*,*) 'xy1 = ', xy1
write(*,*) 'xyn = ', xyn
end program test8

I remember bessel functions from school, so it'd be great to have an
opportunity to use them. Then there's xfinit :

C:\gfortran\clf\bessel_test>type finit.s
.text
.globl _xfinit_
.def _xfinit_; .scl 2; .type 32; .endef
_xfinit_:
finit
ret

Here's where the wheels come off for me. What is this? I notice that
James uses many of the tools that Intel vends. If a person doesn't have a
thousand dollars to spend on the like, does it make any sense to pursue it?
--
Ron Ford

James Van Buskirk

unread,
May 10, 2008, 3:34:10 PM5/10/08
to
"Janne Blomqvist" <f...@bar.invalid> wrote in message
news:slrng215g...@vipunen.hut.fi...

> On 2008-05-06, James Van Buskirk <not_...@comcast.net> wrote:

>> 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).

Now that a new gfortran build is available:

C:\gfortran\james\num_initr1>gfortran -v
Built by Equation Solution (http://www.Equation.com).
Using built-in specs.
Target: x86_64-pc-mingw32
Configured with:
../gcc-4.4-20080509-mingw/configure --host=x86_64-pc-mingw32 --
build=x86_64-unknown-linux-gnu --target=x86_64-pc-mingw32 --prefix=/home/gfortra
n/gcc-home/binary/mingw32/native/x86_64/gcc/4.4-20080509 --with-gmp=/home/gfortr
an/gcc-home/binary/mingw32/native/x86_64/gmp --with-mpfr=/home/gfortran/gcc-home
/binary/mingw32/native/x86_64/mpfr --with-sysroot=/home/gfortran/gcc-home/binary
/mingw32/cross/x86_64/gcc/4.4-20080509 --with-gcc --with-gnu-ld --with-gnu-as
--
disable-shared --disable-nls --disable-tls --enable-languages=c,fortran --enable
-threads=win32 --enable-libgomp --disable-win32-registry
Thread model: win32
gcc version 4.4.0 20080509 (experimental) (GCC)

I have been able to test the bug fixes discussed in this thread.
The MPFR fix has indeed fixed most of the initialization expression
stuff. Only two exceptions for functions that can accept one REAL
argument:

C:\gfortran\test\intrinsics>type erfc_scaled.f90
program erfc_scaled_init
real(4), parameter :: x = 2
real(4), parameter :: y = erfc_scaled(x)
integer(int(3+4*y)) j

write(*,*) kind(j)
end program erfc_scaled_init

C:\gfortran\test\intrinsics>gfortran erfc_scaled.f90 -oerfc_scaled
erfc_scaled.f90:4.21:

integer(int(3+4*y)) j


1
Error: Constant expression required at (1)

C:\gfortran\test\intrinsics>type sizeof.f90
program sizeof_init
real(4), parameter :: x = 2
integer, parameter :: i = sizeof(x)
integer(i) j

write(*,*) kind(j)
end program sizeof_init

C:\gfortran\test\intrinsics>gfortran sizeof.f90 -osizeof
sizeof.f90:3.28:

integer, parameter :: i = sizeof(x)


1
Error: Intrinsic function 'sizeof' at (1) is not permitted in an
initialization
expression

sizeof.f90:4.11:

integer(i) j
1
Error: Parameter 'i' at (1) has not been declared or is a variable, which
does n
ot reduce to a constant expression

But the fix for ordinary expressions going back to the original
post of the thread hasn't yet propagated through mingw to gfortran.

Thanks for the progress on this issue!

0 new messages