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.
27 #include <sys/types.h>
37 #include "readtokens.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. */
48 static bool verbose
= false;
51 static mpz_t
*factor
= NULL
;
52 static size_t nfactors_found
= 0;
53 static size_t nfactors_allocated
= 0;
56 debug (char const *fmt
, ...)
62 vfprintf (stderr
, fmt
, ap
);
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
);
77 emit_ul_factor (unsigned long int i
)
87 factor_using_division (mpz_t t
, unsigned int limit
)
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
);
101 f
= mpz_scan1 (t
, 0);
102 mpz_div_2exp (t
, t
, f
);
111 mpz_tdiv_qr_ui (q
, r
, t
, 3);
112 if (mpz_cmp_ui (r
, 0) != 0)
120 mpz_tdiv_qr_ui (q
, r
, t
, 5);
121 if (mpz_cmp_ui (r
, 0) != 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)
136 if (mpz_cmp_ui (q
, f
) < 0)
140 if (failures
> limit
)
156 factor_using_pollard_rho (mpz_t n
, int a_int
)
164 debug ("[pollard-rho (%d)] ", a_int
);
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);
176 mpz_init_set_ui (P
, 1);
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
);
190 if (mpz_cmp_ui (g
, 1) != 0)
200 if (mpz_cmp_ui (g
, 1) != 0)
206 for (i
= 0; i
< k
; i
++)
208 mpz_mul (x
, x
, x
); mpz_add (x
, x
, a
); mpz_mod (x
, x
, n
);
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))
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
);
243 if (mpz_probab_prime_p (n
, 3))
263 debug (char const *fmt ATTRIBUTE_UNUSED
, ...)
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
[] =
284 #define WHEEL_START (wheel_tab + WHEEL_SIZE)
285 #define WHEEL_END (wheel_tab + (sizeof wheel_tab / sizeof wheel_tab[0]))
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
;
299 /* The exit condition in the following loop is correct because
300 any time it is tested one of these 3 conditions holds:
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. */
313 assert (n_factors
< max_n_factors
);
314 factors
[n_factors
++] = d
;
324 if (n
!= 1 || n0
== 1)
326 assert (n_factors
< max_n_factors
);
327 factors
[n_factors
++] = n
;
333 /* Single-precision factoring */
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
);
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
));
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
);
358 sort_and_print_factors (void)
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
)
373 mpz_out_str (stdout
, 10, *faclist
[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. */
393 /* Arbitrary-precision factoring */
395 print_factors_multi (mpz_t t
)
397 mpz_out_str (stdout
, 10, t
);
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))
415 factor_using_pollard_rho (t
, 1);
420 sort_and_print_factors ();
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. */
432 print_factors (char const *s
)
435 strtol_error err
= xstrtoumax (s
, NULL
, 10, &n
, "");
438 enum { GMP_TURNOVER_POINT
= 100000 };
440 if (err
== LONGINT_OVERFLOW
441 || (err
== LONGINT_OK
&& GMP_TURNOVER_POINT
<= n
))
445 if (gmp_sscanf (s
, "%Zd", t
) == 1)
447 debug ("[%s]", _("using arbitrary-precision arithmetic"));
448 print_factors_multi (t
);
451 err
= LONGINT_INVALID
;
458 debug ("[%s]", _("using single-precision arithmetic"));
459 print_factors_single (n
);
462 case LONGINT_OVERFLOW
:
463 error (0, 0, _("%s is too large"), quote (s
));
467 error (0, 0, _("%s is not a valid positive integer"), quote (s
));
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
},
488 if (status
!= EXIT_SUCCESS
)
489 fprintf (stderr
, _("Try `%s --help' for more information.\n"),
494 Usage: %s [NUMBER]...\n\
497 program_name
, program_name
);
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\
503 fputs (HELP_OPTION_DESCRIPTION
, stdout
);
504 fputs (VERSION_OPTION_DESCRIPTION
, stdout
);
505 emit_bug_reporting_address ();
514 token_buffer tokenbuffer
;
516 init_tokenbuffer (&tokenbuffer
);
520 size_t token_length
= readtoken (stdin
, DELIM
, sizeof (DELIM
) - 1,
522 if (token_length
== (size_t) -1)
524 ok
&= print_factors (tokenbuffer
.buffer
);
526 free (tokenbuffer
.buffer
);
532 main (int argc
, char **argv
)
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)
553 case_GETOPT_HELP_CHAR
;
555 case_GETOPT_VERSION_CHAR (PROGRAM_NAME
, AUTHORS
);
558 usage (EXIT_FAILURE
);
568 for (i
= optind
; i
< argc
; i
++)
569 if (! print_factors (argv
[i
]))
575 exit (ok
? EXIT_SUCCESS
: EXIT_FAILURE
);