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