GCC Bugzilla will be upgraded from version 4.4.9 to 5.0rc3 on Saturday, April 25, starting around 17:00 UTC. The upgrade process should only last a few minutes. Check bug 64968 for details.
Bug 34192 - [4.2, 4.3 Regression] NEAREST can return wrong numbers
[4.2, 4.3 Regression] NEAREST can return wrong numbers
Status: RESOLVED FIXED
Product: gcc
Classification: Unclassified
Component: fortran
4.3.0
: P3 normal
: ---
Assigned To: Not yet assigned to anyone
: rejects-valid, wrong-code
: 34194 (view as bug list)
Depends on:
Blocks: 32834
  Show dependency treegraph
 
Reported: 2007-11-22 13:44 UTC by francois.jacq
Modified: 2007-11-23 21:50 UTC (History)
2 users (show)

See Also:
Host:
Target:
Build:
Known to work: 4.1.3
Known to fail: 4.2.2 4.3.0
Last reconfirmed: 2007-11-22 16:27:20


Attachments
Patch (reverted Rev 117584) plus test case (1.99 KB, text/plain)
2007-11-22 18:41 UTC, Tobias Burnus
Details

Note You need to log in before you can comment on or make changes to this bug.
Description francois.jacq 2007-11-22 13:44:48 UTC
PROGRAM test 
  WRITE(*,*) NEAREST(0.d0,1.d0)
END

>> gfortran test.f90
t.f90:3.21:

  WRITE(*,*) NEAREST(0.d0,1.d0)
                    1
Error: Result of NEAREST underflows its kind at (1)
Comment 1 Tobias Burnus 2007-11-22 15:53:45 UTC
*** Bug 34194 has been marked as a duplicate of this bug. ***
Comment 2 Tobias Burnus 2007-11-22 16:11:36 UTC
I believe there is indeed a bug in gfortran. The standard says that the "the nearest different machine-representable number" is returned.

However, gfortran returns a number which is not representable (denormal or smaller?) and rounded to zero (= same number).

Result (gfortran 4.3.0):
a) The program is wrongly rejected
b) With -fno-range-check, the wrong number is used. For the example the dump shows:
   static real(kind=8) C.859 = 0.0;

The expected result is some denormal number like:
5.E-324, 4.9406564584124654E-324, 2.225073858507201E-308, 4.94065645841246544E-324, 4.940656458412465E-324.
(The ifort result of 2e-308 looks wrong in this regard as there is a nearer number.)

gfortran 4.1.3 returns: 4.940656458412465E-324
gfortran 4.2.2 rejects it and has an ICE with -fno-range-check


Besides fixing this bug, we should also consider mentioning -fno-range-check in the error message. (Though, in this case -fno-range-check would have concealed the real bug.)


The Fortran 2003 standard states:

"NEAREST (X, S)
Description. Returns the nearest different machine-representable number in a given direction.
  X shall be of type real.
  S shall be of type real and not equal to zero.
Result Characteristics. Same as X.
Result Value. The result has a value equal to the machine-representable number distinct from X and nearest to it in the direction of the infinity with the same sign as S.

NOTE 13.16 Unlike other floating-point manipulation functions, NEAREST operates on machine-representable numbers rather than model numbers. On many systems there are machine-representable numbers that lie between adjacent model numbers."
Comment 3 Tobias Burnus 2007-11-22 16:27:20 UTC
Seemingly the problem was introduced by the following patch, which fixed another NEAREST problem.

http://gcc.gnu.org/ml/gcc-cvs/2006-10/msg00249.html
http://gcc.gnu.org/viewcvs?view=rev&revision=117584

2006-10-09  Steven G. Kargl  <kargl@gcc.gnu.org>
[...]
        PR fortran/15441
        PR fortran/29312
[...]
        * simplify.c (gfc_simplify_nearest): Remove explicit subnormalization.
Comment 4 Tobias Burnus 2007-11-22 17:08:01 UTC
Regarding the range check: We need to disable the check for denormal numbers; as overflow etc. cannot occur, this boils down to check only for NaN.


Regarding the result, I think we have a problem with MPFR. The following program prints:

Result (GMP_RNDN): 0
Result (GMP_RNDU): 4.94066e-324

As mpfr_nexttoward should return a representable number, there should be no rounding needed, or should it? By the way, changing GFC_REAL_8_DIGITS did not change the result.

See also:
http://www.mpfr.org/mpfr-2.3.0/mpfr.html#index-mpfr_005fnexttoward-224


#include <gmp.h>
#include <mpfr.h>
#include <stdio.h>
#define GFC_REAL_8_DIGITS 53

