Re: [PATCH] Optimise sqrt reciprocal multiplications

[Kyrill Tkachov <kyrylo dot tkachov at foss dot arm dot com>]

On 04/09/18 17:52, Kyrill Tkachov wrote:
>
> On 04/09/18 15:31, Richard Biener wrote:
>> On Tue, 4 Sep 2018, Kyrill Tkachov wrote:
>>> Hi Richard,
>>>
>>> On 31/08/18 12:07, Richard Biener wrote:
>>>> On Thu, 30 Aug 2018, Kyrill Tkachov wrote:
>>>>> Ping.
>>>>>
>>>>> https://gcc.gnu.org/ml/gcc-patches/2018-08/msg01496.html
>>>>>
>>>>> Thanks,
>>>>> Kyrill
>>>>>
>>>>> On 23/08/18 18:09, Kyrill Tkachov wrote:
>>>>>> Hi Richard,
>>>>>>
>>>>>> On 23/08/18 11:13, Richard Sandiford wrote:
>>>>>>> Kyrill  Tkachov <kyrylo.tkachov@foss.arm.com> writes:
>>>>>>>> Hi all,
>>>>>>>>
>>>>>>>> This patch aims to optimise sequences involving uses of 1.0 / sqrt (a)
>>>>>>>> under -freciprocal-math and -funsafe-math-optimizations.
>>>>>>>> In particular consider:
>>>>>>>>
>>>>>>>> x = 1.0 / sqrt (a);
>>>>>>>> r1 = x * x;  // same as 1.0 / a
>>>>>>>> r2 = a * x; // same as sqrt (a)
>>>>>>>>
>>>>>>>> If x, r1 and r2 are all used further on in the code, this can be
>>>>>>>> transformed into:
>>>>>>>> tmp1 = 1.0 / a
>>>>>>>> tmp2 = sqrt (a)
>>>>>>>> tmp3 = tmp1 * tmp2
>>>>>>>> x = tmp3
>>>>>>>> r1 = tmp1
>>>>>>>> r2 = tmp2
>>>>>>> Nice optimisation :-)  Someone who knows the pass better should review,
>>>>>>> but:
>>>>>> Thanks for the review.
>>>>>>> There seems to be an implicit assumption that this is a win even
>>>>>>> when the r1 and r2 assignments are only conditionally executed.
>>>>>>> That's probably true, but it might be worth saying explicitly.
>>>>>> I'll admit I had not considered that case.
>>>>>> I think it won't make a difference in practice, as the really expensive
>>>>>> operations here
>>>>>> are the sqrt and the division and they are on the executed path in either
>>>>>> case and them
>>>>>> becoming independent should be a benefit of its own.
>>>>>>>> +/* Return TRUE if USE_STMT is a multiplication of DEF by A.  */
>>>>>>>> +
>>>>>>>> +static inline bool
>>>>>>>> +is_mult_by (gimple *use_stmt, tree def, tree a)
>>>>>>>> +{
>>>>>>>> +  if (gimple_code (use_stmt) == GIMPLE_ASSIGN
>>>>>>>> +      && gimple_assign_rhs_code (use_stmt) == MULT_EXPR)
>>>>>>>> +    {
>>>>>>>> +      tree op0 = gimple_assign_rhs1 (use_stmt);
>>>>>>>> +      tree op1 = gimple_assign_rhs2 (use_stmt);
>>>>>>>> +
>>>>>>>> +      return (op0 == def && op1 == a)
>>>>>>>> +          || (op0 == a && op1 == def);
>>>>>>>> +    }
>>>>>>>> +  return 0;
>>>>>>>> +}
>>>>>>> Seems like is_square_of could now be a light-weight wrapper around this.
>>>>>> Indeed, I've done the wrapping now.
>>>>>>>> @@ -652,6 +669,180 @@ execute_cse_reciprocals_1 (gimple_stmt_iterator
>>>>>>>> *def_gsi, tree def)
>>>>>>>>      occ_head = NULL;
>>>>>>>>    }
>>>>>>>>    +/* Transform sequences like
>>>>>>>> +   x = 1.0 / sqrt (a);
>>>>>>>> +   r1 = x * x;
>>>>>>>> +   r2 = a * x;
>>>>>>>> +   into:
>>>>>>>> +   tmp1 = 1.0 / a;
>>>>>>>> +   tmp2 = sqrt (a);
>>>>>>>> +   tmp3 = tmp1 * tmp2;
>>>>>>>> +   x = tmp3;
>>>>>>>> +   r1 = tmp1;
>>>>>>>> +   r2 = tmp2;
>>>>>>>> +   depending on the uses of x, r1, r2.  This removes one multiplication
>>>>>>>> and
>>>>>>>> +   allows the sqrt and division operations to execute in parallel.
>>>>>>>> +   DEF_GSI is the gsi of the initial division by sqrt that defines
>>>>>>>> +   DEF (x in the example abovs). */
>>>>>>>> +
>>>>>>>> +static void
>>>>>>>> +optimize_recip_sqrt (gimple_stmt_iterator *def_gsi, tree def)
>>>>>>>> +{
>>>>>>>> +  use_operand_p use_p;
>>>>>>>> +  imm_use_iterator use_iter;
>>>>>>>> +  gimple *stmt = gsi_stmt (*def_gsi);
>>>>>>>> +  tree x = def;
>>>>>>>> +  tree orig_sqrt_ssa_name = gimple_assign_rhs2 (stmt);
>>>>>>>> +  tree div_rhs1 = gimple_assign_rhs1 (stmt);
>>>>>>>> +
>>>>>>>> +  if (TREE_CODE (orig_sqrt_ssa_name) != SSA_NAME
>>>>>>>> +      || TREE_CODE (div_rhs1) != REAL_CST
>>>>>>>> +      || !real_equal (&TREE_REAL_CST (div_rhs1), &dconst1))
>>>>>>>> +    return;
>>>>>>>> +
>>>>>>>> +  gimple *sqrt_stmt = SSA_NAME_DEF_STMT (orig_sqrt_ssa_name);
>>>>>>>> +  if (!is_gimple_call (sqrt_stmt)
>>>>>>>> +      || !gimple_call_lhs (sqrt_stmt))
>>>>>>>> +    return;
>>>>>>>> +
>>>>>>>> +  gcall *call = as_a <gcall *> (sqrt_stmt);
>>>>>>> Very minor, but:
>>>>>>>
>>>>>>>     gcall *sqrt_stmt
>>>>>>>       = dyn_cast <gcall *> (SSA_NAME_DEF_STMT (orig_sqrt_ssa_name));
>>>>>>>     if (!sqrt_stmt || !gimple_call_lhs (sqrt_stmt))
>>>>>>>       return;
>>>>>>>
>>>>>>> would avoid the need for the separate as_a<>, and would mean that
>>>>>>> we only call gimple_call_* on gcalls.
>>>>>> Ok.
>>>>>>>> +  if (has_other_use)
>>>>>>>> +    {
>>>>>>>> +      /* Using the two temporaries tmp1, tmp2 from above
>>>>>>>> +     the original x is now:
>>>>>>>> +     x = tmp1 * tmp2.  */
>>>>>>>> +      gcc_assert (mult_ssa_name);
>>>>>>>> +      gcc_assert (sqr_ssa_name);
>>>>>>>> +      gimple_stmt_iterator gsi2 = gsi_for_stmt (stmt);
>>>>>>>> +
>>>>>>>> +      tree new_ssa_name
>>>>>>>> +    = make_temp_ssa_name (TREE_TYPE (a), NULL,
>>>>>>>> "recip_sqrt_transformed");
>>>>>>>> +      gimple *new_stmt
>>>>>>>> +    = gimple_build_assign (new_ssa_name, MULT_EXPR,
>>>>>>>> +                   mult_ssa_name, sqr_ssa_name);
>>>>>>>> +      gsi_insert_before (&gsi2, new_stmt, GSI_SAME_STMT);
>>>>>>>> +      gcc_assert (gsi_stmt (gsi2) == stmt);
>>>>>>>> + gimple_assign_set_rhs_from_tree (&gsi2, new_ssa_name);
>>>>>>>> +      fold_stmt (&gsi2);
>>>>>>>> +      update_stmt (stmt);
>>>>>>> In this case we're replacing the statement in its original position,
>>>>>>> so there's no real need to use a temporary.  It seems better to
>>>>>>> change the rhs_code, rhs1 and rhs2 of stmt in-place, with the same
>>>>>>> lhs as before.
>>>>>> Yes, that's cleaner.
>>>>>>>> @@ -762,6 +953,23 @@ pass_cse_reciprocals::execute (function *fun)
>>>>>>>>          if (optimize_bb_for_size_p (bb))
>>>>>>>>            continue;
>>>>>>> Seems unnecessary to skip the new optimisation when optimising for size.
>>>>>>> Like you say, it saves a multiplication overall. Also:
>>>>>> Indeed.
>>>>>>>> +      if (flag_unsafe_math_optimizations)
>>>>>>>> +    {
>>>>>>>> +      for (gimple_stmt_iterator gsi = gsi_after_labels (bb);
>>>>>>>> +           !gsi_end_p (gsi);
>>>>>>>> +           gsi_next (&gsi))
>>>>>>>> +        {
>>>>>>>> +          gimple *stmt = gsi_stmt (gsi);
>>>>>>>> +
>>>>>>>> +          if (gimple_has_lhs (stmt)
>>>>>>>> +          && (def = SINGLE_SSA_TREE_OPERAND (stmt, SSA_OP_DEF)) != NULL
>>>>>>>> +          && FLOAT_TYPE_P (TREE_TYPE (def))
>>>>>>>> +          && TREE_CODE (def) == SSA_NAME
>>>>>>>> +          && is_gimple_assign (stmt)
>>>>>>>> +          && gimple_assign_rhs_code (stmt) == RDIV_EXPR)
>>>>>>>> +        optimize_recip_sqrt (&gsi, def);
>>>>>>>> +        }
>>>>>>>> +    }
>>>>>>> It looks like this could safely be done in one of the existing walks
>>>>>>> (e.g. the execute_cse_reciprocals_1 one, if we do this when optimising
>>>>>>> for size).
>>>>>> You're right. I've moved this into one of the walks above this.
>>>>>>
>>>>>> Bootstrapped and tested on aarch64-none-linux-gnu and
>>>>>> x86_64-unknown-linux-gnu.
>>>>>> CC'ing richi as he's reviewed these kinds of patches in the past.
>>>>>>
>>>>>> Is this ok for trunk?
>>>> I wonder how it interacts with execute_cse_reciprocals_1 given it
>>>> introduces a division 1.0 / a which will be not processed by
>>>> execute_cse_reciprocals_1 given that operates by walking uses of 'a'.
>>>> That may be just a missed optimization of course.
>>> Hmm, I believe right now it doesn't interact with execute_cse_reciprocals_1 as
>>> it's either one or the other.
>>> I've left it as it is for now, but would you like us to call
>>> execute_cse_reciprocals_1 on the potentially transformed division?
>> I think that wouldn't "fix" it given execute_cse_reciprocals_1 works by
>> seeing all reciprocal uses of 'a' so it may have alrady seen two and
>> you introduced the third.  But we can leave "solving" this for
>> future enhacements (I'd avoid doing two passes over the IL just because
>> of said possibility right now).
>>>> +  if (has_other_use)
>>>> +    {
>>>> +      /* Using the two temporaries tmp1, tmp2 from above
>>>> +        the original x is now:
>>>> +        x = tmp1 * tmp2.  */
>>>> +      gcc_assert (mult_ssa_name);
>>>> +      gcc_assert (sqr_ssa_name);
>>>> +
>>>> +      gimple_assign_set_rhs1 (stmt, mult_ssa_name);
>>>> +      gimple_assign_set_rhs2 (stmt, sqr_ssa_name);
>>>> +      gimple_assign_set_rhs_code (stmt, MULT_EXPR);
>>>> +      fold_stmt_inplace (def_gsi);
>>>> +      update_stmt (stmt);
>>>> +    }
>>>>
>>>> so you are leaving the original stmt in place unchanged even if it is
>>>> not used?  Why?  Note that with -fno-call-exceptions this stmt may
>>>> throw, so you should arrange to code-generate 1./a in place of the
>>>> original division to preserve EH behavior.  Watch out because then
>>>> with other uses you have to find a place to insert its computation.
>>> Ok. These are oversights on my part. I've updated the patch to modify the
>>> original division in place. The multiplication is placed after it.
>> If the division throws internally you can't insert after it, you
>> have to insert on the single non-EH edge outgoing from the BB
>> with the division.  Just try a C++ testcase with -fnon-call-exceptions
>> and a try {} catch(...) {} around.
>>
>> Not sure if it's worth to handle so bailing out if
>> stmt_can_throw_internal (stmt) would be an option as well.
>>>> +      if (is_square_of (stmt2, x))
>>>> +       {
>>>> +         if (!sqr_stmts.contains (stmt2))
>>>> +           sqr_stmts.safe_push (stmt2);
>>>> +       }
>>>>
>>>> this is quadratic in the number of square stmts... please consider
>>>> making sqr_stmts a bitmap of SSA defs (so the stmt you have now
>>>> is then SSA_NAME_DEF_STMT (ssa_name (bitmap-element))).
>>> Done. In practice I didn't see there being more than one such use, I expect
>>> them
>>> to be CSE'd, but maybe if they're in different basic blocks...
>> Oh, it just occured to me you could use FOR_EACH_IMM_USE_STMT ()
>> which will only get you each stmt once... ;)  Sorry for the misleading
>> bitmap suggestion.
>>>> You do not seem to restrict placement of the use stmts but insert
>>>> before the 1/sqrt(a) stmt.  That possibly hoists the multiplications
>>>> where they are not needed.  Consider
>>>>
>>>>    x = 1./sqrt(a);
>>>>    if (l)
>>>>      l1 = x * 3.;
>>>>    else if (l2)
>>>>      l1 = x * x;
>>>>    else if (l3)
>>>>      l1 = a * x;
>>>>
>>>> or similar where on the path to x * 3. you now perform two extra
>>>> multiplications.
>>> Ok, I've restricted the optimisation somewhat.
>>> Now if there's an other use it won't perform the transformation unless there
>>> is already
>>> a multiplication present on the main path. That way it won't introduce
>>> multiplications
>>> on paths where there aren't any.
>> OK.
>>>> Your testcases do not cover the case of other uses at all.  Or of
>>>> EH.
>>> gcc.dg/recip_sqrt_mult_1.c should handle the other uses case as tmp is a
>>> global
>>> being written to. I've added more testcases involving uses in different basic
>>> blocks
>>> and a g++.dg testcase with -fnon-call-exceptions.
>>> I'm not sure what functionality to test though apart from that it doesn't ICE
>>> and does
>>> the transformations I expect.
>> That's good enough I guess.  I see you do have a testcase with
>> try/catch so I wonder why foo4 doesn't ICE when you insert after
>> the throwing stmt...  maybe it doesn't throw?  Ah - they probably
>> fail your new mult_on_main_path test?
>
> Yeah. I did manage to reproduce an ICE eventually by removing that check
> and enabling -ftrapping-math.
>
> I think I'll just not enter this at all if stmt_can_throw_internal (stmt) like you suggested.


And here is the version that does that. It also reverts to using a vector for the sqr_stmts for consistency
and uses FOR_EACH_IMM_USE_STMT to iterate over the use statements.

Bootstrapped and tested on aarch64-none-linux-gnu.

Ok?

Thanks,
Kyrill

2018-09-03  Kyrylo Tkachov  <kyrylo.tkachov@arm.com>

    * tree-ssa-math-opts.c (is_mult_by): New function.
    (is_square_of): Use the above.
    (optimize_recip_sqrt): New function.
    (pass_cse_reciprocals::execute): Use the above.

2018-09-03  Kyrylo Tkachov  <kyrylo.tkachov@arm.com>

    * gcc.dg/recip_sqrt_mult_1.c: New test.
    * gcc.dg/recip_sqrt_mult_2.c: Likewise.
    * gcc.dg/recip_sqrt_mult_3.c: Likewise.
    * gcc.dg/recip_sqrt_mult_4.c: Likewise.
    * gcc.dg/recip_sqrt_mult_5.c: Likewise.
    * g++.dg/recip_sqrt_mult_1.C: Likewise.
    * g++.dg/recip_sqrt_mult_2.C: Likewise.

diff --git a/gcc/testsuite/g++.dg/recip_sqrt_mult_1.C b/gcc/testsuite/g++.dg/recip_sqrt_mult_1.C
new file mode 100644
index 0000000000000000000000000000000000000000..11d9c6f758f1529d8ed4cadf85010f6ce379c195
--- /dev/null
+++ b/gcc/testsuite/g++.dg/recip_sqrt_mult_1.C
@@ -0,0 +1,49 @@
+/* { dg-do compile } */
+/* { dg-options "-Ofast -fnon-call-exceptions -fdump-tree-recip" } */
+
+double res, res2, tmp;
+void
+foo1 (double a, double b)
+{
+  try {
+    tmp = 1.0 / __builtin_sqrt (a);
+    res = tmp * tmp;
+    res2 = a * tmp;
+  }
+  catch (...)
+    { ; }
+}
+
+void
+foo4 (double a, double b, int c, int d)
+{
+  try {
+    tmp = 1.0 / __builtin_sqrt (a);
+  }
+  catch (...)
+    {
+      if (c)
+	res = tmp * tmp;
+
+      if (d)
+	res2 = a * tmp;
+    }
+}
+
+void
+foo5 (double a, double b, int c, int d)
+{
+  try {
+    tmp = 1.0 / __builtin_sqrt (a);
+    res = tmp * tmp;
+
+    if (d)
+      res2 = a * tmp;
+  }
+  catch (...)
+    { ; }
+}
+
+/* { dg-final { scan-tree-dump-times "Optimizing reciprocal sqrt multiplications" 2 "recip" } } */
+/* { dg-final { scan-tree-dump-times "Replacing squaring multiplication" 2 "recip" } } */
+/* { dg-final { scan-tree-dump-times "Replacing original division" 2 "recip" } } */
diff --git a/gcc/testsuite/g++.dg/recip_sqrt_mult_2.C b/gcc/testsuite/g++.dg/recip_sqrt_mult_2.C
new file mode 100644
index 0000000000000000000000000000000000000000..cca12caf4d79bfa69933d9b8fce41f38bb5b7a19
--- /dev/null
+++ b/gcc/testsuite/g++.dg/recip_sqrt_mult_2.C
@@ -0,0 +1,49 @@
+/* { dg-do compile } */
+/* { dg-options "-w -Ofast -fnon-call-exceptions -ftrapping-math -fdump-tree-recip" } */
+
+/* Check that the recip_sqrt optimization does not trigger here, causing an
+   ICE due to EH info.  */
+
+
+double res, res2, tmp;
+void
+foo1 (double a, double b)
+{
+  try {
+    tmp = 1.0 / __builtin_sqrt (a);
+    res = tmp * tmp;
+    res2 = a * tmp;
+  }
+  catch (...)
+    { ; }
+}
+
+void
+foo4 (double a, double b, int c, int d)
+{
+  try {
+    tmp = 1.0 / __builtin_sqrt (a);
+  }
+  catch (...)
+    {
+      if (c)
+	res = tmp * tmp;
+
+      if (d)
+	res2 = a * tmp;
+    }
+}
+
+void
+foo5 (double a, double b, int c, int d)
+{
+  try {
+    tmp = 1.0 / __builtin_sqrt (a);
+    res = tmp * tmp;
+
+    if (d)
+      res2 = a * tmp;
+  }
+  catch (...)
+    { ; }
+}
diff --git a/gcc/testsuite/gcc.dg/recip_sqrt_mult_1.c b/gcc/testsuite/gcc.dg/recip_sqrt_mult_1.c
new file mode 100644
index 0000000000000000000000000000000000000000..188390a4ecffca1b21f05c4f19abecfb7ebd188f
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/recip_sqrt_mult_1.c
@@ -0,0 +1,15 @@
+/* { dg-do compile } */
+/* { dg-options "-Ofast -fdump-tree-recip" } */
+
+double res, res2, tmp;
+void
+foo (double a, double b)
+{
+  tmp = 1.0 / __builtin_sqrt (a);
+  res = tmp * tmp;
+  res2 = a * tmp;
+}
+
+/* { dg-final { scan-tree-dump "Optimizing reciprocal sqrt multiplications" "recip" } } */
+/* { dg-final { scan-tree-dump "Replacing squaring multiplication" "recip" } } */
+/* { dg-final { scan-tree-dump "Replacing original division" "recip" } } */
diff --git a/gcc/testsuite/gcc.dg/recip_sqrt_mult_2.c b/gcc/testsuite/gcc.dg/recip_sqrt_mult_2.c
new file mode 100644
index 0000000000000000000000000000000000000000..c5fc3de7b657b1769e76254b4bc874e0595e43ef
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/recip_sqrt_mult_2.c
@@ -0,0 +1,11 @@
+/* { dg-do compile } */
+/* { dg-options "-Ofast -fdump-tree-optimized" } */
+
+float
+foo (float a)
+{
+  float tmp = 1.0f / __builtin_sqrtf (a);
+  return a * tmp;
+}
+
+/* { dg-final { scan-tree-dump-not " / " "optimized" } } */
diff --git a/gcc/testsuite/gcc.dg/recip_sqrt_mult_3.c b/gcc/testsuite/gcc.dg/recip_sqrt_mult_3.c
new file mode 100644
index 0000000000000000000000000000000000000000..e7d185ba7e22cfc7cca72296d5ccc544f24fdb14
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/recip_sqrt_mult_3.c
@@ -0,0 +1,11 @@
+/* { dg-do compile } */
+/* { dg-options "-Ofast -fdump-tree-optimized" } */
+
+double
+foo (double a)
+{
+  double tmp = 1.0f / __builtin_sqrt (a);
+  return tmp * tmp;
+}
+
+/* { dg-final { scan-tree-dump-not "__builtin_sqrt" "optimized" } } */
diff --git a/gcc/testsuite/gcc.dg/recip_sqrt_mult_4.c b/gcc/testsuite/gcc.dg/recip_sqrt_mult_4.c
new file mode 100644
index 0000000000000000000000000000000000000000..e3005f2feb6f4bacbb6eafc0155e196cb866fcdf
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/recip_sqrt_mult_4.c
@@ -0,0 +1,21 @@
+/* { dg-do compile } */
+/* { dg-options "-Ofast -fdump-tree-recip" } */
+
+/* The main path doesn't have any multiplications.
+   Avoid introducing them in the recip pass.  */
+
+double res, res2, tmp;
+void
+foo (double a, double b, int c, int d)
+{
+  tmp = 1.0 / __builtin_sqrt (a);
+  if (c)
+    res = tmp * tmp;
+
+  if (d)
+    res2 = a * tmp;
+}
+
+/* { dg-final { scan-tree-dump-not "Optimizing reciprocal sqrt multiplications" "recip" } } */
+/* { dg-final { scan-tree-dump-not "Replacing squaring multiplication" "recip" } } */
+/* { dg-final { scan-tree-dump-not "Replacing original division" "recip" } } */
diff --git a/gcc/testsuite/gcc.dg/recip_sqrt_mult_5.c b/gcc/testsuite/gcc.dg/recip_sqrt_mult_5.c
new file mode 100644
index 0000000000000000000000000000000000000000..e871f0fcd4feb1687f9815e4babf4d0667a15ea8
--- /dev/null
+++ b/gcc/testsuite/gcc.dg/recip_sqrt_mult_5.c
@@ -0,0 +1,20 @@
+/* { dg-do compile } */
+/* { dg-options "-Ofast -fdump-tree-recip" } */
+
+/* We want to do the recip_sqrt transformations here there is already
+   a multiplication on the main path.  */
+
+double res, res2, tmp;
+void
+foo (double a, double b, int c, int d)
+{
+  tmp = 1.0 / __builtin_sqrt (a);
+  res = tmp * tmp;
+
+  if (d)
+    res2 = a * tmp;
+}
+
+/* { dg-final { scan-tree-dump "Optimizing reciprocal sqrt multiplications" "recip" } } */
+/* { dg-final { scan-tree-dump "Replacing squaring multiplication" "recip" } } */
+/* { dg-final { scan-tree-dump "Replacing original division" "recip" } } */
diff --git a/gcc/tree-ssa-math-opts.c b/gcc/tree-ssa-math-opts.c
index 25378da6f4ab27dbce51e10461efc18cc416a4d7..19bff5c3c3715f3e194e1c091ae01f6c4d6d2ce8 100644
--- a/gcc/tree-ssa-math-opts.c
+++ b/gcc/tree-ssa-math-opts.c
@@ -337,9 +337,9 @@ is_division_by (gimple *use_stmt, tree def)
 	 && gimple_assign_rhs1 (use_stmt) != def;
 }
 
