Rewrite to invoke mpn_sqr when appropriate.

1 files changed, 86 insertions, 48 deletions

diff --git a/mpf/mul.c b/mpf/mul.c
index b36358ce2..bb3ac25a1 100644
--- a/mpf/mul.c
+++ b/mpf/mul.c
@@ -1,6 +1,6 @@
 /* mpf_mul -- Multiply two floats.
 
-Copyright 1993, 1994, 1996, 2001, 2005 Free Software Foundation, Inc.
+Copyright 1993, 1994, 1996, 2001, 2005, 2019 Free Software Foundation, Inc.
 
 This file is part of the GNU MP Library.
 
@@ -33,64 +33,102 @@ see https://www.gnu.org/licenses/.  */
 void
 mpf_mul (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
 {
-  mp_srcptr up, vp;
-  mp_size_t usize, vsize;
   mp_size_t sign_product;
   mp_size_t prec = r->_mp_prec;
+  mp_size_t rsize;
+  mp_limb_t cy_limb;
+  mp_ptr rp, tp;
+  mp_size_t adj;
+  TMP_DECL;
 
-  usize = u->_mp_size;
-  vsize = v->_mp_size;
-  sign_product = usize ^ vsize;
+  if (u == v)
+    {
+      mp_srcptr up;
+      mp_size_t usize;
 
-  usize = ABS (usize);
-  vsize = ABS (vsize);
+      usize = u->_mp_size;
+      sign_product = usize;
 
-  up = u->_mp_d;
-  vp = v->_mp_d;
-  if (usize > prec)
-    {
-      up += usize - prec;
-      usize = prec;
-    }
-  if (vsize > prec)
-    {
-      vp += vsize - prec;
-      vsize = prec;
-    }
+      usize = ABS (usize);
 
-  if (usize == 0 || vsize == 0)
-    {
-      r->_mp_size = 0;
-      r->_mp_exp = 0;		/* ??? */
+      up = u->_mp_d;
+      if (usize > prec)
+	{
+	  up += usize - prec;
+	  usize = prec;
+	}
+
+      if (usize == 0)
+	{
+	  r->_mp_size = 0;
+	  r->_mp_exp = 0;		/* ??? */
+	  return;
+	}
+      else
+	{
+	  TMP_MARK;
+	  rsize = 2 * usize;
+	  tp = TMP_ALLOC_LIMBS (rsize);
+
+	  mpn_sqr (tp, up, usize);
+	  cy_limb = tp[rsize - 1];
+	}
     }
   else
     {
-      mp_size_t rsize;
-      mp_limb_t cy_limb;
-      mp_ptr rp, tp;
-      mp_size_t adj;
-      TMP_DECL;
-
-      TMP_MARK;
-      rsize = usize + vsize;
-      tp = TMP_ALLOC_LIMBS (rsize);
-      cy_limb = (usize >= vsize
-		 ? mpn_mul (tp, up, usize, vp, vsize)
-		 : mpn_mul (tp, vp, vsize, up, usize));
-
-      adj = cy_limb == 0;
-      rsize -= adj;
-      prec++;
-      if (rsize > prec)
+      mp_srcptr up, vp;
+      mp_size_t usize, vsize;
+
+      usize = u->_mp_size;
+      vsize = v->_mp_size;
+      sign_product = usize ^ vsize;
+
+      usize = ABS (usize);
+      vsize = ABS (vsize);
+
+      up = u->_mp_d;
+      vp = v->_mp_d;
+      if (usize > prec)
+	{
+	  up += usize - prec;
+	  usize = prec;
+	}
+      if (vsize > prec)
+	{
+	  vp += vsize - prec;
+	  vsize = prec;
+	}
+
+      if (usize == 0 || vsize == 0)
 	{
-	  tp += rsize - prec;
-	  rsize = prec;
+	  r->_mp_size = 0;
+	  r->_mp_exp = 0;
+	  return;
 	}
-      rp = r->_mp_d;
-      MPN_COPY (rp, tp, rsize);
-      r->_mp_exp = u->_mp_exp + v->_mp_exp - adj;
-      r->_mp_size = sign_product >= 0 ? rsize : -rsize;
+      else
+	{
+	  TMP_MARK;
+	  rsize = usize + vsize;
+	  tp = TMP_ALLOC_LIMBS (rsize);
+	  cy_limb = (usize >= vsize
+		     ? mpn_mul (tp, up, usize, vp, vsize)
+		     : mpn_mul (tp, vp, vsize, up, usize));
 
-      TMP_FREE;
+	}
     }
+
+  adj = cy_limb == 0;
+  rsize -= adj;
+  prec++;
+  if (rsize > prec)
+    {
+      tp += rsize - prec;
+      rsize = prec;
+    }
+  rp = r->_mp_d;
+  MPN_COPY (rp, tp, rsize);
+  r->_mp_exp = u->_mp_exp + v->_mp_exp - adj;
+  r->_mp_size = sign_product >= 0 ? rsize : -rsize;
+
+  TMP_FREE;
 }