| Message ID | CAEFO=4CzZzMDfFRa7_G4LUnAD0+NRqe9KhYQ9h7GPjkk8Uqcaw@mail.gmail.com |
|---|---|
| State | New |
| Headers | show |
| Series | Optimize sin(atan(x)), take 2 | expand |
Ping. On Mon, Sep 3, 2018 at 4:11 PM, Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: > Fixed the issues pointed by the previous discussions. Closes PR86829. > > Adds substitution rules for sin(atan(x)) and cos(atan(x)), being > careful with overflow issues by constructing a assumed convergence > constant (see comment in real.c). > > 2018-09-03 Giuliano Belinassi <giuliano.belinassi@usp.br> > > * match.pd: add simplification rules to sin(atan(x)) and cos(atan(x)). > * real.c: add code for assumed convergence constant to sin(atan(x)). > * real.h: allows the added code from real.c to be called externally. > * tree.c: add code for bulding nodes with the convergence constant. > * tree.h: allows the added code from tree.c to be called externally. > * sinatan-1.c: tests assumed convergence constant. > * sinatan-2.c: tests simplification rule. > * sinatan-3.c: likewise. > > There seems to be no broken tests in trunk that are related to this > modification.
Pinging match.pd and real.c maintainers, as suggested in IRC. Sorry if it is inappropriate On Mon, Sep 17, 2018 at 9:46 AM Giuliano Augusto Faulin Belinassi <giuliano.belinassi@usp.br> wrote: > > Ping. > > On Mon, Sep 3, 2018 at 4:11 PM, Giuliano Augusto Faulin Belinassi > <giuliano.belinassi@usp.br> wrote: > > Fixed the issues pointed by the previous discussions. Closes PR86829. > > > > Adds substitution rules for sin(atan(x)) and cos(atan(x)), being > > careful with overflow issues by constructing a assumed convergence > > constant (see comment in real.c). > > > > 2018-09-03 Giuliano Belinassi <giuliano.belinassi@usp.br> > > > > * match.pd: add simplification rules to sin(atan(x)) and cos(atan(x)). > > * real.c: add code for assumed convergence constant to sin(atan(x)). > > * real.h: allows the added code from real.c to be called externally. > > * tree.c: add code for bulding nodes with the convergence constant. > > * tree.h: allows the added code from tree.c to be called externally. > > * sinatan-1.c: tests assumed convergence constant. > > * sinatan-2.c: tests simplification rule. > > * sinatan-3.c: likewise. > > > > There seems to be no broken tests in trunk that are related to this > > modification.
On 9/3/18 1:11 PM, Giuliano Augusto Faulin Belinassi wrote: > Fixed the issues pointed by the previous discussions. Closes PR86829. > > Adds substitution rules for sin(atan(x)) and cos(atan(x)), being > careful with overflow issues by constructing a assumed convergence > constant (see comment in real.c). > > 2018-09-03 Giuliano Belinassi <giuliano.belinassi@usp.br> > > * match.pd: add simplification rules to sin(atan(x)) and cos(atan(x)). > * real.c: add code for assumed convergence constant to sin(atan(x)). > * real.h: allows the added code from real.c to be called externally. > * tree.c: add code for bulding nodes with the convergence constant. > * tree.h: allows the added code from tree.c to be called externally. > * sinatan-1.c: tests assumed convergence constant. > * sinatan-2.c: tests simplification rule. > * sinatan-3.c: likewise. > > There seems to be no broken tests in trunk that are related to this > modification. Pretty cool. > > > sinatanv2.patch > > Index: gcc/match.pd > =================================================================== > --- gcc/match.pd (revisão 264058) > +++ gcc/match.pd (cópia de trabalho) > @@ -4169,6 +4169,39 @@ > (tans (atans @0)) > @0))) > > + /* Simplify sin(atan(x)) -> x / sqrt(x*x + 1). */ > + (for sins (SIN) > + atans (ATAN) > + sqrts (SQRT) > + (simplify > + (sins (atans:s @0)) > + (if (SCALAR_FLOAT_TYPE_P (type)) > + (switch > + (if (types_match (type, float_type_node)) > + (cond (le (abs @0) { build_sinatan_cst (float_type_node); }) > + (rdiv @0 (sqrts (plus (mult @0 @0) > + {build_one_cst (float_type_node);}))) > + (BUILT_IN_COPYSIGNF { build_one_cst (float_type_node); } @0))) > + (if (types_match (type, double_type_node)) > + (cond (le (abs @0) { build_sinatan_cst (double_type_node); }) > + (rdiv @0 (sqrts (plus (mult @0 @0) > + {build_one_cst (double_type_node);}))) > + (BUILT_IN_COPYSIGN { build_one_cst (double_type_node); } @0))) > + (if (types_match (type, long_double_type_node)) > + (cond (le (abs @0) { build_sinatan_cst (long_double_type_node); }) > + (rdiv @0 (sqrts (plus (mult @0 @0) > + {build_one_cst (long_double_type_node);}))) > + (BUILT_IN_COPYSIGNL { build_one_cst (long_double_type_node); } @0))))))) So you don't want to build the constants as a float, double or long double. Instead you want to build it as "type". I think that should let you simplify this a bit. It turns into (if (SCALAR_FLOAT_TYPE_P (type)) (cond (le (abs @0) {build_sinatan_cst (type); }) (rdiv @0 (sqrts (plus (mult @0 @0) {build_one_cst (type); }))) (BUILT_IN_COPYSIGNL { build_one_cst (type); } @0)))) Or something along those lines I think. Richi is much better at match.pd stuff than I. > + > +/* Simplify cos(atan(x)) -> 1 / sqrt(x*x + 1). */ > + (for coss (COS) > + atans (ATAN) > + sqrts (SQRT) > + (simplify > + (coss (atans:s @0)) > + (rdiv {build_one_cst (type);} > + (sqrts (plus (mult @0 @0) {build_one_cst (type);}))))) Don't we have the same kinds of issues with the x*x in here? As X gets large it will overflow, but the result is going to be approaching zero. So we might be able to use a similar trick here. > + > +/* Build real constant used by sin(atan(x)) optimization. > + The logic here is as follows: It can be mathematically > + shown that sin(atan(x)) = x / sqrt(1 + x*x), but notice > + that the second formula has an x*x, which can overflow > + if x is big enough. However, x / sqrt(1 + x*x) converges > + to 1 for large x. What must be the value of x such that > + when computing x / sqrt (1 + x*x) = 1? > + > + Therefore, we must then solve x / sqrt(1 + x*x) > eps > + for x, where eps is the largest number smaller than 1 > + representable by the target. Hence, solving for eps > + yields that x > eps / sqrt(1 - eps*eps). This eps can > + be easily calculated by calling nextafter. Likewise for > + the negative x. */ Imagine a pause here while I lookup isolation of radicals.... It's been a long time... OK. Sure. I see what you're doing here... > + > +void > +build_sinatan_real (REAL_VALUE_TYPE * r, tree type) > +{ > + REAL_VALUE_TYPE eps; > + mpfr_t mpfr_eps, mpfr_const1, c, divisor; > + const struct real_format * fmt = NULL; > + > + fmt = type ? REAL_MODE_FORMAT (TYPE_MODE (type)) : NULL; I believe we can and should always pass in a value TYPE, so it's never going to be NULL. It also seems like initializing fmt to NULL in the previous statement doesn't make much sense. > + > + mpfr_inits (divisor, mpfr_eps, mpfr_const1, c, NULL); > + mpfr_from_real (mpfr_const1, &dconst1, GMP_RNDN); > + > + real_nextafter (&eps, fmt, &dconst1, &dconst0); > + > + mpfr_from_real (mpfr_eps, &eps, GMP_RNDZ); > + mpfr_mul (divisor, mpfr_eps, mpfr_eps, GMP_RNDU); > + mpfr_sub (divisor, mpfr_const1, divisor, GMP_RNDZ); > + mpfr_sqrt (divisor, divisor, GMP_RNDZ); > + mpfr_div (c, mpfr_eps, divisor, GMP_RNDU); > + > + real_from_mpfr (r, c, fmt, GMP_RNDU); > + > + /* For safety reasons. */ > + times_pten(r, 1); Not sure what you mean for safety reasons. The calculations to produce "c" then convert it into a REAL_VALUE_TYPE all make sense. Just not sure what this line is really meant to do. build_sinatan_cst needs a function comment. Something like /* Build and return the tree constant used by the sin(atan)) optimization. */ Or something like that. This feels really close to being ready. Jeff
On Thu, Oct 4, 2018 at 6:39 AM Jeff Law <law@redhat.com> wrote: > > On 9/3/18 1:11 PM, Giuliano Augusto Faulin Belinassi wrote: > > Fixed the issues pointed by the previous discussions. Closes PR86829. > > > > Adds substitution rules for sin(atan(x)) and cos(atan(x)), being > > careful with overflow issues by constructing a assumed convergence > > constant (see comment in real.c). > > > > 2018-09-03 Giuliano Belinassi <giuliano.belinassi@usp.br> > > > > * match.pd: add simplification rules to sin(atan(x)) and cos(atan(x)). > > * real.c: add code for assumed convergence constant to sin(atan(x)). > > * real.h: allows the added code from real.c to be called externally. > > * tree.c: add code for bulding nodes with the convergence constant. > > * tree.h: allows the added code from tree.c to be called externally. Look at examples to fix up your ChangeLog, you need to list functions you add/modify. I don't think you want the tree.[ch] forwarders, instead use (with { REAL_VALUE_TYPE rcst; build_sinatan_real (&cst, type); tree tcst = build_real (type, cst); } in the pattern then just use (cond (le (abs @0) { tcst; }) or even include the function here in full (well, maybe not, but maybe include some of its comment to make it clear what kind of constant we build here). > > * sinatan-1.c: tests assumed convergence constant. > > * sinatan-2.c: tests simplification rule. > > * sinatan-3.c: likewise. > > > > There seems to be no broken tests in trunk that are related to this > > modification. > Pretty cool. > > > > > > > > sinatanv2.patch > > > > Index: gcc/match.pd > > =================================================================== > > --- gcc/match.pd (revisão 264058) > > +++ gcc/match.pd (cópia de trabalho) > > @@ -4169,6 +4169,39 @@ > > (tans (atans @0)) > > @0))) > > > > + /* Simplify sin(atan(x)) -> x / sqrt(x*x + 1). */ > > + (for sins (SIN) > > + atans (ATAN) > > + sqrts (SQRT) > > + (simplify > > + (sins (atans:s @0)) > > + (if (SCALAR_FLOAT_TYPE_P (type)) > > + (switch > > + (if (types_match (type, float_type_node)) > > + (cond (le (abs @0) { build_sinatan_cst (float_type_node); }) > > + (rdiv @0 (sqrts (plus (mult @0 @0) > > + {build_one_cst (float_type_node);}))) > > + (BUILT_IN_COPYSIGNF { build_one_cst (float_type_node); } @0))) > > + (if (types_match (type, double_type_node)) > > + (cond (le (abs @0) { build_sinatan_cst (double_type_node); }) > > + (rdiv @0 (sqrts (plus (mult @0 @0) > > + {build_one_cst (double_type_node);}))) > > + (BUILT_IN_COPYSIGN { build_one_cst (double_type_node); } @0))) > > + (if (types_match (type, long_double_type_node)) > > + (cond (le (abs @0) { build_sinatan_cst (long_double_type_node); }) > > + (rdiv @0 (sqrts (plus (mult @0 @0) > > + {build_one_cst (long_double_type_node);}))) > > + (BUILT_IN_COPYSIGNL { build_one_cst (long_double_type_node); } @0))))))) > So you don't want to build the constants as a float, double or long > double. Instead you want to build it as "type". I think that should > let you simplify this a bit. It turns into > > (if (SCALAR_FLOAT_TYPE_P (type)) > (cond (le (abs @0) {build_sinatan_cst (type); }) > (rdiv @0 (sqrts (plus (mult @0 @0) > {build_one_cst (type); }))) > (BUILT_IN_COPYSIGNL { build_one_cst (type); } @0)))) > > Or something along those lines I think. Richi is much better at > match.pd stuff than I. Yeah, you'd need to add copysigns (COPYSIGN) to the for iterator and use that instead of BUILT_IN_COPYSIGNL of course. Richard. > > > > > > > > + > > +/* Simplify cos(atan(x)) -> 1 / sqrt(x*x + 1). */ > > + (for coss (COS) > > + atans (ATAN) > > + sqrts (SQRT) > > + (simplify > > + (coss (atans:s @0)) > > + (rdiv {build_one_cst (type);} > > + (sqrts (plus (mult @0 @0) {build_one_cst (type);}))))) > Don't we have the same kinds of issues with the x*x in here? As X gets > large it will overflow, but the result is going to be approaching zero. > So we might be able to use a similar trick here. > > > > + > > +/* Build real constant used by sin(atan(x)) optimization. > > + The logic here is as follows: It can be mathematically > > + shown that sin(atan(x)) = x / sqrt(1 + x*x), but notice > > + that the second formula has an x*x, which can overflow > > + if x is big enough. However, x / sqrt(1 + x*x) converges > > + to 1 for large x. What must be the value of x such that > > + when computing x / sqrt (1 + x*x) = 1? > > + > > + Therefore, we must then solve x / sqrt(1 + x*x) > eps > > + for x, where eps is the largest number smaller than 1 > > + representable by the target. Hence, solving for eps > > + yields that x > eps / sqrt(1 - eps*eps). This eps can > > + be easily calculated by calling nextafter. Likewise for > > + the negative x. */ > Imagine a pause here while I lookup isolation of radicals.... It's been > a long time... OK. Sure. I see what you're doing here... > > > > > + > > +void > > +build_sinatan_real (REAL_VALUE_TYPE * r, tree type) > > +{ > > + REAL_VALUE_TYPE eps; > > + mpfr_t mpfr_eps, mpfr_const1, c, divisor; > > + const struct real_format * fmt = NULL; > > + > > + fmt = type ? REAL_MODE_FORMAT (TYPE_MODE (type)) : NULL; > I believe we can and should always pass in a value TYPE, so it's never > going to be NULL. It also seems like initializing fmt to NULL in the > previous statement doesn't make much sense. > > > > > > + > > + mpfr_inits (divisor, mpfr_eps, mpfr_const1, c, NULL); > > + mpfr_from_real (mpfr_const1, &dconst1, GMP_RNDN); > > + > > + real_nextafter (&eps, fmt, &dconst1, &dconst0); > > + > > + mpfr_from_real (mpfr_eps, &eps, GMP_RNDZ); > > + mpfr_mul (divisor, mpfr_eps, mpfr_eps, GMP_RNDU); > > + mpfr_sub (divisor, mpfr_const1, divisor, GMP_RNDZ); > > + mpfr_sqrt (divisor, divisor, GMP_RNDZ); > > + mpfr_div (c, mpfr_eps, divisor, GMP_RNDU); > > + > > + real_from_mpfr (r, c, fmt, GMP_RNDU); > > + > > + /* For safety reasons. */ > > + times_pten(r, 1); > Not sure what you mean for safety reasons. The calculations to produce > "c" then convert it into a REAL_VALUE_TYPE all make sense. Just not > sure what this line is really meant to do. > > build_sinatan_cst needs a function comment. Something like > > /* Build and return the tree constant used by the sin(atan)) > optimization. */ > > Or something like that. > > This feels really close to being ready. > > Jeff
On Wed, 3 Oct 2018, Jeff Law wrote: >> +/* Simplify cos(atan(x)) -> 1 / sqrt(x*x + 1). */ >> + (for coss (COS) >> + atans (ATAN) >> + sqrts (SQRT) >> + (simplify >> + (coss (atans:s @0)) >> + (rdiv {build_one_cst (type);} >> + (sqrts (plus (mult @0 @0) {build_one_cst (type);}))))) > Don't we have the same kinds of issues with the x*x in here? As X gets > large it will overflow, but the result is going to be approaching zero. > So we might be able to use a similar trick here. (I have not read the patch) The similar trick would say that for X large, this is the same as 1/abs(X) I guess. Note that it may be simpler to generate a call to hypot (C99). A related remark would be: with the precision of double, for x>=cst (about 2^53), atan(x) is constant, within .5 ulp of pi/2 if the math library is super precise (which it probably isn't). Returning 0 for its cos (what happens if x*x gives +Inf) is thus completely fine unless you are using crlibm, but then you wouldn't use flag_unsafe_math_optimizations. The main issue is that if we have -ffast-math, we have -ffinite-math-only, and we are possibly introducing an infinity as intermediate result...
Thank you for the review. I will address all these issues :-). > Imagine a pause here while I lookup isolation of radicals.... It's been > a long time... OK. Sure. I see what you're doing here... Sorry, but I really did not understand your comment. Should I write a shorter comment for that function? > Not sure what you mean for safety reasons. The calculations to produce > "c" then convert it into a REAL_VALUE_TYPE all make sense. Just not > sure what this line is really meant to do. Imagine the following case: Let "c" be the real constant such that it is certain that for every x > "c", 1/sqrt(x*x + 1) = 1. Suppose now that our calculation leads us to a c' < "c" due to a minor imprecision. The logic here is that 10 * c' > "c" and everything will work, thus it is safer. Note however that I cannot prove that 10 * c' > "c", but I would be really surprised if this does not holds. > A related remark would be: with the precision of double, for x>=cst (about > 2^53), atan(x) is constant, within .5 ulp of pi/2 if the math library is > super precise (which it probably isn't). Returning 0 for its cos (what > happens if x*x gives +Inf) is thus completely fine unless you are using > crlibm, but then you wouldn't use flag_unsafe_math_optimizations. The main > issue is that if we have -ffast-math, we have -ffinite-math-only, and we > are possibly introducing an infinity as intermediate result... Thank you. This clarifies the need for a similar constant for the cos(atan(x)). On Thu, Oct 4, 2018 at 9:36 AM Marc Glisse <marc.glisse@inria.fr> wrote: > > On Wed, 3 Oct 2018, Jeff Law wrote: > > >> +/* Simplify cos(atan(x)) -> 1 / sqrt(x*x + 1). */ > >> + (for coss (COS) > >> + atans (ATAN) > >> + sqrts (SQRT) > >> + (simplify > >> + (coss (atans:s @0)) > >> + (rdiv {build_one_cst (type);} > >> + (sqrts (plus (mult @0 @0) {build_one_cst (type);}))))) > > Don't we have the same kinds of issues with the x*x in here? As X gets > > large it will overflow, but the result is going to be approaching zero. > > So we might be able to use a similar trick here. > > (I have not read the patch) > > The similar trick would say that for X large, this is the same as 1/abs(X) > I guess. Note that it may be simpler to generate a call to hypot (C99). > > A related remark would be: with the precision of double, for x>=cst (about > 2^53), atan(x) is constant, within .5 ulp of pi/2 if the math library is > super precise (which it probably isn't). Returning 0 for its cos (what > happens if x*x gives +Inf) is thus completely fine unless you are using > crlibm, but then you wouldn't use flag_unsafe_math_optimizations. The main > issue is that if we have -ffast-math, we have -ffinite-math-only, and we > are possibly introducing an infinity as intermediate result... > > -- > Marc Glisse
On 10/5/18 6:09 AM, Giuliano Augusto Faulin Belinassi wrote: > Thank you for the review. I will address all these issues :-). > >> Imagine a pause here while I lookup isolation of radicals.... It's been >> a long time... OK. Sure. I see what you're doing here... > > Sorry, but I really did not understand your comment. Should I write a > shorter comment for that function? Not at all -- I was just trying to be funny. I was struggling to understand the comment so I tried to recreate the steps necessary to transform the equation myself. And quickly realized that it's been 30 years since I've had to do that kind of algebra :-) > >> Not sure what you mean for safety reasons. The calculations to produce >> "c" then convert it into a REAL_VALUE_TYPE all make sense. Just not >> sure what this line is really meant to do. > > Imagine the following case: > Let "c" be the real constant such that it is certain that for every x >> "c", 1/sqrt(x*x + 1) = 1. > Suppose now that our calculation leads us to a c' < "c" due to a minor > imprecision. > The logic here is that 10 * c' > "c" and everything will work, thus it is safer. > Note however that I cannot prove that 10 * c' > "c", but I would be > really surprised > if this does not holds. Ah. I understand. This probably warrants a slightly better comment. Jeff
Index: gcc/match.pd =================================================================== --- gcc/match.pd (revisão 264058) +++ gcc/match.pd (cópia de trabalho) @@ -4169,6 +4169,39 @@ (tans (atans @0)) @0))) + /* Simplify sin(atan(x)) -> x / sqrt(x*x + 1). */ + (for sins (SIN) + atans (ATAN) + sqrts (SQRT) + (simplify + (sins (atans:s @0)) + (if (SCALAR_FLOAT_TYPE_P (type)) + (switch + (if (types_match (type, float_type_node)) + (cond (le (abs @0) { build_sinatan_cst (float_type_node); }) + (rdiv @0 (sqrts (plus (mult @0 @0) + {build_one_cst (float_type_node);}))) + (BUILT_IN_COPYSIGNF { build_one_cst (float_type_node); } @0))) + (if (types_match (type, double_type_node)) + (cond (le (abs @0) { build_sinatan_cst (double_type_node); }) + (rdiv @0 (sqrts (plus (mult @0 @0) + {build_one_cst (double_type_node);}))) + (BUILT_IN_COPYSIGN { build_one_cst (double_type_node); } @0))) + (if (types_match (type, long_double_type_node)) + (cond (le (abs @0) { build_sinatan_cst (long_double_type_node); }) + (rdiv @0 (sqrts (plus (mult @0 @0) + {build_one_cst (long_double_type_node);}))) + (BUILT_IN_COPYSIGNL { build_one_cst (long_double_type_node); } @0))))))) + +/* Simplify cos(atan(x)) -> 1 / sqrt(x*x + 1). */ + (for coss (COS) + atans (ATAN) + sqrts (SQRT) + (simplify + (coss (atans:s @0)) + (rdiv {build_one_cst (type);} + (sqrts (plus (mult @0 @0) {build_one_cst (type);}))))) + /* cabs(x+0i) or cabs(0+xi) -> abs(x). */ (simplify (CABS (complex:C @0 real_zerop@1)) Index: gcc/real.c =================================================================== --- gcc/real.c (revisão 264058) +++ gcc/real.c (cópia de trabalho) @@ -5279,3 +5279,46 @@ { return HONOR_SIGN_DEPENDENT_ROUNDING (GET_MODE (x)); } + +/* Build real constant used by sin(atan(x)) optimization. + The logic here is as follows: It can be mathematically + shown that sin(atan(x)) = x / sqrt(1 + x*x), but notice + that the second formula has an x*x, which can overflow + if x is big enough. However, x / sqrt(1 + x*x) converges + to 1 for large x. What must be the value of x such that + when computing x / sqrt (1 + x*x) = 1? + + Therefore, we must then solve x / sqrt(1 + x*x) > eps + for x, where eps is the largest number smaller than 1 + representable by the target. Hence, solving for eps + yields that x > eps / sqrt(1 - eps*eps). This eps can + be easily calculated by calling nextafter. Likewise for + the negative x. */ + +void +build_sinatan_real (REAL_VALUE_TYPE * r, tree type) +{ + REAL_VALUE_TYPE eps; + mpfr_t mpfr_eps, mpfr_const1, c, divisor; + const struct real_format * fmt = NULL; + + fmt = type ? REAL_MODE_FORMAT (TYPE_MODE (type)) : NULL; + + mpfr_inits (divisor, mpfr_eps, mpfr_const1, c, NULL); + mpfr_from_real (mpfr_const1, &dconst1, GMP_RNDN); + + real_nextafter (&eps, fmt, &dconst1, &dconst0); + + mpfr_from_real (mpfr_eps, &eps, GMP_RNDZ); + mpfr_mul (divisor, mpfr_eps, mpfr_eps, GMP_RNDU); + mpfr_sub (divisor, mpfr_const1, divisor, GMP_RNDZ); + mpfr_sqrt (divisor, divisor, GMP_RNDZ); + mpfr_div (c, mpfr_eps, divisor, GMP_RNDU); + + real_from_mpfr (r, c, fmt, GMP_RNDU); + + /* For safety reasons. */ + times_pten(r, 1); + + mpfr_clears (divisor, mpfr_eps, mpfr_const1, c, NULL); +} Index: gcc/real.h =================================================================== --- gcc/real.h (revisão 264058) +++ gcc/real.h (cópia de trabalho) @@ -523,4 +523,8 @@ const wide_int_ref &, signop); #endif +/* Build a constant c such that, for every x > c, x / sqrt(1 + x*x) = 1 + in floating point. */ +extern void build_sinatan_real (REAL_VALUE_TYPE *, tree); + #endif /* ! GCC_REAL_H */ Index: gcc/testsuite/gcc.dg/sinatan-1.c =================================================================== --- gcc/testsuite/gcc.dg/sinatan-1.c (nonexistent) +++ gcc/testsuite/gcc.dg/sinatan-1.c (cópia de trabalho) @@ -0,0 +1,136 @@ +/* { dg-do run } */ + +extern float sinf (float); +extern float cosf (float); +extern float atanf (float); +extern float sqrtf (float); +extern float nextafterf (float, float); +extern double sin (double); +extern double cos (double); +extern double atan (double); +extern double sqrt (double); +extern double nextafter (double, double); +extern long double sinl (long double); +extern long double cosl (long double); +extern long double atanl (long double); +extern long double sqrtl (long double); +extern long double nextafterl (long double, long double); + +extern void abort (); + +double __attribute__ ((noinline, optimize("Ofast"))) +sinatan (double x) +{ + return sin (atan (x)); +} + +double __attribute__ ((noinline, optimize("Ofast"))) +cosatan (double x) +{ + return cos (atan (x)); +} + +double __attribute__ ((noinline)) +sinatan_ (double x) +{ + return sin (atan (x)); +} + +float __attribute__ ((noinline, optimize("Ofast"))) +sinatanf(float x) +{ + return sinf (atanf (x)); +} + +float __attribute__ ((noinline, optimize("Ofast"))) +cosatanf(float x) +{ + return cosf (atanf (x)); +} + +float __attribute__ ((noinline)) +sinatanf_ (float x) +{ + return sinf (atanf (x)); +} + +long double __attribute__ ((noinline, optimize("Ofast"))) +sinatanl (long double x) +{ + return sinl (atanl (x)); +} + +long double __attribute__ ((noinline, optimize("Ofast"))) +cosatanl (long double x) +{ + return cosl (atanl (x)); +} + +long double __attribute__ ((noinline)) +sinatanl_ (long double x) +{ + return sinl (atanl (x)); +} + +/* Same logic as implemented in build_sinatan_cst. */ +float +sinatanf_inv (float x) +{ + return (x / sqrtf (1 - x*x))*10; +} + +double +sinatan_inv (double x) +{ + return (x / sqrt (1 - x*x))*10; +} + +long double +sinatanl_inv (long double x) +{ + return (x / sqrtl (1 - x*x))*10; +} + +int +main() +{ + float fc, feps; + double c, eps; + long double lc, leps; + + /* Force move from FPU to memory, otherwise comparison may + fail due to possible more accurate registers (see 387) */ + volatile float fy1, fy2; + volatile double y1, y2; + volatile long double ly1, ly2; + + feps = nextafterf (1.f, 0.f); + eps = nextafter (1. , 0. ); + leps = nextafterl (1.L, 0.L); + + fc = sinatanf_inv (feps); + c = sinatan_inv (eps); + lc = sinatanl_inv (leps); + + fy1 = sinatanf (fc); + fy2 = sinatanf_ (fc); + y1 = sinatan (c); + y2 = sinatan_(c); + ly1 = sinatanl (lc); + ly2 = sinatanl_ (lc); + + if (fy1 != fy2 || y1 != y2 || ly1 != ly2) + abort (); + + fy1 = sinatanf (-fc); + fy2 = sinatanf_ (-fc); + y1 = sinatan (-c); + y2 = sinatan_(-c); + ly1 = sinatanl (-lc); + ly2 = sinatanl_ (-lc); + + if (fy1 != fy2 || y1 != y2 || ly1 != ly2) + abort (); + + return 0; +} Index: gcc/testsuite/gcc.dg/sinatan-2.c =================================================================== --- gcc/testsuite/gcc.dg/sinatan-2.c (nonexistent) +++ gcc/testsuite/gcc.dg/sinatan-2.c (cópia de trabalho) @@ -0,0 +1,59 @@ +/* { dg-do compile } */ +/* { dg-options "-Ofast -fdump-tree-optimized" } */ + +extern float sinf (float); +extern float cosf (float); +extern float atanf (float); +extern double sin (double); +extern double cos (double); +extern double atan (double); +extern long double sinl (long double); +extern long double cosl (long double); +extern long double atanl (long double); + +double __attribute__ ((noinline)) +sinatan_ (double x) +{ + return sin (atan (x)); +} + +double __attribute__ ((noinline)) +cosatan_ (double x) +{ + return cos (atan (x)); +} + +float __attribute__ ((noinline)) +sinatanf_(float x) +{ + return sinf (atanf (x)); +} + +float __attribute__ ((noinline)) +cosatanf_(float x) +{ + return cosf (atanf (x)); +} + +long double __attribute__ ((noinline)) +sinatanl_ (long double x) +{ + return sinl (atanl (x)); +} + +long double __attribute__ ((noinline)) +cosatanl_ (long double x) +{ + return cosl (atanl (x)); +} + +/* There must be no calls to sin, cos, or atan */ +/* {dg-final { scan-tree-dump-not "sin " "optimized" } } */ +/* {dg-final { scan-tree-dump-not "cos " "optimized" } } */ +/* {dg-final { scan-tree-dump-not "atan " "optimized" }} */ +/* {dg-final { scan-tree-dump-not "sinf " "optimized" } } */ +/* {dg-final { scan-tree-dump-not "cosf " "optimized" } } */ +/* {dg-final { scan-tree-dump-not "atanf " "optimized" }} */ +/* {dg-final { scan-tree-dump-not "sinl " "optimized" } } */ +/* {dg-final { scan-tree-dump-not "cosl " "optimized" } } */ +/* {dg-final { scan-tree-dump-not "atanl " "optimized" }} */ Index: gcc/testsuite/gcc.dg/sinatan-3.c =================================================================== --- gcc/testsuite/gcc.dg/sinatan-3.c (nonexistent) +++ gcc/testsuite/gcc.dg/sinatan-3.c (cópia de trabalho) @@ -0,0 +1,65 @@ +/* { dg-do compile } */ +/* { dg-options "-Ofast -fdump-tree-optimized" } */ + +extern float sinf (float); +extern float cosf (float); +extern float atanf (float); +extern double sin (double); +extern double cos (double); +extern double atan (double); +extern long double sinl (long double); +extern long double cosl (long double); +extern long double atanl (long double); + +float __attribute__ ((noinline)) +cosatanf_(float x) +{ + float atg = atanf(x); + return cosf(atg) + atg; +} + +double __attribute__ ((noinline)) +cosatan_(double x) +{ + double atg = atan(x); + return cos(atg) + atg; +} + +long double __attribute__ ((noinline)) +cosatanl_(long double x) +{ + long double atg = atanl(x); + return cosl(atg) + atg; +} + +float __attribute__ ((noinline)) +sinatanf_(float x) +{ + float atg = atanf(x); + return sinf(atg) + atg; +} + +double __attribute__ ((noinline)) +sinatan_(double x) +{ + double atg = atan(x); + return sin(atg) + atg; +} + +long double __attribute__ ((noinline)) +sinatanl_(long double x) +{ + long double atg = atanl(x); + return sinl(atg) + atg; +} + +/* There should be calls to both sin and atan */ +/* { dg-final { scan-tree-dump "cos " "optimized" } } */ +/* { dg-final { scan-tree-dump "sin " "optimized" } } */ +/* { dg-final { scan-tree-dump "atan " "optimized" } } */ +/* { dg-final { scan-tree-dump "cosf " "optimized" } } */ +/* { dg-final { scan-tree-dump "sinf " "optimized" } } */ +/* { dg-final { scan-tree-dump "atanf " "optimized" } } */ +/* { dg-final { scan-tree-dump "cosl " "optimized" } } */ +/* { dg-final { scan-tree-dump "sinl " "optimized" } } */ +/* { dg-final { scan-tree-dump "atanl " "optimized" } } */ Index: gcc/tree.c =================================================================== --- gcc/tree.c (revisão 264058) +++ gcc/tree.c (cópia de trabalho) @@ -2359,7 +2359,16 @@ } } +tree +build_sinatan_cst (tree type) +{ + REAL_VALUE_TYPE cst; + build_sinatan_real (&cst, type); + return build_real (type, cst); +} + + /* Build a BINFO with LEN language slots. */ tree Index: gcc/tree.h =================================================================== --- gcc/tree.h (revisão 264058) +++ gcc/tree.h (cópia de trabalho) @@ -4158,6 +4158,7 @@ extern tree build_minus_one_cst (tree); extern tree build_all_ones_cst (tree); extern tree build_zero_cst (tree); +extern tree build_sinatan_cst (tree); extern tree build_string (int, const char *); extern tree build_poly_int_cst (tree, const poly_wide_int_ref &); extern tree build_tree_list (tree, tree CXX_MEM_STAT_INFO);
Fixed the issues pointed by the previous discussions. Closes PR86829. Adds substitution rules for sin(atan(x)) and cos(atan(x)), being careful with overflow issues by constructing a assumed convergence constant (see comment in real.c). 2018-09-03 Giuliano Belinassi <giuliano.belinassi@usp.br> * match.pd: add simplification rules to sin(atan(x)) and cos(atan(x)). * real.c: add code for assumed convergence constant to sin(atan(x)). * real.h: allows the added code from real.c to be called externally. * tree.c: add code for bulding nodes with the convergence constant. * tree.h: allows the added code from tree.c to be called externally. * sinatan-1.c: tests assumed convergence constant. * sinatan-2.c: tests simplification rule. * sinatan-3.c: likewise. There seems to be no broken tests in trunk that are related to this modification.