-/* Return whether USE_STMT is DEF * DEF.  */
+/* Return TRUE if USE_STMT is a multiplication of DEF by A.  */
 static inline bool
-is_square_of (gimple *use_stmt, tree def)
+is_mult_by (gimple *use_stmt, tree def, tree a)
 {
   if (gimple_code (use_stmt) == GIMPLE_ASSIGN
       && gimple_assign_rhs_code (use_stmt) == MULT_EXPR)
@@ -347,11 +347,19 @@ is_square_of (gimple *use_stmt, tree def)
       tree op0 = gimple_assign_rhs1 (use_stmt);
       tree op1 = gimple_assign_rhs2 (use_stmt);
 
-      return op0 == op1 && op0 == def;
+      return (op0 == def && op1 == a)
+	      || (op0 == a && op1 == def);
     }
   return 0;
 }
 
+/* Return whether USE_STMT is DEF * DEF.  */
+static inline bool
+is_square_of (gimple *use_stmt, tree def)
+{
+  return is_mult_by (use_stmt, def, def);
+}
+
 /* Return whether USE_STMT is a floating-point division by
    DEF * DEF.  */
 static inline bool
@@ -526,6 +534,188 @@ free_bb (struct occurrence *occ)
     }
 }
 
+/* Transform sequences like
+   t = sqrt (a)
+   x = 1.0 / t;
+   r1 = x * x;
+   r2 = a * x;
+   into:
+   t = sqrt (a)
+   r1 = 1.0 / a;
+   r2 = t;
+   x = r1 * r2;
+   depending on the uses of x, r1, r2.  This removes one multiplication and
+   allows the sqrt and division operations to execute in parallel.
+   DEF_GSI is the gsi of the initial division by sqrt that defines
+   DEF (x in the example abovs).  */
+
+static void
+optimize_recip_sqrt (gimple_stmt_iterator *def_gsi, tree def)
+{
+  gimple *use_stmt;
+  imm_use_iterator use_iter;
+  gimple *stmt = gsi_stmt (*def_gsi);
+  tree x = def;
+  tree orig_sqrt_ssa_name = gimple_assign_rhs2 (stmt);
+  tree div_rhs1 = gimple_assign_rhs1 (stmt);
+
+  if (TREE_CODE (orig_sqrt_ssa_name) != SSA_NAME
+      || TREE_CODE (div_rhs1) != REAL_CST
+      || !real_equal (&TREE_REAL_CST (div_rhs1), &dconst1))
+    return;
+
+  gcall *sqrt_stmt
+    = dyn_cast <gcall *> (SSA_NAME_DEF_STMT (orig_sqrt_ssa_name));
+
+  if (!sqrt_stmt || !gimple_call_lhs (sqrt_stmt))
+    return;
+
+  switch (gimple_call_combined_fn (sqrt_stmt))
+    {
+    CASE_CFN_SQRT:
+    CASE_CFN_SQRT_FN:
+      break;
+
+    default:
+      return;
+    }
+  tree a = gimple_call_arg (sqrt_stmt, 0);
+
+  /* We have 'a' and 'x'.  Now analyze the uses of 'x'.  */
+
+  /* Statements that use x in x * x.  */
+  auto_vec<gimple *> sqr_stmts;
+  /* Statements that use x in a * x.  */
+  auto_vec<gimple *> mult_stmts;
+  bool has_other_use = false;
+  bool mult_on_main_path = false;
+
+  FOR_EACH_IMM_USE_STMT (use_stmt, use_iter, x)
+    {
+      if (is_gimple_debug (use_stmt))
+	continue;
+      if (is_square_of (use_stmt, x))
+	{
+	  sqr_stmts.safe_push (use_stmt);
+	  if (gimple_bb (use_stmt) == gimple_bb (stmt))
+	    mult_on_main_path = true;
+	}
+      else if (is_mult_by (use_stmt, x, a))
+	{
+	  mult_stmts.safe_push (use_stmt);
+	  if (gimple_bb (use_stmt) == gimple_bb (stmt))
+	    mult_on_main_path = true;
+	}
+      else
+	has_other_use = true;
+    }
+
+  /* In the x * x and a * x cases we just rewire stmt operands or
+     remove multiplications.  In the has_other_use case we introduce
+     a multiplication so make sure we don't introduce a multiplication
+     on a path where there was none.  */
+  if (has_other_use && !mult_on_main_path)
+    return;
+
+  if (sqr_stmts.is_empty () && mult_stmts.is_empty ())
+    return;
+
+  /* If x = 1.0 / sqrt (a) has uses other than those optimized here we want
+     to be able to compose it from the sqr and mult cases.  */
+  if (has_other_use && (sqr_stmts.is_empty () || mult_stmts.is_empty ()))
+    return;
+
+  if (dump_file)
+    {
+      fprintf (dump_file, "Optimizing reciprocal sqrt multiplications of\n");
+      print_gimple_stmt (dump_file, sqrt_stmt, 0, TDF_NONE);
+      print_gimple_stmt (dump_file, stmt, 0, TDF_NONE);
+      fprintf (dump_file, "\n");
+    }
+
+  bool delete_div = !has_other_use;
+  tree sqr_ssa_name = NULL_TREE;
+  if (!sqr_stmts.is_empty ())
+    {
+      /* r1 = x * x.  Transform the original
+	 x = 1.0 / t
+	 into
+	 tmp1 = 1.0 / a
+	 r1 = tmp1.  */
+
+      sqr_ssa_name
+	= make_temp_ssa_name (TREE_TYPE (a), NULL, "recip_sqrt_sqr");
+
+      if (dump_file)
+	{
+	  fprintf (dump_file, "Replacing original division\n");
+	  print_gimple_stmt (dump_file, stmt, 0, TDF_NONE);
+	  fprintf (dump_file, "with new division\n");
+	}
+      gimple_assign_set_lhs (stmt, sqr_ssa_name);
+      gimple_assign_set_rhs2 (stmt, a);
+      fold_stmt_inplace (def_gsi);
+      update_stmt (stmt);
+
+      if (dump_file)
+	print_gimple_stmt (dump_file, stmt, 0, TDF_NONE);
+
+      delete_div = false;
+      gimple *sqr_stmt;
+      unsigned int i;
+      FOR_EACH_VEC_ELT (sqr_stmts, i, sqr_stmt)
+	{
+	  gimple_stmt_iterator gsi2 = gsi_for_stmt (sqr_stmt);
+	  gimple_assign_set_rhs_from_tree (&gsi2, sqr_ssa_name);
+	  update_stmt (sqr_stmt);
+	}
+    }
+  if (!mult_stmts.is_empty ())
+    {
+      /* r2 = a * x.  Transform this into:
+	 r2 = t (The original sqrt (a)).  */
+      unsigned int i;
+      gimple *mult_stmt = NULL;
+      FOR_EACH_VEC_ELT (mult_stmts, i, mult_stmt)
+	{
+	  gimple_stmt_iterator gsi2 = gsi_for_stmt (mult_stmt);
+
+	  if (dump_file)
+	    {
+	      fprintf (dump_file, "Replacing squaring multiplication\n");
+	      print_gimple_stmt (dump_file, mult_stmt, 0, TDF_NONE);
+	      fprintf (dump_file, "with assignment\n");
+	    }
+	  gimple_assign_set_rhs_from_tree (&gsi2, orig_sqrt_ssa_name);
+	  fold_stmt_inplace (&gsi2);
+	  update_stmt (mult_stmt);
+	  if (dump_file)
+	    print_gimple_stmt (dump_file, mult_stmt, 0, TDF_NONE);
+      }
+    }
+
+  if (has_other_use)
+    {
+      /* Using the two temporaries tmp1, tmp2 from above
+	 the original x is now:
+	 x = tmp1 * tmp2.  */
+      gcc_assert (orig_sqrt_ssa_name);
+      gcc_assert (sqr_ssa_name);
+
+      gimple *new_stmt
+	= gimple_build_assign (x, MULT_EXPR,
+				orig_sqrt_ssa_name, sqr_ssa_name);
+      gsi_insert_after (def_gsi, new_stmt, GSI_NEW_STMT);
+      update_stmt (stmt);
+    }
+  else if (delete_div)
+    {
+      /* Remove the original division.  */
+      gimple_stmt_iterator gsi2 = gsi_for_stmt (stmt);
+      gsi_remove (&gsi2, true);
+      release_defs (stmt);
+    }
+}
 
 /* Look for floating-point divisions among DEF's uses, and try to
    replace them by multiplications with the reciprocal.  Add
@@ -756,7 +946,15 @@ pass_cse_reciprocals::execute (function *fun)
 	      && (def = SINGLE_SSA_TREE_OPERAND (stmt, SSA_OP_DEF)) != NULL
 	      && FLOAT_TYPE_P (TREE_TYPE (def))
 	      && TREE_CODE (def) == SSA_NAME)
-	    execute_cse_reciprocals_1 (&gsi, def);
+	    {
+	      if (flag_unsafe_math_optimizations
+		  && is_gimple_assign (stmt)
+		  && !stmt_can_throw_internal (stmt)
+		  && gimple_assign_rhs_code (stmt) == RDIV_EXPR)
+		optimize_recip_sqrt (&gsi, def);
+	      else
+		execute_cse_reciprocals_1 (&gsi, def);
+	    }
 	}
 
       if (optimize_bb_for_size_p (bb))