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