[Gmp-commit] /var/hg/gmp: 3 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Sun Feb 26 17:29:57 CET 2012
details: /var/hg/gmp/rev/f2f516affc0c
changeset: 14679:f2f516affc0c
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Sun Feb 26 17:25:30 2012 +0100
description:
New function mpz_2fac_ui.
details: /var/hg/gmp/rev/d65a47bed7ad
changeset: 14680:d65a47bed7ad
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Sun Feb 26 17:27:13 2012 +0100
description:
New tests for mpz_2fac_ui.
details: /var/hg/gmp/rev/422fd744313a
changeset: 14681:422fd744313a
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Sun Feb 26 17:29:50 2012 +0100
description:
for -> while
diffstat:
ChangeLog | 7 +++
Makefile.am | 3 +-
gmp-h.in | 3 +
mpz/2fac_ui.c | 98 ++++++++++++++++++++++++++++++++++++++++++++++++++++
mpz/fac_ui.c | 2 +-
tests/mpz/t-fac_ui.c | 25 +++++++++++--
6 files changed, 132 insertions(+), 6 deletions(-)
diffs (214 lines):
diff -r b3d04bbfafa7 -r 422fd744313a ChangeLog
--- a/ChangeLog Fri Feb 24 17:45:54 2012 +0100
+++ b/ChangeLog Sun Feb 26 17:29:50 2012 +0100
@@ -1,3 +1,10 @@
+2012-02-26 Marco Bodrato <bodrato at mail.dm.unipi.it>
+
+ * mpz/2fac_ui.c: New file: implements n!!.
+ * Makefile.am (MPZ_OBJECTS): Add mpz/2fac_ui.
+ * gmp-h.in: Declare mpz_2fac_ui.
+ * tests/mpz/t-fac.c: Test mpz_2fac_ui.
+
2012-02-24 Marco Bodrato <bodrato at mail.dm.unipi.it>
* mpz/invert.c: Use ABSIZ, MPZ_EQUAL_1_P.
diff -r b3d04bbfafa7 -r 422fd744313a Makefile.am
--- a/Makefile.am Fri Feb 24 17:45:54 2012 +0100
+++ b/Makefile.am Sun Feb 26 17:29:50 2012 +0100
@@ -157,7 +157,8 @@
mpz/cong$U.lo mpz/cong_2exp$U.lo mpz/cong_ui$U.lo \
mpz/divexact$U.lo mpz/divegcd$U.lo mpz/dive_ui$U.lo \
mpz/divis$U.lo mpz/divis_ui$U.lo mpz/divis_2exp$U.lo mpz/dump$U.lo \
- mpz/export$U.lo mpz/fac_ui$U.lo mpz/oddfac_1$U.lo mpz/prodlimbs$U.lo \
+ mpz/export$U.lo \
+ mpz/2fac_ui$U.lo mpz/fac_ui$U.lo mpz/oddfac_1$U.lo mpz/prodlimbs$U.lo \
mpz/fdiv_q_ui$U.lo mpz/fdiv_qr$U.lo mpz/fdiv_qr_ui$U.lo \
mpz/fdiv_r$U.lo mpz/fdiv_r_ui$U.lo mpz/fdiv_q$U.lo \
mpz/fdiv_ui$U.lo mpz/fib_ui$U.lo mpz/fib2_ui$U.lo mpz/fits_sint$U.lo \
diff -r b3d04bbfafa7 -r 422fd744313a gmp-h.in
--- a/gmp-h.in Fri Feb 24 17:45:54 2012 +0100
+++ b/gmp-h.in Sun Feb 26 17:29:50 2012 +0100
@@ -811,6 +811,9 @@
#define mpz_fac_ui __gmpz_fac_ui
__GMP_DECLSPEC void mpz_fac_ui __GMP_PROTO ((mpz_ptr, unsigned long int));
+#define mpz_2fac_ui __gmpz_2fac_ui
+__GMP_DECLSPEC void mpz_2fac_ui __GMP_PROTO ((mpz_ptr, unsigned long int));
+
#define mpz_fdiv_q __gmpz_fdiv_q
__GMP_DECLSPEC void mpz_fdiv_q __GMP_PROTO ((mpz_ptr, mpz_srcptr, mpz_srcptr));
diff -r b3d04bbfafa7 -r 422fd744313a mpz/2fac_ui.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpz/2fac_ui.c Sun Feb 26 17:29:50 2012 +0100
@@ -0,0 +1,98 @@
+/* mpz_2fac_ui(RESULT, N) -- Set RESULT to N!!.
+
+Contributed to the GNU project by Marco Bodrato.
+
+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"
+
+#include "fac_ui.h"
+
+/* This is intended for constant THRESHOLDs only, where the compiler
+ can completely fold the result. */
+#define LOG2C(n) \
+ (((n) >= 0x1) + ((n) >= 0x2) + ((n) >= 0x4) + ((n) >= 0x8) + \
+ ((n) >= 0x10) + ((n) >= 0x20) + ((n) >= 0x40) + ((n) >= 0x80) + \
+ ((n) >= 0x100) + ((n) >= 0x200) + ((n) >= 0x400) + ((n) >= 0x800) + \
+ ((n) >= 0x1000) + ((n) >= 0x2000) + ((n) >= 0x4000) + ((n) >= 0x8000))
+
+#define FACTOR_LIST_STORE(P, PR, MAX_PR, VEC, I) \
+ do { \
+ if ((PR) > (MAX_PR)) { \
+ (VEC)[(I)++] = (PR); \
+ (PR) = (P); \
+ } else \
+ (PR) *= (P); \
+ } while (0)
+
+#define FAC_2DSC_THRESHOLD ((FAC_DSC_THRESHOLD << 1) | (FAC_DSC_THRESHOLD & 1))
+#define FACTORS_PER_LIMB (GMP_NUMB_BITS / (LOG2C(FAC_2DSC_THRESHOLD-1)+1))
+
+/* Computes n!!, the 2-multi-factorial of n. (aka double-factorial or semi-factorial)
+ WARNING: it assumes that n fits in a limb!
+ */
+void
+mpz_2fac_ui (mpz_ptr x, unsigned long n)
+{
+ ASSERT (n <= GMP_NUMB_MAX);
+
+ if ((n & 1) == 0) { /* n is even, n = 2k, (2k)!! = k! 2^k */
+ mp_limb_t count;
+
+ popc_limb (count, n);
+ count = n - count;
+ n >>= 1;
+ mpz_oddfac_1 (x, n, 0);
+ mpz_mul_2exp (x, x, count);
+ } else { /* n is odd */
+ static const mp_limb_t tabled[] = { ONE_LIMB_ODD_DOUBLEFACTORIAL_TABLE };
+
+ if (n < 2 * numberof (tabled)) {
+ PTR (x)[0] = tabled[n >> 1];
+ SIZ (x) = 1;
+ } else if (BELOW_THRESHOLD (n, FAC_2DSC_THRESHOLD)) { /* odd basecase, */
+ mp_limb_t *factors, prod, max_prod, j;
+ TMP_SDECL;
+
+ /* FIXME: we might alloc a fixed ammount 1+FAC_2DSC_THRESHOLD/FACTORS_PER_LIMB */
+ TMP_SMARK;
+ factors = TMP_SALLOC_LIMBS (1 + n / (2 * FACTORS_PER_LIMB));
+
+ factors[0] = tabled[numberof (tabled) - 1];
+ j = 1;
+ prod = n;
+
+ max_prod = GMP_NUMB_MAX / FAC_2DSC_THRESHOLD;
+ while ((n -= 2) >= numberof (tabled) * 2 + 1)
+ FACTOR_LIST_STORE (n, prod, max_prod, factors, j);
+
+ factors[j++] = prod;
+ mpz_prodlimbs (x, factors, j);
+
+ TMP_SFREE;
+ } else { /* for the asymptotically fast odd case, let oddfac do the job. */
+ mpz_oddfac_1 (x, n, 1);
+ }
+ }
+}
+
+#undef LOG2C
+#undef FACTORS_PER_LIMB
+#undef FACTOR_LIST_STORE
+#undef FAC_2DSC_THRESHOLD
diff -r b3d04bbfafa7 -r 422fd744313a mpz/fac_ui.c
--- a/mpz/fac_ui.c Fri Feb 24 17:45:54 2012 +0100
+++ b/mpz/fac_ui.c Sun Feb 26 17:29:50 2012 +0100
@@ -81,7 +81,7 @@
#else
max_prod = GMP_NUMB_MAX / (FAC_ODD_THRESHOLD | 1);
#endif
- for (; --n >= numberof (table);)
+ while (--n >= numberof (table))
FACTOR_LIST_STORE (n, prod, max_prod, factors, j);
factors[j++] = prod;
diff -r b3d04bbfafa7 -r 422fd744313a tests/mpz/t-fac_ui.c
--- a/tests/mpz/t-fac_ui.c Fri Feb 24 17:45:54 2012 +0100
+++ b/tests/mpz/t-fac_ui.c Sun Feb 26 17:29:50 2012 +0100
@@ -1,4 +1,4 @@
-/* Exercise mpz_fac_ui.
+/* Exercise mpz_fac_ui and mpz_2fac_ui.
Copyright 2000, 2001, 2002, 2012 Free Software Foundation, Inc.
@@ -35,9 +35,9 @@
int
main (int argc, char *argv[])
{
- unsigned long n;
+ unsigned long n, m;
unsigned long limit = 2222;
- mpz_t f, r;
+ mpz_t df[2], f, r;
tests_start ();
@@ -49,10 +49,12 @@
/* for small limb testing */
limit = MIN (limit, MP_LIMB_T_MAX);
+ mpz_init_set_ui (df[0], 1); /* 0!! = 1 */
+ mpz_init_set_ui (df[1], 1); /* -1!! = 1 */
mpz_init_set_ui (f, 1); /* 0! = 1 */
mpz_init (r);
- for (n = 0; n < limit; n++)
+ for (n = 0, m = 0; n < limit; n++)
{
mpz_fac_ui (r, n);
MPZ_CHECK_FORMAT (r);
@@ -65,9 +67,24 @@
abort ();
}
+ mpz_2fac_ui (r, n);
+ MPZ_CHECK_FORMAT (r);
+
+ if (mpz_cmp (df[m], r) != 0)
+ {
+ printf ("mpz_2fac_ui(%lu) wrong\n", n);
+ printf (" got "); mpz_out_str (stdout, 10, r); printf("\n");
+ printf (" want "); mpz_out_str (stdout, 10, df[m]); printf("\n");
+ abort ();
+ }
+
+ m ^= 1;
+ mpz_mul_ui (df[m], df[m], n+1); /* (n+1)!! = (n-1)!! * (n+1) */
mpz_mul_ui (f, f, n+1); /* (n+1)! = n! * (n+1) */
}
+ mpz_clear (df[0]);
+ mpz_clear (df[1]);
mpz_clear (f);
mpz_clear (r);
More information about the gmp-commit
mailing list