[PATCH] Optimize sin(atan(x)), take 2

[Richard Biener]

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