diff mbox

[tree-ssa-math-opts] Expand pow (x, CONST) using square roots when possible

Message ID 554CD1B3.1060300@foss.arm.com
State New
Headers show

Commit Message

Kyrill Tkachov May 8, 2015, 3:09 p.m. UTC
On 08/05/15 14:56, Kyrill Tkachov wrote:
> On 08/05/15 11:18, Richard Biener wrote:
>> On Fri, May 1, 2015 at 6:02 PM, Kyrill Tkachov
>> <kyrylo.tkachov@foss.arm.com> wrote:
>>> Hi all,
>>>
>>> GCC has some logic to expand calls to pow (x, 0.75), pow (0.25) and pow (x,
>>> (int)k + 0.5)
>>> using square roots. So, for the above examples it would generate sqrt (x) *
>>> sqrt (sqrt (x)),
>>> sqrt (sqrt (x)) and powi (x, k) * sqrt (x) (assuming k > 0. For k < 0 it
>>> will calculate the
>>> reciprocal of that).
>>>
>>> However, the implementation of these optimisations is done on a bit of an
>>> ad-hoc basis with
>>> the 0.25, 0.5, 0.75 cases hardcoded.
>>> Judging by
>>> https://gcc.gnu.org/wiki/summit2010?action=AttachFile&do=get&target=meissner2.pdf
>>> these are the most commonly used exponents (at least in SPEC ;))
>>>
>>> This patch generalises this optimisation into a (hopefully) more robust
>>> algorithm.
>>> In particular, it expands calls to pow (x, CST) by expanding the integer
>>> part of CST
>>> using a powi, like it does already, and then expanding the fractional part
>>> as a product
>>> of repeated applications of a square root if the fractional part can be
>>> expressed
>>> as a multiple of a power of 0.5.
>>>
>>> I try to explain the algorithm in more detail in the comments in the patch
>>> but, for example:
>>>
>>> pow (x, 5.625) is not currently handled, but with this patch will be
>>> expanded
>>> to powi (x, 5) * sqrt (x) * sqrt (sqrt (sqrt (x))) because 5.625 == 5.0 +
>>> 0.5 + 0.5**3
>>>
>>> Negative exponents are handled in either of two ways, depending on the
>>> exponent value:
>>> * Using a simple reciprocal.
>>>     For example:
>>>     pow (x, -5.625) == 1.0 / pow (x, 5.625)
>>>       --> 1.0 / (powi (x, 5) * sqrt (x) * sqrt (sqrt (sqrt (x))))
>>>
>>> * For pow (x, EXP) with negative exponent EXP with integer part INT and
>>> fractional part FRAC:
>>> pow (1.0 - FRAC) / powi (ceil (abs (EXP))).
>>>     For example:
>>>     pow (x, -5.875) == pow (x, 0.125) / powi (X, 6)
>>>       --> sqrt (sqrt (sqrt (x))) / (powi (x, 6))
>>>
>>>
>>> Since hardware square root instructions tend to be expensive, we may want to
>>> reduce the number
>>> of square roots we are willing to calculate. Since we reuse intermediate
>>> square root results,
>>> this boils down to restricting the depth of the square root chains. In all
>>> the examples above
>>> that depth is 3. I've made this maximum depth parametrisable in params.def.
>>> By adjusting that
>>> parameter we can adjust the resolution of this optimisation. So, if it's set
>>> to '4' then we
>>> will synthesize every exponent that is a multiple of 0.5**4 == 0.0625,
>>> including negative
>>> multiples. Currently, GCC will not try to expand negative multiples of
>>> anything else than 0.5
>>>
>>> I have tried to keep the existing functionality intact and activate this
>>> only for
>>> -funsafe-math-optimizations and only when the target has a sqrt instruction.
>>>    An exception to that is pow (x, 0.5) which we prefer to transform to sqrt
>>> even
>>> when a hardware sqrt is not available, presumably because the library
>>> function for
>>> sqrt is usually faster than pow (?).
>> Yes.  It's also a safe transform - which you seem to put under
>> flag_unsafe_math_optimizations only with your patch.
>>
>> It would be clearer to just leave the special-case
>>
>> -  /* Optimize pow(x,0.5) = sqrt(x).  This replacement is always safe
>> -     unless signed zeros must be maintained.  pow(-0,0.5) = +0, while
>> -     sqrt(-0) = -0.  */
>> -  if (sqrtfn
>> -      && REAL_VALUES_EQUAL (c, dconsthalf)
>> -      && !HONOR_SIGNED_ZEROS (mode))
>> -    return build_and_insert_call (gsi, loc, sqrtfn, arg0);
>>
>> in as-is.
> Ok, I'll leave that case explicit.
>
>> You also removed the Os constraint which you should put back in.
>> Basically if !optimize_function_for_speed_p then generate at most
>> two calls to sqrt (iff the HW has a sqrt instruction).
> I tried to move that logic into expand_with_sqrts but
> I'll move it outside it. It seems that this boils down to
> only 0.25, as any other 2xsqrt chain will also involve a
> multiply or a divide which we currently avoid.
>
>> You fail to add a testcase that checks that the optimization applies.
> I'll add one to scan the sincos dump.
> I notice that we don't have a testuite check that the target has
> a hw sqrt instructions. Would you like me to add one? Or can I make
> the testcase aarch64-specific?
>
>> Otherwise the idea looks good though there must be a better way
>> to compute the series than by using real-arithmetic and forcefully
>> trying out all possibilities...
> I get that feeling too. What I need is not only a way
> of figuring out if the fractional part of the exponent can be
> represented in this way, but also compute the depth of the
> sqrt chain and the number of multiplies...
> That being said, the current approach is O(maximum depth) and
> I don't expect the depth to go much beyond 3 or 4 in practice.
>
> Thanks for looking at it!
> I'll respin the patch.

And here it is, with my above comments implemented.
Bootstrapped on x86_64 and tested on aarch64.
Full testing on arm and aarch64 ongoing.

Is this ok if testing comes clean?

Thanks,
Kyrill

