mpz_prodlimbs: removed from mpz/fac_ui.c, now in a new file.

3 files changed, 100 insertions, 75 deletions

diff --git a/mpz/Makefile.am b/mpz/Makefile.am
index 74c2c3444..35ec7a1b1 100644
--- a/mpz/Makefile.am
+++ b/mpz/Makefile.am
@@ -1,6 +1,6 @@
 ## Process this file with automake to generate Makefile.in
 
-# Copyright 1996, 1998, 1999, 2000, 2001, 2002, 2003 Free Software
+# Copyright 1996, 1998, 1999, 2000, 2001, 2002, 2003, 2012 Free Software
 # Foundation, Inc.
 #
 # This file is part of the GNU MP Library.
@@ -45,16 +45,10 @@ libmpz_la_SOURCES = aors.h aors_ui.h fits_s.h mul_i.h \
   lcm.c lcm_ui.c lucnum_ui.c lucnum2_ui.c millerrabin.c \
   mod.c mul.c mul_2exp.c mul_si.c mul_ui.c n_pow_ui.c neg.c nextprime.c \
   out_raw.c out_str.c perfpow.c perfsqr.c popcount.c pow_ui.c powm.c \
-  powm_sec.c powm_ui.c pprime_p.c random.c random2.c \
+  powm_sec.c powm_ui.c pprime_p.c prodlimbs.c random.c random2.c \
   realloc.c realloc2.c remove.c root.c rootrem.c rrandomb.c \
   scan0.c scan1.c set.c set_d.c set_f.c set_q.c set_si.c set_str.c \
   set_ui.c setbit.c size.c sizeinbase.c sqrt.c sqrtrem.c sub.c sub_ui.c \
   swap.c tdiv_ui.c tdiv_q.c tdiv_q_2exp.c tdiv_q_ui.c tdiv_qr.c \
   tdiv_qr_ui.c tdiv_r.c tdiv_r_2exp.c tdiv_r_ui.c tstbit.c ui_pow_ui.c \
   ui_sub.c urandomb.c urandomm.c xor.c
-
-# These are BUILT_SOURCES at the top-level, so normally they're built before
-# recursing into this directory.
-#
-fac_ui.h:
-	cd ..; $(MAKE) $(AM_MAKEFLAGS) mpz/fac_ui.h
diff --git a/mpz/fac_ui.c b/mpz/fac_ui.c
index 58be91b82..b1801509a 100644
--- a/mpz/fac_ui.c
+++ b/mpz/fac_ui.c
@@ -294,73 +294,6 @@ bitwise_primesieve (mp_ptr bit_array, mp_limb_t n)
 #undef BLOCK_SIZE
 
 /*********************************************************/
-/* Section list-prod: product of a list -> mpz_t         */
-/*********************************************************/
-
-/* FIXME: should be tuned */
-#ifndef RECURSIVE_PROD_THRESHOLD
-#define RECURSIVE_PROD_THRESHOLD (MUL_TOOM22_THRESHOLD|(MUL_TOOM22_THRESHOLD>>1))
-#endif
-
-/* Computes the product of the j>1 limbs pointed by factors, puts the
-   result in x.
-   The list in {factors, j} is overwritten.
-*/
-
-static void
-mpz_prodlimbs (mpz_ptr x, mp_limb_t *factors, mp_limb_t j)
-{
-  mp_limb_t i;
-  mp_size_t size;
-  mp_ptr    prod;
-
-  ASSERT (j > 1);
-  ASSERT (RECURSIVE_PROD_THRESHOLD > 3);
-
-  if (BELOW_THRESHOLD (j, RECURSIVE_PROD_THRESHOLD)) {
-    mp_limb_t cy;
-
-    j--;
-    size = 1;
-
-    for (i = 1; i < j; i++)
-      {
-	cy = mpn_mul_1 (factors, factors, size, factors[i]);
-	factors[size] = cy;
-	size += cy != 0;
-      };
-
-    prod = MPZ_REALLOC (x, size + 1);
-
-    cy = mpn_mul_1 (prod, factors, size, factors[i]);
-    prod[size] = cy;
-    SIZ (x) = size + (cy != 0);
-  } else {
-    mpz_t x1, x2;
-    TMP_DECL;
-
-    i = j >> 1;
-    TMP_MARK;
-
-    MPZ_TMP_INIT (x2, j - i);
-
-    PTR (x1) = factors + i;
-    ALLOC (x1) = j - i;
-    mpz_prodlimbs (x2, factors + i, j - i);
-    mpz_prodlimbs (x1, factors, i);
-    size = SIZ (x1) + SIZ (x2);
-    prod = MPZ_REALLOC (x, size);
-    if (SIZ (x1) >= SIZ (x2))
-      i = mpn_mul (prod, PTR(x1), SIZ(x1), PTR(x2), SIZ(x2));
-    else
-      i = mpn_mul (prod, PTR(x2), SIZ(x2), PTR(x1), SIZ(x1));
-    SIZ (x) = size - (i == 0);
-
-    TMP_FREE;
-  }
-}
-
-/*********************************************************/
 /* Section swing: swing factorial                        */
 /*********************************************************/
 
