diff options
| -rw-r--r-- | NEWS | 3 | ||||
| -rw-r--r-- | TODO | 3 | ||||
| -rw-r--r-- | configure.ac | 1 | ||||
| -rw-r--r-- | doc/coreutils.texi | 58 | ||||
| -rw-r--r-- | m4/gmp.m4 | 36 | ||||
| -rw-r--r-- | src/Makefile.am | 3 | ||||
| -rw-r--r-- | src/factor.c | 493 |
7 files changed, 532 insertions, 65 deletions
@@ -19,6 +19,9 @@ GNU coreutils NEWS -*- outline -*- With this new option, after a short read, dd repeatedly calls read, until it fills the incomplete block, reaches EOF, or encounters an error. + factor accepts arbitrarily large numbers and factors them using + Pollard's rho algorithm. + md5sum now accepts the new option, --quiet, to suppress the printing of 'OK' messages. sha1sum, sha224sum, sha384sum, and sha512sum accept it, too. @@ -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 ac93e1ca0..e54479fba 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 76b22e44d..c789a6cc4 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 -@command{factor} 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. +@table @samp +@item --help +Print a short help on standard output, then exit without further +processing. -@example -$ p=`echo '4294967279 * 4294967291'|bc` -$ factor $p -18446743979220271189: 4294967279 4294967291 -@end example +@item --bignum +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: +@item --no-bignum +Forces the use of native operations instead of GNU MP. This causes +@command{factor} to fail for larger inputs. -@example -$ factor 18446744073709551557 - 18446744073709551557: 18446744073709551557 -@end example +@item --version +Print the program version on standard output, then exit without further +processing. +@end 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 000000000..3e6033d47 --- /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 65b20a2a3..f464a70a7 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 08fa2a5c1..c7cd29eb7 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 <phr@ocf.berkeley.edu>. - 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,273 @@ /* 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 /* --no-bignum */ + } 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 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 int 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; + static unsigned int const add[] = {4, 2, 4, 2, 4, 6, 2, 6}; + unsigned int const *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; + size_t n_bits; + 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 to the size of t. */ + n_bits = mpz_sizeinbase (t, 2); + division_limit = MIN (n_bits, 1000); + 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 +332,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 +378,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 +397,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 +405,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 = 0; i < nfactors_found; ++i) + { + faclist[i] = &factor[i]; + } + qsort (faclist, nfactors_found, sizeof *faclist, mpcompare); + + for (i = 0; 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 = 0; 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}, + {"no-bignum", 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 always use arbitrary-precision arithmetic\n\ + --no-bignum always use 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 +574,7 @@ int main (int argc, char **argv) { bool ok; + int c; initialize_main (&argc, &argv); set_program_name (argv[0]); @@ -195,10 +584,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 +619,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); } |