Some mysterious calculations by g95

[Hans Peschke]

Hello g95-folks :)

the following code produces some mysterious results compiled with g95 (with
several versions, see below). I'm running this on a Thinkpad X200 with Ubuntu
9.04 32 Bit.

<code>
program underflow_test
implicit none

integer, parameter :: pr = kind(0.0D0)
real(pr) :: x, hx

x = tiny(x)
hx = 0.5_pr * x

write(*,*) 'x ! = tiny(x) = ', x
write(*,*) 'half_tiny:'
write(*,*) 'hx ! = 0.5_pr * x = ', hx
write(*,*) 'x * 0.5_pr = ', x * 0.5_pr
write(*,*) 'scale(x, -1) = ', scale(x, -1)
write(*,*) 'scale(tiny(x), -1) = ', scale(tiny(x), -1)

write(*,*) 'quarter_tiny:'
write(*,*) 'scale(x, -2) = ', scale(x, -2)
write(*,*) 'scale(hx, -1) = ', scale(hx, -1)
write(*,*) 'scale(x * 0.5_pr, -1) = ', scale(x * 0.5_pr, -1)
write(*,*) 'scale(tiny(x) * 0.5_pr, -1) = ', scale(tiny(x) * 0.5_pr, -1)
write(*,*) 'x * 0.25_pr = ', x * 0.25_pr
write(*,*) 'scale(hx, -2) = ', scale(hx, -2)

end program underflow_test
</code>

<test>
sep@sep-mobile:~/ws/fortran/test$ g95 -v
Using built-in specs.
Target:
Configured with: ../configure --enable-languages=c : (reconfigured) ../configure
--enable-languages=c
Thread model: posix
gcc version 4.1.3 (g95 0.92!) Jun 18 2009
sep@sep-mobile:~/ws/fortran/test$ g95 underflow_test.f90 -o underflow_test
sep@sep-mobile:~/ws/fortran/test$ ./underflow_test
x ! = tiny(x) = 2.2250738585072014E-308
half_tiny:
hx ! = 0.5_pr * x = 1.1125369292536007E-308
x * 0.5_pr = 1.1125369292536007E-308
scale(x, -1) = 0. | ??
scale(tiny(x), -1) = 1.1125369292536007E-308 | ??
quarter_tiny:
scale(x, -2) = +Inf ! very wrong
scale(hx, -1) = 5.562684646268003E-309
scale(x * 0.5_pr, -1) = 5.562684646268003E-309
scale(tiny(x) * 0.5_pr, -1) = 2.781342323134E-309 | wrong
x * 0.25_pr = 5.562684646268003E-309
scale(hx, -2) = 2.781342323134E-309 | wrong
sep@sep-mobile:~/ws/fortran/test$
sep@sep-mobile:~/ws/fortran/test$ # install other g95 version #
sep@sep-mobile:~/ws/fortran/test$ g95 -v
Using built-in specs.
Target:
Configured with: ../configure --enable-languages=c
Thread model: posix
gcc version 4.0.3 (g95 0.92!) Feb 26 2009
sep@sep-mobile:~/ws/fortran/test$ g95 underflow_test.f90 -o underflow_test
sep@sep-mobile:~/ws/fortran/test$ ./underflow_test
x ! = tiny(x) = 2.2250738585072014E-308
half_tiny:
hx ! = 0.5_pr * x = 1.1125369292536007E-308
x * 0.5_pr = 1.1125369292536007E-308
scale(x, -1) = 0. | ??
scale(tiny(x), -1) = 1.1125369292536007E-308 | ??
quarter_tiny:
scale(x, -2) = +Inf ! very wrong
scale(hx, -1) = 5.562684646268003E-309
scale(x * 0.5_pr, -1) = 5.562684646268003E-309
scale(tiny(x) * 0.5_pr, -1) = 2.781342323134E-309 ! loose of precision
x * 0.25_pr = 5.562684646268003E-309
scale(hx, -2) = 2.781342323134E-309 ! loose of precision
sep@sep-mobile:~/ws/fortran/test$
</test>

??: I don't understand this.

gfortran 4.3.3 calculates all correctly. Compiling with -O0 doesn't help.

Do you have any ideas or explenations on this? Confirmations?

Thx and Best regards

Hans Peschke

[Jimmy]

Interesting results which I can confirm by using g95-MinGW.exe on a Vista-32
PC. On first looking at your program I would expect all but the first
result to be zero, so it looks like a g95 bug which should be reported to
Andy. However, might be something subtle here such as conversion from 80
bit FP to 64 bit FP. If no explanation is forthcoming from this Group then
perhaps comp.lang.fortran specialists can help.

