[Bug fortran/109186] New: nearest(huge(x),-1.0_kind(x)) half of correct value

[Bug fortran/109186] New: nearest(huge(x),-1.0_kind(x)) half of correct value

[john.harper at vuw dot ac.nz]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

            Bug ID: 109186
           Summary: nearest(huge(x),-1.0_kind(x)) half of correct value
           Product: gcc
           Version: 12.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: fortran
          Assignee: unassigned at gcc dot gnu.org
          Reporter: john.harper at vuw dot ac.nz
  Target Milestone: ---

With all four real kinds available in gfortran_12 nearest(x,y) with x = huge
and y = -1 this program gave about half the correct result. Ifort gave correct
results for the 3 different kinds it has: sp,dp,qp.

The program:

  implicit none
  integer,parameter::sp = kind(1.0), dp=kind(1d0), ncases = 2, &
       ep = selected_real_kind(18),qp = selected_real_kind(33)
  character(*),parameter:: fmt(4) = &
       "(A,ES"//["16.08E2","26.17E3","30.20E4","44.34E4"]//")"
  character(24),dimension(2):: splist,dplist,eplist,qplist 
  integer n
  real(sp),parameter::sp1 = 1.0_sp
  real(sp):: spx(ncases) = [huge(sp1),nearest(huge(sp1),-sp1)]
  data splist/'huge(sp1)','nearest(huge(sp1),-sp1)'/
  real(dp),parameter::dp1 = 1.0_dp
  real(dp):: dpx(ncases) = [huge(dp1),nearest(huge(dp1),-dp1)]
  data dplist/'huge(dp1)','nearest(huge(dp1),-dp1)'/
  real(ep),parameter::ep1 = 1.0_ep
  real(ep):: epx(ncases) = [huge(ep1),nearest(huge(ep1),-ep1)]
  data eplist/'huge(ep1)','nearest(huge(ep1),-ep1)'/
  real(qp),parameter::qp1 = 1.0_qp
  real(qp):: qpx(ncases) = [huge(qp1),nearest(huge(qp1),-qp1)]
  data qplist/'huge(qp1)','nearest(huge(qp1),-qp1)'/

  print fmt(1),(splist(n),spx(n),n=1,ncases)
  print fmt(2),(dplist(n),dpx(n),n=1,ncases)
  print fmt(3),(eplist(n),epx(n),n=1,ncases)
  print fmt(4),(qplist(n),qpx(n),n=1,ncases)
end program

The output on ubuntu (Linux 86_64)
huge(sp1)                 3.40282347E+38
nearest(huge(sp1),-sp1)   1.70141173E+38
huge(dp1)                 1.79769313486231571E+308
nearest(huge(dp1),-dp1)   8.98846567431157854E+307
huge(ep1)                 1.18973149535723176502E+4932
nearest(huge(ep1),-ep1)   5.94865747678615882511E+4931
huge(qp1)                 1.1897314953572317650857593266280070E+4932
nearest(huge(qp1),-qp1)   5.9486574767861588254287966331400351E+4931

[Bug fortran/109186] nearest(huge(x),-1.0_kind(x)) half of correct value