> Kyrill
>
>> Richard.
>>
>>> Having seen the glibc implementation of a fully IEEE-754-compliant pow
>>> function, I think we
>>> would prefer synthesising the pow call whenever we can for -ffast-math.
>>>
>>> I have seen this optimisation trigger a few times in SPEC2k6, in particular
>>> in 447.dealII
>>> and 481.wrf where it replaced calls to powf (x, -0.25), pow (x, 0.125) and
>>> pow (x, 0.875)
>>> with square roots, multiplies and, in the case of -0.25, divides.
>>> On 481.wrf I saw it remove a total of 22 out of 322 calls to pow
>>>
>>> On 481.wrf on aarch64 I saw about a 1% improvement.
>>> The cycle count on x86_64 was also smaller, but not by a significant amount
>>> (the same calls to
>>> pow were eliminated).
>>>
>>> In general, I think this can shine if multiple expandable calls to pow
>>> appear together.
>>> So, for example for code:
>>> double
>>> baz (double a)
>>> {
>>>     return __builtin_pow (a, -1.25) + __builtin_pow (a, 5.75) - __builtin_pow
>>> (a, 3.375);
>>> }
>>>
>>> we can generate:
>>> baz:
>>>           fsqrt   d3, d0
>>>           fmul    d4, d0, d0
>>>           fmov    d5, 1.0e+0
>>>           fmul    d6, d0, d4
>>>           fsqrt   d2, d3
>>>           fmul    d1, d0, d2
>>>           fsqrt   d0, d2
>>>           fmul    d3, d3, d2
>>>           fdiv    d1, d5, d1
>>>           fmul    d3, d3, d6
>>>           fmul    d2, d2, d0
>>>           fmadd   d0, d4, d3, d1
>>>           fmsub   d0, d6, d2, d0
>>>           ret
>>>
>>> reusing the sqrt results and doing more optimisations rather than the
>>> current:
>>> baz:
>>>           stp     x29, x30, [sp, -48]!
>>>           fmov    d1, -1.25e+0
>>>           add     x29, sp, 0
>>>           stp     d8, d9, [sp, 16]
>>>           fmov    d9, d0
>>>           str     d10, [sp, 32]
>>>           bl      pow
>>>           fmov    d8, d0
>>>           fmov    d0, d9
>>>           fmov    d1, 5.75e+0
>>>           bl      pow
>>>           fmov    d10, d0
>>>           fmov    d0, d9
>>>           fmov    d1, 3.375e+0
>>>           bl      pow
>>>           fadd    d8, d8, d10
>>>           ldr     d10, [sp, 32]
>>>           fsub    d0, d8, d0
>>>           ldp     d8, d9, [sp, 16]
>>>           ldp     x29, x30, [sp], 48
>>>           ret
>>>
>>>
>>> Of course gcc could already do that if the exponents were to fall in the set
>>> {0.25, 0.75, k+0.5}
>>> but with this patch that set can be greatly expanded.
>>>
>>> I suppose if we're really lucky we might even open up new vectorisation
>>> opportunities.
>>> For example:
>>> void
>>> vecfoo (float *a, float *b)
>>> {
>>>     for (int i = 0; i < N; i++)
>>>       a[i] = __builtin_powf (a[i], 1.25f) - __builtin_powf (b[i], 3.625);
>>> }
>>>
>>> will now be vectorisable if we have a hardware vector sqrt instruction.
>>> Though I'm not sure
>>> how often this would appear in real code.
>>>
>>>
>>> I would like advice on some implementation details:
>>> - The new function representable_as_half_series_p is used to check if the
>>> fractional
>>> part of an exponent can be represented as a multiple of a power of 0.5. It
>>> does so
>>> by using REAL_VALUE_TYPE arithmetic and a loop. I wonder whether there is a
>>> smarter
>>> way of doing this, considering that REAL_VALUE_TYPE holds the bit
>>> representation of the
>>> floating point number?
>>>
>>> - Are there any correctness cases that I may have missed? This patch gates
>>> the optimisation
>>> on -funsafe-math-optimizations, but maybe there are some edge cases that I
>>> missed?
>>>
>>> - What should be the default value of the max-pow-sqrt-depth param? In this
>>> patch it's 5 which
>>> on second thought strikes me as a bit aggressive. To catch exponents that
>>> are multiples of
>>> 0.25 we need it to be 2. For multiples of 0.125 it has to be 3 etc... I
>>> suppose it depends on
>>> how likely more fine-grained powers are to appear in real code, how
>>> expensive the C library
>>> implementation of pow is, and how expensive are the sqrt instructions in
>>> hardware.
>>>
>>>
>>> Bootstrapped and tested on x86_64, aarch64, arm (pending)
>>> SPEC2k6 built and ran fine. Can anyone suggest anything other codebase that
>>> might
>>> be of interest to look at?
>>>
>>> Thanks,
>>> Kyrill
>>>
>>> 2015-05-01  Kyrylo Tkachov  <kyrylo.tkachov@arm.com>
>>>
>>>       * params.def (PARAM_MAX_POW_SQRT_DEPTH): New param.
>>>       * tree-ssa-math-opts.c: Include params.h
>>>       (pow_synth_sqrt_info): New struct.
>>>       (representable_as_half_series_p): New function.
>>>       (get_fn_chain): Likewise.
>>>       (print_nested_fn): Likewise.
>>>       (dump_fractional_sqrt_sequence): Likewise.
>>>       (dump_integer_part): Likewise.
>>>       (expand_pow_as_sqrts): Likewise.
>>>       (gimple_expand_builtin_pow): Use above to attempt to expand
>>>       pow as series of square roots.  Removed now unused variables.
>>>
>>> 2015-05-01  Kyrylo Tkachov  <kyrylo.tkachov@arm.com>
>>>
>>>       * gcc.dg/pow-sqrt.x: New file.
>>>       * gcc.dg/pow-sqrt-1.c: New test.
>>>       * gcc.dg/pow-sqrt-2.c: Likewise.
>>>       * gcc.dg/pow-sqrt-3.c: Likewise.

Comments

