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.
Confirmed.
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.
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
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.
Created attachment 7138 [details] More comments
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?
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 :-)
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.
I have a patch, but it requires libm to have ldexp{f,l}.
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
Fixed on trunk (until someone tells me ldexp doesn't exist).