diff --git a/mpz/prodlimbs.c b/mpz/prodlimbs.c
new file mode 100644
index 000000000..0352db5ad
--- /dev/null
+++ b/mpz/prodlimbs.c
@@ -0,0 +1,98 @@
+/* mpz_prodlimps(RESULT, V, LEN) -- Set RESULT to V[0]*V[1]*...*V[LEN-1].
+
+Contributed to the GNU project by Marco Bodrato.
+
+THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.
+IT IS ONLY SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.
+IN FACT, IT IS ALMOST GUARANTEED THAT IT WILL CHANGE OR
+DISAPPEAR IN A FUTURE GNU MP RELEASE.
+
+Copyright 2012 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+/*********************************************************/
+/* Section list-prod: product of a list -> mpz_t         */
+/*********************************************************/
+
+/* FIXME: should be tuned */
+#ifndef RECURSIVE_PROD_THRESHOLD
+#define RECURSIVE_PROD_THRESHOLD (MUL_TOOM22_THRESHOLD|(MUL_TOOM22_THRESHOLD<<1))
+#endif
+
+/* Computes the product of the j>1 limbs pointed by factors, puts the
+ * result in x. It assumes that all limbs are non-zero. Above
+ * Karatsuba's threshold it uses a binary splitting startegy, to gain
+ * speed by the asymptotically fast multiplication algorithms.
+ *
+ * The list in  {factors, j} is overwritten.
+ * Returns the size of the result
+ */
+
+mp_size_t
+mpz_prodlimbs (mpz_ptr x, mp_ptr factors, mp_size_t j)
+{
+  mp_limb_t cy;
+  mp_size_t size, i;
+  mp_ptr    prod;
+
+  ASSERT (j > 1);
+  ASSERT (RECURSIVE_PROD_THRESHOLD > 3);
+
+  if (BELOW_THRESHOLD (j, RECURSIVE_PROD_THRESHOLD)) {
+    j--;
+    size = 1;
+
+    for (i = 1; i < j; i++)
+      {
+	cy = mpn_mul_1 (factors, factors, size, factors[i]);
+	factors[size] = cy;
+	size += cy != 0;
+      };
+
+    prod = MPZ_REALLOC (x, size + 1);
+
+    cy = mpn_mul_1 (prod, factors, size, factors[i]);
+    prod[size] = cy;
+    return SIZ (x) = size + (cy != 0);
+  } else {
+    mpz_t x1, x2;
+    TMP_DECL;
+
+    i = j >> 1;
+    j -= i;
+    TMP_MARK;
+
+    MPZ_TMP_INIT (x2, j);
+
+    PTR (x1) = factors + i;
+    ALLOC (x1) = j;
+    j = mpz_prodlimbs (x2, factors + i, j);
+    i = mpz_prodlimbs (x1, factors, i);
+    size = i + j;
+    prod = MPZ_REALLOC (x, size);
+    if (i >= j)
+      cy = mpn_mul (prod, PTR(x1), i, PTR(x2), j);
+    else
+      cy = mpn_mul (prod, PTR(x2), j, PTR(x1), i);
+    TMP_FREE;
+
+    return SIZ (x) = size - (cy == 0);
+  }
+}