[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[PATCH 1/2] Support arbitrarily-long numbers with GNU MP.
From: |
James Youngman |
Subject: |
[PATCH 1/2] Support arbitrarily-long numbers with GNU MP. |
Date: |
Sun, 27 Jul 2008 22:25:42 +0100 |
2008-07-27 James Youngman <address@hidden>
Support arbitrarily-long numbers with GNU MP.
* m4/gmp.m4: New file; adds cu_GMP, which detects GNU MP.
* configure.ac: Use cu_GMP.
* src/Makefile.am: Link factor against libgmp if available.
* src/factor.c: Use GNU MP if it is available.
(emit_factor, emit_ul_factor, factor_using_division,
factor_using_pollard_rho, extract_factors_multi,
sort_and_print_factors, free_factors): new functions
for the arbitrary-precision implementation, taken from an example
in GNU MP.
(factor_wheel): Renamed; was called factor.
(print_factors_single): Renamed; was called print_factors.
(print_factors): New function, chooses between the single- and
arbitrary-precision algorithms according to availability of GNU MP
and the length of the number to be factored.
(usage, main): New options --use-mp and --nouse-mp.
---
configure.ac | 1 +
m4/gmp.m4 | 20 +++
src/Makefile.am | 3 +
src/factor.c | 497 ++++++++++++++++++++++++++++++++++++++++++++++++++-----
4 files changed, 483 insertions(+), 38 deletions(-)
create mode 100644 m4/gmp.m4
diff --git a/configure.ac b/configure.ac
index ac93e1c..e54479f 100644
--- a/configure.ac
+++ b/configure.ac
@@ -244,6 +244,7 @@ AC_CHECK_DECLS([strsignal, sys_siglist, _sys_siglist,
__sys_siglist], , ,
#include <signal.h>])
cu_LIB_CHECK
+cu_GMP
# Build df only if there's a point to it.
if test $gl_cv_list_mounted_fs = yes && test $gl_cv_fs_space = yes; then
diff --git a/m4/gmp.m4 b/m4/gmp.m4
new file mode 100644
index 0000000..c0d1091
--- /dev/null
+++ b/m4/gmp.m4
@@ -0,0 +1,20 @@
+# Tests for GNU GMP (or any compatible replacement).
+
+dnl Copyright (C) 2008 Free Software Foundation, Inc.
+
+dnl This file is free software; the Free Software Foundation
+dnl gives unlimited permission to copy and/or distribute it,
+dnl with or without modifications, as long as this notice is preserved.
+
+dnl Written by James Youngman.
+
+dnl Check for libgmp. We avoid use of AC_CHECK_LIBS because we don't want to
+dnl add this to $LIBS for all targets.
+AC_DEFUN([cu_GMP],
+[
+ AC_CHECK_LIB(gmp, __gmpz_init,
+ [
+ AC_DEFINE(HAVE_GMP,1,[Define if you have GNU libgmp (or
replacement)])
+ AC_SUBST(LIB_GMP,[-lgmp])
+ ],)
+])
diff --git a/src/Makefile.am b/src/Makefile.am
index 65b20a2..f464a70 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -129,6 +129,9 @@ seq_LDADD = $(LDADD) $(POW_LIB)
# and the `nanosleep' reference in lib/xnanosleep.c.
nanosec_libs = $(LDADD) $(POW_LIB) $(LIB_NANOSLEEP)
+# for various GMP functions
+factor_LDADD = $(LDADD) $(LIB_GMP)
+
sleep_LDADD = $(nanosec_libs)
tail_LDADD = $(nanosec_libs)
diff --git a/src/factor.c b/src/factor.c
index 08fa2a5..9000424 100644
--- a/src/factor.c
+++ b/src/factor.c
@@ -15,12 +15,19 @@
along with this program. If not, see <http://www.gnu.org/licenses/>. */
/* Written by Paul Rubin <address@hidden>.
- Adapted for GNU, fixed to factor UINT_MAX by Jim Meyering. */
+ Adapted for GNU, fixed to factor UINT_MAX by Jim Meyering.
+ Arbitrary-precision code adapted by James Youngman from factorize.c
+ of GNU MP, version 4.2.2.
+*/
#include <config.h>
#include <getopt.h>
#include <stdio.h>
#include <sys/types.h>
+#if HAVE_GMP
+#include <gmp.h>
+#endif
+
#include <assert.h>
#define NDEBUG 1
@@ -40,6 +47,278 @@
/* Token delimiters when reading from a file. */
#define DELIM "\n\t "
+static bool verbose = false;
+
+typedef enum
+ {
+ ALGO_AUTOSELECT, /* default */
+ ALGO_GMP, /* --use-mp */
+ ALGO_SINGLE /* --nouse-mp */
+ } ALGORITHM_CHOICE;
+static ALGORITHM_CHOICE algorithm = ALGO_AUTOSELECT;
+
+
+
+#if HAVE_GMP
+static mpz_t *factors = NULL;
+static size_t nfactors_found = 0u;
+static size_t nfactors_allocated = 0u;
+
+static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
+
+static void
+emit_factor (mpz_t n)
+{
+ if (nfactors_found == nfactors_allocated)
+ factors = x2nrealloc (factors, &nfactors_allocated, sizeof(*factors));
+ mpz_init (factors[nfactors_found]);
+ mpz_set (factors[nfactors_found], n);
+ ++nfactors_found;
+}
+
+static void
+emit_ul_factor (unsigned long i)
+{
+ mpz_t t;
+ mpz_init (t);
+ mpz_set_ui (t, i);
+ emit_factor (t);
+}
+
+
+static void
+factor_using_division (mpz_t t, unsigned int limit)
+{
+ mpz_t q, r;
+ unsigned long int f;
+ int ai;
+ unsigned *addv = add;
+ unsigned int failures;
+
+ if (verbose)
+ {
+ printf ("[trial division (%u)] ", limit);
+ fflush (stdout);
+ }
+
+ mpz_init (q);
+ mpz_init (r);
+
+ f = mpz_scan1 (t, 0);
+ mpz_div_2exp (t, t, f);
+ while (f)
+ {
+ emit_ul_factor (2);
+ --f;
+ }
+
+ for (;;)
+ {
+ mpz_tdiv_qr_ui (q, r, t, 3);
+ if (mpz_cmp_ui (r, 0) != 0)
+ break;
+ mpz_set (t, q);
+ emit_ul_factor (3);
+ }
+
+ for (;;)
+ {
+ mpz_tdiv_qr_ui (q, r, t, 5);
+ if (mpz_cmp_ui (r, 0) != 0)
+ break;
+ mpz_set (t, q);
+ emit_ul_factor (5);
+ }
+
+ failures = 0;
+ f = 7;
+ ai = 0;
+ while (mpz_cmp_ui (t, 1) != 0)
+ {
+ mpz_tdiv_qr_ui (q, r, t, f);
+ if (mpz_cmp_ui (r, 0) != 0)
+ {
+ f += addv[ai];
+ if (mpz_cmp_ui (q, f) < 0)
+ break;
+ ai = (ai + 1) & 7;
+ failures++;
+ if (failures > limit)
+ break;
+ }
+ else
+ {
+ mpz_swap (t, q);
+ emit_ul_factor (f);
+ failures = 0;
+ }
+ }
+
+ mpz_clear (q);
+ mpz_clear (r);
+}
+
+static void
+factor_using_pollard_rho (mpz_t n, int a_int)
+{
+ mpz_t x, x1, y, P;
+ mpz_t a;
+ mpz_t g;
+ mpz_t t1, t2;
+ int k, l, c, i;
+
+ if (verbose)
+ {
+ printf ("[pollard-rho (%d)] ", a_int);
+ fflush (stdout);
+ }
+
+ mpz_init (g);
+ mpz_init (t1);
+ mpz_init (t2);
+
+ mpz_init_set_si (a, a_int);
+ mpz_init_set_si (y, 2);
+ mpz_init_set_si (x, 2);
+ mpz_init_set_si (x1, 2);
+ k = 1;
+ l = 1;
+ mpz_init_set_ui (P, 1);
+ c = 0;
+
+ while (mpz_cmp_ui (n, 1) != 0)
+ {
+S2:
+ mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
+
+ mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
+ c++;
+ if (c == 20)
+ {
+ c = 0;
+ mpz_gcd (g, P, n);
+ if (mpz_cmp_ui (g, 1) != 0)
+ goto S4;
+ mpz_set (y, x);
+ }
+S3:
+ k--;
+ if (k > 0)
+ goto S2;
+
+ mpz_gcd (g, P, n);
+ if (mpz_cmp_ui (g, 1) != 0)
+ goto S4;
+
+ mpz_set (x1, x);
+ k = l;
+ l = 2 * l;
+ for (i = 0; i < k; i++)
+ {
+ mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
+ }
+ mpz_set (y, x);
+ c = 0;
+ goto S2;
+S4:
+ do
+ {
+ mpz_mul (y, y, y); mpz_add (y, y, a); mpz_mod (y, y, n);
+ mpz_sub (t1, x1, y); mpz_gcd (g, t1, n);
+ }
+ while (mpz_cmp_ui (g, 1) == 0);
+
+ mpz_div (n, n, g); /* divide by g, before g is overwritten */
+
+ if (!mpz_probab_prime_p (g, 3))
+ {
+ do
+ {
+ mp_limb_t a_limb;
+ mpn_random (&a_limb, (mp_size_t) 1);
+ a_int = (int) a_limb;
+ }
+ while (a_int == -2 || a_int == 0);
+
+ if (verbose)
+ {
+ printf ("[composite factor--restarting pollard-rho] ");
+ fflush (stdout);
+ }
+ factor_using_pollard_rho (g, a_int);
+ }
+ else
+ {
+ emit_factor (g);
+ }
+ mpz_mod (x, x, n);
+ mpz_mod (x1, x1, n);
+ mpz_mod (y, y, n);
+ if (mpz_probab_prime_p (n, 3))
+ {
+ emit_factor (n);
+ break;
+ }
+ }
+
+ mpz_clear (g);
+ mpz_clear (P);
+ mpz_clear (t2);
+ mpz_clear (t1);
+ mpz_clear (a);
+ mpz_clear (x1);
+ mpz_clear (x);
+ mpz_clear (y);
+}
+
+/* Arbitrary-precision factoring */
+static bool
+extract_factors_multi (const char *s)
+{
+ unsigned int division_limit;
+ mpz_t t;
+
+ mpz_init (t);
+ if (1 != gmp_sscanf (s, "%Zd", t))
+ {
+ error (0, 0, _("%s is not a valid positive integer"), quote (s));
+ return false;
+ }
+
+ printf ("%s:", s);
+
+ if (mpz_sgn (t) == 0)
+ return true; /* 0; no factors */
+
+ /* Set the trial division limit according the size of t. */
+ division_limit = mpz_sizeinbase (t, 2);
+ if (division_limit > 1000)
+ division_limit = 1000 * 1000;
+ else
+ division_limit = division_limit * division_limit;
+
+ factor_using_division (t, division_limit);
+
+ if (mpz_cmp_ui (t, 1) != 0)
+ {
+ if (verbose)
+ {
+ printf ("[is number prime?] ");
+ fflush (stdout);
+ }
+ if (mpz_probab_prime_p (t, 3))
+ {
+ emit_factor (t);
+ }
+ else
+ {
+ factor_using_pollard_rho (t, 1);
+ }
+ }
+ return true;
+}
+#endif
+
/* The maximum number of factors, including -1, for negative numbers. */
#define MAX_N_FACTORS (sizeof (uintmax_t) * CHAR_BIT)
@@ -58,39 +337,10 @@ static const unsigned char wheel_tab[] =
#define WHEEL_START (wheel_tab + WHEEL_SIZE)
#define WHEEL_END (wheel_tab + (sizeof wheel_tab / sizeof wheel_tab[0]))
-void
-usage (int status)
-{
- if (status != EXIT_SUCCESS)
- fprintf (stderr, _("Try `%s --help' for more information.\n"),
- program_name);
- else
- {
- printf (_("\
-Usage: %s [NUMBER]...\n\
- or: %s OPTION\n\
-"),
- program_name, program_name);
- fputs (_("\
-Print the prime factors of each NUMBER.\n\
-\n\
-"), stdout);
- fputs (HELP_OPTION_DESCRIPTION, stdout);
- fputs (VERSION_OPTION_DESCRIPTION, stdout);
- fputs (_("\
-\n\
-Print the prime factors of all specified integer NUMBERs. If no arguments\n\
-are specified on the command line, they are read from standard input.\n\
-"), stdout);
- emit_bug_reporting_address ();
- }
- exit (status);
-}
-
/* FIXME: comment */
static size_t
-factor (uintmax_t n0, size_t max_n_factors, uintmax_t *factors)
+factor_wheel (uintmax_t n0, size_t max_n_factors, uintmax_t *factors)
{
uintmax_t n = n0, d, q;
size_t n_factors = 0;
@@ -133,10 +383,9 @@ factor (uintmax_t n0, size_t max_n_factors, uintmax_t
*factors)
return n_factors;
}
-/* FIXME: comment */
-
+/* Single-precision factoring */
static bool
-print_factors (const char *s)
+print_factors_single (const char *s)
{
uintmax_t factors[MAX_N_FACTORS];
uintmax_t n;
@@ -153,7 +402,7 @@ print_factors (const char *s)
error (0, 0, _("%s is not a valid positive integer"), quote (s));
return false;
}
- n_factors = factor (n, MAX_N_FACTORS, factors);
+ n_factors = factor_wheel (n, MAX_N_FACTORS, factors);
printf ("%s:", umaxtostr (n, buf));
for (i = 0; i < n_factors; i++)
printf (" %s", umaxtostr (factors[i], buf));
@@ -161,6 +410,155 @@ print_factors (const char *s)
return true;
}
+#if HAVE_GMP
+static int
+mpcompare (const void* a, const void *b)
+{
+ return mpz_cmp(**(const mpz_t**)a, **(const mpz_t**)b);
+}
+
+static void
+sort_and_print_factors (void)
+{
+ mpz_t** faclist;
+ size_t i;
+
+ faclist = xcalloc (nfactors_found, sizeof(mpz_t*));
+ for (i=0u; i<nfactors_found; ++i)
+ {
+ faclist[i] = &factors[i];
+ }
+ qsort (faclist, nfactors_found, sizeof(*faclist), mpcompare);
+
+ for (i=0u; i<nfactors_found; ++i)
+ {
+ fputc (' ', stdout);
+ mpz_out_str (stdout, 10, *faclist[i]);
+ }
+ putchar ('\n');
+ free (faclist);
+}
+
+static void
+free_factors (void)
+{
+ size_t i;
+
+ for (i=0u; i<nfactors_found; ++i)
+ {
+ mpz_clear (factors[i]);
+ }
+ /* Don't actually free factors[] because in the case where we are
+ * reasding numbers from stdin, we may be about to use it again.
+ */
+ nfactors_found = 0;
+}
+
+
+/* Emit the factors of the indicated number. If we have the option of using
+ either algorithm, we select on the basis of the length of the number.
+ For longer numbers, we prefer the MP algorithm even if the native algorithm
+ has enough digits, because the algorithm is better. The turnover point
+ depends on the value as well as the length, but since we don't already know
+ if the number presented has small factors, we just switch over at 6 digits.
+*/
+static bool
+print_factors (const char *s)
+{
+ if (ALGO_AUTOSELECT == algorithm)
+ {
+ const size_t digits = strlen(s); /* upper limit on number of digits */
+ algorithm = digits < 5 ? ALGO_SINGLE : ALGO_GMP;
+ }
+ if (ALGO_GMP == algorithm)
+ {
+ if (verbose)
+ printf ("[%s]", _("using arbitrary-precision arithmetic"));
+ if (extract_factors_multi (s))
+ {
+ sort_and_print_factors ();
+ free_factors ();
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+ }
+ else
+ {
+ if (verbose)
+ printf ("[%s]", _("using single-precision arithmetic"));
+ return print_factors_single (s);
+ }
+}
+#else
+static bool
+print_factors (const char *s) /* Single-precision only. */
+{
+ if (ALGO_GMP == algorithm)
+ {
+ error (0, 0, _("arbitrary-precision arithmetic is not available"));
+ return false;
+ }
+ return print_factors_single (s);
+}
+#endif
+
+
+
+
+
+enum
+{
+ VERBOSE_OPTION = CHAR_MAX + 1,
+ USE_MP_OPTION,
+ NO_USE_MP_OPTION
+};
+
+static struct option const long_options[] =
+{
+ {"verbose", no_argument, NULL, VERBOSE_OPTION},
+ {"use-mp", no_argument, NULL, USE_MP_OPTION},
+ {"nouse-mp", no_argument, NULL, NO_USE_MP_OPTION},
+ {GETOPT_HELP_OPTION_DECL},
+ {GETOPT_VERSION_OPTION_DECL},
+ {NULL, 0, NULL, 0}
+};
+
+void
+usage (int status)
+{
+ if (status != EXIT_SUCCESS)
+ fprintf (stderr, _("Try `%s --help' for more information.\n"),
+ program_name);
+ else
+ {
+ printf (_("\
+Usage: %s [NUMBER]...\n\
+ or: %s OPTION\n\
+"),
+ program_name, program_name);
+ fputs (_("\
+Print the prime factors of each NUMBER.\n\
+\n\
+"), stdout);
+ fputs (_("\
+ --use-mp force the use of arbitrary-precision arithmetic\n\
+ --nouse-mp force the use of single-precision arithmetic\n"),
+ stdout);
+ fputs (HELP_OPTION_DESCRIPTION, stdout);
+ fputs (VERSION_OPTION_DESCRIPTION, stdout);
+ fputs (_("\
+\n\
+Print the prime factors of all specified integer NUMBERs. If no arguments\n\
+are specified on the command line, they are read from standard input.\n\
+"), stdout);
+ emit_bug_reporting_address ();
+ }
+ exit (status);
+}
+
static bool
do_stdin (void)
{
@@ -186,6 +584,7 @@ int
main (int argc, char **argv)
{
bool ok;
+ int c;
initialize_main (&argc, &argv);
set_program_name (argv[0]);
@@ -197,8 +596,26 @@ main (int argc, char **argv)
parse_long_options (argc, argv, PROGRAM_NAME, PACKAGE_NAME, VERSION,
usage, AUTHORS, (char const *) NULL);
- if (getopt_long (argc, argv, "", NULL, NULL) != -1)
- usage (EXIT_FAILURE);
+ while ((c = getopt_long (argc, argv, "", long_options, NULL)) != -1)
+ {
+ switch (c)
+ {
+ case VERBOSE_OPTION:
+ verbose = 1;
+ break;
+
+ case USE_MP_OPTION:
+ algorithm = ALGO_GMP;
+ break;
+
+ case NO_USE_MP_OPTION:
+ algorithm = ALGO_SINGLE;
+ break;
+
+ default:
+ usage (EXIT_FAILURE);
+ }
+ }
if (argc <= optind)
ok = do_stdin ();
@@ -210,6 +627,10 @@ main (int argc, char **argv)
if (! print_factors (argv[i]))
ok = false;
}
-
+#if HAVE_GMP
+ free (factors);
+ nfactors_allocated = 0;
+ factors = NULL;
+#endif
exit (ok ? EXIT_SUCCESS : EXIT_FAILURE);
}
--
1.5.6.3
- [PATCH 1/2] Support arbitrarily-long numbers with GNU MP.,
James Youngman <=