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