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