global: convert indentation-TABs to spaces
[coreutils.git] / src / factor.c
blob59fe1f45fd8145a0695f7728df785d5b411f3812
1 /* factor -- print prime factors of n.
2 Copyright (C) 86, 1995-2005, 2007-2009 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.
23 #include <config.h>
24 #include <getopt.h>
25 #include <stdarg.h>
26 #include <stdio.h>
27 #include <sys/types.h>
28 #if HAVE_GMP
29 #include <gmp.h>
30 #endif
32 #include <assert.h>
34 #include "system.h"
35 #include "error.h"
36 #include "quote.h"
37 #include "readtokens.h"
38 #include "xstrtol.h"
40 /* The official name of this program (e.g., no `g' prefix). */
41 #define PROGRAM_NAME "factor"
43 #define AUTHORS proper_name ("Paul Rubin")
45 /* Token delimiters when reading from a file. */
46 #define DELIM "\n\t "
48 static bool verbose = false;
50 #if HAVE_GMP
51 static mpz_t *factor = NULL;
52 static size_t nfactors_found = 0;
53 static size_t nfactors_allocated = 0;
55 static void
56 debug (char const *fmt, ...)
58 if (verbose)
60 va_list ap;
61 va_start (ap, fmt);
62 vfprintf (stderr, fmt, ap);
66 static void
67 emit_factor (mpz_t n)
69 if (nfactors_found == nfactors_allocated)
70 factor = x2nrealloc (factor, &nfactors_allocated, sizeof *factor);
71 mpz_init (factor[nfactors_found]);
72 mpz_set (factor[nfactors_found], n);
73 ++nfactors_found;
76 static void
77 emit_ul_factor (unsigned long int i)
79 mpz_t t;
80 mpz_init (t);
81 mpz_set_ui (t, i);
82 emit_factor (t);
83 mpz_clear (t);
86 static void
87 factor_using_division (mpz_t t, unsigned int limit)
89 mpz_t q, r;
90 unsigned long int f;
91 int ai;
92 static unsigned int const add[] = {4, 2, 4, 2, 4, 6, 2, 6};
93 unsigned int const *addv = add;
94 unsigned int failures;
96 debug ("[trial division (%u)] ", limit);
98 mpz_init (q);
99 mpz_init (r);
101 f = mpz_scan1 (t, 0);
102 mpz_div_2exp (t, t, f);
103 while (f)
105 emit_ul_factor (2);
106 --f;
109 for (;;)
111 mpz_tdiv_qr_ui (q, r, t, 3);
112 if (mpz_cmp_ui (r, 0) != 0)
113 break;
114 mpz_set (t, q);
115 emit_ul_factor (3);
118 for (;;)
120 mpz_tdiv_qr_ui (q, r, t, 5);
121 if (mpz_cmp_ui (r, 0) != 0)
122 break;
123 mpz_set (t, q);
124 emit_ul_factor (5);
127 failures = 0;
128 f = 7;
129 ai = 0;
130 while (mpz_cmp_ui (t, 1) != 0)
132 mpz_tdiv_qr_ui (q, r, t, f);
133 if (mpz_cmp_ui (r, 0) != 0)
135 f += addv[ai];
136 if (mpz_cmp_ui (q, f) < 0)
137 break;
138 ai = (ai + 1) & 7;
139 failures++;
140 if (failures > limit)
141 break;
143 else
145 mpz_swap (t, q);
146 emit_ul_factor (f);
147 failures = 0;
151 mpz_clear (q);
152 mpz_clear (r);
155 static void
156 factor_using_pollard_rho (mpz_t n, int a_int)
158 mpz_t x, x1, y, P;
159 mpz_t a;
160 mpz_t g;
161 mpz_t t1, t2;
162 int k, l, c, i;
164 debug ("[pollard-rho (%d)] ", a_int);
166 mpz_init (g);
167 mpz_init (t1);
168 mpz_init (t2);
170 mpz_init_set_si (a, a_int);
171 mpz_init_set_si (y, 2);
172 mpz_init_set_si (x, 2);
173 mpz_init_set_si (x1, 2);
174 k = 1;
175 l = 1;
176 mpz_init_set_ui (P, 1);
177 c = 0;
179 while (mpz_cmp_ui (n, 1) != 0)
182 mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
184 mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
185 c++;
186 if (c == 20)
188 c = 0;
189 mpz_gcd (g, P, n);
190 if (mpz_cmp_ui (g, 1) != 0)
191 goto S4;
192 mpz_set (y, x);
195 k--;
196 if (k > 0)
197 goto S2;
199 mpz_gcd (g, P, n);
200 if (mpz_cmp_ui (g, 1) != 0)
201 goto S4;
203 mpz_set (x1, x);
204 k = l;
205 l = 2 * l;
206 for (i = 0; i < k; i++)
208 mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
210 mpz_set (y, x);
211 c = 0;
212 goto S2;
216 mpz_mul (y, y, y); mpz_add (y, y, a); mpz_mod (y, y, n);
217 mpz_sub (t1, x1, y); mpz_gcd (g, t1, n);
219 while (mpz_cmp_ui (g, 1) == 0);
221 mpz_div (n, n, g); /* divide by g, before g is overwritten */
223 if (!mpz_probab_prime_p (g, 3))
227 mp_limb_t a_limb;
228 mpn_random (&a_limb, (mp_size_t) 1);
229 a_int = (int) a_limb;
231 while (a_int == -2 || a_int == 0);
233 debug ("[composite factor--restarting pollard-rho] ");
234 factor_using_pollard_rho (g, a_int);
236 else
238 emit_factor (g);
240 mpz_mod (x, x, n);
241 mpz_mod (x1, x1, n);
242 mpz_mod (y, y, n);
243 if (mpz_probab_prime_p (n, 3))
245 emit_factor (n);
246 break;
250 mpz_clear (g);
251 mpz_clear (P);
252 mpz_clear (t2);
253 mpz_clear (t1);
254 mpz_clear (a);
255 mpz_clear (x1);
256 mpz_clear (x);
257 mpz_clear (y);
260 #else
262 static void
263 debug (char const *fmt ATTRIBUTE_UNUSED, ...)
267 #endif
269 /* The maximum number of factors, including -1, for negative numbers. */
270 #define MAX_N_FACTORS (sizeof (uintmax_t) * CHAR_BIT)
272 /* The trial divisor increment wheel. Use it to skip over divisors that
273 are composites of 2, 3, 5, 7, or 11. The part from WHEEL_START up to
274 WHEEL_END is reused periodically, while the "lead in" is used to test
275 for those primes and to jump onto the wheel. For more information, see
276 http://www.utm.edu/research/primes/glossary/WheelFactorization.html */
278 #include "wheel-size.h" /* For the definition of WHEEL_SIZE. */
279 static const unsigned char wheel_tab[] =
281 #include "wheel.h"
284 #define WHEEL_START (wheel_tab + WHEEL_SIZE)
285 #define WHEEL_END (wheel_tab + ARRAY_CARDINALITY (wheel_tab))
287 /* FIXME: comment */
289 static size_t
290 factor_wheel (uintmax_t n0, size_t max_n_factors, uintmax_t *factors)
292 uintmax_t n = n0, d, q;
293 size_t n_factors = 0;
294 unsigned char const *w = wheel_tab;
296 if (n <= 1)
297 return n_factors;
299 /* The exit condition in the following loop is correct because
300 any time it is tested one of these 3 conditions holds:
301 (1) d divides n
302 (2) n is prime
303 (3) n is composite but has no factors less than d.
304 If (1) or (2) obviously the right thing happens.
305 If (3), then since n is composite it is >= d^2. */
307 d = 2;
310 q = n / d;
311 while (n == q * d)
313 assert (n_factors < max_n_factors);
314 factors[n_factors++] = d;
315 n = q;
316 q = n / d;
318 d += *(w++);
319 if (w == WHEEL_END)
320 w = WHEEL_START;
322 while (d <= q);
324 if (n != 1 || n0 == 1)
326 assert (n_factors < max_n_factors);
327 factors[n_factors++] = n;
330 return n_factors;
333 /* Single-precision factoring */
334 static void
335 print_factors_single (uintmax_t n)
337 uintmax_t factors[MAX_N_FACTORS];
338 size_t n_factors = factor_wheel (n, MAX_N_FACTORS, factors);
339 size_t i;
340 char buf[INT_BUFSIZE_BOUND (uintmax_t)];
342 printf ("%s:", umaxtostr (n, buf));
343 for (i = 0; i < n_factors; i++)
344 printf (" %s", umaxtostr (factors[i], buf));
345 putchar ('\n');
348 #if HAVE_GMP
349 static int
350 mpcompare (const void *av, const void *bv)
352 mpz_t *const *a = av;
353 mpz_t *const *b = bv;
354 return mpz_cmp (**a, **b);
357 static void
358 sort_and_print_factors (void)
360 mpz_t **faclist;
361 size_t i;
363 faclist = xcalloc (nfactors_found, sizeof *faclist);
364 for (i = 0; i < nfactors_found; ++i)
366 faclist[i] = &factor[i];
368 qsort (faclist, nfactors_found, sizeof *faclist, mpcompare);
370 for (i = 0; i < nfactors_found; ++i)
372 fputc (' ', stdout);
373 mpz_out_str (stdout, 10, *faclist[i]);
375 putchar ('\n');
376 free (faclist);
379 static void
380 free_factors (void)
382 size_t i;
384 for (i = 0; i < nfactors_found; ++i)
386 mpz_clear (factor[i]);
388 /* Don't actually free factor[] because in the case where we are
389 reading numbers from stdin, we may be about to use it again. */
390 nfactors_found = 0;
393 /* Arbitrary-precision factoring */
394 static void
395 print_factors_multi (mpz_t t)
397 mpz_out_str (stdout, 10, t);
398 putchar (':');
400 if (mpz_sgn (t) != 0)
402 /* Set the trial division limit according to the size of t. */
403 size_t n_bits = mpz_sizeinbase (t, 2);
404 unsigned int division_limit = MIN (n_bits, 1000);
405 division_limit *= division_limit;
407 factor_using_division (t, division_limit);
409 if (mpz_cmp_ui (t, 1) != 0)
411 debug ("[is number prime?] ");
412 if (mpz_probab_prime_p (t, 3))
413 emit_factor (t);
414 else
415 factor_using_pollard_rho (t, 1);
419 mpz_clear (t);
420 sort_and_print_factors ();
421 free_factors ();
423 #endif
426 /* Emit the factors of the indicated number. If we have the option of using
427 either algorithm, we select on the basis of the length of the number.
428 For longer numbers, we prefer the MP algorithm even if the native algorithm
429 has enough digits, because the algorithm is better. The turnover point
430 depends on the value. */
431 static bool
432 print_factors (char const *s)
434 uintmax_t n;
435 strtol_error err = xstrtoumax (s, NULL, 10, &n, "");
437 #if HAVE_GMP
438 enum { GMP_TURNOVER_POINT = 100000 };
440 if (err == LONGINT_OVERFLOW
441 || (err == LONGINT_OK && GMP_TURNOVER_POINT <= n))
443 mpz_t t;
444 mpz_init (t);
445 if (gmp_sscanf (s, "%Zd", t) == 1)
447 debug ("[%s]", _("using arbitrary-precision arithmetic"));
448 print_factors_multi (t);
449 return true;
451 err = LONGINT_INVALID;
453 #endif
455 switch (err)
457 case LONGINT_OK:
458 debug ("[%s]", _("using single-precision arithmetic"));
459 print_factors_single (n);
460 return true;
462 case LONGINT_OVERFLOW:
463 error (0, 0, _("%s is too large"), quote (s));
464 return false;
466 default:
467 error (0, 0, _("%s is not a valid positive integer"), quote (s));
468 return false;
472 enum
474 VERBOSE_OPTION = CHAR_MAX + 1
477 static struct option const long_options[] =
479 {"verbose", no_argument, NULL, VERBOSE_OPTION},
480 {GETOPT_HELP_OPTION_DECL},
481 {GETOPT_VERSION_OPTION_DECL},
482 {NULL, 0, NULL, 0}
485 void
486 usage (int status)
488 if (status != EXIT_SUCCESS)
489 fprintf (stderr, _("Try `%s --help' for more information.\n"),
490 program_name);
491 else
493 printf (_("\
494 Usage: %s [NUMBER]...\n\
495 or: %s OPTION\n\
497 program_name, program_name);
498 fputs (_("\
499 Print the prime factors of each specified integer NUMBER. If none\n\
500 are specified on the command line, read them from standard input.\n\
502 "), stdout);
503 fputs (HELP_OPTION_DESCRIPTION, stdout);
504 fputs (VERSION_OPTION_DESCRIPTION, stdout);
505 emit_bug_reporting_address ();
507 exit (status);
510 static bool
511 do_stdin (void)
513 bool ok = true;
514 token_buffer tokenbuffer;
516 init_tokenbuffer (&tokenbuffer);
518 for (;;)
520 size_t token_length = readtoken (stdin, DELIM, sizeof (DELIM) - 1,
521 &tokenbuffer);
522 if (token_length == (size_t) -1)
523 break;
524 ok &= print_factors (tokenbuffer.buffer);
526 free (tokenbuffer.buffer);
528 return ok;
532 main (int argc, char **argv)
534 bool ok;
535 int c;
537 initialize_main (&argc, &argv);
538 set_program_name (argv[0]);
539 setlocale (LC_ALL, "");
540 bindtextdomain (PACKAGE, LOCALEDIR);
541 textdomain (PACKAGE);
543 atexit (close_stdout);
545 while ((c = getopt_long (argc, argv, "", long_options, NULL)) != -1)
547 switch (c)
549 case VERBOSE_OPTION:
550 verbose = true;
551 break;
553 case_GETOPT_HELP_CHAR;
555 case_GETOPT_VERSION_CHAR (PROGRAM_NAME, AUTHORS);
557 default:
558 usage (EXIT_FAILURE);
562 if (argc <= optind)
563 ok = do_stdin ();
564 else
566 int i;
567 ok = true;
568 for (i = optind; i < argc; i++)
569 if (! print_factors (argv[i]))
570 ok = false;
572 #if HAVE_GMP
573 free (factor);
574 #endif
575 exit (ok ? EXIT_SUCCESS : EXIT_FAILURE);