1 /* Factoring with Pollard's rho method.
2
3 Copyright 1995, 1997-2003, 2005, 2009, 2012, 2015 Free Software
4 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 https://www.gnu.org/licenses/. */
17
18
19 #include <stdlib.h>
20 #include <stdio.h>
21 #include <string.h>
22 #include <inttypes.h>
23
24 #include "gmp.h"
25
26 static unsigned char primes_diff[] = {
27 #define P(a,b,c) a,
28 #include "primes.h"
29 #undef P
30 };
31 #define PRIMES_PTAB_ENTRIES (sizeof(primes_diff) / sizeof(primes_diff[0]))
32
33 int flag_verbose = 0;
34
35 /* Prove primality or run probabilistic tests. */
36 int flag_prove_primality = 1;
37
38 /* Number of Miller-Rabin tests to run when not proving primality. */
39 #define MR_REPS 25
40
41 struct factors
42 {
43 mpz_t *p;
44 unsigned long *e;
45 long nfactors;
46 };
47
48 void factor (mpz_t, struct factors *);
49
50 void
51 factor_init (struct factors *factors)
52 {
53 factors->p = malloc (1);
54 factors->e = malloc (1);
55 factors->nfactors = 0;
56 }
57
58 void
59 factor_clear (struct factors *factors)
60 {
61 int i;
62
63 for (i = 0; i < factors->nfactors; i++)
64 mpz_clear (factors->p[i]);
65
66 free (factors->p);
67 free (factors->e);
68 }
69
70 void
71 factor_insert (struct factors *factors, mpz_t prime)
72 {
73 long nfactors = factors->nfactors;
74 mpz_t *p = factors->p;
75 unsigned long *e = factors->e;
76 long i, j;
77
78 /* Locate position for insert new or increment e. */
79 for (i = nfactors - 1; i >= 0; i--)
80 {
81 if (mpz_cmp (p[i], prime) <= 0)
82 break;
83 }
84
85 if (i < 0 || mpz_cmp (p[i], prime) != 0)
86 {
87 p = realloc (p, (nfactors + 1) * sizeof p[0]);
88 e = realloc (e, (nfactors + 1) * sizeof e[0]);
89
90 mpz_init (p[nfactors]);
91 for (j = nfactors - 1; j > i; j--)
92 {
93 mpz_set (p[j + 1], p[j]);
94 e[j + 1] = e[j];
95 }
96 mpz_set (p[i + 1], prime);
97 e[i + 1] = 1;
98
99 factors->p = p;
100 factors->e = e;
101 factors->nfactors = nfactors + 1;
102 }
103 else
104 {
105 e[i] += 1;
106 }
107 }
108
109 void
110 factor_insert_ui (struct factors *factors, unsigned long prime)
111 {
112 mpz_t pz;
113
114 mpz_init_set_ui (pz, prime);
115 factor_insert (factors, pz);
116 mpz_clear (pz);
117 }
118
119
120 void
121 factor_using_division (mpz_t t, struct factors *factors)
122 {
123 mpz_t q;
124 unsigned long int p;
125 int i;
126
127 if (flag_verbose > 0)
128 {
129 printf ("[trial division] ");
130 }
131
132 mpz_init (q);
133
134 p = mpz_scan1 (t, 0);
135 mpz_fdiv_q_2exp (t, t, p);
136 while (p)
137 {
138 factor_insert_ui (factors, 2);
139 --p;
140 }
141
142 p = 3;
143 for (i = 1; i <= PRIMES_PTAB_ENTRIES;)
144 {
145 if (! mpz_divisible_ui_p (t, p))
146 {
147 p += primes_diff[i++];
148 if (mpz_cmp_ui (t, p * p) < 0)
149 break;
150 }
151 else
152 {
153 mpz_tdiv_q_ui (t, t, p);
154 factor_insert_ui (factors, p);
155 }
156 }
157
158 mpz_clear (q);
159 }
160
161 static int
162 mp_millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
163 mpz_srcptr q, unsigned long int k)
164 {
165 unsigned long int i;
166
167 mpz_powm (y, x, q, n);
168
169 if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
170 return 1;
171
172 for (i = 1; i < k; i++)
173 {
174 mpz_powm_ui (y, y, 2, n);
175 if (mpz_cmp (y, nm1) == 0)
176 return 1;
177 if (mpz_cmp_ui (y, 1) == 0)
178 return 0;
179 }
180 return 0;
181 }
182
183 int
184 mp_prime_p (mpz_t n)
185 {
186 int k, r, is_prime;
187 mpz_t q, a, nm1, tmp;
188 struct factors factors;
189
190 if (mpz_cmp_ui (n, 1) <= 0)
191 return 0;
192
193 /* We have already casted out small primes. */
194 if (mpz_cmp_ui (n, (long) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME) < 0)
195 return 1;
196
197 mpz_inits (q, a, nm1, tmp, NULL);
198
199 /* Precomputation for Miller-Rabin. */
200 mpz_sub_ui (nm1, n, 1);
201
202 /* Find q and k, where q is odd and n = 1 + 2**k * q. */
203 k = mpz_scan1 (nm1, 0);
204 mpz_tdiv_q_2exp (q, nm1, k);
205
206 mpz_set_ui (a, 2);
207
208 /* Perform a Miller-Rabin test, finds most composites quickly. */
209 if (!mp_millerrabin (n, nm1, a, tmp, q, k))
210 {
211 is_prime = 0;
212 goto ret2;
213 }
214
215 if (flag_prove_primality)
216 {
217 /* Factor n-1 for Lucas. */
218 mpz_set (tmp, nm1);
219 factor (tmp, &factors);
220 }
221
222 /* Loop until Lucas proves our number prime, or Miller-Rabin proves our
223 number composite. */
224 for (r = 0; r < PRIMES_PTAB_ENTRIES; r++)
225 {
226 int i;
227
228 if (flag_prove_primality)
229 {
230 is_prime = 1;
231 for (i = 0; i < factors.nfactors && is_prime; i++)
232 {
233 mpz_divexact (tmp, nm1, factors.p[i]);
234 mpz_powm (tmp, a, tmp, n);
235 is_prime = mpz_cmp_ui (tmp, 1) != 0;
236 }
237 }
238 else
239 {
240 /* After enough Miller-Rabin runs, be content. */
241 is_prime = (r == MR_REPS - 1);
242 }
243
244 if (is_prime)
245 goto ret1;
246
247 mpz_add_ui (a, a, primes_diff[r]); /* Establish new base. */
248
249 if (!mp_millerrabin (n, nm1, a, tmp, q, k))
250 {
251 is_prime = 0;
252 goto ret1;
253 }
254 }
255
256 fprintf (stderr, "Lucas prime test failure. This should not happen\n");
257 abort ();
258
259 ret1:
260 if (flag_prove_primality)
261 factor_clear (&factors);
262 ret2:
263 mpz_clears (q, a, nm1, tmp, NULL);
264
265 return is_prime;
266 }
267
268 void
269 factor_using_pollard_rho (mpz_t n, unsigned long a, struct factors *factors)
270 {
271 mpz_t x, z, y, P;
272 mpz_t t, t2;
273 unsigned long long k, l, i;
274
275 if (flag_verbose > 0)
276 {
277 printf ("[pollard-rho (%lu)] ", a);
278 }
279
280 mpz_inits (t, t2, NULL);
281 mpz_init_set_si (y, 2);
282 mpz_init_set_si (x, 2);
283 mpz_init_set_si (z, 2);
284 mpz_init_set_ui (P, 1);
285 k = 1;
286 l = 1;
287
288 while (mpz_cmp_ui (n, 1) != 0)
289 {
290 for (;;)
291 {
292 do
293 {
294 mpz_mul (t, x, x);
295 mpz_mod (x, t, n);
296 mpz_add_ui (x, x, a);
297
298 mpz_sub (t, z, x);
299 mpz_mul (t2, P, t);
300 mpz_mod (P, t2, n);
301
302 if (k % 32 == 1)
303 {
304 mpz_gcd (t, P, n);
305 if (mpz_cmp_ui (t, 1) != 0)
306 goto factor_found;
307 mpz_set (y, x);
308 }
309 }
310 while (--k != 0);
311
312 mpz_set (z, x);
313 k = l;
314 l = 2 * l;
315 for (i = 0; i < k; i++)
316 {
317 mpz_mul (t, x, x);
318 mpz_mod (x, t, n);
319 mpz_add_ui (x, x, a);
320 }
321 mpz_set (y, x);
322 }
323
324 factor_found:
325 do
326 {
327 mpz_mul (t, y, y);
328 mpz_mod (y, t, n);
329 mpz_add_ui (y, y, a);
330
331 mpz_sub (t, z, y);
332 mpz_gcd (t, t, n);
333 }
334 while (mpz_cmp_ui (t, 1) == 0);
335
336 mpz_divexact (n, n, t); /* divide by t, before t is overwritten */
337
338 if (!mp_prime_p (t))
339 {
340 if (flag_verbose > 0)
341 {
342 printf ("[composite factor--restarting pollard-rho] ");
343 }
344 factor_using_pollard_rho (t, a + 1, factors);
345 }
346 else
347 {
348 factor_insert (factors, t);
349 }
350
351 if (mp_prime_p (n))
352 {
353 factor_insert (factors, n);
354 break;
355 }
356
357 mpz_mod (x, x, n);
358 mpz_mod (z, z, n);
359 mpz_mod (y, y, n);
360 }
361
362 mpz_clears (P, t2, t, z, x, y, NULL);
363 }
364
365 void
366 factor (mpz_t t, struct factors *factors)
367 {
368 factor_init (factors);
369
370 if (mpz_sgn (t) != 0)
371 {
372 factor_using_division (t, factors);
373
374 if (mpz_cmp_ui (t, 1) != 0)
375 {
376 if (flag_verbose > 0)
377 {
378 printf ("[is number prime?] ");
379 }
380 if (mp_prime_p (t))
381 factor_insert (factors, t);
382 else
383 factor_using_pollard_rho (t, 1, factors);
384 }
385 }
386 }
387
388 int
389 main (int argc, char *argv[])
390 {
391 mpz_t t;
392 int i, j, k;
393 struct factors factors;
394
395 while (argc > 1)
396 {
397 if (!strcmp (argv[1], "-v"))
398 flag_verbose = 1;
399 else if (!strcmp (argv[1], "-w"))
400 flag_prove_primality = 0;
401 else
402 break;
403
404 argv++;
405 argc--;
406 }
407
408 mpz_init (t);
409 if (argc > 1)
410 {
411 for (i = 1; i < argc; i++)
412 {
413 mpz_set_str (t, argv[i], 0);
414
415 gmp_printf ("%Zd:", t);
416 factor (t, &factors);
417
418 for (j = 0; j < factors.nfactors; j++)
419 for (k = 0; k < factors.e[j]; k++)
420 gmp_printf (" %Zd", factors.p[j]);
421
422 puts ("");
423 factor_clear (&factors);
424 }
425 }
426 else
427 {
428 for (;;)
429 {
430 mpz_inp_str (t, stdin, 0);
431 if (feof (stdin))
432 break;
433
434 gmp_printf ("%Zd:", t);
435 factor (t, &factors);
436
437 for (j = 0; j < factors.nfactors; j++)
438 for (k = 0; k < factors.e[j]; k++)
439 gmp_printf (" %Zd", factors.p[j]);
440
441 puts ("");
442 factor_clear (&factors);
443 }
444 }
445
446 exit (0);
447 }