[PATCH] Fortran: fix compile-time simplification of SET_EXPONENT [PR109511]

[PATCH] Fortran: fix compile-time simplification of SET_EXPONENT [PR109511]

[Harald Anlauf]

[-- Attachment #1: Type: text/plain, Size: 542 bytes --]

[Now with subject added...]

Dear all,

the compile-time simplification of intrinsic SET_EXPONENT was
broken since the early days of gfortran for argument X < 1
(including negative X) and for I < 0.  I identified and fixed
several issues in the implementation.  The testcase explores
argument space comparing compile-time and runtime results and
is checked against Intel.

Regtested on x86_64-pc-linux-gnu.  OK for mainline?

This is not a regression, but can lead to wrong code.
Would it be OK to backport to open branches?

Thanks,
Harald


[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: pr109511.diff --]
[-- Type: text/x-patch, Size: 3902 bytes --]

From fa4cb42870df60deb8888dbd51e2ddc6d6ab9e6a Mon Sep 17 00:00:00 2001
From: Harald Anlauf <anlauf@gmx.de>
Date: Fri, 14 Apr 2023 20:45:19 +0200
Subject: [PATCH] Fortran: fix compile-time simplification of SET_EXPONENT
 [PR109511]

gcc/fortran/ChangeLog:

	PR fortran/109511
	* simplify.cc (gfc_simplify_set_exponent): Fix implementation of
	compile-time simplification of intrinsic SET_EXPONENT for argument
	X < 1 and for I < 0.

gcc/testsuite/ChangeLog:

	PR fortran/109511
	* gfortran.dg/set_exponent_1.f90: New test.
---
 gcc/fortran/simplify.cc                      | 12 +++----
 gcc/testsuite/gfortran.dg/set_exponent_1.f90 | 36 ++++++++++++++++++++
 2 files changed, 42 insertions(+), 6 deletions(-)
 create mode 100644 gcc/testsuite/gfortran.dg/set_exponent_1.f90

diff --git a/gcc/fortran/simplify.cc b/gcc/fortran/simplify.cc
index ecf0e3558df..b65854c1021 100644
--- a/gcc/fortran/simplify.cc
+++ b/gcc/fortran/simplify.cc
@@ -7364,7 +7364,7 @@ gfc_simplify_set_exponent (gfc_expr *x, gfc_expr *i)
 {
   gfc_expr *result;
   mpfr_t exp, absv, log2, pow2, frac;
-  unsigned long exp2;
+  long exp2;

   if (x->expr_type != EXPR_CONSTANT || i->expr_type != EXPR_CONSTANT)
     return NULL;
@@ -7396,19 +7396,19 @@ gfc_simplify_set_exponent (gfc_expr *x, gfc_expr *i)
   mpfr_abs (absv, x->value.real, GFC_RND_MODE);
   mpfr_log2 (log2, absv, GFC_RND_MODE);

-  mpfr_trunc (log2, log2);
+  mpfr_floor (log2, log2);
   mpfr_add_ui (exp, log2, 1, GFC_RND_MODE);

   /* Old exponent value, and fraction.  */
   mpfr_ui_pow (pow2, 2, exp, GFC_RND_MODE);

-  mpfr_div (frac, absv, pow2, GFC_RND_MODE);
+  mpfr_div (frac, x->value.real, pow2, GFC_RND_MODE);

   /* New exponent.  */
-  exp2 = (unsigned long) mpz_get_d (i->value.integer);
-  mpfr_mul_2exp (result->value.real, frac, exp2, GFC_RND_MODE);
+  exp2 = mpz_get_si (i->value.integer);
+  mpfr_mul_2si (result->value.real, frac, exp2, GFC_RND_MODE);

-  mpfr_clears (absv, log2, pow2, frac, NULL);
+  mpfr_clears (absv, log2, exp, pow2, frac, NULL);

   return range_check (result, "SET_EXPONENT");
 }
diff --git a/gcc/testsuite/gfortran.dg/set_exponent_1.f90 b/gcc/testsuite/gfortran.dg/set_exponent_1.f90
new file mode 100644
index 00000000000..4c063e8330b
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/set_exponent_1.f90
@@ -0,0 +1,36 @@
+! { dg-do run }
+! PR fortran/109511
+! Check compile-time simplification of SET_EXPONENT against runtime
+
+program exponent
+  implicit none
+  integer :: i
+  i = 0
+  print *, i, set_exponent(1., 0), set_exponent(1., i)
+  if (set_exponent(1., 0) /= set_exponent(1., i)) stop 1
+  i = 1
+  print *, i, set_exponent(1., 1), set_exponent(1., i)
+  if (set_exponent(1., 1) /= set_exponent(1., i)) stop 2
+  i = 2
+  print *, i, set_exponent(-1.75, 2), set_exponent(-1.75, i)
+  if (set_exponent(-1.75, 2) /= set_exponent(-1.75, i)) stop 3
+  print *, i, set_exponent(0.1875, 2), set_exponent(0.1875, i)
+  if (set_exponent(0.1875, 2) /= set_exponent(0.1875, i)) stop 4
+  i = 3
+  print *, i, set_exponent(0.75, 3), set_exponent(0.75, i)
+  if (set_exponent(0.75, 3) /= set_exponent(0.75, i)) stop 5
+  i = 4
+  print *, i, set_exponent(-2.5, 4), set_exponent(-2.5, i)
+  if (set_exponent(-2.5, 4) /= set_exponent(-2.5, i)) stop 6
+  i = -1
+  print *, i, set_exponent(1., -1), set_exponent(1., i)
+  if (set_exponent(1., -1) /= set_exponent(1., i)) stop 7
+  i = -2
+  print *, i, set_exponent(1.125, -2), set_exponent(1.125, i)
+  if (set_exponent(1.125, -2) /= set_exponent(1.125, i)) stop 8
+  print *, i, set_exponent(-0.25, -2), set_exponent(-0.25, i)
+  if (set_exponent(-0.25, -2) /= set_exponent(-0.25, i)) stop 9
+  i = -3
+  print *, i, set_exponent(0.75, -3), set_exponent(0.75, i)
+  if (set_exponent(0.75, -3) /= set_exponent(0.75, i)) stop 10
+end program exponent
--
2.35.3

Re: [PATCH] Fortran: fix compile-time simplification of SET_EXPONENT [PR109511]

[Steve Kargl]

On Fri, Apr 14, 2023 at 08:59:24PM +0200, Harald Anlauf via Fortran wrote:
> 
> the compile-time simplification of intrinsic SET_EXPONENT was
> broken since the early days of gfortran for argument X < 1
> (including negative X) and for I < 0.  I identified and fixed
> several issues in the implementation.  The testcase explores
> argument space comparing compile-time and runtime results and
> is checked against Intel.
> 
> Regtested on x86_64-pc-linux-gnu.  OK for mainline?
> 

Yes, it is certainly better than the current situation.

> This is not a regression, but can lead to wrong code.
> Would it be OK to backport to open branches?

Sure.  Looks simply and fairly specific.


I was wondering about the difference between set_exponent()
and scale(), and found that set_exponent() talks about IEEE
values while scale() doesn't.  I'm wondering if we should 
add the IEEE special cases to the testsuite.  Of particular
note, I doubt that this is true:

   If X is an IEEE NaN, the result is the same NaN.

program foo
   real x, y
   x = 1
   y = x - x
   x = (x - x) / y
   print '(F4.0,1X,Z8.8)', x, x
   y = set_exponent(x,1)
   print '(F4.0,1X,Z8.8)', y, y
end program foo

 gfcx -o z a.f90 && ./z
 NaN FFC00000
 NaN 7FC00000

Those are not the same NaN.  The second is a qNaN.
The first looks like a qNaN with the sign bit set.

-- 
Steve

Re: [PATCH] Fortran: fix compile-time simplification of SET_EXPONENT [PR109511]

[Steve Kargl]

On Fri, Apr 14, 2023 at 12:33:18PM -0700, Steve Kargl via Fortran wrote:
> 
>    If X is an IEEE NaN, the result is the same NaN.
> 

A better testcase as gfortran will quiet a NaN on assignment.

program foo
   integer i
   equivalence(i,y)
   i = int(z'7FC0BEEF',4)  ! Add payload to NaN.
   print '(F4.0,2(1X,Z8.8))', y, y, set_exponent(y,1)
end program foo

gfortran -o z a.f90 && ./z
N 7FC0BEEF 7FC00000


-- 
Steve

Re: [PATCH] Fortran: fix compile-time simplification of SET_EXPONENT [PR109511]

[Harald Anlauf]

Hi Steve,

On 4/14/23 21:33, Steve Kargl via Gcc-patches wrote:
> On Fri, Apr 14, 2023 at 08:59:24PM +0200, Harald Anlauf via Fortran wrote:
>>
>> the compile-time simplification of intrinsic SET_EXPONENT was
>> broken since the early days of gfortran for argument X < 1
>> (including negative X) and for I < 0.  I identified and fixed
>> several issues in the implementation.  The testcase explores
>> argument space comparing compile-time and runtime results and
>> is checked against Intel.
>>
>> Regtested on x86_64-pc-linux-gnu.  OK for mainline?
>>
>
> Yes, it is certainly better than the current situation.
>
>> This is not a regression, but can lead to wrong code.
>> Would it be OK to backport to open branches?
>
> Sure.  Looks simply and fairly specific.

OK, thanks, will proceed.


> I was wondering about the difference between set_exponent()
> and scale(), and found that set_exponent() talks about IEEE
> values while scale() doesn't.  I'm wondering if we should
> add the IEEE special cases to the testsuite.  Of particular
> note, I doubt that this is true:
>
>     If X is an IEEE NaN, the result is the same NaN.
>
> program foo
>     real x, y
>     x = 1
>     y = x - x
>     x = (x - x) / y
>     print '(F4.0,1X,Z8.8)', x, x
>     y = set_exponent(x,1)
>     print '(F4.0,1X,Z8.8)', y, y
> end program foo
>
>   gfcx -o z a.f90 && ./z
>   NaN FFC00000
>   NaN 7FC00000
>
> Those are not the same NaN.  The second is a qNaN.
> The first looks like a qNaN with the sign bit set.

Until now there was no testing at all of SET_EXPONENT in the testsuite.
It would be really good to have better coverage of compile-time and
runtime behavior of the intrinsics and checking consistency ... ;-)

