To go along with whatever magic we're gonna tack along to the
range-ops sqrt implementation, here is another revision addressing the
VARYING issue you pointed out.
A few things...
Instead of going through trees, I decided to call do_mpfr_arg1
directly. Let's not go the wide int <-> tree rat hole in this one.
The function do_mpfr_arg1 bails on +INF, so I had to handle it manually.
There's a regression in gfortran.dg/ieee/ieee_6.f90, which I'm not
sure how to handle. We are failing because we are calculating
sqrt(-1) and expecting certain IEEE flags set. These flags aren't
set, presumably because we folded sqrt(-1) into a NAN directly:
// All negatives.
if (real_compare (LT_EXPR, &lh_ub, &dconst0))
{
real_nan (&lb, "", 0, TYPE_MODE (type));
ub = lb;
maybe_nan = true;
return;
}
The failing part of the test is:
if (.not. (all(flags .eqv. [.false.,.false.,.true.,.true.,.false.]) &
.or. all(flags .eqv. [.false.,.false.,.true.,.true.,.true.]) &
.or. all(flags .eqv. [.false.,.false.,.true.,.false.,.false.]) &
.or. all(flags .eqv.
[.false.,.false.,.true.,.false.,.true.]))) STOP 5
But we are generating F F F F F. Google has informed me that that 3rd
flag is IEEE_INVALID.
So... is the optimization wrong? Are we not allowed to substitute
that NAN if we know it's gonna happen? Should we also allow F F F F F
in the test? Or something else?
Thanks.
Aldy
On Wed, Nov 16, 2022 at 9:33 PM Jakub Jelinek <ja...@redhat.com> wrote:
>
> On Mon, Nov 14, 2022 at 09:55:29PM +0000, Joseph Myers wrote:
> > On Sun, 13 Nov 2022, Jakub Jelinek via Gcc-patches wrote:
> >
> > > So, I wonder if we don't need to add a target hook where targets will be
> > > able to provide upper bound on error for floating point functions for
> > > different floating point modes and some way to signal unknown
> > > accuracy/can't
> > > be trusted, in which case we would give up or return just the range for
> > > VARYING.
> >
> > Note that the figures given in the glibc manual are purely empirical
> > (largest errors observed for inputs in the glibc testsuite on a system
> > that was then used to update the libm-test-ulps files); they don't
> > constitute any kind of guarantee about either the current implementation
> > or the API, nor are they formally verified, nor do they come from
> > exhaustive testing (though worst cases from exhaustive testing for float
> > may have been added to the glibc testsuite in some cases). (I think the
> > only functions known to give huge errors for some inputs, outside of any
> > IBM long double issues, are the Bessel functions and cpow functions. But
> > even if other functions don't have huge errors, and some
> > architecture-specific implementations might have issues, there are
> > certainly some cases where errors can exceed the 9ulp threshold on what
> > the libm tests will accept in libm-test-ulps files, which are thus
> > considered glibc bugs. (That's 9ulp from the correctly rounded value,
> > computed in ulp of that value. For IBM long double it's 16ulp instead,
> > treating the format as having a fixed 106 bits of precision. Both figures
> > are empirical ones chosen based on what bounds sufficed for most libm
> > functions some years ago; ideally, with better implementations of some
> > functions we could probably bring those numbers down.))
>
> I know I can't get guarantees without formal proofs and even ulps from
> reported errors are better than randomized testing.
> But I think at least for non-glibc we want to be able to get a rough idea
> of the usual error range in ulps.
>
> This is what I came up with so far (link with
> gcc -o ulp-tester{,.c} -O2 -lmpfr -lm
> ), it still doesn't verify that functions are always within the mathematical
> range of results ([-0.0, Inf] for sqrt, [-1.0, 1.0] for sin/cos etc.), guess
> that would be useful and verify the program actually does what is intended.
> One can supply just one argument (number of tests, first 46 aren't really
> random) or two, in the latter case the second should be upward, downward or
> towardzero to use non-default rounding mode.
> The idea is that we'd collect ballpark estimates for roundtonearest and
> then estimates for the other 3 rounding modes, the former would be used
> without -frounding-math, max over all 4 rounding modes for -frounding-math
> as gcc will compute using mpfr always in round to nearest.
>
> Jakub
From 759bcd4b4b6f70fcec045b24fb6874aaca989549 Mon Sep 17 00:00:00 2001
From: Aldy Hernandez <al...@redhat.com>
Date: Sun, 13 Nov 2022 18:39:59 +0100
Subject: [PATCH] [range-ops] Implement sqrt.
gcc/ChangeLog:
* fold-const-call.cc (do_mpfr_arg1): Remove static.
* gimple-range-op.cc (class cfn_sqrt): New.
(gimple_range_op_handler::maybe_builtin_call): Add sqrt case.
* realmpfr.h (do_mpfr_arg1): Add extern.
gcc/testsuite/ChangeLog:
* gcc.dg/tree-ssa/vrp124.c: New test.
---
gcc/fold-const-call.cc | 2 +-
gcc/gimple-range-op.cc | 56 ++++++++++++++++++++++++++
gcc/realmpfr.h | 4 ++
gcc/testsuite/gcc.dg/tree-ssa/vrp124.c | 16 ++++++++
4 files changed, 77 insertions(+), 1 deletion(-)
create mode 100644 gcc/testsuite/gcc.dg/tree-ssa/vrp124.c
diff --git a/gcc/fold-const-call.cc b/gcc/fold-const-call.cc
index 8ceed8f02f9..5c6b852acdc 100644
--- a/gcc/fold-const-call.cc
+++ b/gcc/fold-const-call.cc
@@ -118,7 +118,7 @@ do_mpfr_ckconv (real_value *result, mpfr_srcptr m, bool inexact,
in format FORMAT, given that FUNC is the MPFR implementation of f.
Return true on success. */
-static bool
+bool
do_mpfr_arg1 (real_value *result,
int (*func) (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t),
const real_value *arg, const real_format *format)
diff --git a/gcc/gimple-range-op.cc b/gcc/gimple-range-op.cc
index 7764166d5fb..f1f2f098305 100644
--- a/gcc/gimple-range-op.cc
+++ b/gcc/gimple-range-op.cc
@@ -43,6 +43,8 @@ along with GCC; see the file COPYING3. If not see
#include "range.h"
#include "value-query.h"
#include "gimple-range.h"
+#include "fold-const-call.h"
+#include "realmpfr.h"
// Given stmt S, fill VEC, up to VEC_SIZE elements, with relevant ssa-names
// on the statement. For efficiency, it is an error to not pass in enough
@@ -301,6 +303,54 @@ public:
}
} op_cfn_constant_p;
+// Implement range operator for SQRT.
+class cfn_sqrt : public range_operator_float
+{
+ using range_operator_float::fold_range;
+private:
+ void rv_fold (REAL_VALUE_TYPE &lb, REAL_VALUE_TYPE &ub, bool &maybe_nan,
+ tree type,
+ const REAL_VALUE_TYPE &lh_lb,
+ const REAL_VALUE_TYPE &lh_ub,
+ const REAL_VALUE_TYPE &,
+ const REAL_VALUE_TYPE &,
+ relation_kind) const final override
+ {
+ // All negatives.
+ if (real_compare (LT_EXPR, &lh_ub, &dconst0))
+ {
+ real_nan (&lb, "", 0, TYPE_MODE (type));
+ ub = lb;
+ maybe_nan = true;
+ return;
+ }
+ const real_format *format = REAL_MODE_FORMAT (TYPE_MODE (type));
+ // All positives.
+ if (real_compare (GE_EXPR, &lh_lb, &dconst0))
+ {
+ // ?? Handle +INF manually since do_mpfr_arg1 does not.
+ if (real_isinf (&lh_lb, 0))
+ lb = lh_lb;
+ else if (!do_mpfr_arg1 (&lb, mpfr_sqrt, &lh_lb, format))
+ {
+ lb = dconst0;
+ lb.sign = 1;
+ }
+ maybe_nan = false;
+ }
+ // Both positives and negatives.
+ else
+ {
+ // Range is [-0.0, sqrt(lh_ub)] +-NAN.
+ lb = dconst0;
+ lb.sign = 1;
+ maybe_nan = true;
+ }
+ if (!do_mpfr_arg1 (&ub, mpfr_sqrt, &lh_ub, format))
+ ub = dconstinf;
+ }
+} fop_cfn_sqrt;
+
// Implement range operator for CFN_BUILT_IN_SIGNBIT.
class cfn_signbit : public range_operator_float
{
@@ -907,6 +957,12 @@ gimple_range_op_handler::maybe_builtin_call ()
m_int = &op_cfn_parity;
break;
+ CASE_CFN_SQRT_ALL:
+ m_valid = true;
+ m_op1 = gimple_call_arg (call, 0);
+ m_float = &fop_cfn_sqrt;
+ break;
+
default:
break;
}
diff --git a/gcc/realmpfr.h b/gcc/realmpfr.h
index edc08385fe8..807dd2308d2 100644
--- a/gcc/realmpfr.h
+++ b/gcc/realmpfr.h
@@ -32,4 +32,8 @@ extern void real_from_mpfr (REAL_VALUE_TYPE *, mpfr_srcptr,
const real_format *, mpfr_rnd_t);
extern void mpfr_from_real (mpfr_ptr, const REAL_VALUE_TYPE *, mpfr_rnd_t);
+extern bool do_mpfr_arg1 (real_value *result,
+ int (*func) (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t),
+ const real_value *arg, const real_format *format);
+
#endif /* ! GCC_REALGMP_H */
diff --git a/gcc/testsuite/gcc.dg/tree-ssa/vrp124.c b/gcc/testsuite/gcc.dg/tree-ssa/vrp124.c
new file mode 100644
index 00000000000..ef72d660153
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/tree-ssa/vrp124.c
@@ -0,0 +1,16 @@
+// { dg-do compile }
+// { dg-options "-O2 -fdump-tree-evrp -fdisable-tree-ethread" }
+
+void link_error ();
+
+void foo (float f)
+{
+ float z = __builtin_sqrt (f);
+ if (!__builtin_isnan (z))
+ {
+ if (z < 0.0)
+ link_error ();
+ }
+}
+
+// { dg-final { scan-tree-dump-not "link_error" "evrp" } }
--
2.38.1
