* [PATCH] Implement range-op entry for sin/cos. @ 2023-04-18 13:12 Aldy Hernandez 2023-04-20 12:59 ` Jakub Jelinek 2023-04-27 11:13 ` [PATCH] v2: " Jakub Jelinek 0 siblings, 2 replies; 25+ messages in thread From: Aldy Hernandez @ 2023-04-18 13:12 UTC (permalink / raw) To: GCC patches; +Cc: Andrew MacLeod, Jakub Jelinek, Aldy Hernandez [I don't know why I keep poking at floats. I must really like the pain. Jakub, are you OK with this patch for trunk?] This is the range-op entry for sin/cos. It is meant to serve as an example of what we can do for glibc math functions. It is by no means exhaustive, just a stub to restrict the return range from sin/cos to [-1.0, 1.0] with appropriate smarts of NANs. As can be seen in the testcase, we see sin() as well as __builtin_sin() in the IL, and can resolve the resulting range accordingly. gcc/ChangeLog: * gimple-range-op.cc (class cfn_sincos): New. (gimple_range_op_handler::maybe_builtin_call): Add case for sin/cos. gcc/testsuite/ChangeLog: * gcc.dg/tree-ssa/range-sincos.c: New test. --- gcc/gimple-range-op.cc | 63 ++++++++++++++++++++ gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c | 40 +++++++++++++ 2 files changed, 103 insertions(+) create mode 100644 gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c diff --git a/gcc/gimple-range-op.cc b/gcc/gimple-range-op.cc index 4ca32a7b5d5..36390f2645e 100644 --- a/gcc/gimple-range-op.cc +++ b/gcc/gimple-range-op.cc @@ -402,6 +402,60 @@ public: } } op_cfn_copysign; +class cfn_sincos : public range_operator_float +{ +public: + using range_operator_float::fold_range; + using range_operator_float::op1_range; + virtual bool fold_range (frange &r, tree type, + const frange &lh, const frange &, + relation_trio) const final override + { + if (lh.known_isnan () || lh.known_isinf ()) + { + r.set_nan (type); + return true; + } + r.set (type, dconstm1, dconst1); + if (!lh.maybe_isnan ()) + r.clear_nan (); + return true; + } + virtual bool op1_range (frange &r, tree type, + const frange &lhs, const frange &, + relation_trio) const final override + { + if (!lhs.maybe_isnan ()) + { + // If NAN is not valid result, the input cannot include either + // a NAN nor a +-INF. + REAL_VALUE_TYPE lb = real_min_representable (type); + REAL_VALUE_TYPE ub = real_max_representable (type); + r.set (type, lb, ub, nan_state (false, false)); + return true; + } + // A known NAN means the input is [-INF,-INF][+INF,+INF] U +-NAN, + // which we can't currently represent. + if (lhs.known_isnan ()) + { + r.set_varying (type); + return true; + } + // Results outside of [-1.0, +1.0] are impossible. + REAL_VALUE_TYPE lb = lhs.lower_bound (); + REAL_VALUE_TYPE ub = lhs.upper_bound (); + if (real_less (&lb, &dconstm1) + || real_less (&dconst1, &ub)) + { + r.set_undefined (); + return true; + } + + r.set_varying (type); + return true; + } +} op_cfn_sincos; + // Implement range operator for CFN_BUILT_IN_TOUPPER and CFN_BUILT_IN_TOLOWER. class cfn_toupper_tolower : public range_operator { @@ -880,6 +934,15 @@ gimple_range_op_handler::maybe_builtin_call () m_valid = true; break; + CASE_CFN_SIN: + CASE_CFN_SIN_FN: + CASE_CFN_COS: + CASE_CFN_COS_FN: + m_op1 = gimple_call_arg (call, 0); + m_float = &op_cfn_sincos; + m_valid = true; + break; + case CFN_BUILT_IN_TOUPPER: case CFN_BUILT_IN_TOLOWER: // Only proceed If the argument is compatible with the LHS. diff --git a/gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c b/gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c new file mode 100644 index 00000000000..50fbac1b68f --- /dev/null +++ b/gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c @@ -0,0 +1,40 @@ +// { dg-do compile } +// { dg-options "-O2 -fdump-tree-evrp -fno-thread-jumps" } + +#include <math.h> + +void use(double); +void link_error (); + +void foo(double x) +{ + if (__builtin_isnan (x)) + __builtin_unreachable(); + x = sin (x); + if (x < -1.0 || x > 1.0) + link_error (); + use (x); +} + +void bar(double x) +{ + if (!__builtin_isnan (sin (x))) + { + if (__builtin_isnan (x)) + link_error (); + if (__builtin_isinf (x)) + link_error (); + } +} + +void stool (double x) +{ + double res1 = sin (x); + double res2 = __builtin_sin (x); + if (res1 < -1.0 || res2 < -1.0) + link_error (); + if (res1 > 1.0 || res2 > 1.0) + link_error (); +} + +// { dg-final { scan-tree-dump-not "link_error" "evrp" } } -- 2.39.2 ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-18 13:12 [PATCH] Implement range-op entry for sin/cos Aldy Hernandez @ 2023-04-20 12:59 ` Jakub Jelinek 2023-04-20 13:17 ` Siddhesh Poyarekar 2023-04-21 16:40 ` Jakub Jelinek 2023-04-27 11:13 ` [PATCH] v2: " Jakub Jelinek 1 sibling, 2 replies; 25+ messages in thread From: Jakub Jelinek @ 2023-04-20 12:59 UTC (permalink / raw) To: Aldy Hernandez, Joseph S. Myers; +Cc: GCC patches, Andrew MacLeod On Tue, Apr 18, 2023 at 03:12:50PM +0200, Aldy Hernandez wrote: > [I don't know why I keep poking at floats. I must really like the pain. > Jakub, are you OK with this patch for trunk?] Thanks for working on this. Though expectedly here we are running again into the discussions we had in November about math properties of the functions vs. numeric properties in their implementations, how big maximum error shall we expect for the functions (and whether we should hardcode it for all implementations, or have some more fine-grained list of expected ulp errors for each implementation), whether the implementations at least guarantee the basic mathematical properties of the functions even if they have some errors (say for sin/cos, if they really never return > 1.0 or < -1.0) and the same questions repeated for -frounding-math, what kind of extra errors to expect when using non-default rounding and whether say sin could ever return nextafter (1.0, 2.0) or even larger value say when using non-default rounding mode. https://gcc.gnu.org/pipermail/gcc-patches/2022-November/606466.html was my attempt to get at least some numbers on some targets, I'm afraid for most implementations we aren't going to get any numerical proofs of maximum errors and the like. For sin/cos to check whether the implementation really never returns > 1.0 or < -1.0 perhaps instead of using randomized testing we could exhaustively check some range around 0, M_PI, 3*M_PI, -M_PI, -3*M_PI, and say some much larger multiples of M_PI, say 50 ulps in each direction about those points, and similarly for sin around M_PI/2 etc., in all arounding modes. > This is the range-op entry for sin/cos. It is meant to serve as an > example of what we can do for glibc math functions. It is by no means > exhaustive, just a stub to restrict the return range from sin/cos to > [-1.0, 1.0] with appropriate smarts of NANs. > > As can be seen in the testcase, we see sin() as well as > __builtin_sin() in the IL, and can resolve the resulting range > accordingly. > > gcc/ChangeLog: > > * gimple-range-op.cc (class cfn_sincos): New. > (gimple_range_op_handler::maybe_builtin_call): Add case for sin/cos. > > gcc/testsuite/ChangeLog: > > * gcc.dg/tree-ssa/range-sincos.c: New test. > --- > gcc/gimple-range-op.cc | 63 ++++++++++++++++++++ > gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c | 40 +++++++++++++ > 2 files changed, 103 insertions(+) > create mode 100644 gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c > > diff --git a/gcc/gimple-range-op.cc b/gcc/gimple-range-op.cc > index 4ca32a7b5d5..36390f2645e 100644 > --- a/gcc/gimple-range-op.cc > +++ b/gcc/gimple-range-op.cc > @@ -402,6 +402,60 @@ public: > } > } op_cfn_copysign; > > +class cfn_sincos : public range_operator_float > +{ > +public: > + using range_operator_float::fold_range; > + using range_operator_float::op1_range; > + virtual bool fold_range (frange &r, tree type, > + const frange &lh, const frange &, > + relation_trio) const final override > + { > + if (lh.known_isnan () || lh.known_isinf ()) > + { > + r.set_nan (type); > + return true; > + } > + r.set (type, dconstm1, dconst1); See above, are we sure we can use [-1., 1.] range safely, or should that be [-1.-Nulps, 1.+Nulps] for some kind of expected worse error margin of the implementation? And ditto for -frounding-math, shall we increase that interval in that case, or is [-1., 1.] going to be ok? > + if (!lh.maybe_isnan ()) This condition isn't sufficient, even if lh can't be NAN, but just may be +Inf or -Inf, the result needs to include maybe NAN. > + r.clear_nan (); > + return true; Incrementally, if we decide what to do with the maximum allowed errors in ulps, if lh's range is smaller than 2*M_PI (upper_bound () - lower_bound ()), we could narrow it down further by computing the exact values for the bounds and any local maximums or minimums in between if any and creating a range out of that. > + } > + virtual bool op1_range (frange &r, tree type, > + const frange &lhs, const frange &, > + relation_trio) const final override > + { > + if (!lhs.maybe_isnan ()) > + { > + // If NAN is not valid result, the input cannot include either > + // a NAN nor a +-INF. > + REAL_VALUE_TYPE lb = real_min_representable (type); > + REAL_VALUE_TYPE ub = real_max_representable (type); > + r.set (type, lb, ub, nan_state (false, false)); > + return true; > + } > + // A known NAN means the input is [-INF,-INF][+INF,+INF] U +-NAN, > + // which we can't currently represent. > + if (lhs.known_isnan ()) > + { > + r.set_varying (type); > + return true; > + } > + // Results outside of [-1.0, +1.0] are impossible. > + REAL_VALUE_TYPE lb = lhs.lower_bound (); > + REAL_VALUE_TYPE ub = lhs.upper_bound (); > + if (real_less (&lb, &dconstm1) > + || real_less (&dconst1, &ub)) This condition could fit on one line. > + { > + r.set_undefined (); > + return true; > + } > + > + r.set_varying (type); I'm afraid we can't do really much better for the reverse case, even for a singleton range between -1.0 and +1.0 we'd have infinite number of results. In theory we could perhaps narrow it down from the minimum/maximum representable bounds, but even 1ulp there is much larger than 2*M_PI, so not sure it is worth bothering with it. Jakub ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-20 12:59 ` Jakub Jelinek @ 2023-04-20 13:17 ` Siddhesh Poyarekar 2023-04-20 14:02 ` Jakub Jelinek 2023-04-21 16:40 ` Jakub Jelinek 1 sibling, 1 reply; 25+ messages in thread From: Siddhesh Poyarekar @ 2023-04-20 13:17 UTC (permalink / raw) To: Jakub Jelinek, Aldy Hernandez, Joseph S. Myers Cc: GCC patches, Andrew MacLeod On 2023-04-20 08:59, Jakub Jelinek via Gcc-patches wrote: >> + r.set (type, dconstm1, dconst1); > > See above, are we sure we can use [-1., 1.] range safely, or should that be > [-1.-Nulps, 1.+Nulps] for some kind of expected worse error margin of the > implementation? And ditto for -frounding-math, shall we increase that > interval in that case, or is [-1., 1.] going to be ok? Do any math implementations generate results outside of [-1., 1.]? If yes, then it's a bug in those implementations IMO, not in the range assumption. It feels wrong to cater for what ought to be trivially fixable in libraries if they ever happen to generate such results. Thanks, Sid ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-20 13:17 ` Siddhesh Poyarekar @ 2023-04-20 14:02 ` Jakub Jelinek 2023-04-20 14:20 ` Jakub Jelinek 2023-04-20 15:22 ` Siddhesh Poyarekar 0 siblings, 2 replies; 25+ messages in thread From: Jakub Jelinek @ 2023-04-20 14:02 UTC (permalink / raw) To: Siddhesh Poyarekar Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod [-- Attachment #1: Type: text/plain, Size: 1541 bytes --] On Thu, Apr 20, 2023 at 09:17:10AM -0400, Siddhesh Poyarekar wrote: > On 2023-04-20 08:59, Jakub Jelinek via Gcc-patches wrote: > > > + r.set (type, dconstm1, dconst1); > > > > See above, are we sure we can use [-1., 1.] range safely, or should that be > > [-1.-Nulps, 1.+Nulps] for some kind of expected worse error margin of the > > implementation? And ditto for -frounding-math, shall we increase that > > interval in that case, or is [-1., 1.] going to be ok? > > Do any math implementations generate results outside of [-1., 1.]? If yes, Clearly they do. > then it's a bug in those implementations IMO, not in the range assumption. > It feels wrong to cater for what ought to be trivially fixable in libraries > if they ever happen to generate such results. So, I wrote following test. On x86_64-linux with glibc 2.35, I see for i in FLOAT DOUBLE LDOUBLE FLOAT128; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o /tmp/sincos{,.c} -lm; /tmp/sincos || echo $i $j; done; done Aborted (core dumped) FLOAT UPWARD Aborted (core dumped) FLOAT DOWNWARD On sparc-sun-solaris2.11 I see for i in FLOAT DOUBLE LDOUBLE; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o sincos{,.c} -lm; ./sincos || echo $i $j; done; done Abort (core dumped) DOUBLE UPWARD Abort (core dumped) DOUBLE DOWNWARD Haven't tried anything else. So that shows (but doesn't prove) that maybe [-1., 1.] interval is fine for -fno-rounding-math on those, but not for -frounding-math. Jakub [-- Attachment #2: sincos.c --] [-- Type: text/plain, Size: 1426 bytes --] #define _GNU_SOURCE #include <math.h> #include <fenv.h> #ifdef FLOAT #define TYPE float #define SIN sinf #define COS cosf #ifdef M_PI_2f #define PI2 M_PI_2f #else #define PI2 1.570796326794896619231321691639751442f #endif #define NEXTAFTER nextafterf #elif defined DOUBLE #define TYPE double #define SIN sin #define COS cos #ifdef M_PI_2 #define PI2 M_PI_2 #else #define PI2 1.570796326794896619231321691639751442f #endif #define NEXTAFTER nextafter #elif defined LDOUBLE #define TYPE long double #define SIN sinl #define COS cosl #ifdef M_PI_2l #define PI2 M_PI_2l #else #define PI2 1.570796326794896619231321691639751442f #endif #define NEXTAFTER nextafterl #elif defined FLOAT128 #define TYPE _Float128 #define SIN sinf128 #define COS cosf128 #ifdef #define PI2 M_PI_2f128 #else #define PI2 1.570796326794896619231321691639751442f #endif #define NEXTAFTER nextafterf128 #endif int main () { #ifdef ROUND fesetround (ROUND); #endif for (int i = -1024; i <= 1024; i++) for (int j = -1; j <= 1; j += 2) { TYPE val = ((TYPE) i) * PI2; TYPE inf = j * __builtin_inf (); for (int k = 0; k < 1000; k++) { TYPE res = SIN (val); if (res < (TYPE) -1.0 || res > (TYPE) 1.0) __builtin_abort (); res = COS (val); if (res < (TYPE) -1.0 || res > (TYPE) 1.0) __builtin_abort (); val = NEXTAFTER (val, inf); } } } ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-20 14:02 ` Jakub Jelinek @ 2023-04-20 14:20 ` Jakub Jelinek 2023-04-20 15:22 ` Siddhesh Poyarekar 1 sibling, 0 replies; 25+ messages in thread From: Jakub Jelinek @ 2023-04-20 14:20 UTC (permalink / raw) To: Siddhesh Poyarekar, Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod [-- Attachment #1: Type: text/plain, Size: 1512 bytes --] On Thu, Apr 20, 2023 at 04:02:02PM +0200, Jakub Jelinek via Gcc-patches wrote: > So, I wrote following test. Slightly adjusted to see more info: x86_64-linux glibc 2.35: for i in FLOAT DOUBLE LDOUBLE FLOAT128; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o /tmp/sincos{,.c} -lm; /tmp/sincos || echo $i $j; done; done sin -0x1.2d97c800000000000000p+2 0x1.00000200000000000000p+0 sin -0x1.2d97c800000000000000p+2 0x1.00000200000000000000p+0 sin 0x1.f9cbe200000000000000p+7 0x1.00000200000000000000p+0 FLOAT UPWARD cos -0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0 sin -0x1.f9cbe200000000000000p+7 -0x1.00000200000000000000p+0 sin 0x1.2d97c800000000000000p+2 -0x1.00000200000000000000p+0 sin 0x1.2d97c800000000000000p+2 -0x1.00000200000000000000p+0 cos 0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0 cos 0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0 FLOAT DOWNWARD sparc-sun-solaris2.11 results are too large to post in detail, but are sin -0x1.2d97c7f3321d20000000p+2 0x1.00000000000010000000p+0 ... sin 0x1.f6a7a29553c450000000p+2 0x1.00000000000010000000p+0 DOUBLE UPWARD sin -0x1.f6a7a2955385e0000000p+2 -0x1.00000000000010000000p+0 ... sin 0x1.2d97c7f3325b90000000p+2 -0x1.00000000000010000000p+0 DOUBLE DOWNWARD where all the DOUBLE UPWARD values have 0x1.00000000000010000000p+0 results and all DOUBLE DOWNWARD values -0x1.00000000000010000000p+0. So, I think that is 1ulp in all cases in both directions for -frounding-math. Jakub [-- Attachment #2: sincos.c --] [-- Type: text/plain, Size: 1753 bytes --] #define _GNU_SOURCE #include <math.h> #include <fenv.h> #include <stdio.h> #include <stdlib.h> #ifdef FLOAT #define TYPE float #define SIN sinf #define COS cosf #ifdef M_PI_2f #define PI2 M_PI_2f #else #define PI2 1.570796326794896619231321691639751442f #endif #define PRINT(str) printf ("%s %.20a %.20a\n", str, val, res) #define NEXTAFTER nextafterf #elif defined DOUBLE #define TYPE double #define SIN sin #define COS cos #ifdef M_PI_2 #define PI2 M_PI_2 #else #define PI2 1.570796326794896619231321691639751442f #endif #define NEXTAFTER nextafter #define PRINT(str) printf ("%s %.20a %.20a\n", str, val, res) #elif defined LDOUBLE #define TYPE long double #define SIN sinl #define COS cosl #ifdef M_PI_2l #define PI2 M_PI_2l #else #define PI2 1.570796326794896619231321691639751442f #endif #define NEXTAFTER nextafterl #define PRINT(str) printf ("%s %.20La %.20La\n", str, val, res) #elif defined FLOAT128 #define TYPE _Float128 #define SIN sinf128 #define COS cosf128 #ifdef M_PI_2f128 #define PI2 M_PI_2f128 #else #define PI2 1.570796326794896619231321691639751442f #endif #define NEXTAFTER nextafterf128 #define PRINT(str) __builtin_abort () #endif int main () { int ret = 0; #ifdef ROUND fesetround (ROUND); #endif for (int i = -1024; i <= 1024; i++) for (int j = -1; j <= 1; j += 2) { TYPE val = ((TYPE) i) * PI2; TYPE inf = j * __builtin_inf (); for (int k = 0; k < 1000; k++) { TYPE res = SIN (val); if (res < (TYPE) -1.0 || res > (TYPE) 1.0) { PRINT ("sin"); ret = 1; } res = COS (val); if (res < (TYPE) -1.0 || res > (TYPE) 1.0) { PRINT ("cos"); ret = 1; } val = NEXTAFTER (val, inf); } } return ret; } ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-20 14:02 ` Jakub Jelinek 2023-04-20 14:20 ` Jakub Jelinek @ 2023-04-20 15:22 ` Siddhesh Poyarekar 2023-04-20 15:52 ` Jakub Jelinek 1 sibling, 1 reply; 25+ messages in thread From: Siddhesh Poyarekar @ 2023-04-20 15:22 UTC (permalink / raw) To: Jakub Jelinek Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On 2023-04-20 10:02, Jakub Jelinek wrote: > On x86_64-linux with glibc 2.35, I see > for i in FLOAT DOUBLE LDOUBLE FLOAT128; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o /tmp/sincos{,.c} -lm; /tmp/sincos || echo $i $j; done; done > Aborted (core dumped) > FLOAT UPWARD > Aborted (core dumped) > FLOAT DOWNWARD > On sparc-sun-solaris2.11 I see > for i in FLOAT DOUBLE LDOUBLE; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o sincos{,.c} -lm; ./sincos || echo $i $j; done; done > Abort (core dumped) > DOUBLE UPWARD > Abort (core dumped) > DOUBLE DOWNWARD > Haven't tried anything else. So that shows (but doesn't prove) that > maybe [-1., 1.] interval is fine for -fno-rounding-math on those, but not > for -frounding-math. Would there be a reason to not consider these as bugs? I feel like these should be fixed in glibc, or any math implementation that ends up doing this. I suppose one reason could be the overhead of an additional branch to check for result bounds, but is that serious enough to allow this imprecision? The alternative of output range being defined as [-1.0-ulp, 1.0+ulp] avoids that conversation I guess. Thanks, Sid ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-20 15:22 ` Siddhesh Poyarekar @ 2023-04-20 15:52 ` Jakub Jelinek 2023-04-20 17:57 ` Siddhesh Poyarekar 0 siblings, 1 reply; 25+ messages in thread From: Jakub Jelinek @ 2023-04-20 15:52 UTC (permalink / raw) To: Siddhesh Poyarekar Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On Thu, Apr 20, 2023 at 11:22:24AM -0400, Siddhesh Poyarekar wrote: > On 2023-04-20 10:02, Jakub Jelinek wrote: > > On x86_64-linux with glibc 2.35, I see > > for i in FLOAT DOUBLE LDOUBLE FLOAT128; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o /tmp/sincos{,.c} -lm; /tmp/sincos || echo $i $j; done; done > > Aborted (core dumped) > > FLOAT UPWARD > > Aborted (core dumped) > > FLOAT DOWNWARD > > On sparc-sun-solaris2.11 I see > > for i in FLOAT DOUBLE LDOUBLE; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o sincos{,.c} -lm; ./sincos || echo $i $j; done; done > > Abort (core dumped) > > DOUBLE UPWARD > > Abort (core dumped) > > DOUBLE DOWNWARD > > Haven't tried anything else. So that shows (but doesn't prove) that > > maybe [-1., 1.] interval is fine for -fno-rounding-math on those, but not > > for -frounding-math. > > Would there be a reason to not consider these as bugs? I feel like these > should be fixed in glibc, or any math implementation that ends up doing > this. Why? Unless an implementation guarantees <= 0.5ulps errors, it can be one or more ulps off, why is an error at or near 1.0 or -1.0 error any worse than similar errors for other values? Similarly for other functions which have other ranges, perhaps not with so nice round numbers. Say asin has [-pi/2, pi/2] range, those numbers aren't exactly representable, but is it any worse to round those values to -inf or +inf or worse give something 1-5 ulps further from that interval comparing to other 1-5ulps errors? Jakub ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-20 15:52 ` Jakub Jelinek @ 2023-04-20 17:57 ` Siddhesh Poyarekar 2023-04-21 1:14 ` Siddhesh Poyarekar 2023-04-24 16:03 ` Jakub Jelinek 0 siblings, 2 replies; 25+ messages in thread From: Siddhesh Poyarekar @ 2023-04-20 17:57 UTC (permalink / raw) To: Jakub Jelinek Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On 2023-04-20 11:52, Jakub Jelinek wrote: > Why? Unless an implementation guarantees <= 0.5ulps errors, it can be one > or more ulps off, why is an error at or near 1.0 or -1.0 error any worse > than similar errors for other values? In a general sense, maybe not, but in the sense of breaching the bounds of admissible values, especially when it can be reasonably corrected, it seems worse IMO to let the error slide. > Similarly for other functions which have other ranges, perhaps not with so > nice round numbers. Say asin has [-pi/2, pi/2] range, those numbers aren't > exactly representable, but is it any worse to round those values to -inf or > +inf or worse give something 1-5 ulps further from that interval comparing > to other 1-5ulps errors? I agree the argument in favour of allowing errors breaching the mathematical bounds is certainly stronger for bounds that are not exactly representable. I just feel like the implementation ought to take the additional effort when the bounds are representable and make it easier for the compiler. For bounds that aren't representable, one could get error bounds from libm-test-ulps data in glibc, although I reckon those won't be exhaustive. From a quick peek at the sin/cos data, the arc target seems to be among the worst performers at about 7ulps, although if you include the complex routines we get close to 13 ulps. The very worst imprecision among all math routines (that's gamma) is at 16 ulps for power in glibc tests, so maybe allowing about 25-30 ulps error in bounds might work across the board. But yeah, it's probably going to be guesswork. Thanks, Sid ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-20 17:57 ` Siddhesh Poyarekar @ 2023-04-21 1:14 ` Siddhesh Poyarekar 2023-04-21 6:52 ` Jakub Jelinek 2023-04-24 16:03 ` Jakub Jelinek 1 sibling, 1 reply; 25+ messages in thread From: Siddhesh Poyarekar @ 2023-04-21 1:14 UTC (permalink / raw) To: Jakub Jelinek Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On 2023-04-20 13:57, Siddhesh Poyarekar wrote: > For bounds that aren't representable, one could get error bounds from > libm-test-ulps data in glibc, although I reckon those won't be > exhaustive. From a quick peek at the sin/cos data, the arc target seems > to be among the worst performers at about 7ulps, although if you include > the complex routines we get close to 13 ulps. The very worst > imprecision among all math routines (that's gamma) is at 16 ulps for > power in glibc tests, so maybe allowing about 25-30 ulps error in bounds > might work across the board. I was thinking about this a bit more and it seems like limiting ranges to targets that can generate sane results (i.e. error bounds within, say, 5-6 ulps) and for the rest, avoid emitting the ranges altogether. Emitting a bad range for all architectures seems like a net worse solution again. Thanks, Sid ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-21 1:14 ` Siddhesh Poyarekar @ 2023-04-21 6:52 ` Jakub Jelinek 2023-04-21 11:20 ` Siddhesh Poyarekar 2023-04-25 8:59 ` Aldy Hernandez 0 siblings, 2 replies; 25+ messages in thread From: Jakub Jelinek @ 2023-04-21 6:52 UTC (permalink / raw) To: Siddhesh Poyarekar Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On Thu, Apr 20, 2023 at 09:14:10PM -0400, Siddhesh Poyarekar wrote: > On 2023-04-20 13:57, Siddhesh Poyarekar wrote: > > For bounds that aren't representable, one could get error bounds from > > libm-test-ulps data in glibc, although I reckon those won't be > > exhaustive. From a quick peek at the sin/cos data, the arc target seems > > to be among the worst performers at about 7ulps, although if you include > > the complex routines we get close to 13 ulps. The very worst > > imprecision among all math routines (that's gamma) is at 16 ulps for > > power in glibc tests, so maybe allowing about 25-30 ulps error in bounds > > might work across the board. > > I was thinking about this a bit more and it seems like limiting ranges to > targets that can generate sane results (i.e. error bounds within, say, 5-6 > ulps) and for the rest, avoid emitting the ranges altogether. Emitting a bad > range for all architectures seems like a net worse solution again. Well, at least for basic arithmetics when libm functions aren't involved, there is no point in disabling ranges altogether. And, for libm functions, my plan was to introduce a target hook, which would have combined_fn argument to tell which function is queried, machine_mode to say which floating point format and perhaps a bool whether it is ulps for these basic math boundaries or results somewhere in between, and would return in unsigned int ulps, 0 for 0.5ulps precision. So, we could say for CASE_CFN_SIN: CASE_CFN_COS: in the glibc handler say that ulps is say 3 inside of the ranges and 0 on the boundaries if !flag_rounding_math and 6 and 2 with flag_rounding_math or whatever. And in the generic code don't assume anything if ulps is say 100 or more. The hooks would need to be a union of precision of supported versions of the library through the history, including say libmvec because function calls could be vectorized. And default could be that infinite precision. Back in November I've posted a proglet that can generate some ulps from random number testing, plus on glibc we could pick maximums from ulps files. And if needed, say powerpc*-linux could override the generic glibc version for some subset of functions and call default otherwise (say at least for __ibm128). Jakub ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-21 6:52 ` Jakub Jelinek @ 2023-04-21 11:20 ` Siddhesh Poyarekar 2023-04-25 8:59 ` Aldy Hernandez 1 sibling, 0 replies; 25+ messages in thread From: Siddhesh Poyarekar @ 2023-04-21 11:20 UTC (permalink / raw) To: Jakub Jelinek Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On 2023-04-21 02:52, Jakub Jelinek wrote: > On Thu, Apr 20, 2023 at 09:14:10PM -0400, Siddhesh Poyarekar wrote: >> On 2023-04-20 13:57, Siddhesh Poyarekar wrote: >>> For bounds that aren't representable, one could get error bounds from >>> libm-test-ulps data in glibc, although I reckon those won't be >>> exhaustive. From a quick peek at the sin/cos data, the arc target seems >>> to be among the worst performers at about 7ulps, although if you include >>> the complex routines we get close to 13 ulps. The very worst >>> imprecision among all math routines (that's gamma) is at 16 ulps for >>> power in glibc tests, so maybe allowing about 25-30 ulps error in bounds >>> might work across the board. >> >> I was thinking about this a bit more and it seems like limiting ranges to >> targets that can generate sane results (i.e. error bounds within, say, 5-6 >> ulps) and for the rest, avoid emitting the ranges altogether. Emitting a bad >> range for all architectures seems like a net worse solution again. > > Well, at least for basic arithmetics when libm functions aren't involved, > there is no point in disabling ranges altogether. Oh yeah, I did mean only franges for math function call results. > And, for libm functions, my plan was to introduce a target hook, which > would have combined_fn argument to tell which function is queried, > machine_mode to say which floating point format and perhaps a bool whether > it is ulps for these basic math boundaries or results somewhere in between, > and would return in unsigned int ulps, 0 for 0.5ulps precision. > So, we could say for CASE_CFN_SIN: CASE_CFN_COS: in the glibc handler > say that ulps is say 3 inside of the ranges and 0 on the boundaries if > !flag_rounding_math and 6 and 2 with flag_rounding_math or whatever. > And in the generic code don't assume anything if ulps is say 100 or more. > The hooks would need to be a union of precision of supported versions of > the library through the history, including say libmvec because function > calls could be vectorized. > And default could be that infinite precision. > Back in November I've posted a proglet that can generate some ulps from > random number testing, plus on glibc we could pick maximums from ulps files. > And if needed, say powerpc*-linux could override the generic glibc > version for some subset of functions and call default otherwise (say at > least for __ibm128). Ack, that sounds like a plan. Thanks, Sid ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-21 6:52 ` Jakub Jelinek 2023-04-21 11:20 ` Siddhesh Poyarekar @ 2023-04-25 8:59 ` Aldy Hernandez 1 sibling, 0 replies; 25+ messages in thread From: Aldy Hernandez @ 2023-04-25 8:59 UTC (permalink / raw) To: Jakub Jelinek, Siddhesh Poyarekar Cc: Joseph S. Myers, GCC patches, Andrew MacLeod On 4/21/23 08:52, Jakub Jelinek wrote: > On Thu, Apr 20, 2023 at 09:14:10PM -0400, Siddhesh Poyarekar wrote: >> On 2023-04-20 13:57, Siddhesh Poyarekar wrote: >>> For bounds that aren't representable, one could get error bounds from >>> libm-test-ulps data in glibc, although I reckon those won't be >>> exhaustive. From a quick peek at the sin/cos data, the arc target seems >>> to be among the worst performers at about 7ulps, although if you include >>> the complex routines we get close to 13 ulps. The very worst >>> imprecision among all math routines (that's gamma) is at 16 ulps for >>> power in glibc tests, so maybe allowing about 25-30 ulps error in bounds >>> might work across the board. >> >> I was thinking about this a bit more and it seems like limiting ranges to >> targets that can generate sane results (i.e. error bounds within, say, 5-6 >> ulps) and for the rest, avoid emitting the ranges altogether. Emitting a bad >> range for all architectures seems like a net worse solution again. > > Well, at least for basic arithmetics when libm functions aren't involved, > there is no point in disabling ranges altogether. > And, for libm functions, my plan was to introduce a target hook, which > would have combined_fn argument to tell which function is queried, > machine_mode to say which floating point format and perhaps a bool whether > it is ulps for these basic math boundaries or results somewhere in between, > and would return in unsigned int ulps, 0 for 0.5ulps precision. I wonder if we could export frange_nextafter() to take a final argument for the number of ulps, making it possible to extend in either direction by a number of ulps. And perhaps add a generic function that calls the target hook and extends the range accordingly: extern int TARGET_ULP_ERROR (combined_fn, machine_mode); extern void frange_nextafter (machine_mode, REAL_VALUE_TYPE &value, const REAL_VALUE_TYPE &direction, int ulps); void adjust_frange_for_op (frange &r, combined_fn fn) { ... int ulps = TARGET_ULP_ERROR (fn, mode); frange_nextafter (mode, lb, frange_val_min (type), ulps); frange_nextafter (mode, ub, frange_val_max (type), ulps); ... } bool cfn_sincos::fold_range (...) { ... if (cond) { r.set (type, dconstm1, const1); adjust_frange_for_op (r, m_cfn); } ... } Or if adjusting by ULPS is a common enough idiom, we could promote it to an frange method: void frange::adjust_ulps (int ulps); Would that work, or did you have something else in mind? Aldy ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-20 17:57 ` Siddhesh Poyarekar 2023-04-21 1:14 ` Siddhesh Poyarekar @ 2023-04-24 16:03 ` Jakub Jelinek 2023-04-24 16:05 ` Siddhesh Poyarekar 1 sibling, 1 reply; 25+ messages in thread From: Jakub Jelinek @ 2023-04-24 16:03 UTC (permalink / raw) To: Siddhesh Poyarekar Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod [-- Attachment #1: Type: text/plain, Size: 1349 bytes --] On Thu, Apr 20, 2023 at 01:57:55PM -0400, Siddhesh Poyarekar wrote: > > Similarly for other functions which have other ranges, perhaps not with so > > nice round numbers. Say asin has [-pi/2, pi/2] range, those numbers aren't > > exactly representable, but is it any worse to round those values to -inf or > > +inf or worse give something 1-5 ulps further from that interval comparing > > to other 1-5ulps errors? I've extended my hack^^^test to include sqrt and this time it seems that the boundary in that case holds even for non-standard rounding modes for glibc: for i in FLOAT DOUBLE LDOUBLE FLOAT128; do for j in TONEAREST UPWARD DOWNWARD TOWARDZERO; do gcc -D$i -DROUND=FE_$j -g -O1 -o /tmp/sincossqrt{,.c} -lm; /tmp/sincossqrt || echo $i $j; done; done sin -0x1.2d97c800000000000000p+2 0x1.00000200000000000000p+0 sin -0x1.2d97c800000000000000p+2 0x1.00000200000000000000p+0 sin 0x1.f9cbe200000000000000p+7 0x1.00000200000000000000p+0 FLOAT UPWARD cos -0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0 sin -0x1.f9cbe200000000000000p+7 -0x1.00000200000000000000p+0 sin 0x1.2d97c800000000000000p+2 -0x1.00000200000000000000p+0 sin 0x1.2d97c800000000000000p+2 -0x1.00000200000000000000p+0 cos 0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0 cos 0x1.f9cbe200000000000000p+8 -0x1.00000200000000000000p+0 FLOAT DOWNWARD Jakub [-- Attachment #2: sincossqrt.c --] [-- Type: text/plain, Size: 2124 bytes --] #define _GNU_SOURCE #include <math.h> #include <fenv.h> #include <stdio.h> #include <stdlib.h> #ifdef FLOAT #define TYPE float #define SIN sinf #define COS cosf #define SQRT sqrtf #ifdef M_PI_2f #define PI2 M_PI_2f #else #define PI2 1.570796326794896619231321691639751442f #endif #define PRINT(str) printf ("%s %.20a %.20a\n", str, val, res) #define NEXTAFTER nextafterf #elif defined DOUBLE #define TYPE double #define SIN sin #define COS cos #define SQRT sqrt #ifdef M_PI_2 #define PI2 M_PI_2 #else #define PI2 1.570796326794896619231321691639751442f #endif #define NEXTAFTER nextafter #define PRINT(str) printf ("%s %.20a %.20a\n", str, val, res) #elif defined LDOUBLE #define TYPE long double #define SIN sinl #define COS cosl #define SQRT sqrtl #ifdef M_PI_2l #define PI2 M_PI_2l #else #define PI2 1.570796326794896619231321691639751442f #endif #define NEXTAFTER nextafterl #define PRINT(str) printf ("%s %.20La %.20La\n", str, val, res) #elif defined FLOAT128 #define TYPE _Float128 #define SIN sinf128 #define COS cosf128 #define SQRT sqrtf128 #ifdef M_PI_2f128 #define PI2 M_PI_2f128 #else #define PI2 1.570796326794896619231321691639751442f #endif #define NEXTAFTER nextafterf128 #define PRINT(str) __builtin_abort () #endif int main () { int ret = 0; #ifdef ROUND fesetround (ROUND); #endif for (int i = -1024; i <= 1024; i++) for (int j = -1; j <= 1; j += 2) { TYPE val = ((TYPE) i) * PI2; TYPE inf = j * __builtin_inf (); for (int k = 0; k < 1000; k++) { TYPE res = SIN (val); if (res < (TYPE) -1.0 || res > (TYPE) 1.0) { PRINT ("sin"); ret = 1; } res = COS (val); if (res < (TYPE) -1.0 || res > (TYPE) 1.0) { PRINT ("cos"); ret = 1; } val = NEXTAFTER (val, inf); } } for (int j = -1; j <= 1; j += 2) { TYPE val = 0; TYPE inf = j * __builtin_inf (); if (j < 0) val = -val; for (int k = 0; k < 1000; k++) { TYPE res = SQRT (val); if (res < (TYPE) -0.0) { PRINT ("sqrt"); ret = 1; } } } return ret; } ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-24 16:03 ` Jakub Jelinek @ 2023-04-24 16:05 ` Siddhesh Poyarekar 2023-04-24 16:09 ` Jakub Jelinek 2023-04-24 16:33 ` Jeff Law 0 siblings, 2 replies; 25+ messages in thread From: Siddhesh Poyarekar @ 2023-04-24 16:05 UTC (permalink / raw) To: Jakub Jelinek Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On 2023-04-24 12:03, Jakub Jelinek wrote: > On Thu, Apr 20, 2023 at 01:57:55PM -0400, Siddhesh Poyarekar wrote: >>> Similarly for other functions which have other ranges, perhaps not with so >>> nice round numbers. Say asin has [-pi/2, pi/2] range, those numbers aren't >>> exactly representable, but is it any worse to round those values to -inf or >>> +inf or worse give something 1-5 ulps further from that interval comparing >>> to other 1-5ulps errors? > > I've extended my hack^^^test to include sqrt and this time it seems > that the boundary in that case holds even for non-standard rounding modes > for glibc: IIRC the standard _requires_ sqrt to be correctly rounded. Sid ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-24 16:05 ` Siddhesh Poyarekar @ 2023-04-24 16:09 ` Jakub Jelinek 2023-04-24 16:33 ` Jeff Law 1 sibling, 0 replies; 25+ messages in thread From: Jakub Jelinek @ 2023-04-24 16:09 UTC (permalink / raw) To: Siddhesh Poyarekar Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On Mon, Apr 24, 2023 at 12:05:43PM -0400, Siddhesh Poyarekar wrote: > On 2023-04-24 12:03, Jakub Jelinek wrote: > > On Thu, Apr 20, 2023 at 01:57:55PM -0400, Siddhesh Poyarekar wrote: > > > > Similarly for other functions which have other ranges, perhaps not with so > > > > nice round numbers. Say asin has [-pi/2, pi/2] range, those numbers aren't > > > > exactly representable, but is it any worse to round those values to -inf or > > > > +inf or worse give something 1-5 ulps further from that interval comparing > > > > to other 1-5ulps errors? > > > > I've extended my hack^^^test to include sqrt and this time it seems > > that the boundary in that case holds even for non-standard rounding modes > > for glibc: > > IIRC the standard _requires_ sqrt to be correctly rounded. Well, we still need to make some effort to verify it is the case. BTW, what my hacks didn't check yet and I should eventually is libmvec, because if ranger makes some assumptions on libm calls and then we vectorize it, it will be libmvec rather than libm. Jakub ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-24 16:05 ` Siddhesh Poyarekar 2023-04-24 16:09 ` Jakub Jelinek @ 2023-04-24 16:33 ` Jeff Law 1 sibling, 0 replies; 25+ messages in thread From: Jeff Law @ 2023-04-24 16:33 UTC (permalink / raw) To: Siddhesh Poyarekar, Jakub Jelinek Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On 4/24/23 10:05, Siddhesh Poyarekar wrote: > On 2023-04-24 12:03, Jakub Jelinek wrote: >> On Thu, Apr 20, 2023 at 01:57:55PM -0400, Siddhesh Poyarekar wrote: >>>> Similarly for other functions which have other ranges, perhaps not >>>> with so >>>> nice round numbers. Say asin has [-pi/2, pi/2] range, those numbers >>>> aren't >>>> exactly representable, but is it any worse to round those values to >>>> -inf or >>>> +inf or worse give something 1-5 ulps further from that interval >>>> comparing >>>> to other 1-5ulps errors? >> >> I've extended my hack^^^test to include sqrt and this time it seems >> that the boundary in that case holds even for non-standard rounding modes >> for glibc: > > IIRC the standard _requires_ sqrt to be correctly rounded. Correct. I spent an inordinate amount of time dealing with this in the past ;-) jeff ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-20 12:59 ` Jakub Jelinek 2023-04-20 13:17 ` Siddhesh Poyarekar @ 2023-04-21 16:40 ` Jakub Jelinek 2023-04-21 20:43 ` Mikael Morin 2023-04-25 9:08 ` Aldy Hernandez 1 sibling, 2 replies; 25+ messages in thread From: Jakub Jelinek @ 2023-04-21 16:40 UTC (permalink / raw) To: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On Thu, Apr 20, 2023 at 02:59:35PM +0200, Jakub Jelinek via Gcc-patches wrote: > Thanks for working on this. Though expectedly here we are running > again into the discussions we had in November about math properties of the > functions vs. numeric properties in their implementations, how big maximum > error shall we expect for the functions (and whether we should hardcode > it for all implementations, or have some more fine-grained list of expected > ulp errors for each implementation), whether the implementations at least > guarantee the basic mathematical properties of the functions even if they > have some errors (say for sin/cos, if they really never return > 1.0 or < > -1.0) and the same questions repeated for -frounding-math, what kind of > extra errors to expect when using non-default rounding and whether say sin > could ever return nextafter (1.0, 2.0) or even larger value say when > using non-default rounding mode. > https://gcc.gnu.org/pipermail/gcc-patches/2022-November/606466.html > was my attempt to get at least some numbers on some targets, I'm afraid > for most implementations we aren't going to get any numerical proofs of > maximum errors and the like. For sin/cos to check whether the implementation > really never returns > 1.0 or < -1.0 perhaps instead of using randomized > testing we could exhaustively check some range around 0, M_PI, 3*M_PI, > -M_PI, -3*M_PI, and say some much larger multiples of M_PI, say 50 ulps > in each direction about those points, and similarly for sin around M_PI/2 > etc., in all arounding modes. Appart from the ulps ranges which I plan to introduce a target hook next week... > > + if (!lh.maybe_isnan ()) > > This condition isn't sufficient, even if lh can't be NAN, but just > may be +Inf or -Inf, the result needs to include maybe NAN. I've incorporated my other comments into a patch form. I also found missing undefined_p () checks in two spots, the r.set_undefined (); stuff being misplaced (it was done in the lhs.maybe_isnan () case, which is incorrect, if the lhs may be nan, then even if the finite range is say [-30., -10.] or [1.5., 42.], result shouldn't be invalid, because the result still could be NAN and in that case operand of say +-Inf or NAN would be valid. Actually thinking about it some more, perhaps we should do that check before the maybe_isnan () check and if we find the impossible finite, either use r.set_undefined (); of !maybe_isnan () or handle it like known_isnan () otherwise. Also, I think we should remember if it is SIN or COS, we'll need it both for the ulps case and if we improve it for (ub-lb) < 2*PI ranges. Now, talking about that, I'd very much like to avoid finding if some multiple of PI/2 occurs inside of such ranges, the precision of real.cc is clearly not sufficient for that, but perhaps we could use derivative of sin (cos) or of cos (sin) to see if the function on the boundary is increasing or decreasing and from that on both boundaries and their approximate distance find out if the range needs to be extended to +1 or -1. So, just incremental WIP so far... --- gcc/gimple-range-op.cc.jj 2023-04-21 17:09:48.250367999 +0200 +++ gcc/gimple-range-op.cc 2023-04-21 18:37:26.048325391 +0200 @@ -405,17 +405,20 @@ class cfn_sincos : public range_operator public: using range_operator_float::fold_range; using range_operator_float::op1_range; + cfn_sincos (combined_fn cfn) { m_cfn = cfn; } virtual bool fold_range (frange &r, tree type, const frange &lh, const frange &, relation_trio) const final override { + if (lh.undefined_p ()) + return false; if (lh.known_isnan () || lh.known_isinf ()) { r.set_nan (type); return true; } r.set (type, dconstm1, dconst1); - if (!lh.maybe_isnan ()) + if (!lh.maybe_isnan () && !lh.maybe_isinf ()) r.clear_nan (); return true; } @@ -423,15 +426,9 @@ public: const frange &lhs, const frange &, relation_trio) const final override { - if (!lhs.maybe_isnan ()) - { - // If NAN is not valid result, the input cannot include either - // a NAN nor a +-INF. - REAL_VALUE_TYPE lb = real_min_representable (type); - REAL_VALUE_TYPE ub = real_max_representable (type); - r.set (type, lb, ub, nan_state (false, false)); - return true; - } + if (lhs.undefined_p ()) + return false; + // A known NAN means the input is [-INF,-INF][+INF,+INF] U +-NAN, // which we can't currently represent. if (lhs.known_isnan ()) @@ -439,20 +436,38 @@ public: r.set_varying (type); return true; } + // Results outside of [-1.0, +1.0] are impossible. REAL_VALUE_TYPE lb = lhs.lower_bound (); REAL_VALUE_TYPE ub = lhs.upper_bound (); - if (real_less (&lb, &dconstm1) - || real_less (&dconst1, &ub)) + if (real_less (&lb, &dconstm1) || real_less (&dconst1, &ub)) { - r.set_undefined (); + if (!lhs.maybe_isnan ()) + r.set_undefined (); + else + /* If lhs could be NAN and finite result is impossible, + the range is like lhs.known_isnan () above, + [-INF,-INF][+INF,+INF] U +-NAN. */ + r.set_varying (type); + return true; + } + + if (!lhs.maybe_isnan ()) + { + // If NAN is not valid result, the input cannot include either + // a NAN nor a +-INF. + lb = real_min_representable (type); + ub = real_max_representable (type); + r.set (type, lb, ub, nan_state (false, false)); return true; } r.set_varying (type); return true; } -} op_cfn_sincos; +private: + combined_fn m_cfn; +} op_cfn_sin (CFN_SIN), op_cfn_cos (CFN_COS); // Implement range operator for CFN_BUILT_IN_TOUPPER and CFN_BUILT_IN_TOLOWER. class cfn_toupper_tolower : public range_operator @@ -934,10 +949,15 @@ gimple_range_op_handler::maybe_builtin_c CASE_CFN_SIN: CASE_CFN_SIN_FN: + m_op1 = gimple_call_arg (call, 0); + m_float = &op_cfn_sin; + m_valid = true; + break; + CASE_CFN_COS: CASE_CFN_COS_FN: m_op1 = gimple_call_arg (call, 0); - m_float = &op_cfn_sincos; + m_float = &op_cfn_cos; m_valid = true; break; Jakub ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-21 16:40 ` Jakub Jelinek @ 2023-04-21 20:43 ` Mikael Morin 2023-04-21 20:45 ` Jakub Jelinek 2023-04-25 9:10 ` Aldy Hernandez 2023-04-25 9:08 ` Aldy Hernandez 1 sibling, 2 replies; 25+ messages in thread From: Mikael Morin @ 2023-04-21 20:43 UTC (permalink / raw) To: Jakub Jelinek, Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod Hello, > --- gcc/gimple-range-op.cc.jj 2023-04-21 17:09:48.250367999 +0200 > +++ gcc/gimple-range-op.cc 2023-04-21 18:37:26.048325391 +0200 > @@ -439,20 +436,38 @@ public: > r.set_varying (type); > return true; > } > + > // Results outside of [-1.0, +1.0] are impossible. > REAL_VALUE_TYPE lb = lhs.lower_bound (); > REAL_VALUE_TYPE ub = lhs.upper_bound (); > - if (real_less (&lb, &dconstm1) > - || real_less (&dconst1, &ub)) > + if (real_less (&lb, &dconstm1) || real_less (&dconst1, &ub)) > { Shouldn't lb and ub be swapped in this condition? If I understand correctly, we are looking for ranges like [whatever,x] where x < -1.0 or [y, whatever] where 1.0 < y. Mikael ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-21 20:43 ` Mikael Morin @ 2023-04-21 20:45 ` Jakub Jelinek 2023-04-25 9:10 ` Aldy Hernandez 1 sibling, 0 replies; 25+ messages in thread From: Jakub Jelinek @ 2023-04-21 20:45 UTC (permalink / raw) To: Mikael Morin; +Cc: Aldy Hernandez, Joseph S. Myers, GCC patches, Andrew MacLeod On Fri, Apr 21, 2023 at 10:43:44PM +0200, Mikael Morin wrote: > Hello, > > > --- gcc/gimple-range-op.cc.jj 2023-04-21 17:09:48.250367999 +0200 > > +++ gcc/gimple-range-op.cc 2023-04-21 18:37:26.048325391 +0200 > > @@ -439,20 +436,38 @@ public: > > r.set_varying (type); > > return true; > > } > > + > > // Results outside of [-1.0, +1.0] are impossible. > > REAL_VALUE_TYPE lb = lhs.lower_bound (); > > REAL_VALUE_TYPE ub = lhs.upper_bound (); > > - if (real_less (&lb, &dconstm1) > > - || real_less (&dconst1, &ub)) > > + if (real_less (&lb, &dconstm1) || real_less (&dconst1, &ub)) > > { > > Shouldn't lb and ub be swapped in this condition? Yes, they should. > If I understand correctly, we are looking for ranges like [whatever,x] where > x < -1.0 or [y, whatever] where 1.0 < y. Jakub ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-21 20:43 ` Mikael Morin 2023-04-21 20:45 ` Jakub Jelinek @ 2023-04-25 9:10 ` Aldy Hernandez 1 sibling, 0 replies; 25+ messages in thread From: Aldy Hernandez @ 2023-04-25 9:10 UTC (permalink / raw) To: Mikael Morin, Jakub Jelinek, Joseph S. Myers, GCC patches, Andrew MacLeod On 4/21/23 22:43, Mikael Morin wrote: > Hello, > >> --- gcc/gimple-range-op.cc.jj 2023-04-21 17:09:48.250367999 +0200 >> +++ gcc/gimple-range-op.cc 2023-04-21 18:37:26.048325391 +0200 >> @@ -439,20 +436,38 @@ public: >> r.set_varying (type); >> return true; >> } >> + >> // Results outside of [-1.0, +1.0] are impossible. >> REAL_VALUE_TYPE lb = lhs.lower_bound (); >> REAL_VALUE_TYPE ub = lhs.upper_bound (); >> - if (real_less (&lb, &dconstm1) >> - || real_less (&dconst1, &ub)) >> + if (real_less (&lb, &dconstm1) || real_less (&dconst1, &ub)) >> { > > Shouldn't lb and ub be swapped in this condition? > If I understand correctly, we are looking for ranges like [whatever,x] > where x < -1.0 or [y, whatever] where 1.0 < y. > > Mikael > Thanks. Merged into the latest revision. Aldy ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] Implement range-op entry for sin/cos. 2023-04-21 16:40 ` Jakub Jelinek 2023-04-21 20:43 ` Mikael Morin @ 2023-04-25 9:08 ` Aldy Hernandez 1 sibling, 0 replies; 25+ messages in thread From: Aldy Hernandez @ 2023-04-25 9:08 UTC (permalink / raw) To: Jakub Jelinek, Joseph S. Myers, GCC patches, Andrew MacLeod [-- Attachment #1: Type: text/plain, Size: 7070 bytes --] On 4/21/23 18:40, Jakub Jelinek wrote: > On Thu, Apr 20, 2023 at 02:59:35PM +0200, Jakub Jelinek via Gcc-patches wrote: >> Thanks for working on this. Though expectedly here we are running >> again into the discussions we had in November about math properties of the >> functions vs. numeric properties in their implementations, how big maximum >> error shall we expect for the functions (and whether we should hardcode >> it for all implementations, or have some more fine-grained list of expected >> ulp errors for each implementation), whether the implementations at least >> guarantee the basic mathematical properties of the functions even if they >> have some errors (say for sin/cos, if they really never return > 1.0 or < >> -1.0) and the same questions repeated for -frounding-math, what kind of >> extra errors to expect when using non-default rounding and whether say sin >> could ever return nextafter (1.0, 2.0) or even larger value say when >> using non-default rounding mode. >> https://gcc.gnu.org/pipermail/gcc-patches/2022-November/606466.html >> was my attempt to get at least some numbers on some targets, I'm afraid >> for most implementations we aren't going to get any numerical proofs of >> maximum errors and the like. For sin/cos to check whether the implementation >> really never returns > 1.0 or < -1.0 perhaps instead of using randomized >> testing we could exhaustively check some range around 0, M_PI, 3*M_PI, >> -M_PI, -3*M_PI, and say some much larger multiples of M_PI, say 50 ulps >> in each direction about those points, and similarly for sin around M_PI/2 >> etc., in all arounding modes. > > Appart from the ulps ranges which I plan to introduce a target hook next > week... > >>> + if (!lh.maybe_isnan ()) >> >> This condition isn't sufficient, even if lh can't be NAN, but just >> may be +Inf or -Inf, the result needs to include maybe NAN. > > I've incorporated my other comments into a patch form. > I also found missing undefined_p () checks in two spots, the > r.set_undefined (); stuff being misplaced (it was done in the > lhs.maybe_isnan () case, which is incorrect, if the lhs may be nan, > then even if the finite range is say [-30., -10.] or [1.5., 42.], > result shouldn't be invalid, because the result still could be NAN > and in that case operand of say +-Inf or NAN would be valid. > Actually thinking about it some more, perhaps we should do that check > before the maybe_isnan () check and if we find the impossible finite, > either use r.set_undefined (); of !maybe_isnan () or handle it like > known_isnan () otherwise. > > Also, I think we should remember if it is SIN or COS, we'll need it both > for the ulps case and if we improve it for (ub-lb) < 2*PI ranges. > Now, talking about that, I'd very much like to avoid finding if some > multiple of PI/2 occurs inside of such ranges, the precision of real.cc > is clearly not sufficient for that, but perhaps we could use derivative > of sin (cos) or of cos (sin) to see if the function on the boundary is > increasing or decreasing and from that on both boundaries and their > approximate distance find out if the range needs to be extended to +1 > or -1. Could we do this PI stuff as a follow-up, or should it go in this patch? > > So, just incremental WIP so far... > > --- gcc/gimple-range-op.cc.jj 2023-04-21 17:09:48.250367999 +0200 > +++ gcc/gimple-range-op.cc 2023-04-21 18:37:26.048325391 +0200 > @@ -405,17 +405,20 @@ class cfn_sincos : public range_operator > public: > using range_operator_float::fold_range; > using range_operator_float::op1_range; > + cfn_sincos (combined_fn cfn) { m_cfn = cfn; } > virtual bool fold_range (frange &r, tree type, > const frange &lh, const frange &, > relation_trio) const final override > { > + if (lh.undefined_p ()) > + return false; > if (lh.known_isnan () || lh.known_isinf ()) > { > r.set_nan (type); > return true; > } > r.set (type, dconstm1, dconst1); > - if (!lh.maybe_isnan ()) > + if (!lh.maybe_isnan () && !lh.maybe_isinf ()) > r.clear_nan (); > return true; > } > @@ -423,15 +426,9 @@ public: > const frange &lhs, const frange &, > relation_trio) const final override > { > - if (!lhs.maybe_isnan ()) > - { > - // If NAN is not valid result, the input cannot include either > - // a NAN nor a +-INF. > - REAL_VALUE_TYPE lb = real_min_representable (type); > - REAL_VALUE_TYPE ub = real_max_representable (type); > - r.set (type, lb, ub, nan_state (false, false)); > - return true; > - } > + if (lhs.undefined_p ()) > + return false; > + > // A known NAN means the input is [-INF,-INF][+INF,+INF] U +-NAN, > // which we can't currently represent. > if (lhs.known_isnan ()) > @@ -439,20 +436,38 @@ public: > r.set_varying (type); > return true; > } > + > // Results outside of [-1.0, +1.0] are impossible. > REAL_VALUE_TYPE lb = lhs.lower_bound (); > REAL_VALUE_TYPE ub = lhs.upper_bound (); > - if (real_less (&lb, &dconstm1) > - || real_less (&dconst1, &ub)) > + if (real_less (&lb, &dconstm1) || real_less (&dconst1, &ub)) > { > - r.set_undefined (); > + if (!lhs.maybe_isnan ()) > + r.set_undefined (); > + else > + /* If lhs could be NAN and finite result is impossible, > + the range is like lhs.known_isnan () above, > + [-INF,-INF][+INF,+INF] U +-NAN. */ > + r.set_varying (type); > + return true; > + } > + > + if (!lhs.maybe_isnan ()) > + { > + // If NAN is not valid result, the input cannot include either > + // a NAN nor a +-INF. > + lb = real_min_representable (type); > + ub = real_max_representable (type); > + r.set (type, lb, ub, nan_state (false, false)); I don't like the nan_state(false, false) idiom, and that's my fault. For that matter, the nan_state() default constructor is confusing because it's not clear whether that means a NAN or not a NAN. I'm fixing this in a patch I just posted. > return true; > } > > r.set_varying (type); > return true; > } > -} op_cfn_sincos; > +private: > + combined_fn m_cfn; > +} op_cfn_sin (CFN_SIN), op_cfn_cos (CFN_COS); > > // Implement range operator for CFN_BUILT_IN_TOUPPER and CFN_BUILT_IN_TOLOWER. > class cfn_toupper_tolower : public range_operator > @@ -934,10 +949,15 @@ gimple_range_op_handler::maybe_builtin_c > > CASE_CFN_SIN: > CASE_CFN_SIN_FN: > + m_op1 = gimple_call_arg (call, 0); > + m_float = &op_cfn_sin; > + m_valid = true; > + break; > + > CASE_CFN_COS: > CASE_CFN_COS_FN: > m_op1 = gimple_call_arg (call, 0); > - m_float = &op_cfn_sincos; > + m_float = &op_cfn_cos; > m_valid = true; > break; Thanks. I've merged your work with mine and am retesting. I have also swapped the lb/ub per Mikael's comment downthread. Should we adjust the range by say 50 ulps in either direction, and do your target hook as a follow-up? Aldy [-- Attachment #2: 0002-Implement-range-op-entry-for-sin-cos.patch --] [-- Type: text/x-patch, Size: 4826 bytes --] From 065b8a1c1ec582eb3445ddb7f87ce9adc74394db Mon Sep 17 00:00:00 2001 From: Aldy Hernandez <aldyh@redhat.com> Date: Wed, 29 Mar 2023 12:46:38 +0200 Subject: [PATCH] Implement range-op entry for sin/cos. This is the range-op entry for sin/cos. It is meant to serve as an example of what we can do for glibc math functions. It is by no means exhaustive, just a stub to restrict the return range from sin/cos to [-1.0, 1.0] with appropriate smarts of NANs. As can be seen in the testcase, we see sin() as well as __builtin_sin() in the IL, and can resolve the resulting range accordingly. gcc/ChangeLog: * gimple-range-op.cc (class cfn_sincos): New. (gimple_range_op_handler::maybe_builtin_call): Add case for sin/cos. gcc/testsuite/ChangeLog: * gcc.dg/tree-ssa/range-sincos.c: New test. --- gcc/gimple-range-op.cc | 83 ++++++++++++++++++++ gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c | 40 ++++++++++ 2 files changed, 123 insertions(+) create mode 100644 gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c diff --git a/gcc/gimple-range-op.cc b/gcc/gimple-range-op.cc index f7409e35a99..90305b8614a 100644 --- a/gcc/gimple-range-op.cc +++ b/gcc/gimple-range-op.cc @@ -400,6 +400,75 @@ public: } } op_cfn_copysign; +class cfn_sincos : public range_operator_float +{ +public: + using range_operator_float::fold_range; + using range_operator_float::op1_range; + cfn_sincos (combined_fn cfn) { m_cfn = cfn; } + virtual bool fold_range (frange &r, tree type, + const frange &lh, const frange &, + relation_trio) const final override + { + if (lh.undefined_p ()) + return false; + if (lh.known_isnan () || lh.known_isinf ()) + { + r.set_nan (type); + return true; + } + r.set (type, dconstm1, dconst1); + if (!lh.maybe_isnan () && !lh.maybe_isinf ()) + r.clear_nan (); + return true; + } + virtual bool op1_range (frange &r, tree type, + const frange &lhs, const frange &, + relation_trio) const final override + { + if (lhs.undefined_p ()) + return false; + + // A known NAN means the input is [-INF,-INF][+INF,+INF] U +-NAN, + // which we can't currently represent. + if (lhs.known_isnan ()) + { + r.set_varying (type); + return true; + } + + // Results outside of [-1.0, +1.0] are impossible. + REAL_VALUE_TYPE lb = lhs.lower_bound (); + REAL_VALUE_TYPE ub = lhs.upper_bound (); + if (real_less (&ub, &dconstm1) || real_less (&dconst1, &lb)) + { + if (!lhs.maybe_isnan ()) + r.set_undefined (); + else + /* If lhs could be NAN and finite result is impossible, + the range is like lhs.known_isnan () above, + [-INF,-INF][+INF,+INF] U +-NAN. */ + r.set_varying (type); + return true; + } + + if (!lhs.maybe_isnan ()) + { + // If NAN is not a valid result, the input cannot include + // either a NAN nor a +-INF. + lb = real_min_representable (type); + ub = real_max_representable (type); + r.set (type, lb, ub, nan_state (false)); + return true; + } + + r.set_varying (type); + return true; + } +private: + combined_fn m_cfn; +} op_cfn_sin (CFN_SIN), op_cfn_cos (CFN_COS); + // Implement range operator for CFN_BUILT_IN_TOUPPER and CFN_BUILT_IN_TOLOWER. class cfn_toupper_tolower : public range_operator { @@ -878,6 +947,20 @@ gimple_range_op_handler::maybe_builtin_call () m_valid = true; break; + CASE_CFN_SIN: + CASE_CFN_SIN_FN: + m_op1 = gimple_call_arg (call, 0); + m_float = &op_cfn_sin; + m_valid = true; + break; + + CASE_CFN_COS: + CASE_CFN_COS_FN: + m_op1 = gimple_call_arg (call, 0); + m_float = &op_cfn_cos; + m_valid = true; + break; + case CFN_BUILT_IN_TOUPPER: case CFN_BUILT_IN_TOLOWER: // Only proceed If the argument is compatible with the LHS. diff --git a/gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c b/gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c new file mode 100644 index 00000000000..50fbac1b68f --- /dev/null +++ b/gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c @@ -0,0 +1,40 @@ +// { dg-do compile } +// { dg-options "-O2 -fdump-tree-evrp -fno-thread-jumps" } + +#include <math.h> + +void use(double); +void link_error (); + +void foo(double x) +{ + if (__builtin_isnan (x)) + __builtin_unreachable(); + x = sin (x); + if (x < -1.0 || x > 1.0) + link_error (); + use (x); +} + +void bar(double x) +{ + if (!__builtin_isnan (sin (x))) + { + if (__builtin_isnan (x)) + link_error (); + if (__builtin_isinf (x)) + link_error (); + } +} + +void stool (double x) +{ + double res1 = sin (x); + double res2 = __builtin_sin (x); + if (res1 < -1.0 || res2 < -1.0) + link_error (); + if (res1 > 1.0 || res2 > 1.0) + link_error (); +} + +// { dg-final { scan-tree-dump-not "link_error" "evrp" } } -- 2.40.0 ^ permalink raw reply [flat|nested] 25+ messages in thread
* [PATCH] v2: Implement range-op entry for sin/cos 2023-04-18 13:12 [PATCH] Implement range-op entry for sin/cos Aldy Hernandez 2023-04-20 12:59 ` Jakub Jelinek @ 2023-04-27 11:13 ` Jakub Jelinek 2023-04-27 11:46 ` Aldy Hernandez 1 sibling, 1 reply; 25+ messages in thread From: Jakub Jelinek @ 2023-04-27 11:13 UTC (permalink / raw) To: Aldy Hernandez; +Cc: GCC patches, Andrew MacLeod Hi! On Tue, Apr 18, 2023 at 03:12:50PM +0200, Aldy Hernandez wrote: > [I don't know why I keep poking at floats. I must really like the pain. > Jakub, are you OK with this patch for trunk?] > > This is the range-op entry for sin/cos. It is meant to serve as an > example of what we can do for glibc math functions. It is by no means > exhaustive, just a stub to restrict the return range from sin/cos to > [-1.0, 1.0] with appropriate smarts of NANs. > > As can be seen in the testcase, we see sin() as well as > __builtin_sin() in the IL, and can resolve the resulting range > accordingly. > > gcc/ChangeLog: > > * gimple-range-op.cc (class cfn_sincos): New. > (gimple_range_op_handler::maybe_builtin_call): Add case for sin/cos. > > gcc/testsuite/ChangeLog: > > * gcc.dg/tree-ssa/range-sincos.c: New test. Here is an updated version of the patch on top of the v2: Add targetm.libm_function_max_error patch with all my comments incorporated into your patch (but still no handling of sin/cos ranges shorter than 2*M_PI). So far tested on the new testcase with RUNTESTFLAGS='--target_board=unix\{,-frounding-math} tree-ssa.exp=range-sincos.c' where expectedly without -frounding-math it PASSes and with it fails (as the hooks return 4ulps in that case for the boundaries). Ok for trunk if it passes bootstrap/regtest? I'll defer tweaks to frange_nextafter for ulps or whatever you want to do with it. 2023-04-27 Aldy Hernandez <aldyh@redhat.com> Jakub Jelinek <jakub@redhat.com> * value-range.h (frange_nextafter): Declare. * gimple-range-op.cc (class cfn_sincos): New. (op_cfn_sin, op_cfn_cos): New variables. (gimple_range_op_handler::maybe_builtin_call): Handle CASE_CFN_{SIN,COS}{,_FN}. * gcc.dg/tree-ssa/range-sincos.c: New test. --- gcc/value-range.h.jj 2023-04-27 10:17:46.504485376 +0200 +++ gcc/value-range.h 2023-04-27 12:07:09.891147869 +0200 @@ -1388,4 +1388,7 @@ frange::nan_signbit_p (bool &signbit) co return true; } +void frange_nextafter (enum machine_mode, REAL_VALUE_TYPE &, + const REAL_VALUE_TYPE &); + #endif // GCC_VALUE_RANGE_H --- gcc/gimple-range-op.cc.jj 2023-04-22 10:23:26.942812958 +0200 +++ gcc/gimple-range-op.cc 2023-04-27 11:57:09.865879982 +0200 @@ -400,6 +400,89 @@ public: } } op_cfn_copysign; +class cfn_sincos : public range_operator_float +{ +public: + using range_operator_float::fold_range; + using range_operator_float::op1_range; + cfn_sincos (combined_fn cfn) { m_cfn = cfn; } + virtual bool fold_range (frange &r, tree type, + const frange &lh, const frange &, + relation_trio) const final override + { + if (lh.undefined_p ()) + return false; + if (lh.known_isnan () || lh.known_isinf ()) + { + r.set_nan (type); + return true; + } + unsigned bulps = targetm.libm_function_max_error (m_cfn, TYPE_MODE (type), + true); + if (bulps == ~0U) + r.set_varying (type); + else if (bulps == 0) + r.set (type, dconstm1, dconst1); + else + { + REAL_VALUE_TYPE boundmin, boundmax; + boundmax = dconst1; + while (bulps--) + frange_nextafter (TYPE_MODE (type), boundmax, dconstinf); + real_arithmetic (&boundmin, NEGATE_EXPR, &boundmax, NULL); + r.set (type, boundmin, boundmax); + } + if (!lh.maybe_isnan () && !lh.maybe_isinf ()) + r.clear_nan (); + return true; + } + virtual bool op1_range (frange &r, tree type, + const frange &lhs, const frange &, + relation_trio) const final override + { + if (lhs.undefined_p ()) + return false; + + // A known NAN means the input is [-INF,-INF][+INF,+INF] U +-NAN, + // which we can't currently represent. + if (lhs.known_isnan ()) + { + r.set_varying (type); + return true; + } + + // Results outside of [-1.0, +1.0] are impossible. + REAL_VALUE_TYPE lb = lhs.lower_bound (); + REAL_VALUE_TYPE ub = lhs.upper_bound (); + if (real_less (&ub, &dconstm1) || real_less (&dconst1, &lb)) + { + if (!lhs.maybe_isnan ()) + r.set_undefined (); + else + /* If lhs could be NAN and finite result is impossible, + the range is like lhs.known_isnan () above, + [-INF,-INF][+INF,+INF] U +-NAN. */ + r.set_varying (type); + return true; + } + + if (!lhs.maybe_isnan ()) + { + // If NAN is not valid result, the input cannot include either + // a NAN nor a +-INF. + lb = real_min_representable (type); + ub = real_max_representable (type); + r.set (type, lb, ub, nan_state (false, false)); + return true; + } + + r.set_varying (type); + return true; + } +private: + combined_fn m_cfn; +} op_cfn_sin (CFN_SIN), op_cfn_cos (CFN_COS); + // Implement range operator for CFN_BUILT_IN_TOUPPER and CFN_BUILT_IN_TOLOWER. class cfn_toupper_tolower : public range_operator { @@ -878,6 +961,20 @@ gimple_range_op_handler::maybe_builtin_c m_valid = true; break; + CASE_CFN_SIN: + CASE_CFN_SIN_FN: + m_op1 = gimple_call_arg (call, 0); + m_float = &op_cfn_sin; + m_valid = true; + break; + + CASE_CFN_COS: + CASE_CFN_COS_FN: + m_op1 = gimple_call_arg (call, 0); + m_float = &op_cfn_cos; + m_valid = true; + break; + case CFN_BUILT_IN_TOUPPER: case CFN_BUILT_IN_TOLOWER: // Only proceed If the argument is compatible with the LHS. --- gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c.jj 2023-04-27 10:27:08.135348612 +0200 +++ gcc/testsuite/gcc.dg/tree-ssa/range-sincos.c 2023-04-27 12:14:50.788437504 +0200 @@ -0,0 +1,43 @@ +// { dg-do compile } +// { dg-options "-O2 -fdump-tree-evrp -fno-thread-jumps" } + +#include <math.h> + +void use (double); +void link_error (); + +void +foo (double x) +{ + if (__builtin_isnan (x)) + __builtin_unreachable (); + x = sin (x); + if (x < -1.0 || x > 1.0) + link_error (); + use (x); +} + +void +bar (double x) +{ + if (!__builtin_isnan (sin (x))) + { + if (__builtin_isnan (x)) + link_error (); + if (__builtin_isinf (x)) + link_error (); + } +} + +void +stool (double x) +{ + double res1 = sin (x); + double res2 = __builtin_sin (x); + if (res1 < -1.0 || res2 < -1.0) + link_error (); + if (res1 > 1.0 || res2 > 1.0) + link_error (); +} + +// { dg-final { scan-tree-dump-not "link_error" "evrp" { target { { *-*-linux* } && { glibc } } } } } Jakub ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] v2: Implement range-op entry for sin/cos 2023-04-27 11:13 ` [PATCH] v2: " Jakub Jelinek @ 2023-04-27 11:46 ` Aldy Hernandez 2023-04-27 11:53 ` Jakub Jelinek 0 siblings, 1 reply; 25+ messages in thread From: Aldy Hernandez @ 2023-04-27 11:46 UTC (permalink / raw) To: Jakub Jelinek; +Cc: GCC patches, Andrew MacLeod On 4/27/23 13:13, Jakub Jelinek wrote: > + unsigned bulps = targetm.libm_function_max_error (m_cfn, TYPE_MODE (type), > + true); > + if (bulps == ~0U) > + r.set_varying (type); > + else if (bulps == 0) > + r.set (type, dconstm1, dconst1); > + else > + { > + REAL_VALUE_TYPE boundmin, boundmax; > + boundmax = dconst1; > + while (bulps--) > + frange_nextafter (TYPE_MODE (type), boundmax, dconstinf); > + real_arithmetic (&boundmin, NEGATE_EXPR, &boundmax, NULL); > + r.set (type, boundmin, boundmax); > + } This seems like something we'll do over and over for other operations, right? If so, could you abstract it into a separate function? Thanks. Aldy ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] v2: Implement range-op entry for sin/cos 2023-04-27 11:46 ` Aldy Hernandez @ 2023-04-27 11:53 ` Jakub Jelinek 2023-04-27 12:03 ` Aldy Hernandez 0 siblings, 1 reply; 25+ messages in thread From: Jakub Jelinek @ 2023-04-27 11:53 UTC (permalink / raw) To: Aldy Hernandez; +Cc: GCC patches, Andrew MacLeod On Thu, Apr 27, 2023 at 01:46:19PM +0200, Aldy Hernandez wrote: > On 4/27/23 13:13, Jakub Jelinek wrote: > > > + unsigned bulps = targetm.libm_function_max_error (m_cfn, TYPE_MODE (type), > > + true); > > + if (bulps == ~0U) > > + r.set_varying (type); > > + else if (bulps == 0) > > + r.set (type, dconstm1, dconst1); > > + else > > + { > > + REAL_VALUE_TYPE boundmin, boundmax; > > + boundmax = dconst1; > > + while (bulps--) > > + frange_nextafter (TYPE_MODE (type), boundmax, dconstinf); > > + real_arithmetic (&boundmin, NEGATE_EXPR, &boundmax, NULL); > > + r.set (type, boundmin, boundmax); > > + } > > This seems like something we'll do over and over for other operations, > right? If so, could you abstract it into a separate function? Not easily. E.g. take the difference between sin/cos, where the above grows the interval on both sides by bulps, vs. what I'm working on right now (sqrt), where it grows in one direction only (from dconstm0 toward dconstninf). There could be other functions which only grow the positive bound. Jakub ^ permalink raw reply [flat|nested] 25+ messages in thread
* Re: [PATCH] v2: Implement range-op entry for sin/cos 2023-04-27 11:53 ` Jakub Jelinek @ 2023-04-27 12:03 ` Aldy Hernandez 0 siblings, 0 replies; 25+ messages in thread From: Aldy Hernandez @ 2023-04-27 12:03 UTC (permalink / raw) To: Jakub Jelinek; +Cc: GCC patches, Andrew MacLeod On 4/27/23 13:53, Jakub Jelinek wrote: > On Thu, Apr 27, 2023 at 01:46:19PM +0200, Aldy Hernandez wrote: >> On 4/27/23 13:13, Jakub Jelinek wrote: >> >>> + unsigned bulps = targetm.libm_function_max_error (m_cfn, TYPE_MODE (type), >>> + true); >>> + if (bulps == ~0U) >>> + r.set_varying (type); >>> + else if (bulps == 0) >>> + r.set (type, dconstm1, dconst1); >>> + else >>> + { >>> + REAL_VALUE_TYPE boundmin, boundmax; >>> + boundmax = dconst1; >>> + while (bulps--) >>> + frange_nextafter (TYPE_MODE (type), boundmax, dconstinf); >>> + real_arithmetic (&boundmin, NEGATE_EXPR, &boundmax, NULL); >>> + r.set (type, boundmin, boundmax); >>> + } >> >> This seems like something we'll do over and over for other operations, >> right? If so, could you abstract it into a separate function? > > Not easily. E.g. take the difference between sin/cos, where the above > grows the interval on both sides by bulps, vs. what I'm working on right now > (sqrt), where it grows in one direction only (from dconstm0 toward > dconstninf). There could be other functions which only grow the positive > bound. Fair enough. Then OK. Thanks for working on this. Aldy ^ permalink raw reply [flat|nested] 25+ messages in thread
end of thread, other threads:[~2023-04-27 12:03 UTC | newest] Thread overview: 25+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2023-04-18 13:12 [PATCH] Implement range-op entry for sin/cos Aldy Hernandez 2023-04-20 12:59 ` Jakub Jelinek 2023-04-20 13:17 ` Siddhesh Poyarekar 2023-04-20 14:02 ` Jakub Jelinek 2023-04-20 14:20 ` Jakub Jelinek 2023-04-20 15:22 ` Siddhesh Poyarekar 2023-04-20 15:52 ` Jakub Jelinek 2023-04-20 17:57 ` Siddhesh Poyarekar 2023-04-21 1:14 ` Siddhesh Poyarekar 2023-04-21 6:52 ` Jakub Jelinek 2023-04-21 11:20 ` Siddhesh Poyarekar 2023-04-25 8:59 ` Aldy Hernandez 2023-04-24 16:03 ` Jakub Jelinek 2023-04-24 16:05 ` Siddhesh Poyarekar 2023-04-24 16:09 ` Jakub Jelinek 2023-04-24 16:33 ` Jeff Law 2023-04-21 16:40 ` Jakub Jelinek 2023-04-21 20:43 ` Mikael Morin 2023-04-21 20:45 ` Jakub Jelinek 2023-04-25 9:10 ` Aldy Hernandez 2023-04-25 9:08 ` Aldy Hernandez 2023-04-27 11:13 ` [PATCH] v2: " Jakub Jelinek 2023-04-27 11:46 ` Aldy Hernandez 2023-04-27 11:53 ` Jakub Jelinek 2023-04-27 12:03 ` Aldy Hernandez
This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).