Jimmy.

[Udo Grabowski]

Jimmy wrote:

>> sep@sep-mobile:~/ws/fortran/test$ ./underflow_test
>> x ! = tiny(x) = 2.2250738585072014E-308
>> half_tiny:
>> hx ! = 0.5_pr * x = 1.1125369292536007E-308
>> x * 0.5_pr = 1.1125369292536007E-308
>> scale(x, -1) = 0. | ??
>> scale(tiny(x), -1) = 1.1125369292536007E-308 | ??
>> quarter_tiny:
>> scale(x, -2) = +Inf ! very wrong
>> scale(hx, -1) = 5.562684646268003E-309
>> scale(x * 0.5_pr, -1) = 5.562684646268003E-309
>> scale(tiny(x) * 0.5_pr, -1) = 2.781342323134E-309 | wrong
>> x * 0.25_pr = 5.562684646268003E-309
>> scale(hx, -2) = 2.781342323134E-309 | wrong

Here are the results of the opensolaris jury (06.2009 with SunStudio 12.1):

x ! = tiny(x) = 2.2250738585072013E-308
half_tiny:
hx ! = 0.5_pr * x = 1.1125369292536006E-308
x * 0.5_pr = 1.1125369292536006E-308
scale(x, -1) = 1.1125369292536006E-308
scale(tiny(x), -1) = 1.1125369292536006E-308
quarter_tiny:
scale(x, -2) = 5.562684646268003E-309

scale(hx, -1) = 5.562684646268003E-309
scale(x * 0.5_pr, -1) = 5.562684646268003E-309

scale(tiny(x) * 0.5_pr, -1) = 5.562684646268003E-309

x * 0.25_pr = 5.562684646268003E-309

scale(hx, -2) = 2.781342323134002E-309

The handling of 80 vs. 64 bit precision seems to be commonly done wrong.
When setting flag -fns (floating point nonstandard(!)) this changes to the
standard (!) behaviour (almost):

x ! = tiny(x) = 2.2250738585072013E-308
half_tiny:
hx ! = 0.5_pr * x = 0.0E+0
x * 0.5_pr = 0.0E+0
scale(x, -1) = 0.0E+0
scale(tiny(x), -1) = 1.1125369292536006E-308
quarter_tiny:
scale(x, -2) = 0.0E+0
scale(hx, -1) = 0.0E+0
scale(x * 0.5_pr, -1) = 0.0E+0
scale(tiny(x) * 0.5_pr, -1) = 5.562684646268003E-309
x * 0.25_pr = 0.0E+0
scale(hx, -2) = 0.0E+0

So g95 is not alone at least ...
--
Dr. Udo Grabowski email: udo.gr...@imk.fzk.de
Institut f. Meteorologie und Klimaforschung ASF,Forschungszentrum Karlsruhe
Postfach 3640, 76021 Karlsruhe, Germany Tel: (+49) 7247 82-6026
http://www.fzk.de/imk/asf/sat/grabowski/ Fax: " -7026

[jfh]

> Dr. Udo Grabowski email: udo.grabow...@imk.fzk.de

> Institut f. Meteorologie und Klimaforschung ASF,Forschungszentrum Karlsruhe
> Postfach 3640, 76021 Karlsruhe, Germany Tel: (+49) 7247 82-6026http://www.fzk.de/imk/asf/sat/grabowski/ Fax: " -7026

The result of scale(x,i) should be x*b**i, where b is the base of the
floating-point model of x (usually 2), provided this result is within
range; if not, the Fortran 95 or 2003 standard says the result is
processor dependent. Tiny(x) is of course the smallest positive number
within range, so don't be surprised that different processors give
different answers when the "expected" result is even smaller.

-- John Harper

[sep]

@Jimmy: already sent to Andy.

I think, there are some distinct problems to consider:

1) Either a Compiler does gradual-underflow unless zero or flushes-to-
zero. Imo this policy has to hold for every function/operator.
2) Obviously scale() produces different results depending on the first
argument is either a variable holding a value or the result of a
function with the same value, which is not the way it should be from
the developer-side of view.
3) Scale() produces inconsistent results (especially the Inf+) and is
not monotone in the second argument.
4) The processor isn't the main point i think, we get compiler
dependence results, which is certainly not possible to prevent in all
situations, but it is possible to achieve in mathematical
calculations, because there are not so many correct results.

I don't understand the problem with the conversion from 80 to 64 FP
Bits.

thx and Happy coding :)

Hans Peschke