1 /* factor -- print prime factors of n.
2 Copyright (C) 86, 1995-2005, 2007-2008 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
17 /* Written by Paul Rubin <phr@ocf.berkeley.edu>.
18 Adapted for GNU, fixed to factor UINT_MAX by Jim Meyering.
19 Arbitrary-precision code adapted by James Youngman from Torbjörn
20 Granlund's factorize.c, from GNU MP version 4.2.2.
27 #include <sys/types.h>
39 #include "readtokens.h"
42 /* The official name of this program (e.g., no `g' prefix). */
43 #define PROGRAM_NAME "factor"
45 #define AUTHORS proper_name ("Paul Rubin")
47 /* Token delimiters when reading from a file. */
50 static bool verbose = false;
54 ALGO_AUTOSELECT, /* default */
55 ALGO_GMP, /* --bignum */
56 ALGO_SINGLE /* --no-bignum */
58 static ALGORITHM_CHOICE algorithm = ALGO_AUTOSELECT;
62 static mpz_t *factor = NULL;
63 static size_t nfactors_found = 0;
64 static size_t nfactors_allocated = 0;
67 debug (char const *fmt, ...)
73 vfprintf (stderr, fmt, ap);
80 if (nfactors_found == nfactors_allocated)
81 factor = x2nrealloc (factor, &nfactors_allocated, sizeof *factor);
82 mpz_init (factor[nfactors_found]);
83 mpz_set (factor[nfactors_found], n);
88 emit_ul_factor (unsigned long int i)
97 factor_using_division (mpz_t t, unsigned int limit)
102 static unsigned int const add[] = {4, 2, 4, 2, 4, 6, 2, 6};
103 unsigned int const *addv = add;
104 unsigned int failures;
106 debug ("[trial division (%u)] ", limit);
111 f = mpz_scan1 (t, 0);
112 mpz_div_2exp (t, t, f);
121 mpz_tdiv_qr_ui (q, r, t, 3);
122 if (mpz_cmp_ui (r, 0) != 0)
130 mpz_tdiv_qr_ui (q, r, t, 5);
131 if (mpz_cmp_ui (r, 0) != 0)
140 while (mpz_cmp_ui (t, 1) != 0)
142 mpz_tdiv_qr_ui (q, r, t, f);
143 if (mpz_cmp_ui (r, 0) != 0)
146 if (mpz_cmp_ui (q, f) < 0)
150 if (failures > limit)
166 factor_using_pollard_rho (mpz_t n, int a_int)
174 debug ("[pollard-rho (%d)] ", a_int);
180 mpz_init_set_si (a, a_int);
181 mpz_init_set_si (y, 2);
182 mpz_init_set_si (x, 2);
183 mpz_init_set_si (x1, 2);
186 mpz_init_set_ui (P, 1);
189 while (mpz_cmp_ui (n, 1) != 0)
192 mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
194 mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
200 if (mpz_cmp_ui (g, 1) != 0)
210 if (mpz_cmp_ui (g, 1) != 0)
216 for (i = 0; i < k; i++)
218 mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
226 mpz_mul (y, y, y); mpz_add (y, y, a); mpz_mod (y, y, n);
227 mpz_sub (t1, x1, y); mpz_gcd (g, t1, n);
229 while (mpz_cmp_ui (g, 1) == 0);
231 mpz_div (n, n, g); /* divide by g, before g is overwritten */
233 if (!mpz_probab_prime_p (g, 3))
238 mpn_random (&a_limb, (mp_size_t) 1);
239 a_int = (int) a_limb;
241 while (a_int == -2 || a_int == 0);
243 debug ("[composite factor--restarting pollard-rho] ");
244 factor_using_pollard_rho (g, a_int);
253 if (mpz_probab_prime_p (n, 3))
270 /* Arbitrary-precision factoring */
272 extract_factors_multi (char const *s)
274 unsigned int division_limit;
279 if (1 != gmp_sscanf (s, "%Zd", t))
281 error (0, 0, _("%s is not a valid positive integer"), quote (s));
287 if (mpz_sgn (t) == 0)
290 return true; /* 0; no factors */
293 /* Set the trial division limit according to the size of t. */
294 n_bits = mpz_sizeinbase (t, 2);
295 division_limit = MIN (n_bits, 1000);
296 division_limit *= division_limit;
298 factor_using_division (t, division_limit);
300 if (mpz_cmp_ui (t, 1) != 0)
302 debug ("[is number prime?] ");
303 if (mpz_probab_prime_p (t, 3))
309 factor_using_pollard_rho (t, 1);
317 /* The maximum number of factors, including -1, for negative numbers. */
318 #define MAX_N_FACTORS (sizeof (uintmax_t) * CHAR_BIT)
320 /* The trial divisor increment wheel. Use it to skip over divisors that
321 are composites of 2, 3, 5, 7, or 11. The part from WHEEL_START up to
322 WHEEL_END is reused periodically, while the "lead in" is used to test
323 for those primes and to jump onto the wheel. For more information, see
324 http://www.utm.edu/research/primes/glossary/WheelFactorization.html */
326 #include "wheel-size.h" /* For the definition of WHEEL_SIZE. */
327 static const unsigned char wheel_tab[] =
332 #define WHEEL_START (wheel_tab + WHEEL_SIZE)
333 #define WHEEL_END (wheel_tab + (sizeof wheel_tab / sizeof wheel_tab[0]))
338 factor_wheel (uintmax_t n0, size_t max_n_factors, uintmax_t *factors)
340 uintmax_t n = n0, d, q;
341 size_t n_factors = 0;
342 unsigned char const *w = wheel_tab;
347 /* The exit condition in the following loop is correct because
348 any time it is tested one of these 3 conditions holds:
351 (3) n is composite but has no factors less than d.
352 If (1) or (2) obviously the right thing happens.
353 If (3), then since n is composite it is >= d^2. */
361 assert (n_factors < max_n_factors);
362 factors[n_factors++] = d;
372 if (n != 1 || n0 == 1)
374 assert (n_factors < max_n_factors);
375 factors[n_factors++] = n;
381 /* Single-precision factoring */
383 print_factors_single (char const *s)
385 uintmax_t factors[MAX_N_FACTORS];
389 char buf[INT_BUFSIZE_BOUND (uintmax_t)];
392 if ((err = xstrtoumax (s, NULL, 10, &n, "")) != LONGINT_OK)
394 if (err == LONGINT_OVERFLOW)
395 error (0, 0, _("%s is too large"), quote (s));
397 error (0, 0, _("%s is not a valid positive integer"), quote (s));
400 n_factors = factor_wheel (n, MAX_N_FACTORS, factors);
401 printf ("%s:", umaxtostr (n, buf));
402 for (i = 0; i < n_factors; i++)
403 printf (" %s", umaxtostr (factors[i], buf));
410 mpcompare (const void *av, const void *bv)
412 mpz_t *const *a = av;
413 mpz_t *const *b = bv;
414 return mpz_cmp (**a, **b);
418 sort_and_print_factors (void)
423 faclist = xcalloc (nfactors_found, sizeof *faclist);
424 for (i = 0; i < nfactors_found; ++i)
426 faclist[i] = &factor[i];
428 qsort (faclist, nfactors_found, sizeof *faclist, mpcompare);
430 for (i = 0; i < nfactors_found; ++i)
433 mpz_out_str (stdout, 10, *faclist[i]);
444 for (i = 0; i < nfactors_found; ++i)
446 mpz_clear (factor[i]);
448 /* Don't actually free factor[] because in the case where we are
449 reading numbers from stdin, we may be about to use it again. */
454 /* Emit the factors of the indicated number. If we have the option of using
455 either algorithm, we select on the basis of the length of the number.
456 For longer numbers, we prefer the MP algorithm even if the native algorithm
457 has enough digits, because the algorithm is better. The turnover point
458 depends on the value as well as the length, but since we don't already know
459 if the number presented has small factors, we just switch over at 6 digits.
462 print_factors (char const *s)
464 if (ALGO_AUTOSELECT == algorithm)
466 const size_t digits = strlen (s); /* upper limit on number of digits */
467 algorithm = digits < 6 ? ALGO_SINGLE : ALGO_GMP;
469 if (ALGO_GMP == algorithm)
471 debug ("[%s]", _("using arbitrary-precision arithmetic"));
472 if (extract_factors_multi (s))
474 sort_and_print_factors ();
485 debug ("[%s]", _("using single-precision arithmetic"));
486 return print_factors_single (s);
491 print_factors (const char *s) /* Single-precision only. */
493 if (ALGO_GMP == algorithm)
495 error (0, 0, _("arbitrary-precision arithmetic is not available"));
498 return print_factors_single (s);
504 VERBOSE_OPTION = CHAR_MAX + 1,
509 static struct option const long_options[] =
511 {"verbose", no_argument, NULL, VERBOSE_OPTION},
512 {"bignum", no_argument, NULL, USE_BIGNUM},
513 {"no-bignum", no_argument, NULL, NO_USE_BIGNUM},
514 {GETOPT_HELP_OPTION_DECL},
515 {GETOPT_VERSION_OPTION_DECL},
522 if (status != EXIT_SUCCESS)
523 fprintf (stderr, _("Try `%s --help' for more information.\n"),
528 Usage: %s [NUMBER]...\n\
531 program_name, program_name);
533 Print the prime factors of each NUMBER.\n\
537 --bignum always use arbitrary-precision arithmetic\n\
538 --no-bignum always use single-precision arithmetic\n"),
540 fputs (HELP_OPTION_DESCRIPTION, stdout);
541 fputs (VERSION_OPTION_DESCRIPTION, stdout);
544 Print the prime factors of all specified integer NUMBERs. If no arguments\n\
545 are specified on the command line, they are read from standard input.\n\
547 emit_bug_reporting_address ();
556 token_buffer tokenbuffer;
558 init_tokenbuffer (&tokenbuffer);
562 size_t token_length = readtoken (stdin, DELIM, sizeof (DELIM) - 1,
564 if (token_length == (size_t) -1)
566 ok &= print_factors (tokenbuffer.buffer);
568 free (tokenbuffer.buffer);
574 main (int argc, char **argv)
579 initialize_main (&argc, &argv);
580 set_program_name (argv[0]);
581 setlocale (LC_ALL, "");
582 bindtextdomain (PACKAGE, LOCALEDIR);
583 textdomain (PACKAGE);
585 atexit (close_stdout);
587 while ((c = getopt_long (argc, argv, "", long_options, NULL)) != -1)
596 algorithm = ALGO_GMP;
600 algorithm = ALGO_SINGLE;
603 case_GETOPT_HELP_CHAR;
605 case_GETOPT_VERSION_CHAR (PROGRAM_NAME, AUTHORS);
608 usage (EXIT_FAILURE);
618 for (i = optind; i < argc; i++)
619 if (! print_factors (argv[i]))
625 exit (ok ? EXIT_SUCCESS : EXIT_FAILURE);