Richard Biener May 13, 2015, 11:27 a.m. UTC | #1
On Fri, May 8, 2015 at 5:09 PM, Kyrill Tkachov
<kyrylo.tkachov@foss.arm.com> wrote:
>
> On 08/05/15 14:56, Kyrill Tkachov wrote:
>>
>> On 08/05/15 11:18, Richard Biener wrote:
>>>
>>> On Fri, May 1, 2015 at 6:02 PM, Kyrill Tkachov
>>> <kyrylo.tkachov@foss.arm.com> wrote:
>>>>
>>>> Hi all,
>>>>
>>>> GCC has some logic to expand calls to pow (x, 0.75), pow (0.25) and pow
>>>> (x,
>>>> (int)k + 0.5)
>>>> using square roots. So, for the above examples it would generate sqrt
>>>> (x) *
>>>> sqrt (sqrt (x)),
>>>> sqrt (sqrt (x)) and powi (x, k) * sqrt (x) (assuming k > 0. For k < 0 it
>>>> will calculate the
>>>> reciprocal of that).
>>>>
>>>> However, the implementation of these optimisations is done on a bit of
>>>> an
>>>> ad-hoc basis with
>>>> the 0.25, 0.5, 0.75 cases hardcoded.
>>>> Judging by
>>>>
>>>> https://gcc.gnu.org/wiki/summit2010?action=AttachFile&do=get&target=meissner2.pdf
>>>> these are the most commonly used exponents (at least in SPEC ;))
>>>>
>>>> This patch generalises this optimisation into a (hopefully) more robust
>>>> algorithm.
>>>> In particular, it expands calls to pow (x, CST) by expanding the integer
>>>> part of CST
>>>> using a powi, like it does already, and then expanding the fractional
>>>> part
>>>> as a product
>>>> of repeated applications of a square root if the fractional part can be
>>>> expressed
>>>> as a multiple of a power of 0.5.
>>>>
>>>> I try to explain the algorithm in more detail in the comments in the
>>>> patch
>>>> but, for example:
>>>>
>>>> pow (x, 5.625) is not currently handled, but with this patch will be
>>>> expanded
>>>> to powi (x, 5) * sqrt (x) * sqrt (sqrt (sqrt (x))) because 5.625 == 5.0
>>>> +
>>>> 0.5 + 0.5**3
>>>>
>>>> Negative exponents are handled in either of two ways, depending on the
>>>> exponent value:
>>>> * Using a simple reciprocal.
>>>>     For example:
>>>>     pow (x, -5.625) == 1.0 / pow (x, 5.625)
>>>>       --> 1.0 / (powi (x, 5) * sqrt (x) * sqrt (sqrt (sqrt (x))))
>>>>
>>>> * For pow (x, EXP) with negative exponent EXP with integer part INT and
>>>> fractional part FRAC:
>>>> pow (1.0 - FRAC) / powi (ceil (abs (EXP))).
>>>>     For example:
>>>>     pow (x, -5.875) == pow (x, 0.125) / powi (X, 6)
>>>>       --> sqrt (sqrt (sqrt (x))) / (powi (x, 6))
>>>>
>>>>
>>>> Since hardware square root instructions tend to be expensive, we may
>>>> want to
>>>> reduce the number
>>>> of square roots we are willing to calculate. Since we reuse intermediate
>>>> square root results,
>>>> this boils down to restricting the depth of the square root chains. In
>>>> all
>>>> the examples above
>>>> that depth is 3. I've made this maximum depth parametrisable in
>>>> params.def.
>>>> By adjusting that
>>>> parameter we can adjust the resolution of this optimisation. So, if it's
>>>> set
>>>> to '4' then we
>>>> will synthesize every exponent that is a multiple of 0.5**4 == 0.0625,
>>>> including negative
>>>> multiples. Currently, GCC will not try to expand negative multiples of
>>>> anything else than 0.5
>>>>
>>>> I have tried to keep the existing functionality intact and activate this
>>>> only for
>>>> -funsafe-math-optimizations and only when the target has a sqrt
>>>> instruction.
>>>>    An exception to that is pow (x, 0.5) which we prefer to transform to
>>>> sqrt
>>>> even
>>>> when a hardware sqrt is not available, presumably because the library
>>>> function for
>>>> sqrt is usually faster than pow (?).
>>>
>>> Yes.  It's also a safe transform - which you seem to put under
>>> flag_unsafe_math_optimizations only with your patch.
>>>
>>> It would be clearer to just leave the special-case
>>>
>>> -  /* Optimize pow(x,0.5) = sqrt(x).  This replacement is always safe
>>> -     unless signed zeros must be maintained.  pow(-0,0.5) = +0, while
>>> -     sqrt(-0) = -0.  */
>>> -  if (sqrtfn
>>> -      && REAL_VALUES_EQUAL (c, dconsthalf)
>>> -      && !HONOR_SIGNED_ZEROS (mode))
>>> -    return build_and_insert_call (gsi, loc, sqrtfn, arg0);
>>>
>>> in as-is.
>>
>> Ok, I'll leave that case explicit.
>>
>>> You also removed the Os constraint which you should put back in.
>>> Basically if !optimize_function_for_speed_p then generate at most
>>> two calls to sqrt (iff the HW has a sqrt instruction).
>>
>> I tried to move that logic into expand_with_sqrts but
>> I'll move it outside it. It seems that this boils down to
>> only 0.25, as any other 2xsqrt chain will also involve a
>> multiply or a divide which we currently avoid.
>>
>>> You fail to add a testcase that checks that the optimization applies.
>>
>> I'll add one to scan the sincos dump.
>> I notice that we don't have a testuite check that the target has
>> a hw sqrt instructions. Would you like me to add one? Or can I make
>> the testcase aarch64-specific?

Would be great to have a testsuite check for this.

>>
>>> Otherwise the idea looks good though there must be a better way
>>> to compute the series than by using real-arithmetic and forcefully
>>> trying out all possibilities...
>>
>> I get that feeling too. What I need is not only a way
>> of figuring out if the fractional part of the exponent can be
>> represented in this way, but also compute the depth of the
>> sqrt chain and the number of multiplies...
>> That being said, the current approach is O(maximum depth) and
>> I don't expect the depth to go much beyond 3 or 4 in practice.
>>
>> Thanks for looking at it!
>> I'll respin the patch.
>
>
> And here it is, with my above comments implemented.
> Bootstrapped on x86_64 and tested on aarch64.
> Full testing on arm and aarch64 ongoing.
>
> Is this ok if testing comes clean?

Ok.

Thanks,
Richard.

