Bug 15441 - RRSPACING broken for denormals
Summary: RRSPACING broken for denormals
Alias: None
Product: gcc
Classification: Unclassified
Component: fortran (show other bugs)
Version: 4.0.0
: P2 enhancement
Target Milestone: 4.2.0
Assignee: kargl
Keywords: wrong-code
Depends on:
Reported: 2004-05-14 17:30 UTC by Tobias Schlüter
Modified: 2006-10-09 20:57 UTC (History)
2 users (show)

See Also:
Known to work:
Known to fail:
Last reconfirmed: 2006-10-06 00:33:10

More comments (1.02 KB, patch)
2004-09-15 09:39 UTC, Feng Wang
Details | Diff

Note You need to log in before you can comment on or make changes to this bug.
Description Tobias Schlüter 2004-05-14 17:30:38 UTC
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.
Comment 1 Andrew Pinski 2004-05-14 17:32:55 UTC
Comment 2 Feng Wang 2004-09-07 09:56:07 UTC
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.
Comment 3 Tobias Schlüter 2004-09-13 17:36:41 UTC
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
Comment 4 Tobias Schlüter 2004-09-14 12:59:46 UTC
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
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.
Comment 5 Feng Wang 2004-09-15 09:39:10 UTC
Created attachment 7138 [details]
More comments
Comment 6 Feng Wang 2004-09-15 09:40:01 UTC
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 

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?
Comment 7 Tobias Schlüter 2004-09-15 10:24:35 UTC
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 :-)
Comment 8 kargl 2006-10-04 17:20:58 UTC
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'
    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.
       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.
Comment 9 kargl 2006-10-06 00:33:10 UTC
I have a patch, but it requires libm to have ldexp{f,l}.
Comment 10 kargl 2006-10-09 20:55:40 UTC
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
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.


Comment 11 kargl 2006-10-09 20:57:22 UTC
Fixed on trunk (until someone tells me ldexp doesn't exist).