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.belina...@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
