[Bug fortran/29383] Fortran 2003/F95[TR15580:1999]: Floating point exception (IEEE) support

["kargl at gcc dot gnu.org" <gcc-bugzilla at gcc dot gnu dot org>]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=29383

kargl at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |kargl at gcc dot gnu.org

--- Comment #9 from kargl at gcc dot gnu.org 2012-10-09 17:00:42 UTC ---
(In reply to comment #8)

> In response to your other email (that it is very very hard), can you explain
> a bit why if you've got the time/inclination? My first thought was that it 
> would be a fairly simple pass through/reimplementation of the fpclassify 
> functionality of c.

There are a few reason (and this is just my opinion).

1) gcc runs on numerous cpu architectures (i386, x86_64,
   arm, sparc, powerpc, etc).  Some (many? all?) do not
   have hardware support for IEEE754.  An implementation
   needs to be able to work on most (all?) of these cpus.

2) gcc runs on numerous operating systems.  The operating
   systems may or may not have sufficient support to allow
   an efficient IEEE754 implementation (a software
   implementation of IEEE754 will be slow).

3) Using C interoperability to access the fenv.h facility
   seems appealing, but gcc would need to deal with systems
   that do not have fenv.h.

4) gcc performs constant folding with mpfr, and in gfortran
   this constant folding is performed in round-to-nearest mode.
   gfortran would need to be update to delay constant folding
   until the rounding mode has been established.

   program foo
      use ieee_arithmetic
      real x
      call ieee_set_rounding_mode(ieee_down)
      x = 0.3 + 0.1999  ! Should be down in round-down?
      print '(Z0)', x
   end program foo

5) Finally, for some intrinsic functions, it is not possible
   to map to a libm routine without now putting a wrapper around
   the libm routine.  For example, consider FRACTION().  The standard
   states:

       Result Value.  The result has the value X*b**(-e), where b
       and e are as defined in 13.4 for the model representation of
       X.  If X has the value zero, the result has the value zero.
       If X is an IEEE infinity, the result is that infinity. If X
       is an IEEE NaN, the result is that NaN.

   The word "that" in "the result is that NaN" is problematic. 
   There are something like 2**p bit patterns that much NaN.
   Unfortunately, mapping FRACTION to libm's frexp gives the
   wrong NaN.  See PR 48979 for a longer explanation.

One can side step 1), 2), and 3) by unilaterally returning
false for ieee_support_standard(), which of course indicates
that the IEEE754 isn't standard.  4) and 5) are much harder
to fix.

Of course, I could be wrong.