[Gmp-commit] /var/hg/gmp: oddfac: small speed-up for small argument.

[mercurial at gmplib.org]

details:   /var/hg/gmp/rev/07e7c66b0c8a
changeset: 14690:07e7c66b0c8a
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Feb 27 08:25:17 2012 +0100
description:
oddfac: small speed-up for small argument.

diffstat:

 ChangeLog      |   2 ++
 mpz/oddfac_1.c |  38 +++++++++++++++++++++++---------------
 2 files changed, 25 insertions(+), 15 deletions(-)

diffs (84 lines):

diff -r bfe92c3a6b72 -r 07e7c66b0c8a ChangeLog
--- a/ChangeLog	Mon Feb 27 07:36:24 2012 +0100
+++ b/ChangeLog	Mon Feb 27 08:25:17 2012 +0100
@@ -17,6 +17,8 @@
 	* doc/gmp.texi: Document mpz_2fac_ui.
 	* mpz/Makefile.am (libmpz_la_SOURCES): Add 2fac_ui.c.
 
+	* mpz/oddfac_1.c (mpz_oddfac_1): Use umul_ppmm when size = 2. 
+
 2012-02-26  Niels Möller  <nisse at lysator.liu.se>
 
 	* mini-gmp: New subdirectory. For use by GMP bootstrap, and as a
diff -r bfe92c3a6b72 -r 07e7c66b0c8a mpz/oddfac_1.c
--- a/mpz/oddfac_1.c	Mon Feb 27 07:36:24 2012 +0100
+++ b/mpz/oddfac_1.c	Mon Feb 27 08:25:17 2012 +0100
@@ -488,6 +488,11 @@
    if n is odd, n > FAC_DSC_THRESHOLD and flag = 1, x is set to n!!.
 
    If n is too small, flag is ignored, and an ASSERT can be triggered.
+
+   TODO: FAC_DSC_THRESHOLD is used here with two different roles:
+    - to decide when prime factorisation is needed,
+    - to stop the recursion, once sieving is done.
+   Maybe two thresholds can do a better job.
  */
 void
 mpz_oddfac_1 (mpz_ptr x, mp_limb_t n, unsigned flag)
@@ -512,8 +517,6 @@
       {
 	static const mp_limb_t tabled[] = { ONE_LIMB_ODD_DOUBLEFACTORIAL_TABLE };
 	mp_limb_t tn;
-	mp_size_t j;
-	TMP_SDECL;
 
 #if TUNE_PROGRAM_BUILD
 	ASSERT (FAC_DSC_THRESHOLD_LIMIT >= FAC_DSC_THRESHOLD);
@@ -523,14 +526,16 @@
 	for (tn = n; ABOVE_THRESHOLD (tn, FAC_DSC_THRESHOLD); s++)
 	  tn >>= 1;
 
-	j = 0;
-
-	TMP_SMARK;
-	factors = TMP_SALLOC_LIMBS (1 + tn / FACTORS_PER_LIMB);
-	ASSERT (tn >= FACTORS_PER_LIMB);
-
 	if (tn >= numberof (tabled) * 2 + 1) {
 	  mp_limb_t prod, max_prod, i;
+	  mp_size_t j;
+	  TMP_SDECL;
+
+	  j = 0;
+
+	  TMP_SMARK;
+	  factors = TMP_SALLOC_LIMBS (1 + tn / FACTORS_PER_LIMB);
+	  ASSERT (tn >= FACTORS_PER_LIMB);
 
 	  prod = 1;
 #if TUNE_PROGRAM_BUILD
@@ -551,14 +556,17 @@
 	  } while (tn >= numberof (tabled) * 2 + 1);
 
 	  factors[j++] = prod;
+	  ASSERT (numberof (tablef) > numberof (tabled));
+	  factors[j++] = tabled[(tn - 1) >> 1];
+	  factors[j++] = tablef[tn >> 1];
+	  mpz_prodlimbs (x, factors, j);
+
+	  TMP_SFREE;
+	} else {
+	  MPZ_REALLOC (x, 2);
+	  umul_ppmm (PTR (x)[1], PTR (x)[0], tabled[(tn - 1) >> 1], tablef[tn >> 1]);
+	  SIZ (x) = 2;
 	}
-
-	ASSERT (numberof (tablef) > numberof (tabled));
-	factors[j++] = tabled[(tn - 1) >> 1];
-	factors[j++] = tablef[tn >> 1];
-	mpz_prodlimbs (x, factors, j);
-
-	TMP_SFREE;
       }
 
       if (s != 0)