The following test
integer :: i, n
real :: x, y(10)
x =1.0
n = 10
do i = 1, n
y(i) = 1.0/x
x = 2.0*x
end do
print *, y
end
when compiled with -O -ffast-math -mrecip, gives
0.99999994 0.49999997 0.24999999 0.12499999
6.24999963E-02 3.12499981E-02 1.56249991E-02 7.81249953E-03 3.90624977E-03
1.95312488E-03
the inverse of an integer power of 2 is not an integer power of 2 with -mrecip
on i686-apple-darwin9 and x86_64-unknown-linux-gnu.
Since x0*(2.0-x*x0) does not have round-off errors when x=2.0**i and
x0=2.0**(-i), I assume that rcp* don't return an integer power of 2 if the
input is an integer power of 2.
In addition the result of a Newton-Raphson iteration is not the same from above
of from below:
real :: x, y
x = 1.0
y = nearest(x,x)
print *, y*(2.0-x*y)
y = nearest(x,-x)
print *, y*(2.0-x*y)
end
gives
1.00000000
0.99999994
This result is quite unfortunate since the integer powers of 2 are the only
floating point numbers having an exact inverse. This numerical error probably
not be fixed in an effective way, but should probably be documented.
As a side note, I stumbled on the problem while trying to run the aermod.f90
polyhedron test with -mrecip. I have been chasing a resulting bus error at
execution for a while until I found the culprit as line 35369 of the subroutine
NUMRISE:
IF ( FLOAT(NNP/NP).EQ.FLOAT(NNP)/XNP ) THEN
for which the comparison fails even if NP=1 and XNP=1.0 with -mrecip. It
follows that the variable NN, initialized within this IF block, is not
initialized, thus in some case leading to a variable DELN negative or bigger
than 1 at line 35514, leading to an access out of bounds. Note also that there
is no point to discuss (at least with me) the way the program is written: I am
well aware of the dangers of testing floating points for equality!
A way to limit this kind of problems while letting people use -mrecip if it
speeds up their code, could be (if possible) to use the exact division within
IF expressions.
--
Summary: 1.0 is not the inverse of 1.0 with -mrecip on x86
Product: gcc
Version: 4.3.0
Status: UNCONFIRMED
Severity: normal
Priority: P3
Component: target
AssignedTo: unassigned at gcc dot gnu dot org
ReportedBy: dominiq at lps dot ens dot fr
GCC build triplet: i686-apple-darwin9
GCC host triplet: i686-apple-darwin9
GCC target triplet: i686-apple-darwin9
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34702