> Thanks,
> Kyrill
>
>
>> Kyrill
>>
>>> Richard.
>>>
>>>> Having seen the glibc implementation of a fully IEEE-754-compliant pow
>>>> function, I think we
>>>> would prefer synthesising the pow call whenever we can for -ffast-math.
>>>>
>>>> I have seen this optimisation trigger a few times in SPEC2k6, in
>>>> particular
>>>> in 447.dealII
>>>> and 481.wrf where it replaced calls to powf (x, -0.25), pow (x, 0.125)
>>>> and
>>>> pow (x, 0.875)
>>>> with square roots, multiplies and, in the case of -0.25, divides.
>>>> On 481.wrf I saw it remove a total of 22 out of 322 calls to pow
>>>>
>>>> On 481.wrf on aarch64 I saw about a 1% improvement.
>>>> The cycle count on x86_64 was also smaller, but not by a significant
>>>> amount
>>>> (the same calls to
>>>> pow were eliminated).
>>>>
>>>> In general, I think this can shine if multiple expandable calls to pow
>>>> appear together.
>>>> So, for example for code:
>>>> double
>>>> baz (double a)
>>>> {
>>>>     return __builtin_pow (a, -1.25) + __builtin_pow (a, 5.75) -
>>>> __builtin_pow
>>>> (a, 3.375);
>>>> }
>>>>
>>>> we can generate:
>>>> baz:
>>>>           fsqrt   d3, d0
>>>>           fmul    d4, d0, d0
>>>>           fmov    d5, 1.0e+0
>>>>           fmul    d6, d0, d4
>>>>           fsqrt   d2, d3
>>>>           fmul    d1, d0, d2
>>>>           fsqrt   d0, d2
>>>>           fmul    d3, d3, d2
>>>>           fdiv    d1, d5, d1
>>>>           fmul    d3, d3, d6
>>>>           fmul    d2, d2, d0
>>>>           fmadd   d0, d4, d3, d1
>>>>           fmsub   d0, d6, d2, d0
>>>>           ret
>>>>
>>>> reusing the sqrt results and doing more optimisations rather than the
>>>> current:
>>>> baz:
>>>>           stp     x29, x30, [sp, -48]!
>>>>           fmov    d1, -1.25e+0
>>>>           add     x29, sp, 0
>>>>           stp     d8, d9, [sp, 16]
>>>>           fmov    d9, d0
>>>>           str     d10, [sp, 32]
>>>>           bl      pow
>>>>           fmov    d8, d0
>>>>           fmov    d0, d9
>>>>           fmov    d1, 5.75e+0
>>>>           bl      pow
>>>>           fmov    d10, d0
>>>>           fmov    d0, d9
>>>>           fmov    d1, 3.375e+0
>>>>           bl      pow
>>>>           fadd    d8, d8, d10
>>>>           ldr     d10, [sp, 32]
>>>>           fsub    d0, d8, d0
>>>>           ldp     d8, d9, [sp, 16]
>>>>           ldp     x29, x30, [sp], 48
>>>>           ret
>>>>
>>>>
>>>> Of course gcc could already do that if the exponents were to fall in the
>>>> set
>>>> {0.25, 0.75, k+0.5}
>>>> but with this patch that set can be greatly expanded.
>>>>
>>>> I suppose if we're really lucky we might even open up new vectorisation
>>>> opportunities.
>>>> For example:
>>>> void
>>>> vecfoo (float *a, float *b)
>>>> {
>>>>     for (int i = 0; i < N; i++)
>>>>       a[i] = __builtin_powf (a[i], 1.25f) - __builtin_powf (b[i],
>>>> 3.625);
>>>> }
>>>>
>>>> will now be vectorisable if we have a hardware vector sqrt instruction.
>>>> Though I'm not sure
>>>> how often this would appear in real code.
>>>>
>>>>
>>>> I would like advice on some implementation details:
>>>> - The new function representable_as_half_series_p is used to check if
>>>> the
>>>> fractional
>>>> part of an exponent can be represented as a multiple of a power of 0.5.
>>>> It
>>>> does so
>>>> by using REAL_VALUE_TYPE arithmetic and a loop. I wonder whether there
>>>> is a
>>>> smarter
>>>> way of doing this, considering that REAL_VALUE_TYPE holds the bit
>>>> representation of the
>>>> floating point number?
>>>>
>>>> - Are there any correctness cases that I may have missed? This patch
>>>> gates
>>>> the optimisation
>>>> on -funsafe-math-optimizations, but maybe there are some edge cases that
>>>> I
>>>> missed?
>>>>
>>>> - What should be the default value of the max-pow-sqrt-depth param? In
>>>> this
>>>> patch it's 5 which
>>>> on second thought strikes me as a bit aggressive. To catch exponents
>>>> that
>>>> are multiples of
>>>> 0.25 we need it to be 2. For multiples of 0.125 it has to be 3 etc... I
>>>> suppose it depends on
>>>> how likely more fine-grained powers are to appear in real code, how
>>>> expensive the C library
>>>> implementation of pow is, and how expensive are the sqrt instructions in
>>>> hardware.
>>>>
>>>>
>>>> Bootstrapped and tested on x86_64, aarch64, arm (pending)
>>>> SPEC2k6 built and ran fine. Can anyone suggest anything other codebase
>>>> that
>>>> might
>>>> be of interest to look at?
>>>>
>>>> Thanks,
>>>> Kyrill
>>>>
>>>> 2015-05-01  Kyrylo Tkachov  <kyrylo.tkachov@arm.com>
>>>>
>>>>       * params.def (PARAM_MAX_POW_SQRT_DEPTH): New param.
>>>>       * tree-ssa-math-opts.c: Include params.h
>>>>       (pow_synth_sqrt_info): New struct.
>>>>       (representable_as_half_series_p): New function.
>>>>       (get_fn_chain): Likewise.
>>>>       (print_nested_fn): Likewise.
>>>>       (dump_fractional_sqrt_sequence): Likewise.
>>>>       (dump_integer_part): Likewise.
>>>>       (expand_pow_as_sqrts): Likewise.
>>>>       (gimple_expand_builtin_pow): Use above to attempt to expand
>>>>       pow as series of square roots.  Removed now unused variables.
>>>>
>>>> 2015-05-01  Kyrylo Tkachov  <kyrylo.tkachov@arm.com>
>>>>
>>>>       * gcc.dg/pow-sqrt.x: New file.
>>>>       * gcc.dg/pow-sqrt-1.c: New test.
>>>>       * gcc.dg/pow-sqrt-2.c: Likewise.
>>>>       * gcc.dg/pow-sqrt-3.c: Likewise.
>
>
diff mbox

Patch

commit 1c55f44e3294b17ba113032c2d05bb9e09c8fc69
Author: Kyrylo Tkachov <kyrylo.tkachov@arm.com>
Date:   Wed Apr 29 13:27:07 2015 +0100

    [tree-ssa-math-opts] Expand pow (x, CONST) using square roots when possible

diff --git a/gcc/params.def b/gcc/params.def
index 48b39a2..3e4ba3a 100644
--- a/gcc/params.def
+++ b/gcc/params.def
@@ -262,6 +262,14 @@  DEFPARAM(PARAM_MAX_HOIST_DEPTH,
 	 "Maximum depth of search in the dominator tree for expressions to hoist",
 	 30, 0, 0)
 
