src/factor.c

   1 /* factor -- print prime factors of n.
   2    Copyright (C) 86, 1995-2005, 2007-2009 Free Software Foundation, Inc.
   3
   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.
   8
   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.
  13
  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/>.  */
  16
  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.
  21 */
  22
  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
  31
  32 #include <assert.h>
  33
  34 #include "system.h"
  35 #include "error.h"
  36 #include "quote.h"
  37 #include "readtokens.h"
  38 #include "xstrtol.h"
  39
  40 /* The official name of this program (e.g., no `g' prefix).  */
  41 #define PROGRAM_NAME "factor"
  42
  43 #define AUTHORS proper_name ("Paul Rubin")
  44
  45 /* Token delimiters when reading from a file.  */
  46 #define DELIM "\n\t "
  47
  48 static bool verbose = false;
  49
  50 #if HAVE_GMP
  51 static mpz_t *factor = NULL;
  52 static size_t nfactors_found = 0;
  53 static size_t nfactors_allocated = 0;
  54
  55 static void
  56 debug (char const *fmt, ...)
  57 {
  58   if (verbose)
  59     {
  60       va_list ap;
  61       va_start (ap, fmt);
  62       vfprintf (stderr, fmt, ap);
  63     }
  64 }
  65
  66 static void
  67 emit_factor (mpz_t n)
  68 {
  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;
  74 }
  75
  76 static void
  77 emit_ul_factor (unsigned long int i)
  78 {
  79   mpz_t t;
  80   mpz_init (t);
  81   mpz_set_ui (t, i);
  82   emit_factor (t);
  83   mpz_clear (t);
  84 }
  85
  86 static void
  87 factor_using_division (mpz_t t, unsigned int limit)
  88 {
  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;
  95
  96   debug ("[trial division (%u)] ", limit);
  97
  98   mpz_init (q);
  99   mpz_init (r);
 100
 101   f = mpz_scan1 (t, 0);
 102   mpz_div_2exp (t, t, f);
 103   while (f)
 104     {
 105       emit_ul_factor (2);
 106       --f;
 107     }
 108
 109   for (;;)
 110     {
 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);
 116     }
 117
 118   for (;;)
 119     {
 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);
 125     }
 126
 127   failures = 0;
 128   f = 7;
 129   ai = 0;
 130   while (mpz_cmp_ui (t, 1) != 0)
 131     {
 132       mpz_tdiv_qr_ui (q, r, t, f);
 133       if (mpz_cmp_ui (r, 0) != 0)
 134         {
 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;
 142         }
 143       else
 144         {
 145           mpz_swap (t, q);
 146           emit_ul_factor (f);
 147           failures = 0;
 148         }
 149     }
 150
 151   mpz_clear (q);
 152   mpz_clear (r);
 153 }
 154
 155 static void
 156 factor_using_pollard_rho (mpz_t n, int a_int)
 157 {
 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;
 163
 164   debug ("[pollard-rho (%d)] ", a_int);
 165
 166   mpz_init (g);
 167   mpz_init (t1);
 168   mpz_init (t2);
 169
 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;
 178
 179   while (mpz_cmp_ui (n, 1) != 0)
 180     {
 181 S2:
 182       mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
 183
 184       mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
 185       c++;
 186       if (c == 20)
 187         {
 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);
 193         }
 194
 195       k--;
 196       if (k > 0)
 197         goto S2;
 198
 199       mpz_gcd (g, P, n);
 200       if (mpz_cmp_ui (g, 1) != 0)
 201         goto S4;
 202
 203       mpz_set (x1, x);
 204       k = l;
 205       l = 2 * l;
 206       for (i = 0; i < k; i++)
 207         {
 208           mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
 209         }
 210       mpz_set (y, x);
 211       c = 0;
 212       goto S2;
 213 S4:
 214       do
 215         {
 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);
 218         }
 219       while (mpz_cmp_ui (g, 1) == 0);
 220
 221       mpz_div (n, n, g);        /* divide by g, before g is overwritten */
 222
 223       if (!mpz_probab_prime_p (g, 3))
 224         {
 225           do
 226             {
 227               mp_limb_t a_limb;
 228               mpn_random (&a_limb, (mp_size_t) 1);
 229               a_int = (int) a_limb;
 230             }
 231           while (a_int == -2 || a_int == 0);
 232
 233           debug ("[composite factor--restarting pollard-rho] ");
 234           factor_using_pollard_rho (g, a_int);
 235         }
 236       else
 237         {
 238           emit_factor (g);
 239         }
 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))
 244         {
 245           emit_factor (n);
 246           break;
 247         }
 248     }
 249
 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);
 258 }
 259
 260 #else
 261
 262 static void
 263 debug (char const *fmt ATTRIBUTE_UNUSED, ...)
 264 {
 265 }
 266
 267 #endif
 268
 269 /* The maximum number of factors, including -1, for negative numbers.  */
 270 #define MAX_N_FACTORS (sizeof (uintmax_t) * CHAR_BIT)
 271
 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  */
 277
 278 #include "wheel-size.h"  /* For the definition of WHEEL_SIZE.  */
 279 static const unsigned char wheel_tab[] =
 280   {
 281 #include "wheel.h"
 282   };
 283
 284 #define WHEEL_START (wheel_tab + WHEEL_SIZE)
 285 #define WHEEL_END (wheel_tab + (sizeof wheel_tab / sizeof wheel_tab[0]))
 286
 287 /* FIXME: comment */
 288
 289 static size_t
 290 factor_wheel (uintmax_t n0, size_t max_n_factors, uintmax_t *factors)
 291 {
 292   uintmax_t n = n0, d, q;
 293   size_t n_factors = 0;
 294   unsigned char const *w = wheel_tab;
 295
 296   if (n <= 1)
 297     return n_factors;
 298
 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. */
 306
 307   d = 2;
 308   do
 309     {
 310       q = n / d;
 311       while (n == q * d)
 312         {
 313           assert (n_factors < max_n_factors);
 314           factors[n_factors++] = d;
 315           n = q;
 316           q = n / d;
 317         }
 318       d += *(w++);
 319       if (w == WHEEL_END)
 320         w = WHEEL_START;
 321     }
 322   while (d <= q);
 323
 324   if (n != 1 || n0 == 1)
 325     {
 326       assert (n_factors < max_n_factors);
 327       factors[n_factors++] = n;
 328     }
 329
 330   return n_factors;
 331 }
 332
 333 /* Single-precision factoring */
 334 static void
 335 print_factors_single (uintmax_t n)
 336 {
 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)];
 341
 342   printf ("%s:", umaxtostr (n, buf));
 343   for (i = 0; i < n_factors; i++)
 344     printf (" %s", umaxtostr (factors[i], buf));
 345   putchar ('\n');
 346 }
 347
 348 #if HAVE_GMP
 349 static int
 350 mpcompare (const void *av, const void *bv)
 351 {
 352   mpz_t *const *a = av;
 353   mpz_t *const *b = bv;
 354   return mpz_cmp (**a, **b);
 355 }
 356
 357 static void
 358 sort_and_print_factors (void)
 359 {
 360   mpz_t **faclist;
 361   size_t i;
 362
 363   faclist = xcalloc (nfactors_found, sizeof *faclist);
 364   for (i = 0; i < nfactors_found; ++i)
 365     {
 366       faclist[i] = &factor[i];
 367     }
 368   qsort (faclist, nfactors_found, sizeof *faclist, mpcompare);
 369
 370   for (i = 0; i < nfactors_found; ++i)
 371     {
 372       fputc (' ', stdout);
 373       mpz_out_str (stdout, 10, *faclist[i]);
 374     }
 375   putchar ('\n');
 376   free (faclist);
 377 }
 378
 379 static void
 380 free_factors (void)
 381 {
 382   size_t i;
 383
 384   for (i = 0; i < nfactors_found; ++i)
 385     {
 386       mpz_clear (factor[i]);
 387     }
 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;
 391 }
 392
 393 /* Arbitrary-precision factoring */
 394 static void
 395 print_factors_multi (mpz_t t)
 396 {
 397   mpz_out_str (stdout, 10, t);
 398   putchar (':');
 399
 400   if (mpz_sgn (t) != 0)
 401     {
 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;
 406
 407       factor_using_division (t, division_limit);
 408
 409       if (mpz_cmp_ui (t, 1) != 0)
 410         {
 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);
 416         }
 417     }
 418
 419   mpz_clear (t);
 420   sort_and_print_factors ();
 421   free_factors ();
 422 }
 423 #endif
 424
 425
 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)
 433 {
 434   uintmax_t n;
 435   strtol_error err = xstrtoumax (s, NULL, 10, &n, "");
 436
 437 #if HAVE_GMP
 438   enum { GMP_TURNOVER_POINT = 100000 };
 439
 440   if (err == LONGINT_OVERFLOW
 441       || (err == LONGINT_OK && GMP_TURNOVER_POINT <= n))
 442     {
 443       mpz_t t;
 444       mpz_init (t);
 445       if (gmp_sscanf (s, "%Zd", t) == 1)
 446         {
 447           debug ("[%s]", _("using arbitrary-precision arithmetic"));
 448           print_factors_multi (t);
 449           return true;
 450         }
 451       err = LONGINT_INVALID;
 452     }
 453 #endif
 454
 455   switch (err)
 456     {
 457     case LONGINT_OK:
 458       debug ("[%s]", _("using single-precision arithmetic"));
 459       print_factors_single (n);
 460       return true;
 461
 462     case LONGINT_OVERFLOW:
 463       error (0, 0, _("%s is too large"), quote (s));
 464       return false;
 465
 466     default:
 467       error (0, 0, _("%s is not a valid positive integer"), quote (s));
 468       return false;
 469     }
 470 }
 471
 472 enum
 473 {
 474   VERBOSE_OPTION = CHAR_MAX + 1
 475 };
 476
 477 static struct option const long_options[] =
 478 {
 479   {"verbose", no_argument, NULL, VERBOSE_OPTION},
 480   {GETOPT_HELP_OPTION_DECL},
 481   {GETOPT_VERSION_OPTION_DECL},
 482   {NULL, 0, NULL, 0}
 483 };
 484
 485 void
 486 usage (int status)
 487 {
 488   if (status != EXIT_SUCCESS)
 489     fprintf (stderr, _("Try `%s --help' for more information.\n"),
 490              program_name);
 491   else
 492     {
 493       printf (_("\
 494 Usage: %s [NUMBER]...\n\
 495   or:  %s OPTION\n\
 496 "),
 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\
 501 \n\
 502 "), stdout);
 503       fputs (HELP_OPTION_DESCRIPTION, stdout);
 504       fputs (VERSION_OPTION_DESCRIPTION, stdout);
 505       emit_bug_reporting_address ();
 506     }
 507   exit (status);
 508 }
 509
 510 static bool
 511 do_stdin (void)
 512 {
 513   bool ok = true;
 514   token_buffer tokenbuffer;
 515
 516   init_tokenbuffer (&tokenbuffer);
 517
 518   for (;;)
 519     {
 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);
 525     }
 526   free (tokenbuffer.buffer);
 527
 528   return ok;
 529 }
 530
 531 int
 532 main (int argc, char **argv)
 533 {
 534   bool ok;
 535   int c;
 536
 537   initialize_main (&argc, &argv);
 538   set_program_name (argv[0]);
 539   setlocale (LC_ALL, "");
 540   bindtextdomain (PACKAGE, LOCALEDIR);
 541   textdomain (PACKAGE);
 542
 543   atexit (close_stdout);
 544
 545   while ((c = getopt_long (argc, argv, "", long_options, NULL)) != -1)
 546     {
 547       switch (c)
 548         {
 549         case VERBOSE_OPTION:
 550           verbose = true;
 551           break;
 552
 553         case_GETOPT_HELP_CHAR;
 554
 555         case_GETOPT_VERSION_CHAR (PROGRAM_NAME, AUTHORS);
 556
 557         default:
 558           usage (EXIT_FAILURE);
 559         }
 560     }
 561
 562   if (argc <= optind)
 563     ok = do_stdin ();
 564   else
 565     {
 566       int i;
 567       ok = true;
 568       for (i = optind; i < argc; i++)
 569         if (! print_factors (argv[i]))
 570           ok = false;
 571     }
 572 #if HAVE_GMP
 573   free (factor);
 574 #endif
 575   exit (ok ? EXIT_SUCCESS : EXIT_FAILURE);
 576 }