1/* Factoring with Pollard's rho method.
2
3Copyright 1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2009
4Free Software Foundation, Inc.
5
6This program is free software; you can redistribute it and/or modify it under
7the terms of the GNU General Public License as published by the Free Software
8Foundation; either version 3 of the License, or (at your option) any later
9version.
10
11This program is distributed in the hope that it will be useful, but WITHOUT ANY
12WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
13PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14
15You should have received a copy of the GNU General Public License along with
16this 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
25int flag_verbose = 0;
26
27static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
28
29void
30factor_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
104void
105factor_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
137void
138factor_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, 25))
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, 25))
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
272void
273factor (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, 25))
300	mpz_out_str (stdout, 10, t);
301      else
302	factor_using_pollard_rho (t, 1L, p);
303    }
304}
305
306int
307main (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}
376