+
+/* When synthesizing expnonentiation by a real constant operations using square
+   roots, this controls how deep sqrt chains we are willing to generate.  */
+DEFPARAM(PARAM_MAX_POW_SQRT_DEPTH,
+	 "max-pow-sqrt-depth",
+	 "Maximum depth of sqrt chains to use when synthesizing exponentiation by a real constant",
+	 5, 1, 32)
+
 /* This parameter limits the number of insns in a loop that will be unrolled,
    and by how much the loop is unrolled.
 
diff --git a/gcc/testsuite/gcc.dg/pow-sqrt-1.c b/gcc/testsuite/gcc.dg/pow-sqrt-1.c
new file mode 100644
index 0000000..0793b6f
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/pow-sqrt-1.c
@@ -0,0 +1,6 @@ 
+/* { dg-do run } */
+/* { dg-options "-O2 -ffast-math --param max-pow-sqrt-depth=5" } */
+
+#define EXPN (-6 * (0.5*0.5*0.5*0.5))
+
+#include "pow-sqrt.x"
diff --git a/gcc/testsuite/gcc.dg/pow-sqrt-2.c b/gcc/testsuite/gcc.dg/pow-sqrt-2.c
new file mode 100644
index 0000000..b2fada4
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/pow-sqrt-2.c
@@ -0,0 +1,5 @@ 
+/* { dg-do run } */
+/* { dg-options "-O2 -ffast-math --param max-pow-sqrt-depth=5" } */
+
+#define EXPN (-5.875)
+#include "pow-sqrt.x"
diff --git a/gcc/testsuite/gcc.dg/pow-sqrt-3.c b/gcc/testsuite/gcc.dg/pow-sqrt-3.c
new file mode 100644
index 0000000..18c7231
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/pow-sqrt-3.c
@@ -0,0 +1,5 @@ 
+/* { dg-do run } */
+/* { dg-options "-O2 -ffast-math --param max-pow-sqrt-depth=3" } */
+
+#define EXPN (1.25)
+#include "pow-sqrt.x"
diff --git a/gcc/testsuite/gcc.dg/pow-sqrt.x b/gcc/testsuite/gcc.dg/pow-sqrt.x
new file mode 100644
index 0000000..bd744c6
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/pow-sqrt.x
@@ -0,0 +1,30 @@ 
+
+extern void abort (void);
+
+
+__attribute__((noinline)) double
+real_pow (double x, double pow_exp)
+{
+  return __builtin_pow (x, pow_exp);
+}
+
+#define EPS (0.000000000000000000001)
+
+#define SYNTH_POW(X, Y) __builtin_pow (X, Y)
+volatile double arg;
+
+int
+main (void)
+{
+  double i_arg = 0.1;
+
+  for (arg = i_arg; arg < 100.0; arg += 1.0)
+    {
+      double synth_res = SYNTH_POW (arg, EXPN);
+      double real_res = real_pow (arg, EXPN);
+
+      if (__builtin_abs (SYNTH_POW (arg, EXPN) - real_pow (arg, EXPN)) > EPS)
+	abort ();
+    }
+  return 0;
+}
diff --git a/gcc/testsuite/gcc.target/aarch64/pow-sqrt-synth-1.c b/gcc/testsuite/gcc.target/aarch64/pow-sqrt-synth-1.c
new file mode 100644
index 0000000..52514fb
--- /dev/null
+++ b/gcc/testsuite/gcc.target/aarch64/pow-sqrt-synth-1.c
@@ -0,0 +1,38 @@ 
+/* { dg-do compile } */
+/* { dg-options "-fdump-tree-sincos -Ofast --param max-pow-sqrt-depth=8" } */
+
+
+double
+foo (double a)
+{
+  return __builtin_pow (a, -5.875);
+}
+
+double
+foof (double a)
+{
+  return __builtin_pow (a, 0.75f);
+}
+
+double
+bar (double a)
+{
+  return __builtin_pow (a, 1.0 + 0.00390625);
+}
+
+double
+baz (double a)
+{
+  return __builtin_pow (a, -1.25) + __builtin_pow (a, 5.75) - __builtin_pow (a, 3.375);
+}
+
+#define N 256
+void
+vecfoo (double *a)
+{
+  for (int i = 0; i < N; i++)
+    a[i] = __builtin_pow (a[i], 1.25);
+}
+
+/* { dg-final { scan-tree-dump-times "synthesizing" 7 "sincos" } } */
+/* { dg-final { cleanup-tree-dump "sincos" } } */
\ No newline at end of file
diff --git a/gcc/tree-ssa-math-opts.c b/gcc/tree-ssa-math-opts.c
index c22a677..11288b1 100644
--- a/gcc/tree-ssa-math-opts.c
+++ b/gcc/tree-ssa-math-opts.c
@@ -143,6 +143,7 @@  along with GCC; see the file COPYING3.  If not see
 #include "target.h"
 #include "gimple-pretty-print.h"
 #include "builtins.h"
+#include "params.h"
 
 /* FIXME: RTL headers have to be included here for optabs.  */
 #include "rtl.h"		/* Because optabs.h wants enum rtx_code.  */
@@ -1148,6 +1149,357 @@  build_and_insert_cast (gimple_stmt_iterator *gsi, location_t loc,
   return result;
 }
 
