[Gmp-commit] /var/hg/gmp: 4 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Mon Apr 30 17:38:27 CEST 2012
details: /var/hg/gmp/rev/f08862bcddd1
changeset: 14916:f08862bcddd1
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Apr 30 17:20:12 2012 +0200
description:
Added to mercurial the new files...
details: /var/hg/gmp/rev/ea7367b71680
changeset: 14917:ea7367b71680
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Apr 30 17:24:57 2012 +0200
description:
Changelog
details: /var/hg/gmp/rev/b458c58fe985
changeset: 14918:b458c58fe985
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Apr 30 17:25:55 2012 +0200
description:
doc/gmp.texi: Document mpz_mfac_uiui.
details: /var/hg/gmp/rev/f0e71d6941c6
changeset: 14919:f0e71d6941c6
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Apr 30 17:38:21 2012 +0200
description:
doc/gmp.texi: Collapse factorial functions.
diffstat:
ChangeLog | 18 +++--
doc/gmp.texi | 25 ++++----
mpz/mfac_uiui.c | 119 ++++++++++++++++++++++++++++++++++++++++++
tests/mpz/t-mfac_uiui.c | 136 ++++++++++++++++++++++++++++++++++++++++++++++++
4 files changed, 277 insertions(+), 21 deletions(-)
diffs (truncated from 334 to 300 lines):
diff -r 8703c179e3b3 -r f0e71d6941c6 ChangeLog
--- a/ChangeLog Mon Apr 30 16:25:05 2012 +0200
+++ b/ChangeLog Mon Apr 30 17:38:21 2012 +0200
@@ -10,6 +10,16 @@
* mpz/oddfac_1.c: Likewise.
* mpz/primorial_ui.c: Likewise.
+ * mpz/mfac_uiui.c: New file.
+ * Makefile.am: Compile it.
+ * mpz/Makefile.am (libmpz_la_SOURCES): Add mpz_mfac_uiui.c
+ * gmp-h.in (mpz_mfac_uiui): Declare.
+
+ * tests/mpz/t-mfac_uiui.c: New file.
+ * tests/mpz/Makefile.am: Run it.
+
+ * doc/gmp.texi: Document mpz_mfac_uiui, collapsing with other factorial functions.
+
2012-04-29 Torbjorn Granlund <tege at gmplib.org>
* mpn/arm/aors_n.asm: Tune for more stable performance.
@@ -30,14 +40,6 @@
* mpn/arm/v6/addmul_1.asm: Rewrite for stable speed, smaller size.
* mpn/arm/v6/mul_1.asm: Likewise.
- * mpz/mfac_uiui.c: New file.
- * Makefile.am: Compile it.
- * mpz/Makefile.am (libmpz_la_SOURCES): Add mpz_mfac_uiui.c
- * gmp-h.in (mpz_mfac_uiui): Declare.
-
- * tests/mpz/t-mfac_uiui.c: New file.
- * tests/mpz/Makefile.am: Run it.
-
2012-04-27 Torbjorn Granlund <tege at gmplib.org>
* configure.in: Search arm/v6t2 for arm7.
diff -r 8703c179e3b3 -r f0e71d6941c6 doc/gmp.texi
--- a/doc/gmp.texi Mon Apr 30 16:25:05 2012 +0200
+++ b/doc/gmp.texi Mon Apr 30 17:38:21 2012 +0200
@@ -3631,20 +3631,19 @@
removed.
@end deftypefun
- at deftypefun void mpz_fac_ui (mpz_t @var{rop}, unsigned long int @var{op})
+ at deftypefun void mpz_fac_ui (mpz_t @var{rop}, unsigned long int @var{n})
+ at deftypefunx void mpz_2fac_ui (mpz_t @var{rop}, unsigned long int @var{n})
+ at deftypefunx void mpz_mfac_uiui (mpz_t @var{rop}, unsigned long int @var{n}, unsigned long int @var{m})
@cindex Factorial functions
-Set @var{rop} to @var{op}!, the factorial of @var{op}.
- at end deftypefun
-
- at deftypefun void mpz_2fac_ui (mpz_t @var{rop}, unsigned long int @var{op})
- at cindex Factorial functions
-Set @var{rop} to @var{op}!!, the double-factorial of @var{op}.
- at end deftypefun
-
- at deftypefun void mpz_primorial_ui (mpz_t @var{rop}, unsigned long int @var{op})
- at cindex Factorial functions
-Set @var{rop} to the primorial of @var{op}, i.e. the product of all positive
-prime numbers smaller then or equal to @var{op}.
+Set @var{rop} to the factorial of @var{n}: @code{mpz_fac_ui} computes the plain factorial @var{n}!,
+ at code{mpz_2fac_ui} computes the double-factorial @var{n}!!, and @code{mpz_mfac_uiui} the
+ at var{m}-multi-factorial @m{n!^{(m)}, @var{n}!^(@var{m})}.
+ at end deftypefun
+
+ at deftypefun void mpz_primorial_ui (mpz_t @var{rop}, unsigned long int @var{n})
+ at cindex Primorial functions
+Set @var{rop} to the primorial of @var{n}, i.e. the product of all positive
+prime numbers smaller then or equal to @var{n}.
@end deftypefun
@deftypefun void mpz_bin_ui (mpz_t @var{rop}, mpz_t @var{n}, unsigned long int @var{k})
diff -r 8703c179e3b3 -r f0e71d6941c6 mpz/mfac_uiui.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpz/mfac_uiui.c Mon Apr 30 17:38:21 2012 +0200
@@ -0,0 +1,119 @@
+/* mpz_mfac_uiui(RESULT, N, M) -- Set RESULT to N!^(M) = N(N-M)(N-2M)...
+
+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"
+
+/*************************************************************/
+/* Section macros: common macros, for swing/fac/bin (&sieve) */
+/*************************************************************/
+
+#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)
+
+/*********************************************************/
+/* Section oder factorials: */
+/*********************************************************/
+
+/* mpz_mfac_uiui (x, n, m) computes x = n!^(m) = n*(n-m)*(n-2m)*... */
+
+void
+mpz_mfac_uiui (mpz_ptr x, unsigned long n, unsigned long m)
+{
+ ASSERT (n <= GMP_NUMB_MAX);
+ ASSERT (m != 0);
+
+ if (n < 3 || n - 3 < m - 1) { /* (n < 3 || n - 1 <= m || m == 0) */
+ PTR (x)[0] = n + (n == 0);
+ SIZ (x) = 1;
+ } else { /* m < n - 1 < GMP_NUMB_MAX */
+ mp_limb_t g, sn;
+ mpz_t t;
+
+ g = mpn_gcd_1 (&n, 1, m);
+ if (g != 1) { n/=g; m/=g; }
+
+ if (m <= 2) { /* fac or 2fac */
+ if (m == 1) {
+ if (g > 2) {
+ mpz_init (t);
+ mpz_fac_ui (t, n);
+ sn = n;
+ } else if (g == 2) {
+ mpz_2fac_ui (x, n << 1);
+ g = 1;
+ } else
+ mpz_fac_ui (x, n);
+ } else { /* m == 2 */
+ if (g != 1) {
+ mpz_init (t);
+ mpz_2fac_ui (t, n);
+ sn = n / 2 + 1;
+ } else
+ mpz_2fac_ui (x, n);
+ }
+ } else { /* m >= 3, gcd(n,m) = 1 */
+ mp_limb_t *factors;
+ mp_limb_t prod, max_prod, j;
+ TMP_DECL;
+
+ sn = n / m + 1;
+
+ j = 0;
+ prod = n;
+ n -= m;
+ max_prod = GMP_NUMB_MAX / n;
+
+ TMP_MARK;
+ factors = TMP_SALLOC_LIMBS (sn / log_n_max (n) + 1);
+
+ for (; n > m; n -= m)
+ FACTOR_LIST_STORE (n, prod, max_prod, factors, j);
+
+ factors[j++] = n;
+ factors[j++] = prod;
+
+ if (g > 1) {
+ mpz_init (t);
+ mpz_prodlimbs (t, factors, j);
+ } else
+ mpz_prodlimbs (x, factors, j);
+
+ TMP_FREE;
+ }
+
+ if (g > 1) {
+ mpz_t p;
+
+ mpz_init (p);
+ mpz_ui_pow_ui (p, g, sn); /* g^sn */
+ mpz_mul (x, p, t);
+ mpz_clear (p);
+ mpz_clear (t);
+ }
+ }
+}
diff -r 8703c179e3b3 -r f0e71d6941c6 tests/mpz/t-mfac_uiui.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/mpz/t-mfac_uiui.c Mon Apr 30 17:38:21 2012 +0200
@@ -0,0 +1,136 @@
+/* Exercise mpz_mfac_uiui.
+
+Copyright 2000, 2001, 2002, 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 <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "tests.h"
+
+
+/* Usage: t-mfac_uiui [x|num]
+
+ With no arguments testing goes up to the initial value of "limit" below.
+ With a number argument tests are carried that far, or with a literal "x"
+ tests are continued without limit (this being meant only for development
+ purposes). */
+
+#define MULTIFAC_WHEEL (2*3*11)
+#define MULTIFAC_WHEEL2 (5*13)
+
+int
+main (int argc, char *argv[])
+{
+ mpz_t ref[MULTIFAC_WHEEL], ref2[MULTIFAC_WHEEL2], res;
+ unsigned long n, j, m, m2;
+ unsigned long limit = 2222, step = 1;
+
+ tests_start ();
+
+ if (argc > 1 && argv[1][0] == 'x')
+ limit = ULONG_MAX;
+ else if (argc > 1)
+ limit = atoi (argv[1]);
+
+ /* for small limb testing */
+ limit = MIN (limit, MP_LIMB_T_MAX);
+
+ for (m = 0; m < MULTIFAC_WHEEL; m++)
+ mpz_init_set_ui(ref [m],1);
+ for (m2 = 0; m2 < MULTIFAC_WHEEL2; m2++)
+ mpz_init_set_ui(ref2 [m2],1);
+
+ mpz_init (res);
+
+ m = 0;
+ m2 = 0;
+ for (n = 0; n <= limit;)
+ {
+ mpz_mfac_uiui (res, n, MULTIFAC_WHEEL);
+ MPZ_CHECK_FORMAT (res);
+ if (mpz_cmp (ref[m], res) != 0)
+ {
+ printf ("mpz_mfac_uiui(%lu,&i) wrong\n", n, MULTIFAC_WHEEL);
+ printf (" got "); mpz_out_str (stdout, 10, res); printf("\n");
+ printf (" want "); mpz_out_str (stdout, 10, ref[m]); printf("\n");
+ abort ();
+ }
+ mpz_mfac_uiui (res, n, MULTIFAC_WHEEL2);
+ MPZ_CHECK_FORMAT (res);
+ if (mpz_cmp (ref2[m2], res) != 0)
+ {
+ printf ("mpz_mfac_uiui(%lu,&i) wrong\n", n, MULTIFAC_WHEEL2);
+ printf (" got "); mpz_out_str (stdout, 10, res); printf("\n");
+ printf (" want "); mpz_out_str (stdout, 10, ref2[m2]); printf("\n");
+ abort ();
+ }
+ if (n + step <= limit)
+ for (j = 0; j < step; j++) {
+ n++; m++; m2++;
+ if (m >= MULTIFAC_WHEEL) m -= MULTIFAC_WHEEL;
+ if (m2 >= MULTIFAC_WHEEL2) m2 -= MULTIFAC_WHEEL2;
+ mpz_mul_ui (ref[m], ref[m], n); /* Compute a reference, with current library */
+ mpz_mul_ui (ref2[m2], ref2[m2], n); /* Compute a reference, with current library */
+ }
+ else n += step;
+ }
+ mpz_fac_ui (ref[0], n);
+ mpz_mfac_uiui (res, n, 1);
+ MPZ_CHECK_FORMAT (res);
+ if (mpz_cmp (ref[0], res) != 0)
+ {
+ printf ("mpz_mfac_uiui(%lu,1) wrong\n", n);
+ printf (" got "); mpz_out_str (stdout, 10, res); printf("\n");
+ printf (" want "); mpz_out_str (stdout, 10, ref[0]); printf("\n");
+ abort ();
+ }
More information about the gmp-commit
mailing list