demos/factorize.c

   1 /* Factoring with Pollard's rho method.
   2
   3 Copyright 1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2009
   4 Free Software Foundation, Inc.
   5
   6 This program is free software; you can redistribute it and/or modify it under
   7 the terms of the GNU General Public License as published by the Free Software
   8 Foundation; either version 3 of the License, or (at your option) any later
   9 version.
  10
  11 This program is distributed in the hope that it will be useful, but WITHOUT ANY
  12 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
  13 PARTICULAR PURPOSE.  See the GNU General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License along with
  16 this program.  If not, see http://www.gnu.org/licenses/.  */
  17
  18
  19 #include <stdlib.h>
  20 #include <stdio.h>
  21 #include <string.h>
  22
  23 #include "gmp.h"
  24
  25 int flag_verbose = 0;
  26
  27 static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
  28
  29 void
  30 factor_using_division (mpz_t t, unsigned int limit)
  31 {
  32   mpz_t q, r;
  33   unsigned long int f;
  34   int ai;
  35   unsigned *addv = add;
  36   unsigned int failures;
  37
  38   if (flag_verbose > 0)
  39     {
  40       printf ("[trial division (%u)] ", limit);
  41       fflush (stdout);
  42     }
  43
  44   mpz_init (q);
  45   mpz_init (r);
  46
  47   f = mpz_scan1 (t, 0);
  48   mpz_div_2exp (t, t, f);
  49   while (f)
  50     {
  51       printf ("2 ");
  52       fflush (stdout);
  53       --f;
  54     }
  55
  56   for (;;)
  57     {
  58       mpz_tdiv_qr_ui (q, r, t, 3);
  59       if (mpz_cmp_ui (r, 0) != 0)
  60         break;
  61       mpz_set (t, q);
  62       printf ("3 ");
  63       fflush (stdout);
  64     }
  65
  66   for (;;)
  67     {
  68       mpz_tdiv_qr_ui (q, r, t, 5);
  69       if (mpz_cmp_ui (r, 0) != 0)
  70         break;
  71       mpz_set (t, q);
  72       printf ("5 ");
  73       fflush (stdout);
  74     }
  75
  76   failures = 0;
  77   f = 7;
  78   ai = 0;
  79   while (mpz_cmp_ui (t, 1) != 0)
  80     {
  81       mpz_tdiv_qr_ui (q, r, t, f);
  82       if (mpz_cmp_ui (r, 0) != 0)
  83         {
  84           f += addv[ai];
  85           if (mpz_cmp_ui (q, f) < 0)
  86             break;
  87           ai = (ai + 1) & 7;
  88           failures++;
  89           if (failures > limit)
  90             break;
  91         }
  92       else
  93         {
  94           mpz_swap (t, q);
  95           printf ("%lu ", f);
  96           fflush (stdout);
  97           failures = 0;
  98         }
  99     }
 100
 101   mpz_clears (q, r, NULL);
 102 }
 103
 104 void
 105 factor_using_division_2kp (mpz_t t, unsigned int limit, unsigned long p)
 106 {
 107   mpz_t r;
 108   mpz_t f;
 109   unsigned int k;
 110
 111   if (flag_verbose > 0)
 112     {
 113       printf ("[trial division (%u)] ", limit);
 114       fflush (stdout);
 115     }
 116
 117   mpz_init (r);
 118   mpz_init_set_ui (f, 2 * p);
 119   mpz_add_ui (f, f, 1);
 120   for (k = 1; k < limit; k++)
 121     {
 122       mpz_tdiv_r (r, t, f);
 123       while (mpz_cmp_ui (r, 0) == 0)
 124         {
 125           mpz_tdiv_q (t, t, f);
 126           mpz_tdiv_r (r, t, f);
 127           mpz_out_str (stdout, 10, f);
 128           fflush (stdout);
 129           fputc (' ', stdout);
 130         }
 131       mpz_add_ui (f, f, 2 * p);
 132     }
 133
 134   mpz_clears (f, r, NULL);
 135 }
 136
 137 void
 138 factor_using_pollard_rho (mpz_t n, unsigned long a, unsigned long p)
 139 {
 140   mpz_t x, x1, y, P;
 141   mpz_t t1, t2;
 142   unsigned long long k, l, i;
 143
 144   if (flag_verbose > 0)
 145     {
 146       printf ("[pollard-rho (%lu)] ", a);
 147       fflush (stdout);
 148     }
 149
 150   mpz_inits (t1, t2, NULL);
 151   mpz_init_set_si (y, 2);
 152   mpz_init_set_si (x, 2);
 153   mpz_init_set_si (x1, 2);
 154   mpz_init_set_ui (P, 1);
 155   k = 1;
 156   l = 1;
 157
 158   while (mpz_cmp_ui (n, 1) != 0)
 159     {
 160       for (;;)
 161         {
 162           do
 163             {
 164               if (p != 0)
 165                 {
 166                   mpz_powm_ui (x, x, p, n);
 167                   mpz_add_ui (x, x, a);
 168                 }
 169               else
 170                 {
 171                   mpz_mul (t1, x, x);
 172                   mpz_mod (x, t1, n);
 173                   mpz_add_ui (x, x, a);
 174                 }
 175
 176               mpz_sub (t1, x1, x);
 177               mpz_mul (t2, P, t1);
 178               mpz_mod (P, t2, n);
 179
 180               if (k % 32 == 1)
 181                 {
 182                   mpz_gcd (t1, P, n);
 183                   if (mpz_cmp_ui (t1, 1) != 0)
 184                     goto factor_found;
 185                   mpz_set (y, x);
 186                 }
 187             }
 188           while (--k != 0);
 189
 190           mpz_gcd (t1, P, n);
 191           if (mpz_cmp_ui (t1, 1) != 0)
 192             goto factor_found;
 193
 194           mpz_set (x1, x);
 195           k = l;
 196           l = 2 * l;
 197           for (i = 0; i < k; i++)
 198             {
 199               if (p != 0)
 200                 {
 201                   mpz_powm_ui (x, x, p, n);
 202                   mpz_add_ui (x, x, a);
 203                 }
 204               else
 205                 {
 206                   mpz_mul (t1, x, x);
 207                   mpz_mod (x, t1, n);
 208                   mpz_add_ui (x, x, a);
 209                 }
 210             }
 211           mpz_set (y, x);
 212         }
 213
 214     factor_found:
 215       do
 216         {
 217           if (p != 0)
 218             {
 219               mpz_powm_ui (y, y, p, n); mpz_add_ui (y, y, a);
 220             }
 221           else
 222             {
 223               mpz_mul (t1, y, y);
 224               mpz_mod (y, t1, n);
 225               mpz_add_ui (y, y, a);
 226             }
 227           mpz_sub (t1, x1, y);
 228           mpz_gcd (t1, t1, n);
 229         }
 230       while (mpz_cmp_ui (t1, 1) == 0);
 231
 232       mpz_divexact (n, n, t1);  /* divide by t1, before t1 is overwritten */
 233
 234       if (!mpz_probab_prime_p (t1, 10))
 235         {
 236           do
 237             {
 238               mp_limb_t a_limb;
 239               mpn_random (&a_limb, (mp_size_t) 1);
 240               a = a_limb;
 241             }
 242           while (a == 0);
 243
 244           if (flag_verbose > 0)
 245             {
 246               printf ("[composite factor--restarting pollard-rho] ");
 247               fflush (stdout);
 248             }
 249           factor_using_pollard_rho (t1, a, p);
 250         }
 251       else
 252         {
 253           mpz_out_str (stdout, 10, t1);
 254           fflush (stdout);
 255           fputc (' ', stdout);
 256         }
 257       mpz_mod (x, x, n);
 258       mpz_mod (x1, x1, n);
 259       mpz_mod (y, y, n);
 260       if (mpz_probab_prime_p (n, 10))
 261         {
 262           mpz_out_str (stdout, 10, n);
 263           fflush (stdout);
 264           fputc (' ', stdout);
 265           break;
 266         }
 267     }
 268
 269   mpz_clears (P, t2, t1, x1, x, y, NULL);
 270 }
 271
 272 void
 273 factor (mpz_t t, unsigned long p)
 274 {
 275   unsigned int division_limit;
 276
 277   if (mpz_sgn (t) == 0)
 278     return;
 279
 280   /* Set the trial division limit according the size of t.  */
 281   division_limit = mpz_sizeinbase (t, 2);
 282   if (division_limit > 1000)
 283     division_limit = 1000 * 1000;
 284   else
 285     division_limit = division_limit * division_limit;
 286
 287   if (p != 0)
 288     factor_using_division_2kp (t, division_limit / 10, p);
 289   else
 290     factor_using_division (t, division_limit);
 291
 292   if (mpz_cmp_ui (t, 1) != 0)
 293     {
 294       if (flag_verbose > 0)
 295         {
 296           printf ("[is number prime?] ");
 297           fflush (stdout);
 298         }
 299       if (mpz_probab_prime_p (t, 10))
 300         mpz_out_str (stdout, 10, t);
 301       else
 302         factor_using_pollard_rho (t, 1L, p);
 303     }
 304 }
 305
 306 int
 307 main (int argc, char *argv[])
 308 {
 309   mpz_t t;
 310   unsigned long p;
 311   int i;
 312
 313   if (argc > 1 && !strcmp (argv[1], "-v"))
 314     {
 315       flag_verbose = 1;
 316       argv++;
 317       argc--;
 318     }
 319   if (argc > 1 && !strcmp (argv[1], "-q"))
 320     {
 321       flag_verbose = -1;
 322       argv++;
 323       argc--;
 324     }
 325
 326   mpz_init (t);
 327   if (argc > 1)
 328     {
 329       p = 0;
 330       for (i = 1; i < argc; i++)
 331         {
 332           if (!strncmp (argv[i], "-Mp", 3))
 333             {
 334               p = atoi (argv[i] + 3);
 335               mpz_set_ui (t, 1);
 336               mpz_mul_2exp (t, t, p);
 337               mpz_sub_ui (t, t, 1);
 338             }
 339           else if (!strncmp (argv[i], "-2kp", 4))
 340             {
 341               p = atoi (argv[i] + 4);
 342               continue;
 343             }
 344           else
 345             {
 346               mpz_set_str (t, argv[i], 0);
 347             }
 348
 349           if (mpz_cmp_ui (t, 0) == 0)
 350             puts ("-");
 351           else
 352             {
 353               factor (t, p);
 354               puts ("");
 355             }
 356         }
 357     }
 358   else
 359     {
 360       for (;;)
 361         {
 362           mpz_inp_str (t, stdin, 0);
 363           if (feof (stdin))
 364             break;
 365           if (flag_verbose >= 0)
 366             {
 367               mpz_out_str (stdout, 10, t); printf (" = ");
 368             }
 369           factor (t, 0);
 370           puts ("");
 371         }
 372     }
 373
 374   exit (0);
 375 }