+struct pow_synth_sqrt_info
+{
+  bool *factors;
+  unsigned int deepest;
+  unsigned int num_mults;
+};
+
+/* Return true iff the real value C can be represented as a
+   sum of powers of 0.5 up to N.  That is:
+   C == SUM<i from 1..N> (a[i]*(0.5**i)) where a[i] is either 0 or 1.
+   Record in INFO the various parameters of the synthesis algorithm such
+   as the factors a[i], the maximum 0.5 power and the number of
+   multiplications that will be required.  */
+
+bool
+representable_as_half_series_p (REAL_VALUE_TYPE c, unsigned n,
+				 struct pow_synth_sqrt_info *info)
+{
+  REAL_VALUE_TYPE factor = dconsthalf;
+  REAL_VALUE_TYPE remainder = c;
+
+  info->deepest = 0;
+  info->num_mults = 0;
+  memset (info->factors, 0, n * sizeof (bool));
+
+  for (unsigned i = 0; i < n; i++)
+    {
+      REAL_VALUE_TYPE res;
+
+      /* If something inexact happened bail out now.  */
+      if (REAL_ARITHMETIC (res, MINUS_EXPR, remainder, factor))
+	return false;
+
+      /* We have hit zero.  The number is representable as a sum
+         of powers of 0.5.  */
+      if (REAL_VALUES_EQUAL (res, dconst0))
+	{
+	  info->factors[i] = true;
+	  info->deepest = i + 1;
+	  return true;
+	}
+      else if (!REAL_VALUE_NEGATIVE (res))
+	{
+	  remainder = res;
+	  info->factors[i] = true;
+	  info->num_mults++;
+	}
+      else
+	info->factors[i] = false;
+
+      REAL_ARITHMETIC (factor, MULT_EXPR, factor, dconsthalf);
+    }
+  return false;
+}
+
+/* Return the tree corresponding to FN being applied
+   to ARG N times at GSI and LOC.
+   Look up previous results from CACHE if need be.
+   cache[0] should contain just plain ARG i.e. FN applied to ARG 0 times.  */
+
+static tree
+get_fn_chain (tree arg, unsigned int n, gimple_stmt_iterator *gsi,
+	      tree fn, location_t loc, tree *cache)
+{
+  tree res = cache[n];
+  if (!res)
+    {
+      tree prev = get_fn_chain (arg, n - 1, gsi, fn, loc, cache);
+      res = build_and_insert_call (gsi, loc, fn, prev);
+      cache[n] = res;
+    }
+
+  return res;
+}
+
+/* Print to STREAM the repeated application of function FNAME to ARG
+   N times.  So, for FNAME = "foo", ARG = "x", N = 2 it would print:
+   "foo (foo (x))".  */
+
+static void
+print_nested_fn (FILE* stream, const char *fname, const char* arg,
+		 unsigned int n)
+{
+  if (n == 0)
+    fprintf (stream, "%s", arg);
+  else
+    {
+      fprintf (stream, "%s (", fname);
+      print_nested_fn (stream, fname, arg, n - 1);
+      fprintf (stream, ")");
+    }
+}
+
+/* Print to STREAM the fractional sequence of sqrt chains
+   applied to ARG, described by INFO.  Used for the dump file.  */
+
+static void
+dump_fractional_sqrt_sequence (FILE *stream, const char *arg,
+			        struct pow_synth_sqrt_info *info)
+{
+  for (unsigned int i = 0; i < info->deepest; i++)
+    {
+      bool is_set = info->factors[i];
+      if (is_set)
+	{
+	  print_nested_fn (stream, "sqrt", arg, i + 1);
+	  if (i != info->deepest - 1)
+	    fprintf (stream, " * ");
+	}
+    }
+}
+
+/* Print to STREAM a representation of raising ARG to an integer
+   power N.  Used for the dump file.  */
+
+static void
+dump_integer_part (FILE *stream, const char* arg, HOST_WIDE_INT n)
+{
+  if (n > 1)
+    fprintf (stream, "powi (%s, " HOST_WIDE_INT_PRINT_DEC ")", arg, n);
+  else if (n == 1)
+    fprintf (stream, "%s", arg);
+}
+
+/* Attempt to synthesize a POW[F] (ARG0, ARG1) call using chains of
+   square roots.  Place at GSI and LOC.  Limit the maximum depth
+   of the sqrt chains to MAX_DEPTH.  Return the tree holding the
+   result of the expanded sequence or NULL_TREE if the expansion failed.
+
+   This routine assumes that ARG1 is a real number with a fractional part
+   (the integer exponent case will have been handled earlier in
+   gimple_expand_builtin_pow).
+
+   For ARG1 > 0.0:
+   * For ARG1 composed of a whole part WHOLE_PART and a fractional part
+     FRAC_PART i.e. WHOLE_PART == floor (ARG1) and
+                    FRAC_PART == ARG1 - WHOLE_PART:
+     Produce POWI (ARG0, WHOLE_PART) * POW (ARG0, FRAC_PART) where
+     POW (ARG0, FRAC_PART) is expanded as a product of square root chains
+     if it can be expressed as such, that is if FRAC_PART satisfies:
+     FRAC_PART == <SUM from i = 1 until MAX_DEPTH> (a[i] * (0.5**i))
+     where integer a[i] is either 0 or 1.
+
+     Example:
+     POW (x, 3.625) == POWI (x, 3) * POW (x, 0.625)
+       --> POWI (x, 3) * SQRT (x) * SQRT (SQRT (SQRT (x)))
+
+   For ARG1 < 0.0 there are two approaches:
+   * (A) Expand to 1.0 / POW (ARG0, -ARG1) where POW (ARG0, -ARG1)
+         is calculated as above.
+
+     Example:
+     POW (x, -5.625) == 1.0 / POW (x, 5.625)
+       -->  1.0 / (POWI (x, 5) * SQRT (x) * SQRT (SQRT (SQRT (x))))
+
+   * (B) : WHOLE_PART := - ceil (abs (ARG1))
+           FRAC_PART  := ARG1 - WHOLE_PART
+     and expand to POW (x, FRAC_PART) / POWI (x, WHOLE_PART).
+     Example:
+     POW (x, -5.875) == POW (x, 0.125) / POWI (X, 6)
+       --> SQRT (SQRT (SQRT (x))) / (POWI (x, 6))
+
+   For ARG1 < 0.0 we choose between (A) and (B) depending on
+   how many multiplications we'd have to do.
+   So, for the example in (B): POW (x, -5.875), if we were to
+   follow algorithm (A) we would produce:
+   1.0 / POWI (X, 5) * SQRT (X) * SQRT (SQRT (X)) * SQRT (SQRT (SQRT (X)))
+   which contains more multiplications than approach (B).
+
+   Hopefully, this approach will eliminate potentially expensive POW library
+   calls when unsafe floating point math is enabled and allow the compiler to
+   further optimise the multiplies, square roots and divides produced by this
+   function.  */
+
+static tree
+expand_pow_as_sqrts (gimple_stmt_iterator *gsi, location_t loc,
+		     tree arg0, tree arg1, HOST_WIDE_INT max_depth)
+{
+  tree type = TREE_TYPE (arg0);
+  machine_mode mode = TYPE_MODE (type);
+  tree sqrtfn = mathfn_built_in (type, BUILT_IN_SQRT);
+  bool one_over = true;
+
+  if (!sqrtfn)
+    return NULL_TREE;
+
+  if (TREE_CODE (arg1) != REAL_CST)
+    return NULL_TREE;
+
+  REAL_VALUE_TYPE exp_init = TREE_REAL_CST (arg1);
+
+  gcc_assert (max_depth > 0);
+  tree *cache = XALLOCAVEC (tree, max_depth + 1);
+
+  struct pow_synth_sqrt_info synth_info;
+  synth_info.factors = XALLOCAVEC (bool, max_depth + 1);
+  synth_info.deepest = 0;
+  synth_info.num_mults = 0;
+
+  bool neg_exp = REAL_VALUE_NEGATIVE (exp_init);
+  REAL_VALUE_TYPE exp = real_value_abs (&exp_init);
+
+  /* The whole and fractional parts of exp.  */
+  REAL_VALUE_TYPE whole_part;
+  REAL_VALUE_TYPE frac_part;
+
+  real_floor (&whole_part, mode, &exp);
+  REAL_ARITHMETIC (frac_part, MINUS_EXPR, exp, whole_part);
+
+
+  REAL_VALUE_TYPE ceil_whole = dconst0;
+  REAL_VALUE_TYPE ceil_fract = dconst0;
+
+  if (neg_exp)
+    {
+      real_ceil (&ceil_whole, mode, &exp);
+      REAL_ARITHMETIC (ceil_fract, MINUS_EXPR, ceil_whole, exp);
+    }
+
+  if (!representable_as_half_series_p (frac_part, max_depth, &synth_info))
+    return NULL_TREE;
+
+  /* Check whether it's more profitable to not use 1.0 / ...  */
+  if (neg_exp)
+    {
+      struct pow_synth_sqrt_info alt_synth_info;
+      alt_synth_info.factors = XALLOCAVEC (bool, max_depth + 1);
+      alt_synth_info.deepest = 0;
+      alt_synth_info.num_mults = 0;
+
+      if (representable_as_half_series_p (ceil_fract, max_depth,
+					   &alt_synth_info)
+	  && alt_synth_info.deepest <= synth_info.deepest
+	  && alt_synth_info.num_mults < synth_info.num_mults)
+	{
+	  whole_part = ceil_whole;
+	  frac_part = ceil_fract;
+	  synth_info.deepest = alt_synth_info.deepest;
+	  synth_info.num_mults = alt_synth_info.num_mults;
+	  memcpy (synth_info.factors, alt_synth_info.factors,
+		  (max_depth + 1) * sizeof (bool));
+	  one_over = false;
+	}
+    }
+
+  HOST_WIDE_INT n = real_to_integer (&whole_part);
+  REAL_VALUE_TYPE cint;
+  real_from_integer (&cint, VOIDmode, n, SIGNED);
+
+  if (!real_identical (&whole_part, &cint))
+    return NULL_TREE;
+
+  if (powi_cost (n) + synth_info.num_mults > POWI_MAX_MULTS)
+    return NULL_TREE;
+
+  memset (cache, 0, (max_depth + 1) * sizeof (tree));
+
+  tree integer_res = n == 0 ? build_real (type, dconst1) : arg0;
+
+  /* Calculate the integer part of the exponent.  */
+  if (n > 1)
+    {
+      integer_res = gimple_expand_builtin_powi (gsi, loc, arg0, n);
+      if (!integer_res)
+	return NULL_TREE;
+    }
+
+  if (dump_file)
+    {
+      char string[64];
+
+      real_to_decimal (string, &exp_init, sizeof (string), 0, 1);
+      fprintf (dump_file, "synthesizing pow (x, %s) as:\n", string);
+
+      if (neg_exp)
+	{
+	  if (one_over)
+	    {
+	      fprintf (dump_file, "1.0 / (");
+	      dump_integer_part (dump_file, "x", n);
+	      if (n > 0)
+	        fprintf (dump_file, " * ");
+	      dump_fractional_sqrt_sequence (dump_file, "x", &synth_info);
+	      fprintf (dump_file, ")");
+	    }
+	  else
+	    {
+	      dump_fractional_sqrt_sequence (dump_file, "x", &synth_info);
+	      fprintf (dump_file, " / (");
+	      dump_integer_part (dump_file, "x", n);
+	      fprintf (dump_file, ")");
+	    }
+	}
+      else
+	{
+	  dump_fractional_sqrt_sequence (dump_file, "x", &synth_info);
+	  if (n > 0)
+	    fprintf (dump_file, " * ");
+	  dump_integer_part (dump_file, "x", n);
+	}
+
+      fprintf (dump_file, "\ndeepest sqrt chain: %d\n", synth_info.deepest);
+    }
+
+
+  tree fract_res = NULL_TREE;
+  cache[0] = arg0;
+
+  /* Calculate the fractional part of the exponent.  */
+  for (unsigned i = 0; i < synth_info.deepest; i++)
+    {
+      if (synth_info.factors[i])
+	{
+	  tree sqrt_chain = get_fn_chain (arg0, i + 1, gsi, sqrtfn, loc, cache);
+
+	  if (!fract_res)
+	      fract_res = sqrt_chain;
+
+	  else
+	    fract_res = build_and_insert_binop (gsi, loc, "powroot", MULT_EXPR,
+					   fract_res, sqrt_chain);
+	}
+    }
+
+  tree res = NULL_TREE;
+
+  if (neg_exp)
+    {
+      if (one_over)
+	{
+	  if (n > 0)
+	    res = build_and_insert_binop (gsi, loc, "powroot", MULT_EXPR,
+					   fract_res, integer_res);
+	  else
+	    res = fract_res;
+
+	  res = build_and_insert_binop (gsi, loc, "powrootrecip", RDIV_EXPR,
+					  build_real (type, dconst1), res);
+	}
+      else
+	{
+	  res = build_and_insert_binop (gsi, loc, "powroot", RDIV_EXPR,
+					 fract_res, integer_res);
+	}
+    }
+  else
+    res = build_and_insert_binop (gsi, loc, "powroot", MULT_EXPR,
+				   fract_res, integer_res);
+  return res;
+}
+
 /* ARG0 and ARG1 are the two arguments to a pow builtin call in GSI
    with location info LOC.  If possible, create an equivalent and
    less expensive sequence of statements prior to GSI, and return an
@@ -1157,13 +1509,17 @@  static tree
 gimple_expand_builtin_pow (gimple_stmt_iterator *gsi, location_t loc, 
 			   tree arg0, tree arg1)
 {
-  REAL_VALUE_TYPE c, cint, dconst1_4, dconst3_4, dconst1_3, dconst1_6;
+  REAL_VALUE_TYPE c, cint, dconst1_3, dconst1_4, dconst1_6;
   REAL_VALUE_TYPE c2, dconst3;
   HOST_WIDE_INT n;
-  tree type, sqrtfn, cbrtfn, sqrt_arg0, sqrt_sqrt, result, cbrt_x, powi_cbrt_x;
+  tree type, sqrtfn, cbrtfn, sqrt_arg0, result, cbrt_x, powi_cbrt_x;
   machine_mode mode;
+  bool speed_p = optimize_bb_for_speed_p (gsi_bb (*gsi));
   bool hw_sqrt_exists, c_is_int, c2_is_int;
 
+  dconst1_4 = dconst1;
+  SET_REAL_EXP (&dconst1_4, REAL_EXP (&dconst1_4) - 2);
+
   /* If the exponent isn't a constant, there's nothing of interest
      to be done.  */
   if (TREE_CODE (arg1) != REAL_CST)
