[Bug fortran/15441] New: [gfortran] RRSPACING broken for denormals

[Bug fortran/15441] New: [gfortran] RRSPACING broken for denormals

[tobi at gcc dot gnu dot org]

Quoting from http://gcc.gnu.org/ml/gcc-patches/2004-05/msg00851.html:

The existing rrspacing code to handle denormals is broken. This occurs when x <
TINY(x), ie. expn == 0, frac != 0 I suspect the right shift should be a logical
(unsigned) shift, not an arithmetic (signed) one.

-- 
           Summary: [gfortran] RRSPACING broken for denormals
           Product: gcc
           Version: 3.5.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P2
         Component: fortran
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: tobi at gcc dot gnu dot org
                CC: gcc-bugs at gcc dot gnu dot org,pbrook at gcc dot gnu
                    dot org


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

[Bug fortran/15441] [gfortran] RRSPACING broken for denormals

[pinskia at gcc dot gnu dot org]

------- Additional Comments From pinskia at gcc dot gnu dot org  2004-05-14 17:32 -------
Confirmed.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |NEW
     Ever Confirmed|                            |1
           Keywords|                            |wrong-code
   Last reconfirmed|0000-00-00 00:00:00         |2004-05-14 17:32:55
               date|                            |


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

[Bug fortran/15441] RRSPACING broken for denormals

[wf_cs at yahoo dot com]

------- Additional Comments From wf_cs at yahoo dot com  2004-09-07 09:56 -------
Hi, Tobi
  I think the code considered the denormal reals. Can you file a concrete test? 
In fact, the original code just didn't deal with zero. After applying your 
patch, I could not find another test to confirm this.

-- 


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

[Bug fortran/15441] RRSPACING broken for denormals

[Tobias dot Schlueter at physik dot uni-muenchen dot de]

------- Additional Comments From Tobias dot Schlueter at physik dot uni-muenchen dot de  2004-09-13 17:36 -------
Subject: Re:  RRSPACING broken for denormals

wf_cs at yahoo dot com wrote:
> ------- Additional Comments From wf_cs at yahoo dot com  2004-09-07 09:56 -------
> Hi, Tobi
>   I think the code considered the denormal reals. Can you file a concrete test? 
> In fact, the original code just didn't deal with zero. After applying your 
> patch, I could not find another test to confirm this.
> 
I have to admit that I don't know a test. I added this because Paul said the
code wouldn't work for denormals, but didn't verify this myself. I will verify
this myself and will close the bug, if appropriate.

- Tobi


-- 


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

[Bug fortran/15441] RRSPACING broken for denormals

[tobi at gcc dot gnu dot org]

------- Additional Comments From tobi at gcc dot gnu dot org  2004-09-14 12:59 -------
I tried this:
x = 1.
i = 0
do while (x>0.)
   x = x/2.; i = i+1
   write(*,10) i, x, rrspacing(x)
10 format("2**(-",i3,") rrspacing(",e8.2,") =",e8.2)
end do
end
the output reads:
2**(-  1) rrspacing(0.50E+00) =0.84E+07
2**(-  2) rrspacing(0.25E+00) =0.84E+07
2**(-  3) rrspacing(0.13E+00) =0.84E+07
...
2**(-147) rrspacing(0.56E-44) =0.84E+07
2**(-148) rrspacing(0.28E-44) =0.84E+07
2**(-149) rrspacing(0.14E-44) =0.84E+07
2**(-150) rrspacing(0.00E+00) =0.00E+00

Now if I understand things correctly and got my facts right, the relative
spacing should become larger once we hit denormals, because we don't have as
many digits left for recording the mantissa. More precisely, IIRC numbers with
exponent in [-149,-126] are denormals, and their relative spacing should
increase by a factor of two at every iteration of the above program, for
RRSPACING to be also meaningful in the case of denormals.

I understand that denormals are not part of the floating point model specified
in the standard, and therefore the behavior of RRSPACING is not specified for
denormals, therefore I'm changing the severity to 'enhancement', but I'm not
convinced that we're doing The Right Thing.

BTW the program compiled with ifc gives the same output as with gfortran.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
           Severity|normal                      |enhancement


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

[Bug fortran/15441] RRSPACING broken for denormals

