[PATCH] Fix up norm2 simplification (PR middle-end/88074)

[PATCH] Fix up norm2 simplification (PR middle-end/88074)

[Jakub Jelinek]

Hi!

As discussed recently on the mailing list, the norm2 simplification doesn't
work if we limit mpfr emin/emax to some values derived from maximum floating
exponents (and precision for denormals).

The following patch adjusts the computation, so that it is scaled down if
needed.  In particular, if any value in the array is so large that **2 will
overflow on it or will be very close to it, the scale down is set to
2**(max_exponent/2+4), and if the result during computation gets close to
overflowing, it is scaled down a little bit too.  The scaling is always done
using powers of two, operands by that and the sum by **2 of that, and at the
end it multiplies the sqrt back.  I had to change
simplify_transformation_to_array, so that post_op is done immediately after
finishing ops corresponding to that, so that there can be just one global
variable for the scale.  From my understanding of e.g. the libgfortran norm2
code where sqrt is done basically in this spot I hope it isn't possible that
the same *dest is updated multiple times with dest increments/decrements in
between.

Bootstrapped/regtested on x86_64-linux and i686-linux (together with
Richard's patch), ok for trunk?

2019-02-12  Jakub Jelinek  <jakub@redhat.com>

	PR middle-end/88074
	* simplify.c (simplify_transformation_to_array): Run post_op
	immediately after processing corresponding row, rather than at the
	end.
	(norm2_scale): New variable.
	(add_squared): Rename to ...
	(norm2_add_squared): ... this.  Scale down operand and/or result
	if needed.
	(do_sqrt): Rename to ...
	(norm2_do_sqrt): ... this.  Handle the result == e case.  Scale up
	result and clear norm2_scale.
	(gfc_simplify_norm2): Clear norm2_scale.  Change add_squared to
	norm2_add_squared and &do_sqrt to norm2_do_sqrt.  Scale up result
	and clear norm2_scale again.

--- gcc/fortran/simplify.c.jj	2019-01-10 11:43:12.452409482 +0100
+++ gcc/fortran/simplify.c	2019-02-12 19:54:03.726526824 +0100
@@ -636,6 +636,9 @@ simplify_transformation_to_array (gfc_ex
 	if (*src)
 	  *dest = op (*dest, gfc_copy_expr (*src));
 
+      if (post_op)
+	*dest = post_op (*dest, *dest);
+
       count[0]++;
       base += sstride[0];
       dest += dstride[0];
@@ -671,10 +674,7 @@ simplify_transformation_to_array (gfc_ex
   result_ctor = gfc_constructor_first (result->value.constructor);
   for (i = 0; i < resultsize; ++i)
     {
-      if (post_op)
-	result_ctor->expr = post_op (result_ctor->expr, resultvec[i]);
-      else
-	result_ctor->expr = resultvec[i];
+      result_ctor->expr = resultvec[i];
       result_ctor = gfc_constructor_next (result_ctor);
     }
 
@@ -6048,9 +6048,10 @@ gfc_simplify_idnint (gfc_expr *e)
   return simplify_nint ("IDNINT", e, NULL);
 }
 
+static int norm2_scale;
 
 static gfc_expr *
-add_squared (gfc_expr *result, gfc_expr *e)
+norm2_add_squared (gfc_expr *result, gfc_expr *e)
 {
   mpfr_t tmp;
 
@@ -6059,8 +6060,45 @@ add_squared (gfc_expr *result, gfc_expr
 	      && result->expr_type == EXPR_CONSTANT);
 
   gfc_set_model_kind (result->ts.kind);
+  int index = gfc_validate_kind (BT_REAL, result->ts.kind, false);
+  mpfr_exp_t exp;
+  if (mpfr_regular_p (result->value.real))
+    {
+      exp = mpfr_get_exp (result->value.real);
+      /* If result is getting close to overflowing, scale down.  */
+      if (exp >= gfc_real_kinds[index].max_exponent - 4
+	  && norm2_scale <= gfc_real_kinds[index].max_exponent - 2)
+	{
+	  norm2_scale += 2;
+	  mpfr_div_ui (result->value.real, result->value.real, 16,
+		       GFC_RND_MODE);
+	}
+    }
+
   mpfr_init (tmp);
-  mpfr_pow_ui (tmp, e->value.real, 2, GFC_RND_MODE);
+  if (mpfr_regular_p (e->value.real))
+    {
+      exp = mpfr_get_exp (e->value.real);
+      /* If e**2 would overflow or close to overflowing, scale down.  */
+      if (exp - norm2_scale >= gfc_real_kinds[index].max_exponent / 2 - 2)
+	{
+	  int new_scale = gfc_real_kinds[index].max_exponent / 2 + 4;
+	  mpfr_set_ui (tmp, 1, GFC_RND_MODE);
+	  mpfr_set_exp (tmp, new_scale - norm2_scale);
+	  mpfr_div (result->value.real, result->value.real, tmp, GFC_RND_MODE);
+	  mpfr_div (result->value.real, result->value.real, tmp, GFC_RND_MODE);
+	  norm2_scale = new_scale;
+	}
+    }
+  if (norm2_scale)
+    {
+      mpfr_set_ui (tmp, 1, GFC_RND_MODE);
+      mpfr_set_exp (tmp, norm2_scale);
+      mpfr_div (tmp, e->value.real, tmp, GFC_RND_MODE);
+    }
+  else
+    mpfr_set (tmp, e->value.real, GFC_RND_MODE);
+  mpfr_pow_ui (tmp, tmp, 2, GFC_RND_MODE);
   mpfr_add (result->value.real, result->value.real, tmp,
 	    GFC_RND_MODE);
   mpfr_clear (tmp);
@@ -6070,14 +6108,26 @@ add_squared (gfc_expr *result, gfc_expr
 
 
 static gfc_expr *
-do_sqrt (gfc_expr *result, gfc_expr *e)
+norm2_do_sqrt (gfc_expr *result, gfc_expr *e)
 {
   gcc_assert (e->ts.type == BT_REAL && e->expr_type == EXPR_CONSTANT);
   gcc_assert (result->ts.type == BT_REAL
 	      && result->expr_type == EXPR_CONSTANT);
 
-  mpfr_set (result->value.real, e->value.real, GFC_RND_MODE);
+  if (result != e)
+    mpfr_set (result->value.real, e->value.real, GFC_RND_MODE);
   mpfr_sqrt (result->value.real, result->value.real, GFC_RND_MODE);
+  if (norm2_scale && mpfr_regular_p (result->value.real))
+    {
+      mpfr_t tmp;
+      mpfr_init (tmp);
+      mpfr_set_ui (tmp, 1, GFC_RND_MODE);
+      mpfr_set_exp (tmp, norm2_scale);
+      mpfr_mul (result->value.real, result->value.real, tmp, GFC_RND_MODE);
+      mpfr_clear (tmp);
+    }
+  norm2_scale = 0;
+
   return result;
 }
 
@@ -6100,15 +6150,27 @@ gfc_simplify_norm2 (gfc_expr *e, gfc_exp
   if (size_zero)
     return result;
 
+  norm2_scale = 0;
   if (!dim || e->rank == 1)
     {
       result = simplify_transformation_to_scalar (result, e, NULL,
-						  add_squared);
+						  norm2_add_squared);
       mpfr_sqrt (result->value.real, result->value.real, GFC_RND_MODE);
+      if (norm2_scale && mpfr_regular_p (result->value.real))
+	{
+	  mpfr_t tmp;
+	  mpfr_init (tmp);
+	  mpfr_set_ui (tmp, 1, GFC_RND_MODE);
+	  mpfr_set_exp (tmp, norm2_scale);
+	  mpfr_mul (result->value.real, result->value.real, tmp, GFC_RND_MODE);
+	  mpfr_clear (tmp);
+	}
+      norm2_scale = 0;
     }
   else
     result = simplify_transformation_to_array (result, e, dim, NULL,
-					       add_squared, &do_sqrt);
+					       norm2_add_squared,
+					       norm2_do_sqrt);
 
   return result;
 }

	Jakub

Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074)

[Thomas Koenig]

Hi Jakub,

> Bootstrapped/regtested on x86_64-linux and i686-linux (together with
> Richard's patch), ok for trunk?

A couple of points:

with this patch, we will have an algorithm that will evaluate NORM2
in a different way than before, possibly leading to regressions
in some cases where we previously generated good results and
would return +/- Inf now.

Also, with this patch, we would use a very different algorithm
at run-time to compile-time, potentially leading to differences
far exceeding a few ulps. I am not sure I like that idea.

I'll look at the literature and generate some test cases, and
come back to you.

 From curiosity, did you chose this algorithm from literature?

Regards

	Thomas

Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074)

[Jakub Jelinek]

On Wed, Feb 13, 2019 at 07:17:53PM +0100, Thomas Koenig wrote:
> > Bootstrapped/regtested on x86_64-linux and i686-linux (together with
> > Richard's patch), ok for trunk?
> 
> A couple of points:
> 
> with this patch, we will have an algorithm that will evaluate NORM2
> in a different way than before, possibly leading to regressions
> in some cases where we previously generated good results and
> would return +/- Inf now.

I'm not aware of such a case, do you have examples what do you have in mind?
The adjustments are only done by powers of two, the only reason I'm not
using just mpfr_set_exp is to deal with the case where there are extreme
differences in exponents, say 1.0e-8100L*1.0e-8100L+1.0e+16382L*1.0e+16382L,
where because 1.0e+16382L*1.0e+16382L would overflow the very tiny temporary
result needs to be divided by a very large power of 2 and so I want to let
mpfr do the rounding.

> Also, with this patch, we would use a very different algorithm
> at run-time to compile-time, potentially leading to differences
> far exceeding a few ulps. I am not sure I like that idea.

The runtime code is already so much different that I don't see how this
holds.  And it is actually much worse, exactly because it doesn't scale just
by powers of two and thus every scaling has a rounding error.

> I'll look at the literature and generate some test cases, and
> come back to you.
> 
> From curiosity, did you chose this algorithm from literature?

No.  But, first of all, the patch should change nothing at all unless
either at least one extremely large number (whose **2 is not representable
in the kind or close to that) appears, or the array is so huge that the
result is attacking the maximum.  And, for the cases where it does change
something, unless the differences in exponents between numbers are extreme,
there should be again no rounding difference, the exponent of the numbers
will simply change and nothing else.

	Jakub

Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074)

[Steve Kargl]

On Wed, Feb 13, 2019 at 07:30:53PM +0100, Jakub Jelinek wrote:
> On Wed, Feb 13, 2019 at 07:17:53PM +0100, Thomas Koenig wrote:
> > > Bootstrapped/regtested on x86_64-linux and i686-linux (together with
> > > Richard's patch), ok for trunk?
> > 
> > A couple of points:
> > 
> > with this patch, we will have an algorithm that will evaluate NORM2
> > in a different way than before, possibly leading to regressions
> > in some cases where we previously generated good results and
> > would return +/- Inf now.
> 
> I'm not aware of such a case, do you have examples what do you have in mind?
> The adjustments are only done by powers of two, the only reason I'm not
> using just mpfr_set_exp is to deal with the case where there are extreme
> differences in exponents, say 1.0e-8100L*1.0e-8100L+1.0e+16382L*1.0e+16382L,
> where because 1.0e+16382L*1.0e+16382L would overflow the very tiny temporary
> result needs to be divided by a very large power of 2 and so I want to let
> mpfr do the rounding.
> 
> > Also, with this patch, we would use a very different algorithm
> > at run-time to compile-time, potentially leading to differences
> > far exceeding a few ulps. I am not sure I like that idea.
> 
> The runtime code is already so much different that I don't see how this
> holds.  And it is actually much worse, exactly because it doesn't scale just
> by powers of two and thus every scaling has a rounding error.
> 
> > I'll look at the literature and generate some test cases, and
> > come back to you.
> > 
> > From curiosity, did you chose this algorithm from literature?
> 
> No.  But, first of all, the patch should change nothing at all unless
> either at least one extremely large number (whose **2 is not representable
> in the kind or close to that) appears, or the array is so huge that the
> result is attacking the maximum.  And, for the cases where it does change
> something, unless the differences in exponents between numbers are extreme,
> there should be again no rounding difference, the exponent of the numbers
> will simply change and nothing else.
> 

I think we can agree that gfortran's norm2 needs some work.

To get Richard's emin/emax patch into the tree so adequate
testing can be done prior to the 9.1 release, I'll suggest
that we XFAIL norm2_3.f90 while we work out the details for
norm2.

Yesteraday, I looked at Jakub's patch, and thought the approach
was fine.  I haven't had time to look at the library side.
Perhaps, we can adopt Jakub's approach for the library as well.

-- 
Steve

Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074)

[Jakub Jelinek]

On Wed, Feb 13, 2019 at 10:51:07AM -0800, Steve Kargl wrote:
> To get Richard's emin/emax patch into the tree so adequate
> testing can be done prior to the 9.1 release, I'll suggest
> that we XFAIL norm2_3.f90 while we work out the details for
> norm2.
> 
> Yesteraday, I looked at Jakub's patch, and thought the approach
> was fine.  I haven't had time to look at the library side.
> Perhaps, we can adopt Jakub's approach for the library as well.

Yeah.  I guess fastest would be to precompute as constants the
two boundaries when we might want to scale something (of course different
for each kind), then just compare the operand and/or result against
those and only if they are above that do frexp* to determine exponent,
(re)compute the scale variable (scalbn* ?) and then just multiply or divide
by that (if division is much slower, we could have two scale variables, one
for scaling down and one for scaling up).

	Jakub

Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074)