I think you have much more experience in that area.  (Hint!)

Cheers,
Harald

Re: [PATCH] Fortran: fix compile-time simplification of SET_EXPONENT [PR109511]

[Harald Anlauf]

Hi Steve,

On 4/14/23 21:33, Steve Kargl via Gcc-patches wrote:
> On Fri, Apr 14, 2023 at 08:59:24PM +0200, Harald Anlauf via Fortran wrote:
>>
>> the compile-time simplification of intrinsic SET_EXPONENT was
>> broken since the early days of gfortran for argument X < 1
>> (including negative X) and for I < 0.  I identified and fixed
>> several issues in the implementation.  The testcase explores
>> argument space comparing compile-time and runtime results and
>> is checked against Intel.
>>
>> Regtested on x86_64-pc-linux-gnu.  OK for mainline?
>>
> 
> Yes, it is certainly better than the current situation.
> 
>> This is not a regression, but can lead to wrong code.
>> Would it be OK to backport to open branches?
> 
> Sure.  Looks simply and fairly specific.

OK, thanks, will proceed.


> I was wondering about the difference between set_exponent()
> and scale(), and found that set_exponent() talks about IEEE
> values while scale() doesn't.  I'm wondering if we should
> add the IEEE special cases to the testsuite.  Of particular
> note, I doubt that this is true:
> 
>     If X is an IEEE NaN, the result is the same NaN.
> 
> program foo
>     real x, y
>     x = 1
>     y = x - x
>     x = (x - x) / y
>     print '(F4.0,1X,Z8.8)', x, x
>     y = set_exponent(x,1)
>     print '(F4.0,1X,Z8.8)', y, y
> end program foo
> 
>   gfcx -o z a.f90 && ./z
>   NaN FFC00000
>   NaN 7FC00000
> 
> Those are not the same NaN.  The second is a qNaN.
> The first looks like a qNaN with the sign bit set.

