src/factor.c

   1 /* factor -- print prime factors of n.
   2    Copyright (C) 86, 1995-2005, 2007-2008 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
  33 #include <assert.h>
  34 #define NDEBUG 1
  35
  36 #include "system.h"
  37 #include "error.h"
  38 #include "quote.h"
  39 #include "readtokens.h"
  40 #include "xstrtol.h"
  41
  42 /* The official name of this program (e.g., no `g' prefix).  */
  43 #define PROGRAM_NAME "factor"
  44
  45 #define AUTHORS proper_name ("Paul Rubin")
  46
  47 /* Token delimiters when reading from a file.  */
  48 #define DELIM "\n\t "
  49
  50 static bool verbose = false;
  51
  52 typedef enum
  53   {
  54     ALGO_AUTOSELECT,            /* default */
  55     ALGO_GMP,                   /* --bignum */
  56     ALGO_SINGLE                 /* --no-bignum */
  57   } ALGORITHM_CHOICE;
  58 static ALGORITHM_CHOICE algorithm = ALGO_AUTOSELECT;
  59
  60
  61 #if HAVE_GMP
  62 static mpz_t *factor = NULL;
  63 static size_t nfactors_found = 0;
  64 static size_t nfactors_allocated = 0;
  65
  66 static void
  67 debug (char const *fmt, ...)
  68 {
  69   if (verbose)
  70     {
  71       va_list ap;
  72       va_start (ap, fmt);
  73       vfprintf (stderr, fmt, ap);
  74     }
  75 }
  76
  77 static void
  78 emit_factor (mpz_t n)
  79 {
  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);
  84   ++nfactors_found;
  85 }
  86
  87 static void
  88 emit_ul_factor (unsigned long int i)
  89 {
  90   mpz_t t;
  91   mpz_init (t);
  92   mpz_set_ui (t, i);
  93   emit_factor (t);
  94 }
  95
  96 static void
  97 factor_using_division (mpz_t t, unsigned int limit)
  98 {
  99   mpz_t q, r;
 100   unsigned long int f;
 101   int ai;
 102   static unsigned int const add[] = {4, 2, 4, 2, 4, 6, 2, 6};
 103   unsigned int const *addv = add;
 104   unsigned int failures;
 105
 106   debug ("[trial division (%u)] ", limit);
 107
 108   mpz_init (q);
 109   mpz_init (r);
 110
 111   f = mpz_scan1 (t, 0);
 112   mpz_div_2exp (t, t, f);
 113   while (f)
 114     {
 115       emit_ul_factor (2);
 116       --f;
 117     }
 118
 119   for (;;)
 120     {
 121       mpz_tdiv_qr_ui (q, r, t, 3);
 122       if (mpz_cmp_ui (r, 0) != 0)
 123         break;
 124       mpz_set (t, q);
 125       emit_ul_factor (3);
 126     }
 127
 128   for (;;)
 129     {
 130       mpz_tdiv_qr_ui (q, r, t, 5);
 131       if (mpz_cmp_ui (r, 0) != 0)
 132         break;
 133       mpz_set (t, q);
 134       emit_ul_factor (5);
 135     }
 136
 137   failures = 0;
 138   f = 7;
 139   ai = 0;
 140   while (mpz_cmp_ui (t, 1) != 0)
 141     {
 142       mpz_tdiv_qr_ui (q, r, t, f);
 143       if (mpz_cmp_ui (r, 0) != 0)
 144         {
 145           f += addv[ai];
 146           if (mpz_cmp_ui (q, f) < 0)
 147             break;
 148           ai = (ai + 1) & 7;
 149           failures++;
 150           if (failures > limit)
 151             break;
 152         }
 153       else
 154         {
 155           mpz_swap (t, q);
 156           emit_ul_factor (f);
 157           failures = 0;
 158         }
 159     }
 160
 161   mpz_clear (q);
 162   mpz_clear (r);
 163 }
 164
 165 static void
 166 factor_using_pollard_rho (mpz_t n, int a_int)
 167 {
 168   mpz_t x, x1, y, P;
 169   mpz_t a;
 170   mpz_t g;
 171   mpz_t t1, t2;
 172   int k, l, c, i;
 173
 174   debug ("[pollard-rho (%d)] ", a_int);
 175
 176   mpz_init (g);
 177   mpz_init (t1);
 178   mpz_init (t2);
 179
 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);
 184   k = 1;
 185   l = 1;
 186   mpz_init_set_ui (P, 1);
 187   c = 0;
 188
 189   while (mpz_cmp_ui (n, 1) != 0)
 190     {
 191 S2:
 192       mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
 193
 194       mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
 195       c++;
 196       if (c == 20)
 197         {
 198           c = 0;
 199           mpz_gcd (g, P, n);
 200           if (mpz_cmp_ui (g, 1) != 0)
 201             goto S4;
 202           mpz_set (y, x);
 203         }
 204
 205       k--;
 206       if (k > 0)
 207         goto S2;
 208
 209       mpz_gcd (g, P, n);
 210       if (mpz_cmp_ui (g, 1) != 0)
 211         goto S4;
 212
 213       mpz_set (x1, x);
 214       k = l;
 215       l = 2 * l;
 216       for (i = 0; i < k; i++)
 217         {
 218           mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
 219         }
 220       mpz_set (y, x);
 221       c = 0;
 222       goto S2;
 223 S4:
 224       do
 225         {
 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);
 228         }
 229       while (mpz_cmp_ui (g, 1) == 0);
 230
 231       mpz_div (n, n, g);        /* divide by g, before g is overwritten */
 232
 233       if (!mpz_probab_prime_p (g, 3))
 234         {
 235           do
 236             {
 237               mp_limb_t a_limb;
 238               mpn_random (&a_limb, (mp_size_t) 1);
 239               a_int = (int) a_limb;
 240             }
 241           while (a_int == -2 || a_int == 0);
 242
 243           debug ("[composite factor--restarting pollard-rho] ");
 244           factor_using_pollard_rho (g, a_int);
 245         }
 246       else
 247         {
 248           emit_factor (g);
 249         }
 250       mpz_mod (x, x, n);
 251       mpz_mod (x1, x1, n);
 252       mpz_mod (y, y, n);
 253       if (mpz_probab_prime_p (n, 3))
 254         {
 255           emit_factor (n);
 256           break;
 257         }
 258     }
 259
 260   mpz_clear (g);
 261   mpz_clear (P);
 262   mpz_clear (t2);
 263   mpz_clear (t1);
 264   mpz_clear (a);
 265   mpz_clear (x1);
 266   mpz_clear (x);
 267   mpz_clear (y);
 268 }
 269
 270 /* Arbitrary-precision factoring */
 271 static bool
 272 extract_factors_multi (char const *s)
 273 {
 274   unsigned int division_limit;
 275   size_t n_bits;
 276   mpz_t t;
 277
 278   mpz_init (t);
 279   if (1 != gmp_sscanf (s, "%Zd", t))
 280     {
 281       error (0, 0, _("%s is not a valid positive integer"), quote (s));
 282       return false;
 283     }
 284
 285   printf ("%s:", s);
 286
 287   if (mpz_sgn (t) == 0)
 288     {
 289       mpz_clear (t);
 290       return true;              /* 0; no factors */
 291     }
 292
 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;
 297
 298   factor_using_division (t, division_limit);
 299
 300   if (mpz_cmp_ui (t, 1) != 0)
 301     {
 302       debug ("[is number prime?] ");
 303       if (mpz_probab_prime_p (t, 3))
 304         {
 305           emit_factor (t);
 306         }
 307       else
 308         {
 309           factor_using_pollard_rho (t, 1);
 310         }
 311     }
 312   mpz_clear (t);
 313   return true;
 314 }
 315 #endif
 316
 317 /* The maximum number of factors, including -1, for negative numbers.  */
 318 #define MAX_N_FACTORS (sizeof (uintmax_t) * CHAR_BIT)
 319
 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  */
 325
 326 #include "wheel-size.h"  /* For the definition of WHEEL_SIZE.  */
 327 static const unsigned char wheel_tab[] =
 328   {
 329 #include "wheel.h"
 330   };
 331
 332 #define WHEEL_START (wheel_tab + WHEEL_SIZE)
 333 #define WHEEL_END (wheel_tab + (sizeof wheel_tab / sizeof wheel_tab[0]))
 334
 335 /* FIXME: comment */
 336
 337 static size_t
 338 factor_wheel (uintmax_t n0, size_t max_n_factors, uintmax_t *factors)
 339 {
 340   uintmax_t n = n0, d, q;
 341   size_t n_factors = 0;
 342   unsigned char const *w = wheel_tab;
 343
 344   if (n <= 1)
 345     return n_factors;
 346
 347   /* The exit condition in the following loop is correct because
 348      any time it is tested one of these 3 conditions holds:
 349       (1) d divides n
 350       (2) n is prime
 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. */
 354
 355   d = 2;
 356   do
 357     {
 358       q = n / d;
 359       while (n == q * d)
 360         {
 361           assert (n_factors < max_n_factors);
 362           factors[n_factors++] = d;
 363           n = q;
 364           q = n / d;
 365         }
 366       d += *(w++);
 367       if (w == WHEEL_END)
 368         w = WHEEL_START;
 369     }
 370   while (d <= q);
 371
 372   if (n != 1 || n0 == 1)
 373     {
 374       assert (n_factors < max_n_factors);
 375       factors[n_factors++] = n;
 376     }
 377
 378   return n_factors;
 379 }
 380
 381 /* Single-precision factoring */
 382 static bool
 383 print_factors_single (char const *s)
 384 {
 385   uintmax_t factors[MAX_N_FACTORS];
 386   uintmax_t n;
 387   size_t n_factors;
 388   size_t i;
 389   char buf[INT_BUFSIZE_BOUND (uintmax_t)];
 390   strtol_error err;
 391
 392   if ((err = xstrtoumax (s, NULL, 10, &n, "")) != LONGINT_OK)
 393     {
 394       if (err == LONGINT_OVERFLOW)
 395         error (0, 0, _("%s is too large"), quote (s));
 396       else
 397         error (0, 0, _("%s is not a valid positive integer"), quote (s));
 398       return false;
 399     }
 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));
 404   putchar ('\n');
 405   return true;
 406 }
 407
 408 #if HAVE_GMP
 409 static int
 410 mpcompare (const void *av, const void *bv)
 411 {
 412   mpz_t *const *a = av;
 413   mpz_t *const *b = bv;
 414   return mpz_cmp (**a, **b);
 415 }
 416
 417 static void
 418 sort_and_print_factors (void)
 419 {
 420   mpz_t **faclist;
 421   size_t i;
 422
 423   faclist = xcalloc (nfactors_found, sizeof *faclist);
 424   for (i = 0; i < nfactors_found; ++i)
 425     {
 426       faclist[i] = &factor[i];
 427     }
 428   qsort (faclist, nfactors_found, sizeof *faclist, mpcompare);
 429
 430   for (i = 0; i < nfactors_found; ++i)
 431     {
 432       fputc (' ', stdout);
 433       mpz_out_str (stdout, 10, *faclist[i]);
 434     }
 435   putchar ('\n');
 436   free (faclist);
 437 }
 438
 439 static void
 440 free_factors (void)
 441 {
 442   size_t i;
 443
 444   for (i = 0; i < nfactors_found; ++i)
 445     {
 446       mpz_clear (factor[i]);
 447     }
 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.  */
 450   nfactors_found = 0;
 451 }
 452
 453
 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.
 460 */
 461 static bool
 462 print_factors (char const *s)
 463 {
 464   if (ALGO_AUTOSELECT == algorithm)
 465     {
 466       const size_t digits = strlen (s); /* upper limit on number of digits */
 467       algorithm = digits < 6 ? ALGO_SINGLE : ALGO_GMP;
 468     }
 469   if (ALGO_GMP == algorithm)
 470     {
 471       debug ("[%s]", _("using arbitrary-precision arithmetic"));
 472       if (extract_factors_multi (s))
 473         {
 474           sort_and_print_factors ();
 475           free_factors ();
 476           return true;
 477         }
 478       else
 479         {
 480           return false;
 481         }
 482     }
 483   else
 484     {
 485       debug ("[%s]", _("using single-precision arithmetic"));
 486       return print_factors_single (s);
 487     }
 488 }
 489 #else
 490 static bool
 491 print_factors (const char *s)   /* Single-precision only. */
 492 {
 493   if (ALGO_GMP == algorithm)
 494     {
 495       error (0, 0, _("arbitrary-precision arithmetic is not available"));
 496       return false;
 497     }
 498   return print_factors_single (s);
 499 }
 500 #endif
 501
 502 enum
 503 {
 504   VERBOSE_OPTION = CHAR_MAX + 1,
 505   USE_BIGNUM,
 506   NO_USE_BIGNUM
 507 };
 508
 509 static struct option const long_options[] =
 510 {
 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},
 516   {NULL, 0, NULL, 0}
 517 };
 518
 519 void
 520 usage (int status)
 521 {
 522   if (status != EXIT_SUCCESS)
 523     fprintf (stderr, _("Try `%s --help' for more information.\n"),
 524              program_name);
 525   else
 526     {
 527       printf (_("\
 528 Usage: %s [NUMBER]...\n\
 529   or:  %s OPTION\n\
 530 "),
 531               program_name, program_name);
 532       fputs (_("\
 533 Print the prime factors of each NUMBER.\n\
 534 \n\
 535 "), stdout);
 536       fputs (_("\
 537       --bignum     always use arbitrary-precision arithmetic\n\
 538       --no-bignum  always use single-precision arithmetic\n"),
 539                stdout);
 540       fputs (HELP_OPTION_DESCRIPTION, stdout);
 541       fputs (VERSION_OPTION_DESCRIPTION, stdout);
 542       fputs (_("\
 543 \n\
 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\
 546 "), stdout);
 547       emit_bug_reporting_address ();
 548     }
 549   exit (status);
 550 }
 551
 552 static bool
 553 do_stdin (void)
 554 {
 555   bool ok = true;
 556   token_buffer tokenbuffer;
 557
 558   init_tokenbuffer (&tokenbuffer);
 559
 560   for (;;)
 561     {
 562       size_t token_length = readtoken (stdin, DELIM, sizeof (DELIM) - 1,
 563                                        &tokenbuffer);
 564       if (token_length == (size_t) -1)
 565         break;
 566       ok &= print_factors (tokenbuffer.buffer);
 567     }
 568   free (tokenbuffer.buffer);
 569
 570   return ok;
 571 }
 572
 573 int
 574 main (int argc, char **argv)
 575 {
 576   bool ok;
 577   int c;
 578
 579   initialize_main (&argc, &argv);
 580   set_program_name (argv[0]);
 581   setlocale (LC_ALL, "");
 582   bindtextdomain (PACKAGE, LOCALEDIR);
 583   textdomain (PACKAGE);
 584
 585   atexit (close_stdout);
 586
 587   while ((c = getopt_long (argc, argv, "", long_options, NULL)) != -1)
 588     {
 589       switch (c)
 590         {
 591         case VERBOSE_OPTION:
 592           verbose = true;
 593           break;
 594
 595         case USE_BIGNUM:
 596           algorithm = ALGO_GMP;
 597           break;
 598
 599         case NO_USE_BIGNUM:
 600           algorithm = ALGO_SINGLE;
 601           break;
 602
 603         case_GETOPT_HELP_CHAR;
 604
 605         case_GETOPT_VERSION_CHAR (PROGRAM_NAME, AUTHORS);
 606
 607         default:
 608           usage (EXIT_FAILURE);
 609         }
 610     }
 611
 612   if (argc <= optind)
 613     ok = do_stdin ();
 614   else
 615     {
 616       int i;
 617       ok = true;
 618       for (i = optind; i < argc; i++)
 619         if (! print_factors (argv[i]))
 620           ok = false;
 621     }
 622 #if HAVE_GMP
 623   free (factor);
 624 #endif
 625   exit (ok ? EXIT_SUCCESS : EXIT_FAILURE);
 626 }