[jvdelisle at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

Jerry DeLisle <jvdelisle at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jvdelisle at gcc dot gnu.org

--- Comment #1 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> ---
Can you clarify what you expect the correct results should be that you expect?

[Bug fortran/109186] nearest(huge(x),-1.0_kind(x)) half of correct value

[john.harper at vuw dot ac.nz]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

--- Comment #2 from john.harper at vuw dot ac.nz ---
The results I expected were that

nearest(huge(sp1),-sp1) would be near huge(sp1)/(1+epsilon(sp1)),
nearest(huge(dp1),-dp1) would be near huge(dp1)/(1+epsilon(dp1)),
nearest(huge(ep1),-ep1) would be near huge(ep1)/(1+epsilon(ep1)),
nearest(huge(qp1),-qp1) would be near huge(qp1)/(1+epsilon(qp1)).

The results from gfortran were different from those by a factor near 2.
On compiling with ifort I got results close to what I had expected, which
is why I thought the bug was not in my program. The g95 compiler (now 
about 10 years old) gave results identical to those of ifort except
for ep = selected_real_kind(18) which ifort treats as if it were 
selected_real_kind(33) because ifort has only three real kinds.

By the way I thought I was using gfortran 12.0. Compiling with -v said
GNU Fortran (Ubuntu 12.1.0-2ubuntu1~22.04) version 12.1.0 
(x86_64-linux-gnu)
I wish to apologise for that error in my bug report. I got the same output
on another machine with GNU Fortran (GCC) version 12.2.0 
(x86_64-pc-linux-gnu)

The factor 2 discrepancies become clearer in this version of my program,
which has 4 more lines printing huge/nearest for the 4 real kinds.

program nearesttest3
   implicit none
   integer,parameter::sp = kind(1.0), dp=kind(1d0), ncases = 2, &
        ep = selected_real_kind(18),qp = selected_real_kind(33)
   character(*),parameter:: fmt(4) = &
        "(A,ES"//["16.08E2","26.17E3","30.20E4","44.34E4"]//")"
   character(24),dimension(2):: splist,dplist,eplist,qplist
   integer n
   real(sp),parameter::sp1 = 1.0_sp
   real(sp):: spx(ncases) = [huge(sp1),nearest(huge(sp1),-sp1)]
   data splist/'huge(sp1)','nearest(huge(sp1),-sp1)'/
   real(dp),parameter::dp1 = 1.0_dp
   real(dp):: dpx(ncases) = [huge(dp1),nearest(huge(dp1),-dp1)]
   data dplist/'huge(dp1)','nearest(huge(dp1),-dp1)'/
   real(ep),parameter::ep1 = 1.0_ep
   real(ep):: epx(ncases) = [huge(ep1),nearest(huge(ep1),-ep1)]
   data eplist/'huge(ep1)','nearest(huge(ep1),-ep1)'/
   real(qp),parameter::qp1 = 1.0_qp
   real(qp):: qpx(ncases) = [huge(qp1),nearest(huge(qp1),-qp1)]
   data qplist/'huge(qp1)','nearest(huge(qp1),-qp1)'/

   print fmt(1),(splist(n),spx(n),n=1,ncases)
   print fmt(1),'huge/nearest = ',huge(sp1)/nearest(huge(sp1),-sp1)
   print fmt(2),(dplist(n),dpx(n),n=1,ncases)
   print fmt(2),'huge/nearest = ',huge(dp1)/nearest(huge(dp1),-dp1)
   print fmt(3),(eplist(n),epx(n),n=1,ncases)
   print fmt(3),'huge/nearest = ',huge(ep1)/nearest(huge(ep1),-ep1)
   print fmt(4),(qplist(n),qpx(n),n=1,ncases)
   print fmt(4),'huge/nearest = ',huge(qp1)/nearest(huge(qp1),-qp1)
end program nearesttest3

Its gfortran output:

huge(sp1)                 3.40282347E+38
nearest(huge(sp1),-sp1)   1.70141173E+38
huge/nearest =   2.00000000E+00
huge(dp1)                 1.79769313486231571E+308
nearest(huge(dp1),-dp1)   8.98846567431157854E+307
huge/nearest =   2.00000000000000000E+000
huge(ep1)                 1.18973149535723176502E+4932
nearest(huge(ep1),-ep1)   5.94865747678615882511E+4931
huge/nearest =   2.00000000000000000000E+0000
huge(qp1)                 1.1897314953572317650857593266280070E+4932
nearest(huge(qp1),-qp1)   5.9486574767861588254287966331400351E+4931
huge/nearest =   2.0000000000000000000000000000000000E+0000

Its g95 output:

uge(sp1)                 3.40282347E+38
nearest(huge(sp1),-sp1)   3.40282326E+38
huge/nearest =   1.00000012E+00
huge(dp1)                 1.79769313486231571E+308
nearest(huge(dp1),-dp1)   1.79769313486231551E+308
huge/nearest =   1.00000000000000022E+000
huge(ep1)                 1.18973149535723176502E+4932
nearest(huge(ep1),-ep1)   1.18973149535723176496E+4932
huge/nearest =   1.00000000000000000011E+0000
huge(qp1)                 1.1897314953572317650857593266280070E+4932
nearest(huge(qp1),-qp1)   1.1897314953572317650857593266280069E+4932
huge/nearest =   1.0000000000000000000000000000000002E+0000

On Sun, 19 Mar 2023, jvdelisle at gcc dot gnu.org wrote:

> Date: Sun, 19 Mar 2023 05:06:47 +0000
> From: jvdelisle at gcc dot gnu.org <gcc-bugzilla@gcc.gnu.org>
> To: John Harper <john.harper@vuw.ac.nz>
> Subject: [Bug fortran/109186] nearest(huge(x),
>     -1.0_kind(x)) half of correct value
> Resent-Date: Sun, 19 Mar 2023 18:06:59 +1300 (NZDT)
> Resent-From: <john.harper@vuw.ac.nz>
> 
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186
>
> Jerry DeLisle <jvdelisle at gcc dot gnu.org> changed:
>
>           What    |Removed                     |Added
> ----------------------------------------------------------------------------
>                 CC|                            |jvdelisle at gcc dot gnu.org
>
> --- Comment #1 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> ---
> Can you clarify what you expect the correct results should be that you expect?
>
> -- 
> You are receiving this mail because:
> You reported the bug.
>


-- John Harper, School of Mathematics and Statistics
Victoria Univ. of Wellington, PO Box 600, Wellington 6140, New Zealand.
e-mail john.harper@vuw.ac.nz phone +64(0) 4 463 5276

[Bug fortran/109186] nearest(huge(x),-1.0_kind(x)) half of correct value

[anlauf at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

anlauf at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |anlauf at gcc dot gnu.org
             Status|UNCONFIRMED                 |NEW
     Ever confirmed|0                           |1
   Last reconfirmed|                            |2023-03-19
           Keywords|                            |wrong-code

--- Comment #3 from anlauf at gcc dot gnu.org ---
Confirmed.

This is a compile-time issue.  When the first argument of NEAREST is a
variable and not a constant expression, I get the correct result:

program p
  real, parameter :: x = huge(1.0)
  real :: y, z
  y = x
  z = nearest (y, -1.0)
  print *, nearest (x, -1.0), z
end

This prints:

   1.70141173E+38   3.40282326E+38

Seems to affect all versions at least down to 7.

[Bug fortran/109186] nearest(huge(x),-1.0_kind(x)) half of correct value

[anlauf at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

--- Comment #4 from anlauf at gcc dot gnu.org ---
The following patch fixes this:

diff --git a/gcc/fortran/simplify.cc b/gcc/fortran/simplify.cc
index 20ea38e0007..35ae637a483 100644
--- a/gcc/fortran/simplify.cc
+++ b/gcc/fortran/simplify.cc
@@ -6114,7 +6112,7 @@ gfc_simplify_nearest (gfc_expr *x, gfc_expr *s)
   kind = gfc_validate_kind (BT_REAL, x->ts.kind, 0);
   mpfr_set_emin ((mpfr_exp_t) gfc_real_kinds[kind].min_exponent -
                mpfr_get_prec(result->value.real) + 1);
-  mpfr_set_emax ((mpfr_exp_t) gfc_real_kinds[kind].max_exponent - 1);
+  mpfr_set_emax ((mpfr_exp_t) gfc_real_kinds[kind].max_exponent);
   mpfr_check_range (result->value.real, 0, MPFR_RNDU);

   if (mpfr_sgn (s->value.real) > 0)

The above code was apparently introduced with r0-84566-gb6f63e898498e6

[Bug fortran/109186] nearest(huge(x),-1.0_kind(x)) half of correct value

[anlauf at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

anlauf at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           Assignee|unassigned at gcc dot gnu.org      |anlauf at gcc dot gnu.org
             Status|NEW                         |ASSIGNED

--- Comment #5 from anlauf at gcc dot gnu.org ---
Submitted: https://gcc.gnu.org/pipermail/fortran/2023-March/059099.html

[Bug fortran/109186] nearest(huge(x),-1.0_kind(x)) half of correct value

[cvs-commit at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

--- Comment #6 from CVS Commits <cvs-commit at gcc dot gnu.org> ---
The master branch has been updated by Harald Anlauf <anlauf@gcc.gnu.org>:

https://gcc.gnu.org/g:4410a08b80cc40342eeaa5b6af824cd4352b218c

commit r13-6762-g4410a08b80cc40342eeaa5b6af824cd4352b218c
Author: Harald Anlauf <anlauf@gmx.de>
Date:   Sun Mar 19 21:29:46 2023 +0100

    Fortran: simplification of NEAREST for large argument [PR109186]

    gcc/fortran/ChangeLog:

            PR fortran/109186
            * simplify.cc (gfc_simplify_nearest): Fix off-by-one error in
setting
            up real kind-specific maximum exponent for mpfr.

    gcc/testsuite/ChangeLog:

            PR fortran/109186
            * gfortran.dg/nearest_6.f90: New test.

[Bug fortran/109186] nearest(huge(x),-1.0_kind(x)) half of correct value

[cvs-commit at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

--- Comment #7 from CVS Commits <cvs-commit at gcc dot gnu.org> ---
The releases/gcc-12 branch has been updated by Harald Anlauf
<anlauf@gcc.gnu.org>:

https://gcc.gnu.org/g:f1342c21b2acbdd1daa0f2519154f8ac911a800f

commit r12-9310-gf1342c21b2acbdd1daa0f2519154f8ac911a800f
Author: Harald Anlauf <anlauf@gmx.de>
Date:   Sun Mar 19 21:29:46 2023 +0100

    Fortran: simplification of NEAREST for large argument [PR109186]

    gcc/fortran/ChangeLog:

            PR fortran/109186
            * simplify.cc (gfc_simplify_nearest): Fix off-by-one error in
setting
            up real kind-specific maximum exponent for mpfr.

    gcc/testsuite/ChangeLog:

            PR fortran/109186
            * gfortran.dg/nearest_6.f90: New test.

    (cherry picked from commit 4410a08b80cc40342eeaa5b6af824cd4352b218c)

[Bug fortran/109186] nearest(huge(x),-1.0_kind(x)) half of correct value

[cvs-commit at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

--- Comment #8 from CVS Commits <cvs-commit at gcc dot gnu.org> ---
The releases/gcc-11 branch has been updated by Harald Anlauf
<anlauf@gcc.gnu.org>:

https://gcc.gnu.org/g:4b722050fe8e7f06a3e653f290a39ccdc4d03174

commit r11-10595-g4b722050fe8e7f06a3e653f290a39ccdc4d03174
Author: Harald Anlauf <anlauf@gmx.de>
Date:   Sun Mar 19 21:29:46 2023 +0100

    Fortran: simplification of NEAREST for large argument [PR109186]

    gcc/fortran/ChangeLog:

            PR fortran/109186
            * simplify.c (gfc_simplify_nearest): Fix off-by-one error in
setting
            up real kind-specific maximum exponent for mpfr.

    gcc/testsuite/ChangeLog:

            PR fortran/109186
            * gfortran.dg/nearest_6.f90: New test.

    (cherry picked from commit 4410a08b80cc40342eeaa5b6af824cd4352b218c)

[Bug fortran/109186] nearest(huge(x),-1.0_kind(x)) half of correct value

[cvs-commit at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

--- Comment #9 from CVS Commits <cvs-commit at gcc dot gnu.org> ---
The releases/gcc-10 branch has been updated by Harald Anlauf
<anlauf@gcc.gnu.org>:

https://gcc.gnu.org/g:b9f6e7e234bab6ef8b5425c3e8b88bf8dfbc5089

commit r10-11265-gb9f6e7e234bab6ef8b5425c3e8b88bf8dfbc5089
Author: Harald Anlauf <anlauf@gmx.de>
Date:   Sun Mar 19 21:29:46 2023 +0100

    Fortran: simplification of NEAREST for large argument [PR109186]

    gcc/fortran/ChangeLog:

            PR fortran/109186
            * simplify.c (gfc_simplify_nearest): Fix off-by-one error in
setting
            up real kind-specific maximum exponent for mpfr.

    gcc/testsuite/ChangeLog:

            PR fortran/109186
            * gfortran.dg/nearest_6.f90: New test.

    (cherry picked from commit 4410a08b80cc40342eeaa5b6af824cd4352b218c)

[Bug fortran/109186] nearest(huge(x),-1.0_kind(x)) half of correct value

[anlauf at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=109186

anlauf at gcc dot gnu.org changed:

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

--- Comment #10 from anlauf at gcc dot gnu.org ---
Fixed on all open branches.  Closing.

Thanks for the report!