Until now there was no testing at all of SET_EXPONENT in the testsuite.
It would be really good to have better coverage of compile-time and
runtime behavior of the intrinsics and checking consistency ... ;-)

I think you have much more experience in that area.  (Hint!)

Cheers,
Harald

Re: [PATCH] Fortran: fix compile-time simplification of SET_EXPONENT [PR109511]

[Steve Kargl]

On Fri, Apr 14, 2023 at 09:49:22PM +0200, Harald Anlauf wrote:
> 
> On 4/14/23 21:33, Steve Kargl via Gcc-patches wrote:
> > I was wondering about the difference between set_exponent()
> > and scale(), and found that set_exponent() talks about IEEE
> > values while scale() doesn't.  I'm wondering if we should
> > add the IEEE special cases to the testsuite.  Of particular
> > note, I doubt that this is true:
> > 
> >     If X is an IEEE NaN, the result is the same NaN.
> > 
> > program foo
> >     real x, y
> >     x = 1
> >     y = x - x
> >     x = (x - x) / y
> >     print '(F4.0,1X,Z8.8)', x, x
> >     y = set_exponent(x,1)
> >     print '(F4.0,1X,Z8.8)', y, y
> > end program foo
> > 
> >   gfcx -o z a.f90 && ./z
> >   NaN FFC00000
> >   NaN 7FC00000
> > 
> > Those are not the same NaN.  The second is a qNaN.
> > The first looks like a qNaN with the sign bit set.
> 
> Until now there was no testing at all of SET_EXPONENT in the testsuite.
> It would be really good to have better coverage of compile-time and
> runtime behavior of the intrinsics and checking consistency ... ;-)
> 
> I think you have much more experience in that area.  (Hint!)
> 

I might have some experience. :-)  Unfortunately, I have no time
(for at least 4-5 years) to jump down that rabbit hole (as I would
try to fix things ;-).

However, writing testcases to exercise the intrinsic subprograms
is something that the LURKERS here in the mailing list might
comtemplate.  Any takers?  

Note returning the 'same Nan" is not special to set_exponent().
At least, fraction() and rrspacing() have "If X is an IEEE NaN,
the result is that NaN."


Oh, and thanks for your relentless assault on bugs.

-- 
Steve

Re: [PATCH] Fortran: fix compile-time simplification of SET_EXPONENT [PR109511]

[Bernhard Reutner-Fischer]

On Fri, 14 Apr 2023 20:59:24 +0200
Harald Anlauf via Fortran <fortran@gcc.gnu.org> wrote:

> -  mpfr_clears (absv, log2, pow2, frac, NULL);
> +  mpfr_clears (absv, log2, exp, pow2, frac, NULL);

Thanks for taking care of that leak.
Consider that hunk dropped from
https://inbox.sourceware.org/fortran/D9328809-A533-4D6B-A474-DE4C582D9C02@gmail.com/

although i had it in the order of allocation to ease reading..
cheers,