Hi *, there is a nasty bug with EXPONENT(): EXPONENT(1.0) evaluates to 0 (which is wrong), while real :: a=1.0; EXPONENT(a) evaluates to 1 (which is correct). For more details and a further example see the attached code sample. Cheers, -ha
Created attachment 11838 [details] Example code
Confirmed. I have a patch, but it may depend on MPFR 2.2.0.
(CC myself) I see the same problem for real(8) and for real(10) as for real = real(4), except for print *, scale (fraction (a), exponent (a)) / a where I get "NaN" with real(10) instead of the correct "1.0", which one gets for real and real(8).
Subject: Re: EXPONENT() broken for real constants On Mon, Jul 17, 2006 at 05:54:46PM -0000, tobias dot burnus at physik dot fu-berlin dot de wrote: > > I see the same problem for real(8) and for real(10) as for real = real(4), > except for > print *, scale (fraction (a), exponent (a)) / a > where I get "NaN" with real(10) instead of the correct "1.0", which one gets > for real and real(8). > Can you try this patch? It required MPFR 2.2.0. Index: gcc4x/gcc/fortran/simplify.c =================================================================== --- gcc4x/gcc/fortran/simplify.c (revision 115505) +++ gcc4x/gcc/fortran/simplify.c (working copy) @@ -1049,12 +1049,10 @@ gfc_simplify_exp (gfc_expr * x) return range_check (result, "EXP"); } -/* FIXME: MPFR should be able to do this better */ gfc_expr * gfc_simplify_exponent (gfc_expr * x) { int i; - mpfr_t tmp; gfc_expr *result; if (x->expr_type != EXPR_CONSTANT) @@ -1071,20 +1069,9 @@ gfc_simplify_exponent (gfc_expr * x) return result; } - mpfr_init (tmp); - - mpfr_abs (tmp, x->value.real, GFC_RND_MODE); - mpfr_log2 (tmp, tmp, GFC_RND_MODE); - - gfc_mpfr_to_mpz (result->value.integer, tmp); - - /* The model number for tiny(x) is b**(emin - 1) where b is the base and emin - is the smallest exponent value. So, we need to add 1 if x is tiny(x). */ - i = gfc_validate_kind (x->ts.type, x->ts.kind, false); - if (mpfr_cmp (x->value.real, gfc_real_kinds[i].tiny) == 0) - mpz_add_ui (result->value.integer,result->value.integer, 1); - - mpfr_clear (tmp); + /* Requires MPFR 2.2.0 or newer. */ + i = (int) mpfr_get_exp (x->value.real); + mpz_set_si (result->value.integer, i); return range_check (result, "EXPONENT"); }
> Can you try this patch? It required MPFR 2.2.0. Tested on openSUSE/AMD64 with current gfortran trunk + patch and mpfr-2.2.0-9. It works ok with kind = 4 and kind = 8. However, with kind = 10 I get a NaN for print *, scale (fraction (a), exponent (a)) / a My test case: ---------------- integer, parameter :: k = 10 real(k), parameter :: one = 1.0_k real(k) :: a a = one print *, "This ratio should be 1:" print *, scale (fraction (a), exponent (a)) / a end
Subject: Re: EXPONENT() broken for real constants On Sun, Sep 10, 2006 at 10:43:21PM -0000, tobias dot burnus at physik dot fu-berlin dot de wrote: > > Comment #5 from tobias dot burnus at physik dot fu-berlin dot de 2006-09-10 22:43 ------- > > Can you try this patch? It required MPFR 2.2.0. > Tested on openSUSE/AMD64 with current gfortran trunk + patch and mpfr-2.2.0-9. > > It works ok with kind = 4 and kind = 8. However, with kind = 10 I get a > NaN for > print *, scale (fraction (a), exponent (a)) / a > > My test case: > ---------------- > integer, parameter :: k = 10 > real(k), parameter :: one = 1.0_k > real(k) :: a > a = one > print *, "This ratio should be 1:" > print *, scale (fraction (a), exponent (a)) / a > end > > gfortran is calling the wrong libm function. -fdump-tree-original on FreeBSD shows { real10 D.1254; D.1254 = scalbn(_gfortran_fraction_r10(a10),_gfortran_exponent_r10(a10))/a10; _gfortran_transfer_real (&dt_parm.2, &D.1254, 10); } That should be scalbnl(). I'll see if I can track down the bug.
Subject: Re: EXPONENT() broken for real constants On Mon, Sep 11, 2006 at 07:19:41PM -0000, sgk at troutmask dot apl dot washington dot edu wrote: > > gfortran is calling the wrong libm function. -fdump-tree-original > on FreeBSD shows > > { > real10 D.1254; > > D.1254 = scalbn(_gfortran_fraction_r10(a10),_gfortran_exponent_r10(a10))/a10; > _gfortran_transfer_real (&dt_parm.2, &D.1254, 10); > } > > That should be scalbnl(). I'll see if I can track down the bug. > I should also note that the parsing is correct as shown by -fdump-parse-tree ASSIGN MAIN__:a10 1.00000000000000000000_10 ASSIGN MAIN__:a10 (/ __scale_10[[((__fraction_10[[((MAIN__:a10))]]) (__exponent_10[[((MAIN__:a10))]]))]] MAIN__:a10) So it's the translation of __scale_10 to scalbnl.
Subject: Re: EXPONENT() broken for real constants On Mon, Sep 11, 2006 at 07:33:35PM -0000, sgk at troutmask dot apl dot washington dot edu wrote: > > I should also note that the parsing is correct as shown by -fdump-parse-tree > > > ASSIGN MAIN__:a10 1.00000000000000000000_10 > ASSIGN MAIN__:a10 (/ __scale_10[[((__fraction_10[[((MAIN__:a10))]]) > (__exponent_10[[((MAIN__:a10))]]))]] MAIN__:a10) > > So it's the translation of __scale_10 to scalbnl. > I have a patch for this problem.
Subject: Bug 28276 Author: kargl Date: Wed Sep 27 20:15:22 2006 New Revision: 117257 URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=117257 Log: * configure.in: Check for GMP 4.1 or newer. Check for MPFR 2.2.0 or newer. * configure: Regenerated. * doc/install.texi: Document required versions of GMP and MPFR. * fortran/arith.c: Conditionally include arctangent2(). (gfc_check_real_range): Use mpfr_subnormalize in preference to local hack. * fortran/trans-intrinsic.c (gfc_get_intrinsic_lib_fndecl): Append l for long double functions. * fortran/simplify.c: Wrap Copyright to new line. (gfc_simplify_atan2): Use mpfr_atan2 in preference to arctangent2(). (gfc_simplify_log): Ditto. PR fortran/28276 * fortran/simplify.c (gfc_simplify_exponent): Use mpfr_get_exp in preference to broken local hack. PR fortran/27021 * fortran/simplify.c (gfc_simplify_nearest): Use mpfr_nexttoward and mpfr_subnormalize to handle numbers near zero in preference to broken local hack. PR fortran/28276 * testsuite/gfortran.dg/exponent_1.f90: New test. PR fortran/27021 * testsuite/gfortran.dg/nearest_1.f90: New test. Added: trunk/gcc/testsuite/gfortran.dg/exponent_1.f90 trunk/gcc/testsuite/gfortran.dg/nearest_1.f90 Modified: trunk/ChangeLog trunk/configure trunk/configure.in trunk/gcc/ChangeLog trunk/gcc/doc/install.texi trunk/gcc/fortran/ChangeLog trunk/gcc/fortran/arith.c trunk/gcc/fortran/simplify.c trunk/gcc/fortran/trans-intrinsic.c trunk/gcc/testsuite/ChangeLog
Fixed on trunk. Won't be fixed on 4.1.