[Thomas Koenig]

Hi Jakub,

I checked the patch together with Richard's (by which I assume you
mean https://gcc.gnu.org/bugzilla/attachment.cgi?id=45052 ), and
thinks looked good.

So, the Fortran part of this is OK.

However, we should really also do power-of-two scaling
for the runtime method.

Also, we seem to have a lot of issues with IEEE flags
when calculating NORM2, this would also need to be
addressed.

Regards

	Thomas

Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074)

[Steve Kargl]

On Sat, Feb 16, 2019 at 05:23:58PM +0100, Thomas Koenig wrote:
> 
> Also, we seem to have a lot of issues with IEEE flags
> when calculating NORM2, this would also need to be
> addressed.
> 

Which IEEE flags and are you referring using the
Fortran modules or -ffpe-trap?

-- 
Steve

Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074)

[Thomas Koenig]

Hi Steve,

> On Sat, Feb 16, 2019 at 05:23:58PM +0100, Thomas Koenig wrote:
>>
>> Also, we seem to have a lot of issues with IEEE flags
>> when calculating NORM2, this would also need to be
>> addressed.
>>
> 
> Which IEEE flags and are you referring using the
> Fortran modules or -ffpe-trap?

I checked out the software from that paper (thanks for the link), and
gfortran ran up quite a lot of errors running the test suite even with
the author's preferred version.  That's what I meant.