[wf_cs at yahoo dot com]

------- Additional Comments From wf_cs at yahoo dot com  2004-09-15 09:40 -------
I understand this problem this way. In fortran standard, 

RRSPACING (X) = |X * POW (2, -e)| * POW (2, p) = FRACTION (X) * POW (2, p)
 
So the result's exponenet is p. And if X is normalized, X's fraction part
is the result's fraction. If X is denormalized, to get the X's fraction we
shift X's fraction part to left until the first '1' is removed.

That's why if X's fraction has only '1', we will get the same result as you 
test.

I agree that the behavior of RRSPACING is not specified for denormals. But I 
think we do the Right Thing as I explained.
I attach a pathch to give more comments about the code and add some folding. 
How do you think?

-- 


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

[Bug fortran/15441] RRSPACING broken for denormals

[tobi at gcc dot gnu dot org]

------- Additional Comments From tobi at gcc dot gnu dot org  2004-09-15 10:24 -------
I looked at the standard again, and it says ""[RRSPACING] Returns the reciprocal
of the relative spacing of model numbers near the argument". I.e., it is not
only specified for model numbers (as I wrongly assumed), but for a denormal it
should return the spacing of a nearby normalized number, and so indeed, you're
right (we could also return RRSPACING(0), and agree with the standard, but I
don't think that's what people would expect).

> I attach a pathch to give more comments about the code and add some folding. 
> How do you think?
> 

Post it to the mailing list, and let the maintainers decide :-)

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
         Resolution|                            |INVALID


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

[Bug fortran/15441] RRSPACING broken for denormals

[kargl at gcc dot gnu dot org]

