src/factor.c

   1 /* factor -- print prime factors of n.
   2    Copyright (C) 1986, 1995-2005, 2007-2011 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       va_end (ap);
  64     }
  65 }
  66
  67 static void
  68 emit_factor (mpz_t n)
  69 {
  70   if (nfactors_found == nfactors_allocated)
  71     factor = X2NREALLOC (factor, &nfactors_allocated);
  72   mpz_init (factor[nfactors_found]);
  73   mpz_set (factor[nfactors_found], n);
  74   ++nfactors_found;
  75 }
  76
  77 static void
  78 emit_ul_factor (unsigned long int i)
  79 {
  80   mpz_t t;
  81   mpz_init (t);
  82   mpz_set_ui (t, i);
  83   emit_factor (t);
  84   mpz_clear (t);
  85 }
  86
  87 static void
  88 factor_using_division (mpz_t t, unsigned int limit)
  89 {
  90   mpz_t q, r;
  91   unsigned long int f;
  92   int ai;
  93   static unsigned int const add[] = {4, 2, 4, 2, 4, 6, 2, 6};
  94   unsigned int const *addv = add;
  95   unsigned int failures;
  96
  97   debug ("[trial division (%u)] ", limit);
  98
  99   mpz_init (q);
 100   mpz_init (r);
 101
 102   f = mpz_scan1 (t, 0);
 103   mpz_div_2exp (t, t, f);
 104   while (f)
 105     {
 106       emit_ul_factor (2);
 107       --f;
 108     }
 109
 110   while (true)
 111     {
 112       mpz_tdiv_qr_ui (q, r, t, 3);
 113       if (mpz_cmp_ui (r, 0) != 0)
 114         break;
 115       mpz_set (t, q);
 116       emit_ul_factor (3);
 117     }
 118
 119   while (true)
 120     {
 121       mpz_tdiv_qr_ui (q, r, t, 5);
 122       if (mpz_cmp_ui (r, 0) != 0)
 123         break;
 124       mpz_set (t, q);
 125       emit_ul_factor (5);
 126     }
 127
 128   failures = 0;
 129   f = 7;
 130   ai = 0;
 131   while (mpz_cmp_ui (t, 1) != 0)
 132     {
 133       mpz_tdiv_qr_ui (q, r, t, f);
 134       if (mpz_cmp_ui (r, 0) != 0)
 135         {
 136           f += addv[ai];
 137           if (mpz_cmp_ui (q, f) < 0)
 138             break;
 139           ai = (ai + 1) & 7;
 140           failures++;
 141           if (failures > limit)
 142             break;
 143         }
 144       else
 145         {
 146           mpz_swap (t, q);
 147           emit_ul_factor (f);
 148           failures = 0;
 149         }
 150     }
 151
 152   mpz_clear (q);
 153   mpz_clear (r);
 154 }
 155
 156 static void
 157 factor_using_pollard_rho (mpz_t n, int a_int)
 158 {
 159   mpz_t x, x1, y, P;
 160   mpz_t a;
 161   mpz_t g;
 162   mpz_t t1, t2;
 163   int k, l, c;
 164
 165   debug ("[pollard-rho (%d)] ", a_int);
 166
 167   mpz_init (g);
 168   mpz_init (t1);
 169   mpz_init (t2);
 170
 171   mpz_init_set_si (a, a_int);
 172   mpz_init_set_si (y, 2);
 173   mpz_init_set_si (x, 2);
 174   mpz_init_set_si (x1, 2);
 175   k = 1;
 176   l = 1;
 177   mpz_init_set_ui (P, 1);
 178   c = 0;
 179
 180   while (mpz_cmp_ui (n, 1) != 0)
 181     {
 182 S2:
 183       mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
 184
 185       mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
 186       c++;
 187       if (c == 20)
 188         {
 189           c = 0;
 190           mpz_gcd (g, P, n);
 191           if (mpz_cmp_ui (g, 1) != 0)
 192             goto S4;
 193           mpz_set (y, x);
 194         }
 195
 196       k--;
 197       if (k > 0)
 198         goto S2;
 199
 200       mpz_gcd (g, P, n);
 201       if (mpz_cmp_ui (g, 1) != 0)
 202         goto S4;
 203
 204       mpz_set (x1, x);
 205       k = l;
 206       l = 2 * l;
 207       unsigned int i;
 208       for (i = 0; i < k; i++)
 209         {
 210           mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
 211         }
 212       mpz_set (y, x);
 213       c = 0;
 214       goto S2;
 215 S4:
 216       do
 217         {
 218           mpz_mul (y, y, y); mpz_add (y, y, a); mpz_mod (y, y, n);
 219           mpz_sub (t1, x1, y); mpz_gcd (g, t1, n);
 220         }
 221       while (mpz_cmp_ui (g, 1) == 0);
 222
 223       mpz_div (n, n, g);        /* divide by g, before g is overwritten */
 224
 225       if (!mpz_probab_prime_p (g, 3))
 226         {
 227           do
 228             {
 229               mp_limb_t a_limb;
 230               mpn_random (&a_limb, (mp_size_t) 1);
 231               a_int = (int) a_limb;
 232             }
 233           while (a_int == -2 || a_int == 0);
 234
 235           debug ("[composite factor--restarting pollard-rho] ");
 236           factor_using_pollard_rho (g, a_int);
 237         }
 238       else
 239         {
 240           emit_factor (g);
 241         }
 242       mpz_mod (x, x, n);
 243       mpz_mod (x1, x1, n);
 244       mpz_mod (y, y, n);
 245       if (mpz_probab_prime_p (n, 3))
 246         {
 247           emit_factor (n);
 248           break;
 249         }
 250     }
 251
 252   mpz_clear (g);
 253   mpz_clear (P);
 254   mpz_clear (t2);
 255   mpz_clear (t1);
 256   mpz_clear (a);
 257   mpz_clear (x1);
 258   mpz_clear (x);
 259   mpz_clear (y);
 260 }
 261
 262 #else
 263
 264 static void
 265 debug (char const *fmt ATTRIBUTE_UNUSED, ...)
 266 {
 267 }
 268
 269 #endif
 270
 271 /* The maximum number of factors, including -1, for negative numbers.  */
 272 #define MAX_N_FACTORS (sizeof (uintmax_t) * CHAR_BIT)
 273
 274 /* The trial divisor increment wheel.  Use it to skip over divisors that
 275    are composites of 2, 3, 5, 7, or 11.  The part from WHEEL_START up to
 276    WHEEL_END is reused periodically, while the "lead in" is used to test
 277    for those primes and to jump onto the wheel.  For more information, see
 278    http://www.utm.edu/research/primes/glossary/WheelFactorization.html  */
 279
 280 #include "wheel-size.h"  /* For the definition of WHEEL_SIZE.  */
 281 static const unsigned char wheel_tab[] =
 282   {
 283 #include "wheel.h"
 284   };
 285
 286 #define WHEEL_START (wheel_tab + WHEEL_SIZE)
 287 #define WHEEL_END (wheel_tab + ARRAY_CARDINALITY (wheel_tab))
 288
 289 /* FIXME: comment */
 290
 291 static size_t
 292 factor_wheel (uintmax_t n0, size_t max_n_factors, uintmax_t *factors)
 293 {
 294   uintmax_t n = n0, d, q;
 295   size_t n_factors = 0;
 296   unsigned char const *w = wheel_tab;
 297
 298   if (n <= 1)
 299     return n_factors;
 300
 301   /* The exit condition in the following loop is correct because
 302      any time it is tested one of these 3 conditions holds:
 303       (1) d divides n
 304       (2) n is prime
 305       (3) n is composite but has no factors less than d.
 306      If (1) or (2) obviously the right thing happens.
 307      If (3), then since n is composite it is >= d^2. */
 308
 309   d = 2;
 310   do
 311     {
 312       q = n / d;
 313       while (n == q * d)
 314         {
 315           assert (n_factors < max_n_factors);
 316           factors[n_factors++] = d;
 317           n = q;
 318           q = n / d;
 319         }
 320       d += *(w++);
 321       if (w == WHEEL_END)
 322         w = WHEEL_START;
 323     }
 324   while (d <= q);
 325
 326   if (n != 1 || n0 == 1)
 327     {
 328       assert (n_factors < max_n_factors);
 329       factors[n_factors++] = n;
 330     }
 331
 332   return n_factors;
 333 }
 334
 335 /* Single-precision factoring */
 336 static void
 337 print_factors_single (uintmax_t n)
 338 {
 339   uintmax_t factors[MAX_N_FACTORS];
 340   size_t n_factors = factor_wheel (n, MAX_N_FACTORS, factors);
 341   size_t i;
 342   char buf[INT_BUFSIZE_BOUND (uintmax_t)];
 343
 344   printf ("%s:", umaxtostr (n, buf));
 345   for (i = 0; i < n_factors; i++)
 346     printf (" %s", umaxtostr (factors[i], buf));
 347   putchar ('\n');
 348 }
 349
 350 #if HAVE_GMP
 351 static int
 352 mpcompare (const void *av, const void *bv)
 353 {
 354   mpz_t *const *a = av;
 355   mpz_t *const *b = bv;
 356   return mpz_cmp (**a, **b);
 357 }
 358
 359 static void
 360 sort_and_print_factors (void)
 361 {
 362   mpz_t **faclist;
 363   size_t i;
 364
 365   faclist = xcalloc (nfactors_found, sizeof *faclist);
 366   for (i = 0; i < nfactors_found; ++i)
 367     {
 368       faclist[i] = &factor[i];
 369     }
 370   qsort (faclist, nfactors_found, sizeof *faclist, mpcompare);
 371
 372   for (i = 0; i < nfactors_found; ++i)
 373     {
 374       fputc (' ', stdout);
 375       mpz_out_str (stdout, 10, *faclist[i]);
 376     }
 377   putchar ('\n');
 378   free (faclist);
 379 }
 380
 381 static void
 382 free_factors (void)
 383 {
 384   size_t i;
 385
 386   for (i = 0; i < nfactors_found; ++i)
 387     {
 388       mpz_clear (factor[i]);
 389     }
 390   /* Don't actually free factor[] because in the case where we are
 391      reading numbers from stdin, we may be about to use it again.  */
 392   nfactors_found = 0;
 393 }
 394
 395 /* Arbitrary-precision factoring */
 396 static void
 397 print_factors_multi (mpz_t t)
 398 {
 399   mpz_out_str (stdout, 10, t);
 400   putchar (':');
 401
 402   if (mpz_sgn (t) != 0)
 403     {
 404       /* Set the trial division limit according to the size of t.  */
 405       size_t n_bits = mpz_sizeinbase (t, 2);
 406       unsigned int division_limit = MIN (n_bits, 1000);
 407       division_limit *= division_limit;
 408
 409       factor_using_division (t, division_limit);
 410
 411       if (mpz_cmp_ui (t, 1) != 0)
 412         {
 413           debug ("[is number prime?] ");
 414           if (mpz_probab_prime_p (t, 3))
 415             emit_factor (t);
 416           else
 417             factor_using_pollard_rho (t, 1);
 418         }
 419     }
 420
 421   mpz_clear (t);
 422   sort_and_print_factors ();
 423   free_factors ();
 424 }
 425 #endif
 426
 427
 428 /* Emit the factors of the indicated number.  If we have the option of using
 429    either algorithm, we select on the basis of the length of the number.
 430    For longer numbers, we prefer the MP algorithm even if the native algorithm
 431    has enough digits, because the algorithm is better.  The turnover point
 432    depends on the value.  */
 433 static bool
 434 print_factors (char const *s)
 435 {
 436   uintmax_t n;
 437   strtol_error err = xstrtoumax (s, NULL, 10, &n, "");
 438
 439 #if HAVE_GMP
 440   enum { GMP_TURNOVER_POINT = 100000 };
 441
 442   if (err == LONGINT_OVERFLOW
 443       || (err == LONGINT_OK && GMP_TURNOVER_POINT <= n))
 444     {
 445       mpz_t t;
 446       mpz_init (t);
 447       if (gmp_sscanf (s, "%Zd", t) == 1)
 448         {
 449           debug ("[%s]", _("using arbitrary-precision arithmetic"));
 450           print_factors_multi (t);
 451           return true;
 452         }
 453       err = LONGINT_INVALID;
 454     }
 455 #endif
 456
 457   switch (err)
 458     {
 459     case LONGINT_OK:
 460       debug ("[%s]", _("using single-precision arithmetic"));
 461       print_factors_single (n);
 462       return true;
 463
 464     case LONGINT_OVERFLOW:
 465       error (0, 0, _("%s is too large"), quote (s));
 466       return false;
 467
 468     default:
 469       error (0, 0, _("%s is not a valid positive integer"), quote (s));
 470       return false;
 471     }
 472 }
 473
 474 enum
 475 {
 476   VERBOSE_OPTION = CHAR_MAX + 1
 477 };
 478
 479 static struct option const long_options[] =
 480 {
 481   {"verbose", no_argument, NULL, VERBOSE_OPTION},
 482   {GETOPT_HELP_OPTION_DECL},
 483   {GETOPT_VERSION_OPTION_DECL},
 484   {NULL, 0, NULL, 0}
 485 };
 486
 487 void
 488 usage (int status)
 489 {
 490   if (status != EXIT_SUCCESS)
 491     fprintf (stderr, _("Try `%s --help' for more information.\n"),
 492              program_name);
 493   else
 494     {
 495       printf (_("\
 496 Usage: %s [NUMBER]...\n\
 497   or:  %s OPTION\n\
 498 "),
 499               program_name, program_name);
 500       fputs (_("\
 501 Print the prime factors of each specified integer NUMBER.  If none\n\
 502 are specified on the command line, read them from standard input.\n\
 503 \n\
 504 "), stdout);
 505       fputs (HELP_OPTION_DESCRIPTION, stdout);
 506       fputs (VERSION_OPTION_DESCRIPTION, stdout);
 507       emit_ancillary_info ();
 508     }
 509   exit (status);
 510 }
 511
 512 static bool
 513 do_stdin (void)
 514 {
 515   bool ok = true;
 516   token_buffer tokenbuffer;
 517
 518   init_tokenbuffer (&tokenbuffer);
 519
 520   while (true)
 521     {
 522       size_t token_length = readtoken (stdin, DELIM, sizeof (DELIM) - 1,
 523                                        &tokenbuffer);
 524       if (token_length == (size_t) -1)
 525         break;
 526       ok &= print_factors (tokenbuffer.buffer);
 527     }
 528   free (tokenbuffer.buffer);
 529
 530   return ok;
 531 }
 532
 533 int
 534 main (int argc, char **argv)
 535 {
 536   bool ok;
 537   int c;
 538
 539   initialize_main (&argc, &argv);
 540   set_program_name (argv[0]);
 541   setlocale (LC_ALL, "");
 542   bindtextdomain (PACKAGE, LOCALEDIR);
 543   textdomain (PACKAGE);
 544
 545   atexit (close_stdout);
 546
 547   while ((c = getopt_long (argc, argv, "", long_options, NULL)) != -1)
 548     {
 549       switch (c)
 550         {
 551         case VERBOSE_OPTION:
 552           verbose = true;
 553           break;
 554
 555         case_GETOPT_HELP_CHAR;
 556
 557         case_GETOPT_VERSION_CHAR (PROGRAM_NAME, AUTHORS);
 558
 559         default:
 560           usage (EXIT_FAILURE);
 561         }
 562     }
 563
 564   if (argc <= optind)
 565     ok = do_stdin ();
 566   else
 567     {
 568       int i;
 569       ok = true;
 570       for (i = optind; i < argc; i++)
 571         if (! print_factors (argv[i]))
 572           ok = false;
 573     }
 574 #if HAVE_GMP
 575   free (factor);
 576 #endif
 577   exit (ok ? EXIT_SUCCESS : EXIT_FAILURE);
 578 }