Regards

	Thomas

[PATCH] Minimum version of mpfr? (was Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074))

[David Malcolm]

On Wed, 2019-02-13 at 00:09 +0100, Jakub Jelinek wrote:
> Hi!
>
> As discussed recently on the mailing list, the norm2 simplification
> doesn't
> work if we limit mpfr emin/emax to some values derived from maximum
> floating
> exponents (and precision for denormals).
>
> The following patch adjusts the computation, so that it is scaled
> down if
> needed.  In particular, if any value in the array is so large that
> **2 will
> overflow on it or will be very close to it, the scale down is set to
> 2**(max_exponent/2+4), and if the result during computation gets
> close to
> overflowing, it is scaled down a little bit too.  The scaling is
> always done
> using powers of two, operands by that and the sum by **2 of that, and
> at the
> end it multiplies the sqrt back.  I had to change
> simplify_transformation_to_array, so that post_op is done immediately
> after
> finishing ops corresponding to that, so that there can be just one
> global
> variable for the scale.  From my understanding of e.g. the
> libgfortran norm2
> code where sqrt is done basically in this spot I hope it isn't
> possible that
> the same *dest is updated multiple times with dest
> increments/decrements in
> between.
>
> Bootstrapped/regtested on x86_64-linux and i686-linux (together with
> Richard's patch), ok for trunk?
>
> 2019-02-12  Jakub Jelinek  <jakub@redhat.com>
>
> 	PR middle-end/88074
> 	* simplify.c (simplify_transformation_to_array): Run post_op
> 	immediately after processing corresponding row, rather than at
> the
> 	end.
> 	(norm2_scale): New variable.
> 	(add_squared): Rename to ...
> 	(norm2_add_squared): ... this.  Scale down operand and/or
> result
> 	if needed.
> 	(do_sqrt): Rename to ...
> 	(norm2_do_sqrt): ... this.  Handle the result == e case.  Scale
> up
> 	result and clear norm2_scale.
> 	(gfc_simplify_norm2): Clear norm2_scale.  Change add_squared to
> 	norm2_add_squared and &do_sqrt to norm2_do_sqrt.  Scale up
> result
> 	and clear norm2_scale again.
>

[...]

>  static gfc_expr *
> -add_squared (gfc_expr *result, gfc_expr *e)
> +norm2_add_squared (gfc_expr *result, gfc_expr *e)
>  {
>    mpfr_t tmp;
>
> @@ -6059,8 +6060,45 @@ add_squared (gfc_expr *result, gfc_expr
>  	      && result->expr_type == EXPR_CONSTANT);
>
>    gfc_set_model_kind (result->ts.kind);
> +  int index = gfc_validate_kind (BT_REAL, result->ts.kind, false);
> +  mpfr_exp_t exp;
> +  if (mpfr_regular_p (result->value.real))
> +    {

I've started seeing build failures in my testing setup here, apparently
introduced with this patch:

../../../src/gcc/fortran/simplify.c: In function ‘gfc_expr* norm2_add_squared(gfc_expr*, gfc_expr*)’:
../../../src/gcc/fortran/simplify.c:6064:3: error: ‘mpfr_exp_t’ was not declared in this scope; did you mean ‘mpfr_expm1’?
 6064 |   mpfr_exp_t exp;
      |   ^~~~~~~~~~
      |   mpfr_expm1

I'm using mpfr-2.4.2.tar.bz2, which is the minimum version stated at:
  https://gcc.gnu.org/install/prerequisites.html
which says "MPFR Library version 2.4.2 (or later)"; that
version doesn't seem to have that type.

However, I see on:
  https://www.mpfr.org/mpfr-3.0.1/mpfr.html#API-Compatibility
that "The official type for exponent values changed from mp_exp_t to
mpfr_exp_t in MPFR 3.0."

So should the patch be tweaked to use the old type name, like in the
following patch (caveat: not tested yet), or do we need to update
the minimum version of mpfr, and if so, to what version?

[...]

Dave

gcc/fortran/ChangeLog:
	PR middle-end/88074
	* simplify.c (norm2_add_squared): Use mp_exp_t rather than
	mpfr_exp_t.
---
 gcc/fortran/simplify.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
index fa6396b..a1df735 100644
--- a/gcc/fortran/simplify.c
+++ b/gcc/fortran/simplify.c
@@ -6061,7 +6061,7 @@ norm2_add_squared (gfc_expr *result, gfc_expr *e)
 
   gfc_set_model_kind (result->ts.kind);
   int index = gfc_validate_kind (BT_REAL, result->ts.kind, false);
-  mpfr_exp_t exp;
+  mp_exp_t exp;
   if (mpfr_regular_p (result->value.real))
     {
       exp = mpfr_get_exp (result->value.real);
-- 
1.8.5.3

Re: [PATCH] Minimum version of mpfr? (was Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074))

[Jakub Jelinek]

On Thu, Feb 21, 2019 at 04:50:27PM -0500, David Malcolm wrote:
>  gcc/fortran/simplify.c | 2 +-
>  1 file changed, 1 insertion(+), 1 deletion(-)
> 
> diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
> index fa6396b..a1df735 100644
> --- a/gcc/fortran/simplify.c
> +++ b/gcc/fortran/simplify.c
> @@ -6061,7 +6061,7 @@ norm2_add_squared (gfc_expr *result, gfc_expr *e)
>  
>    gfc_set_model_kind (result->ts.kind);
>    int index = gfc_validate_kind (BT_REAL, result->ts.kind, false);
> -  mpfr_exp_t exp;
> +  mp_exp_t exp;
>    if (mpfr_regular_p (result->value.real))
>      {
>        exp = mpfr_get_exp (result->value.real);

Ah, sorry for that.
I guess either we can use something like:
#if MPFR_VERSION < MPFR_VERSION_NUM(3,0,0)
  mp_exp_t exp;
#else
  mpfr_exp_t exp;
#endif

or perhaps just
  long exp;
or
  int exp;
?  We don't have pointers to that passed somewhere, so all we care is that
we cover the range (but, because of the previous limiting of the exponents
in toplev.c, it must fit into integer, because that is the type of the
min_exp/max_exp passed to mpfr_set_e{min,max}.

	Jakub

[PATCH] Minimum version of mpfr? (was Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074)), take 2

[Jakub Jelinek]

On Thu, Feb 21, 2019 at 04:50:27PM -0500, David Malcolm wrote:
> gcc/fortran/ChangeLog:
> 	PR middle-end/88074
> 	* simplify.c (norm2_add_squared): Use mp_exp_t rather than
> 	mpfr_exp_t.

I have discussed this with Richard on IRC earlier today, there is another
issue that mpfr_regular_p is only 3.0 and later.  And he prefers that
mp_exp_t over say
#if MPFR_VERSION < MPFR_VERSION_NUM(3,0,0)
typedef mp_exp_t mpfr_exp_t;
#endif

So, here is the full patch I'll bootstrap/regtest soon:

2019-02-22  David Malcolm  <dmalcolm@redhat.com>
	    Jakub Jelinek  <jakub@redhat.com>

	PR middle-end/88074
	* simplify.c (norm2_do_sqrt, gfc_simplify_norm2): Use
	mpfr_number_p && !mpfr_zero_p instead of mpfr_regular_p.
	(norm2_add_squared): Likewise.  Use mp_exp_t rather than mpfr_exp_t.

--- gcc/fortran/simplify.c.jj	2019-02-21 22:20:07.272388243 +0100
+++ gcc/fortran/simplify.c	2019-02-22 11:07:36.093374586 +0100
@@ -6061,8 +6061,8 @@ norm2_add_squared (gfc_expr *result, gfc
 
   gfc_set_model_kind (result->ts.kind);
   int index = gfc_validate_kind (BT_REAL, result->ts.kind, false);
-  mpfr_exp_t exp;
-  if (mpfr_regular_p (result->value.real))
+  mp_exp_t exp;
+  if (mpfr_number_p (result->value.real) && !mpfr_zero_p (result->value.real))
     {
       exp = mpfr_get_exp (result->value.real);
       /* If result is getting close to overflowing, scale down.  */
@@ -6076,7 +6076,7 @@ norm2_add_squared (gfc_expr *result, gfc
     }
 
   mpfr_init (tmp);
-  if (mpfr_regular_p (e->value.real))
+  if (mpfr_number_p (e->value.real) && !mpfr_zero_p (e->value.real))
     {
       exp = mpfr_get_exp (e->value.real);
       /* If e**2 would overflow or close to overflowing, scale down.  */
@@ -6117,7 +6117,9 @@ norm2_do_sqrt (gfc_expr *result, gfc_exp
   if (result != e)
     mpfr_set (result->value.real, e->value.real, GFC_RND_MODE);
   mpfr_sqrt (result->value.real, result->value.real, GFC_RND_MODE);
-  if (norm2_scale && mpfr_regular_p (result->value.real))
+  if (norm2_scale
+      && mpfr_number_p (result->value.real)
+      && !mpfr_zero_p (result->value.real))
     {
       mpfr_t tmp;
       mpfr_init (tmp);
@@ -6156,7 +6158,9 @@ gfc_simplify_norm2 (gfc_expr *e, gfc_exp
       result = simplify_transformation_to_scalar (result, e, NULL,
 						  norm2_add_squared);
       mpfr_sqrt (result->value.real, result->value.real, GFC_RND_MODE);
-      if (norm2_scale && mpfr_regular_p (result->value.real))
+      if (norm2_scale
+	  && mpfr_number_p (result->value.real)
+	  && !mpfr_zero_p (result->value.real))
 	{
 	  mpfr_t tmp;
 	  mpfr_init (tmp);

	Jakub

Re: [PATCH] Minimum version of mpfr? (was Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074)), take 2

[Richard Biener]

On Fri, 22 Feb 2019, Jakub Jelinek wrote:

> On Thu, Feb 21, 2019 at 04:50:27PM -0500, David Malcolm wrote:
> > gcc/fortran/ChangeLog:
> > 	PR middle-end/88074
> > 	* simplify.c (norm2_add_squared): Use mp_exp_t rather than
> > 	mpfr_exp_t.
> 
> I have discussed this with Richard on IRC earlier today, there is another
> issue that mpfr_regular_p is only 3.0 and later.  And he prefers that
> mp_exp_t over say
> #if MPFR_VERSION < MPFR_VERSION_NUM(3,0,0)
> typedef mp_exp_t mpfr_exp_t;
> #endif
> 
> So, here is the full patch I'll bootstrap/regtest soon:

OK if it passes.

Apart from autotesters where some of ours run very old outdated systems...
for SLES 12 which I think is the oldest relevant distro I'd like to
see newer GCCs to work on we have mpfr 3.1.2 and gmp 5.1.3 just in
case we might finally consider bumping the minimum required versions
at some point.  Though I didn't yet see any technical reason to do so,
we have more "issues" supporting old host GCC versions (SLE12 uses
GCC 4.8 so that would be a considerable bump already).

Richard.

> 2019-02-22  David Malcolm  <dmalcolm@redhat.com>
> 	    Jakub Jelinek  <jakub@redhat.com>
> 
> 	PR middle-end/88074
> 	* simplify.c (norm2_do_sqrt, gfc_simplify_norm2): Use
> 	mpfr_number_p && !mpfr_zero_p instead of mpfr_regular_p.
> 	(norm2_add_squared): Likewise.  Use mp_exp_t rather than mpfr_exp_t.
> 
> --- gcc/fortran/simplify.c.jj	2019-02-21 22:20:07.272388243 +0100
> +++ gcc/fortran/simplify.c	2019-02-22 11:07:36.093374586 +0100
> @@ -6061,8 +6061,8 @@ norm2_add_squared (gfc_expr *result, gfc
>  
>    gfc_set_model_kind (result->ts.kind);
>    int index = gfc_validate_kind (BT_REAL, result->ts.kind, false);
> -  mpfr_exp_t exp;
> -  if (mpfr_regular_p (result->value.real))
> +  mp_exp_t exp;
> +  if (mpfr_number_p (result->value.real) && !mpfr_zero_p (result->value.real))
>      {
>        exp = mpfr_get_exp (result->value.real);
>        /* If result is getting close to overflowing, scale down.  */
> @@ -6076,7 +6076,7 @@ norm2_add_squared (gfc_expr *result, gfc
>      }
>  
>    mpfr_init (tmp);
> -  if (mpfr_regular_p (e->value.real))
> +  if (mpfr_number_p (e->value.real) && !mpfr_zero_p (e->value.real))
>      {
>        exp = mpfr_get_exp (e->value.real);
>        /* If e**2 would overflow or close to overflowing, scale down.  */
> @@ -6117,7 +6117,9 @@ norm2_do_sqrt (gfc_expr *result, gfc_exp
>    if (result != e)
>      mpfr_set (result->value.real, e->value.real, GFC_RND_MODE);
>    mpfr_sqrt (result->value.real, result->value.real, GFC_RND_MODE);
> -  if (norm2_scale && mpfr_regular_p (result->value.real))
> +  if (norm2_scale
> +      && mpfr_number_p (result->value.real)
> +      && !mpfr_zero_p (result->value.real))
>      {
>        mpfr_t tmp;
>        mpfr_init (tmp);
> @@ -6156,7 +6158,9 @@ gfc_simplify_norm2 (gfc_expr *e, gfc_exp
>        result = simplify_transformation_to_scalar (result, e, NULL,
>  						  norm2_add_squared);
>        mpfr_sqrt (result->value.real, result->value.real, GFC_RND_MODE);
> -      if (norm2_scale && mpfr_regular_p (result->value.real))
> +      if (norm2_scale
> +	  && mpfr_number_p (result->value.real)
> +	  && !mpfr_zero_p (result->value.real))
>  	{
>  	  mpfr_t tmp;
>  	  mpfr_init (tmp);
> 
> 	Jakub
> 
> 

-- 
Richard Biener <rguenther@suse.de>
SUSE LINUX GmbH, GF: Felix Imendoerffer, Jane Smithard, Graham Norton, HRB 21284 (AG Nuernberg)

Re: [PATCH] Minimum version of mpfr? (was Re: [PATCH] Fix up norm2 simplification (PR middle-end/88074)), take 2

[David Malcolm]

On Fri, 2019-02-22 at 14:19 +0100, Jakub Jelinek wrote:
> On Thu, Feb 21, 2019 at 04:50:27PM -0500, David Malcolm wrote:
> > gcc/fortran/ChangeLog:
> > 	PR middle-end/88074
> > 	* simplify.c (norm2_add_squared): Use mp_exp_t rather than
> > 	mpfr_exp_t.
> 
> I have discussed this with Richard on IRC earlier today, there is
> another
> issue that mpfr_regular_p is only 3.0 and later.  And he prefers that
> mp_exp_t over say
> #if MPFR_VERSION < MPFR_VERSION_NUM(3,0,0)
> typedef mp_exp_t mpfr_exp_t;
> #endif
> 
> So, here is the full patch I'll bootstrap/regtest soon:

FWIW, your patch fixed the bootstrap for me (and the regression tests
look sane:
  gfortran.sum : total: 49370 PASS: 49164 XFAIL: 126 UNSUPPORTED: 80
though sadly I've blown away my known-good baseline for comparison)

Dave

> 2019-02-22  David Malcolm  <dmalcolm@redhat.com>
> 	    Jakub Jelinek  <jakub@redhat.com>
> 
> 	PR middle-end/88074
> 	* simplify.c (norm2_do_sqrt, gfc_simplify_norm2): Use
> 	mpfr_number_p && !mpfr_zero_p instead of mpfr_regular_p.
> 	(norm2_add_squared): Likewise.  Use mp_exp_t rather than
> mpfr_exp_t.
> 
> --- gcc/fortran/simplify.c.jj	2019-02-21 22:20:07.272388243
> +0100
> +++ gcc/fortran/simplify.c	2019-02-22 11:07:36.093374586 +0100
> @@ -6061,8 +6061,8 @@ norm2_add_squared (gfc_expr *result, gfc
>  
>    gfc_set_model_kind (result->ts.kind);
>    int index = gfc_validate_kind (BT_REAL, result->ts.kind, false);
> -  mpfr_exp_t exp;
> -  if (mpfr_regular_p (result->value.real))
> +  mp_exp_t exp;
> +  if (mpfr_number_p (result->value.real) && !mpfr_zero_p (result-
> >value.real))
>      {
>        exp = mpfr_get_exp (result->value.real);
>        /* If result is getting close to overflowing, scale down.  */
> @@ -6076,7 +6076,7 @@ norm2_add_squared (gfc_expr *result, gfc
>      }
>  
>    mpfr_init (tmp);
> -  if (mpfr_regular_p (e->value.real))
> +  if (mpfr_number_p (e->value.real) && !mpfr_zero_p (e->value.real))
>      {
>        exp = mpfr_get_exp (e->value.real);
>        /* If e**2 would overflow or close to overflowing, scale
> down.  */
> @@ -6117,7 +6117,9 @@ norm2_do_sqrt (gfc_expr *result, gfc_exp
>    if (result != e)
>      mpfr_set (result->value.real, e->value.real, GFC_RND_MODE);
>    mpfr_sqrt (result->value.real, result->value.real, GFC_RND_MODE);
> -  if (norm2_scale && mpfr_regular_p (result->value.real))
> +  if (norm2_scale
> +      && mpfr_number_p (result->value.real)
> +      && !mpfr_zero_p (result->value.real))
>      {
>        mpfr_t tmp;
>        mpfr_init (tmp);
> @@ -6156,7 +6158,9 @@ gfc_simplify_norm2 (gfc_expr *e, gfc_exp
>        result = simplify_transformation_to_scalar (result, e, NULL,
>  						  norm2_add_squared)
> ;
>        mpfr_sqrt (result->value.real, result->value.real,
> GFC_RND_MODE);
> -      if (norm2_scale && mpfr_regular_p (result->value.real))
> +      if (norm2_scale
> +	  && mpfr_number_p (result->value.real)
> +	  && !mpfr_zero_p (result->value.real))
>  	{
>  	  mpfr_t tmp;
>  	  mpfr_init (tmp);
> 
> 	Jakub