------- Comment #11 from kargl at gcc dot gnu dot org  2006-10-09 20:57 -------
Fixed on trunk (until someone tells me ldexp doesn't exist).


-- 

kargl at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|ASSIGNED                    |RESOLVED
         Resolution|                            |FIXED
   Target Milestone|---                         |4.2.0


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

[Bug fortran/15441] RRSPACING broken for denormals

[kargl at gcc dot gnu dot org]

------- Comment #10 from kargl at gcc dot gnu dot org  2006-10-09 20:55 -------
Subject: Bug 15441

Author: kargl
Date: Mon Oct  9 20:55:29 2006
New Revision: 117584

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=117584
Log:
2006-10-06  Steven G. Kargl  <kargl@gcc.gnu.org>

        * gfortran.h: Define GFC_MPFR_TOO_OLD via mpfr version info.
        * arith.c (arctangent, gfc_check_real_range): Use it.   
        * simplify.c (gfc_simplify_atan2, gfc_simplify_exponent,
        gfc_simplify_log, gfc_simplify_nearest): Use it.

        PR fortran/15441
        PR fortran/29312
        * iresolve.c (gfc_resolve_rrspacing): Give rrspacing library
        routine hidden precision argument.
        (gfc_resolve_spacing): Give spacing library routine hidden
        precision, emin - 1, and tiny(x) arguments.
        * simplify.c (gfc_simplify_nearest): Remove explicit subnormalization.
        (gfc_simplify_rrspacing): Implement formula from Fortran 95 standard.
        (gfc_simplify_spacing): Implement formula from Fortran 2003 standard.
        * trans-intrinsic.c (gfc_intrinsic_map_t) Declare rrspacing and
        spacing via LIBF_FUNCTION
        (prepare_arg_info, call_builtin_clz, gfc_conv_intrinsic_spacing,
        gfc_conv_intrinsic_rrspacing): Remove functions.
        (gfc_conv_intrinsic_function): Remove calls to
        gfc_conv_intrinsic_spacing and gfc_conv_intrinsic_rrspacing.
        * f95-lang.c (gfc_init_builtin_functions): Remove __builtin_clz,
        __builtin_clzl and __builtin_clzll


2006-10-06  Steven G. Kargl  <kargl@gcc.gnu.org>

        PR fortran/15441
        PR fortran/29312
        * configure.ac: Add HAVE_LDEXPF, HAVE_LDEXP, and HAVE_LDEXPL
        * m4/spacing.m4: New file.  Use new HAVE_* defines.
        * m4/rrspacing.m4: Ditto.
        * Makefile.am: Handle new files.
        * configure: Regenerated.
        * Makefile.in: Ditto.
        * config.h.in: Ditto.
        * generated/spacing_r4.c: Generated.
        * generated/spacing_r8.c: Ditto.
        * generated/spacing_r10.c: Ditto.
        * generated/spacing_r16.c: Ditto.
        * generated/rrspacing_r4.c: Ditto.
        * generated/rrspacing_r8.c: Ditto.
        * generated/rrspacing_r10.c: Ditto.
        * generated/rrspacing_r16.c: Ditto.


Added:
    trunk/libgfortran/generated/rrspacing_r10.c
    trunk/libgfortran/generated/rrspacing_r16.c
    trunk/libgfortran/generated/rrspacing_r4.c
    trunk/libgfortran/generated/rrspacing_r8.c
    trunk/libgfortran/generated/spacing_r10.c
    trunk/libgfortran/generated/spacing_r16.c
    trunk/libgfortran/generated/spacing_r4.c
    trunk/libgfortran/generated/spacing_r8.c
    trunk/libgfortran/m4/rrspacing.m4
    trunk/libgfortran/m4/spacing.m4
Modified:
    trunk/gcc/fortran/ChangeLog
    trunk/gcc/fortran/arith.c
    trunk/gcc/fortran/f95-lang.c
    trunk/gcc/fortran/gfortran.h
    trunk/gcc/fortran/iresolve.c
    trunk/gcc/fortran/simplify.c
    trunk/gcc/fortran/trans-intrinsic.c
    trunk/libgfortran/ChangeLog
    trunk/libgfortran/Makefile.am
    trunk/libgfortran/Makefile.in
    trunk/libgfortran/config.h.in
    trunk/libgfortran/configure
    trunk/libgfortran/configure.ac


-- 


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

[Bug fortran/15441] RRSPACING broken for denormals

[kargl at gcc dot gnu dot org]

------- Comment #9 from kargl at gcc dot gnu dot org  2006-10-06 00:33 -------
I have a patch, but it requires libm to have ldexp{f,l}.


-- 

kargl at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
         AssignedTo|unassigned at gcc dot gnu   |kargl at gcc dot gnu dot org
                   |dot org                     |
             Status|REOPENED                    |ASSIGNED
   Last reconfirmed|2004-05-14 17:32:55         |2006-10-06 00:33:10
               date|                            |


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

[Bug fortran/15441] RRSPACING broken for denormals

[kargl at gcc dot gnu dot org]

------- Comment #8 from kargl at gcc dot gnu dot org  2006-10-04 17:20 -------
gfortran still gets rrspacing of subnormal wrong.

troutmask:sgk[215] ./z
 s               =  4.4081038E-39  <-- Subnormal named constant
 x               =  4.4081038E-39  <-- Subnormal variable
 rrspacing(s)    =  1.2582912E+07  <-- Constant folding
 my_rrspacing(x) =  1.2582912E+07  <-- Fortran code
 rrspacing(x)    =            NaN  <-- gfortran inline

program h
  real, parameter :: s = 0.375 * tiny(1.)
  real x
  x = s
  print *, 's               =', s, ' <-- Subnormal named constant'
  print *, 'x               =', x, ' <-- Subnormal variable'
  print *, 'rrspacing(s)    =', rrspacing(s), ' <-- Constant folding'
  print *, 'my_rrspacing(x) =', my_rrspacing(x), ' <-- Fortran code'
  print *, 'rrspacing(x)    =', rrspacing(x),   ' <-- gfortran inline'
  contains
    function my_rrspacing(x)
    implicit none
    real my_rrspacing
    real, intent(in) :: x
    real y
    y = abs(x)
    if (y == 0.) then
       my_rrspacing = 0.
    else
       my_rrspacing = scale(scale(y, - exponent(y)), digits(y))
    end if
    end function my_rrspacing
end program h

The F95 standard says that rrspacing(x) = | x * b**(-e) | 2**p where b, e,
and p are the radix, exponent of x, and prcision.

trans-intrinsics.c (gfc_conv_intrinsic_rrspacing) appears to mutilate the 
mathematical expression.


-- 

kargl at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|RESOLVED                    |REOPENED
         Resolution|INVALID                     |


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