int main()
{
  mpfr_t x;
  mpfr_t y;
  int base2prec;
  double result;

  mpfr_set_default_prec (GFC_REAL_8_DIGITS);
  mpfr_init(x);
  mpfr_init(y);
  mpfr_set_d (x, 0.0, GMP_RNDN); /* x = 0.0         */
  mpfr_set_inf (y, 1);           /* y = INF         */
  mpfr_nexttoward (x, y);        /* x = x + epsilon */
  printf("Result (GMP_RNDN): %g\n", mpfr_get_d(x, GMP_RNDN));
  printf("Result (GMP_RNDU): %g\n", mpfr_get_d(x, GMP_RNDU));
  return 0;
}
Comment 5 Tobias Burnus 2007-11-22 18:05:57 UTC
Hmm, this does not seem to be easily possible in MPFR.
  http://websympa.loria.fr/wwsympa/arc/mpfr/2007-11/msg00026.html

Possible implementation (as of 2005) by Tobias Schlüter:
  http://gcc.gnu.org/ml/fortran/2005-04/msg00389.html

Steve writes to my MPFR email of today:
| This is a point I had not fully appreciated when I patched
| gfc_simplify_nearest.  Try restoring the code that I removed. 
| I don't remember why I removed the code, although I suspect
| that range_check function should properly handle subnormal
| numbers.
Comment 6 Tobias Burnus 2007-11-22 18:41:28 UTC
Created attachment 14609 [details]
Patch (reverted Rev 117584) plus test case

Initial implementation by reverting Rev. 117584 as suggested by Steve.

Problems:
- nearest(x,-1) for negative denormal x does not give a smaller number
- returned values are not the smallest possible:
  4.9406564584124654E-324 vs. 1.11253692925360069E-308
  1.4012985E-45 vs. 5.87747175E-39
- NaN, INF, -INF test cases missing.

Probably makes sense to check this in first and worry about the rest later.
Comment 7 kargl 2007-11-22 19:00:55 UTC
(In reply to comment #6)
> Created an attachment (id=14609) [edit]
> Patch (reverted Rev 117584) plus test case
> 
> Initial implementation by reverting Rev. 117584 as suggested by Steve.
> 
> Problems:
> - nearest(x,-1) for negative denormal x does not give a smaller number
> - returned values are not the smallest possible:
>   4.9406564584124654E-324 vs. 1.11253692925360069E-308
>   1.4012985E-45 vs. 5.87747175E-39

To fix this problem, I believe that you need to set emin to emin - prec + 1

Something like

+  mpfr_set_emin ((mp_exp_t) (gfc_real_kinds[sgn].min_exponent - 
                              gfc_real_kinds[sgn].digit + 1));

It might be a -1.
 
IIRC, mpfr always stores a normalized significand, so you need to
play with emin to properly account for gradual underflow.

Comment 8 patchapp@dberlin.org 2007-11-23 20:05:52 UTC
Subject: Bug number PR 34192

A patch for this bug has been added to the patch tracker.
The mailing list url for the patch is http://gcc.gnu.org/ml/gcc-patches/2007-11/msg01240.html
Comment 9 Tobias Burnus 2007-11-23 21:03:58 UTC
Subject: Bug 34192

Author: burnus
Date: Fri Nov 23 21:03:48 2007
New Revision: 130383

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=130383
Log:
2007-11-23  Tobias Burnus  <burnus@net-b.de>
            Steven G. Kargl  <kargl@gcc.gnu.org>

        PR fortran/34192
        * simplify.c (gfc_simplify_nearest): Fix NEAREST for
        subnormal numbers.

2007-11-23  Tobias Burnus  <burnus@net-b.de>

        PR fortran/34192
        * gfortran.dg/nearest_2.f90: New.


Added:
    trunk/gcc/testsuite/gfortran.dg/nearest_2.f90
Modified:
    trunk/gcc/fortran/ChangeLog
    trunk/gcc/fortran/simplify.c
    trunk/gcc/testsuite/ChangeLog

Comment 10 Tobias Burnus 2007-11-23 21:50:54 UTC
FIXED on the trunk (4.3.0). I do not indent to backport it to GCC 4.2.x; but I can be convince otherwise.
Comment 11 Tobias Burnus 2007-11-24 13:18:36 UTC
Subject: Bug 34192

Author: burnus
Date: Sat Nov 24 13:18:27 2007
New Revision: 130396

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=130396
Log:
2007-11-24  Tobias Burnus  <burnus@net-b.de>

       PR fortran/34192
       * gfortran.dg/nearest_2.f90: Add INF/NAN tests.


Modified:
    trunk/gcc/testsuite/ChangeLog
    trunk/gcc/testsuite/gfortran.dg/nearest_2.f90