@@ -1179,7 +1535,7 @@  gimple_expand_builtin_pow (gimple_stmt_iterator *gsi, location_t loc,
   if (c_is_int
       && ((n >= -1 && n <= 2)
 	  || (flag_unsafe_math_optimizations
-	      && optimize_bb_for_speed_p (gsi_bb (*gsi))
+	      && speed_p
 	      && powi_cost (n) <= POWI_MAX_MULTS)))
     return gimple_expand_builtin_powi (gsi, loc, arg0, n);
 
@@ -1196,49 +1552,8 @@  gimple_expand_builtin_pow (gimple_stmt_iterator *gsi, location_t loc,
       && !HONOR_SIGNED_ZEROS (mode))
     return build_and_insert_call (gsi, loc, sqrtfn, arg0);
 
-  /* Optimize pow(x,0.25) = sqrt(sqrt(x)).  Assume on most machines that
-     a builtin sqrt instruction is smaller than a call to pow with 0.25,
-     so do this optimization even if -Os.  Don't do this optimization
-     if we don't have a hardware sqrt insn.  */
-  dconst1_4 = dconst1;
-  SET_REAL_EXP (&dconst1_4, REAL_EXP (&dconst1_4) - 2);
   hw_sqrt_exists = optab_handler (sqrt_optab, mode) != CODE_FOR_nothing;
 
-  if (flag_unsafe_math_optimizations
-      && sqrtfn
-      && REAL_VALUES_EQUAL (c, dconst1_4)
-      && hw_sqrt_exists)
-    {
-      /* sqrt(x)  */
-      sqrt_arg0 = build_and_insert_call (gsi, loc, sqrtfn, arg0);
-
-      /* sqrt(sqrt(x))  */
-      return build_and_insert_call (gsi, loc, sqrtfn, sqrt_arg0);
-    }
-      
-  /* Optimize pow(x,0.75) = sqrt(x) * sqrt(sqrt(x)) unless we are
-     optimizing for space.  Don't do this optimization if we don't have
-     a hardware sqrt insn.  */
-  real_from_integer (&dconst3_4, VOIDmode, 3, SIGNED);
-  SET_REAL_EXP (&dconst3_4, REAL_EXP (&dconst3_4) - 2);
-
-  if (flag_unsafe_math_optimizations
-      && sqrtfn
-      && optimize_function_for_speed_p (cfun)
-      && REAL_VALUES_EQUAL (c, dconst3_4)
-      && hw_sqrt_exists)
-    {
-      /* sqrt(x)  */
-      sqrt_arg0 = build_and_insert_call (gsi, loc, sqrtfn, arg0);
-
-      /* sqrt(sqrt(x))  */
-      sqrt_sqrt = build_and_insert_call (gsi, loc, sqrtfn, sqrt_arg0);
-
-      /* sqrt(x) * sqrt(sqrt(x))  */
-      return build_and_insert_binop (gsi, loc, "powroot", MULT_EXPR,
-				     sqrt_arg0, sqrt_sqrt);
-    }
-
   /* Optimize pow(x,1./3.) = cbrt(x).  This requires unsafe math
      optimizations since 1./3. is not exactly representable.  If x
      is negative and finite, the correct value of pow(x,1./3.) is
@@ -1263,7 +1578,7 @@  gimple_expand_builtin_pow (gimple_stmt_iterator *gsi, location_t loc,
       && sqrtfn
       && cbrtfn
       && (gimple_val_nonnegative_real_p (arg0) || !HONOR_NANS (mode))
-      && optimize_function_for_speed_p (cfun)
+      && speed_p
       && hw_sqrt_exists
       && REAL_VALUES_EQUAL (c, dconst1_6))
     {
@@ -1274,54 +1589,31 @@  gimple_expand_builtin_pow (gimple_stmt_iterator *gsi, location_t loc,
       return build_and_insert_call (gsi, loc, cbrtfn, sqrt_arg0);
     }
 
-  /* Optimize pow(x,c), where n = 2c for some nonzero integer n
-     and c not an integer, into
-
-       sqrt(x) * powi(x, n/2),                n > 0;
-       1.0 / (sqrt(x) * powi(x, abs(n/2))),   n < 0.
-
-     Do not calculate the powi factor when n/2 = 0.  */
-  real_arithmetic (&c2, MULT_EXPR, &c, &dconst2);
-  n = real_to_integer (&c2);
-  real_from_integer (&cint, VOIDmode, n, SIGNED);
-  c2_is_int = real_identical (&c2, &cint);
 
+  /* Attempt to expand the POW as a product of square root chains.
+     Expand the 0.25 case even when otpimising for size.  */
   if (flag_unsafe_math_optimizations
       && sqrtfn
-      && c2_is_int
-      && !c_is_int
-      && optimize_function_for_speed_p (cfun))
+      && hw_sqrt_exists
+      && (speed_p || REAL_VALUES_EQUAL (c, dconst1_4))
+      && !HONOR_SIGNED_ZEROS (mode))
     {
-      tree powi_x_ndiv2 = NULL_TREE;
-
-      /* Attempt to fold powi(arg0, abs(n/2)) into multiplies.  If not
-         possible or profitable, give up.  Skip the degenerate case when
-         n is 1 or -1, where the result is always 1.  */
-      if (absu_hwi (n) != 1)
-	{
-	  powi_x_ndiv2 = gimple_expand_builtin_powi (gsi, loc, arg0,
-						     abs_hwi (n / 2));
-	  if (!powi_x_ndiv2)
-	    return NULL_TREE;
-	}
+      unsigned int max_depth = speed_p
+				? PARAM_VALUE (PARAM_MAX_POW_SQRT_DEPTH)
+				: 2;
 
-      /* Calculate sqrt(x).  When n is not 1 or -1, multiply it by the
-	 result of the optimal multiply sequence just calculated.  */
-      sqrt_arg0 = build_and_insert_call (gsi, loc, sqrtfn, arg0);
+      tree expand_with_sqrts
+	= expand_pow_as_sqrts (gsi, loc, arg0, arg1, max_depth);
 
-      if (absu_hwi (n) == 1)
-	result = sqrt_arg0;
-      else
-	result = build_and_insert_binop (gsi, loc, "powroot", MULT_EXPR,
-					 sqrt_arg0, powi_x_ndiv2);
-
-      /* If n is negative, reciprocate the result.  */
-      if (n < 0)
-	result = build_and_insert_binop (gsi, loc, "powroot", RDIV_EXPR,
-					 build_real (type, dconst1), result);
-      return result;
+      if (expand_with_sqrts)
+	return expand_with_sqrts;
     }
 
+  real_arithmetic (&c2, MULT_EXPR, &c, &dconst2);
+  n = real_to_integer (&c2);
+  real_from_integer (&cint, VOIDmode, n, SIGNED);
+  c2_is_int = real_identical (&c2, &cint);
+
   /* Optimize pow(x,c), where 3c = n for some nonzero integer n, into
 
      powi(x, n/3) * powi(cbrt(x), n%3),                    n > 0;