2008-07-27 James Youngman <[EMAIL PROTECTED]>
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 --bignum and --nobignum.
* coreutils.texi (factor invocation): Document new command-line
options for the MP implementation and update the performance
numbers to take into account the asymptotically faster algorithm.
* TODO: Removed item about factoring large primes (it's done).
* m4/gmp.m4: Add support for --without-gmp.
---
TODO | 3 -
configure.ac | 1 +
doc/coreutils.texi | 58 ++++---
m4/gmp.m4 | 36 ++++
src/Makefile.am | 3 +
src/factor.c | 497 +++++++++++++++++++++++++++++++++++++++++++++++-----
6 files changed, 533 insertions(+), 65 deletions(-)
create mode 100644 m4/gmp.m4
diff --git a/TODO b/TODO
index e7f0508..c7bb820 100644
--- a/TODO
+++ b/TODO
@@ -127,9 +127,6 @@ Changes expected to go in, someday.
an implicit --NO-dereference-command-line-symlink-to-dir meaning.
Pointed out by Karl Berry.
- A more efficient version of factor, and possibly one that
- accepts inputs of size 2^64 and larger.
-
dd: consider adding an option to suppress `bytes/block read/written'
output to stderr. Suggested here:
http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=165045
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/doc/coreutils.texi b/doc/coreutils.texi
index 76b22e4..54798fd 100644
--- a/doc/coreutils.texi
+++ b/doc/coreutils.texi
@@ -14359,36 +14359,52 @@ factor @var{option}
If no @var{number} is specified on the command line, @command{factor} reads
numbers from standard input, delimited by newlines, tabs, or spaces.
-The only options are @option{--help} and @option{--version}. @xref{Common
-options}.
+The @command{factor} command supports only a small number of options:
-The algorithm it uses is not very sophisticated, so for some inputs
[EMAIL PROTECTED] runs for a long time. The hardest numbers to factor are
-the products of large primes. Factoring the product of the two largest 32-bit
-prime numbers takes about 80 seconds of CPU time on a 1.6 GHz Athlon.
[EMAIL PROTECTED] @samp
[EMAIL PROTECTED] --help
+Print a short help on standard output, then exit without further
+processing.
[EMAIL PROTECTED]
-$ p=`echo '4294967279 * 4294967291'|bc`
-$ factor $p
-18446743979220271189: 4294967279 4294967291
[EMAIL PROTECTED] example
[EMAIL PROTECTED] --use-mp
+Forces the use of the GNU MP library. By default, @command{factor}
+selects between using GNU MP and using native operations on the basis
+of the length of the number to be factored.
-Similarly, it takes about 80 seconds for GNU factor (from coreutils-5.1.2)
-to ``factor'' the largest 64-bit prime:
[EMAIL PROTECTED] --nouse-mp
+Forces the use of native operations instead of GNU MP. This causes
[EMAIL PROTECTED] to fail for larger inputs.
[EMAIL PROTECTED]
-$ factor 18446744073709551557
- 18446744073709551557: 18446744073709551557
[EMAIL PROTECTED] example
[EMAIL PROTECTED] --version
+Print the program version on standard output, then exit without further
+processing.
[EMAIL PROTECTED] table
-In contrast, @command{factor} factors the largest 64-bit number in just
-over a tenth of a second:
+Factoring the product of the eighth and ninth Mersenne primes
+takes about 30 milliseconds of CPU time on a 2.2 GHz Athlon.
@example
-$ factor `echo '2^64-1'|bc`
-18446744073709551615: 3 5 17 257 641 65537 6700417
+M8=`echo 2^31-1|bc` ; M9=`echo 2^61-1|bc`
+/usr/bin/time -f '%U' factor $(echo "$M8 * $M9" | bc)
+4951760154835678088235319297: 2147483647 2305843009213693951
+0.03
@end example
+Similarly, factoring the eighth Fermat number @math{2^{256}+1} takes
+about 20 seconds on the same machine.
+
+Factoring large prime numbers is, in general, hard. The Pollard Rho
+algorithm used by @command{factor} is particularly effective for
+numbers with relatively small factors. If you wish to factor large
+numbers which do not have small factors (for example, numbers which
+are the product of two large primes), other methods are far better.
+
+If @command{factor} is built without using GNU MP, only
+single-precision arithmetic is available, and so large numbers
+(typically @math{2^{64}} and above) will not be supported. The
single-precision
+code uses an algorithm which is designed for factoring smaller
+numbers.
+
@exitstatus
diff --git a/m4/gmp.m4 b/m4/gmp.m4
new file mode 100644
index 0000000..3e6033d
--- /dev/null
+++ b/m4/gmp.m4
@@ -0,0 +1,36 @@
+# 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],
+[
+ LIB_GMP=
+ AC_SUBST([LIB_GMP])
+
+ AC_ARG_WITH([gmp],
+ AS_HELP_STRING([--without-gmp],
+ [do not use the GNU MP library for arbitrary precision
+ calculation (default: use it if available)]),
+ [cu_use_gmp=$withval],
+ [cu_use_gmp=auto])
+
+ if test $cu_use_gmp != no; then
+ cu_saved_libs=$LIBS
+ AC_SEARCH_LIBS([__gmpz_init], [gmp],
+ [test "$ac_cv_search___gmpz_init" = "none required" ||
+ {
+ LIB_GMP=$ac_cv_search___gmpz_init
+ AC_DEFINE([HAVE_GMP], 1,
+ [Define if you have GNU libgmp (or replacement)])
+ }])
+ LIBS=$cu_saved_libs
+ fi
+])
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..707ccb0 100644
--- a/src/factor.c
+++ b/src/factor.c
@@ -15,19 +15,26 @@
along with this program. If not, see <http://www.gnu.org/licenses/>. */
/* Written by Paul Rubin <[EMAIL PROTECTED]>.
- 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 <stdarg.h>
#include <stdio.h>
#include <sys/types.h>
+#if HAVE_GMP
+#include <gmp.h>
+#endif
+
#include <assert.h>
#define NDEBUG 1
#include "system.h"
#include "error.h"
-#include "long-options.h"
#include "quote.h"
#include "readtokens.h"
#include "xstrtol.h"
@@ -40,6 +47,277 @@
/* Token delimiters when reading from a file. */
#define DELIM "\n\t "
+static bool verbose = false;
+
+typedef enum
+ {
+ ALGO_AUTOSELECT, /* default */
+ ALGO_GMP, /* --bignum */
+ ALGO_SINGLE /* --nobignum */
+ } ALGORITHM_CHOICE;
+static ALGORITHM_CHOICE algorithm = ALGO_AUTOSELECT;
+
+
+
+#if HAVE_GMP
+static mpz_t *factor = NULL;
+static size_t nfactors_found = 0;
+static size_t nfactors_allocated = 0;
+
+static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
+
+static void
+debug (char const *fmt, ...)
+{
+ if (verbose)
+ {
+ va_list ap;
+ va_start (ap, fmt);
+ vfprintf (stderr, fmt, ap);
+ }
+}
+
+
+static void
+emit_factor (mpz_t n)
+{
+ if (nfactors_found == nfactors_allocated)
+ factor = x2nrealloc (factor, &nfactors_allocated, sizeof *factor);
+ mpz_init (factor[nfactors_found]);
+ mpz_set (factor[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;
+
+ debug ("[trial division (%u)] ", limit);
+
+ 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;
+
+ debug ("[pollard-rho (%d)] ", a_int);
+
+ 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);
+ }
+
+ 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);
+
+ debug ("[composite factor--restarting pollard-rho] ");
+ 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 (char const *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)
+ {
+ mpz_clear (t);
+ 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)
+ {
+ debug ("[is number prime?] ");
+ if (mpz_probab_prime_p (t, 3))
+ {
+ emit_factor (t);
+ }
+ else
+ {
+ factor_using_pollard_rho (t, 1);
+ }
+ }
+ mpz_clear (t);
+ return true;
+}
+#endif
+
/* The maximum number of factors, including -1, for negative numbers. */
#define MAX_N_FACTORS (sizeof (uintmax_t) * CHAR_BIT)
@@ -58,39 +336,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 +382,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 (char const *s)
{
uintmax_t factors[MAX_N_FACTORS];
uintmax_t n;
@@ -153,7 +401,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 +409,150 @@ print_factors (const char *s)
return true;
}
+#if HAVE_GMP
+static int
+mpcompare (const void* av, const void *bv)
+{
+ mpz_t *const *a = av;
+ mpz_t *const *b = bv;
+ return mpz_cmp (**a, **b);
+}
+
+static void
+sort_and_print_factors (void)
+{
+ mpz_t** faclist;
+ size_t i;
+
+ faclist = xcalloc (nfactors_found, sizeof *faclist);
+ for (i=0u; i < nfactors_found; ++i)
+ {
+ faclist[i] = &factor[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 (factor[i]);
+ }
+ /* Don't actually free factor[] because in the case where we are
+ reading 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 (char const *s)
+{
+ if (ALGO_AUTOSELECT == algorithm)
+ {
+ const size_t digits = strlen(s); /* upper limit on number of digits */
+ algorithm = digits < 6 ? ALGO_SINGLE : ALGO_GMP;
+ }
+ if (ALGO_GMP == algorithm)
+ {
+ debug ("[%s]", _("using arbitrary-precision arithmetic"));
+ if (extract_factors_multi (s))
+ {
+ sort_and_print_factors ();
+ free_factors ();
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+ }
+ else
+ {
+ debug ("[%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_BIGNUM,
+ NO_USE_BIGNUM
+};
+
+static struct option const long_options[] =
+{
+ {"verbose", no_argument, NULL, VERBOSE_OPTION},
+ {"bignum", no_argument, NULL, USE_BIGNUM},
+ {"nobignum", no_argument, NULL, NO_USE_BIGNUM},
+ {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 (_("\
+ --bignum force the use of arbitrary-precision arithmetic\n\
+ --nobignum 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 +578,7 @@ int
main (int argc, char **argv)
{
bool ok;
+ int c;
initialize_main (&argc, &argv);
set_program_name (argv[0]);
@@ -195,10 +588,30 @@ main (int argc, char **argv)
atexit (close_stdout);
- 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 = true;
+ break;
+
+ case USE_BIGNUM:
+ algorithm = ALGO_GMP;
+ break;
+
+ case NO_USE_BIGNUM:
+ algorithm = ALGO_SINGLE;
+ break;
+
+ case_GETOPT_HELP_CHAR;
+
+ case_GETOPT_VERSION_CHAR (PROGRAM_NAME, AUTHORS);
+
+ default:
+ usage (EXIT_FAILURE);
+ }
+ }
if (argc <= optind)
ok = do_stdin ();
@@ -210,6 +623,8 @@ main (int argc, char **argv)
if (! print_factors (argv[i]))
ok = false;
}
-
+#if HAVE_GMP
+ free (factor);
+#endif
exit (ok ? EXIT_SUCCESS : EXIT_FAILURE);
}
--
1.5.6.3
_______________________________________________
Bug-coreutils mailing list
Bug-coreutils@gnu.org
http://lists.gnu.org/mailman/listinfo/bug-coreutils
