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