Line data Source code
1 : /* The 'P+1' algorithm.
2 :
3 : Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012
4 : Paul Zimmermann and Alexander Kruppa.
5 :
6 : This file is part of the ECM Library.
7 :
8 : The ECM Library is free software; you can redistribute it and/or modify
9 : it under the terms of the GNU Lesser General Public License as published by
10 : the Free Software Foundation; either version 3 of the License, or (at your
11 : option) any later version.
12 :
13 : The ECM Library is distributed in the hope that it will be useful, but
14 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 : License for more details.
17 :
18 : You should have received a copy of the GNU Lesser General Public License
19 : along with the ECM Library; see the file COPYING.LIB. If not, see
20 : http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 : 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22 :
23 : /* References:
24 :
25 : A p+1 Method of Factoring, H. C. Williams, Mathematics of Computation,
26 : volume 39, number 159, pages 225-234, 1982.
27 :
28 : Evaluating recurrences of form X_{m+n} = f(X_m, X_n, X_{m-n}) via
29 : Lucas chains, Peter L. Montgomery, December 1983, revised January 1992. */
30 :
31 : #include <math.h>
32 : #include <stdlib.h>
33 : #include "ecm-impl.h"
34 : #include "getprime_r.h"
35 :
36 : #ifdef HAVE_LIMITS_H
37 : # include <limits.h>
38 : #else
39 : # ifndef ULONG_MAX
40 : # define ULONG_MAX __GMP_ULONG_MAX
41 : # endif
42 : #endif
43 :
44 :
45 : /******************************************************************************
46 : * *
47 : * Stage 1 *
48 : * *
49 : ******************************************************************************/
50 :
51 : /* prime powers are accumulated up to about n^L1 */
52 : #define L1 1
53 :
54 : /* P1 <- V_e(P0), using P, Q as auxiliary variables,
55 : where V_{2k}(P0) = V_k(P0)^2 - 2
56 : V_{2k-1}(P0) = V_k(P0)*V_{k-1}(P0) - P0.
57 : (More generally V_{m+n} = V_m * V_n - V_{m-n}.)
58 : Warning: P1 and P0 may be equal.
59 : Assume e >= 1.
60 : */
61 : static void
62 4248 : pp1_mul (mpres_t P1, mpres_t P0, mpz_t e, mpmod_t n, mpres_t P, mpres_t Q)
63 : {
64 : mp_size_t size_e, i;
65 :
66 : ASSERT (mpz_cmp_ui (e, 1) >= 0);
67 :
68 4248 : if (mpz_cmp_ui (e, 1) == 0)
69 : {
70 31 : mpres_set (P1, P0, n);
71 31 : return;
72 : }
73 :
74 : /* now e >= 2 */
75 4217 : mpz_sub_ui (e, e, 1);
76 4217 : mpres_sqr (P, P0, n);
77 4217 : mpres_sub_ui (P, P, 2, n); /* P = V_2(P0) = P0^2-2 */
78 4217 : mpres_set (Q, P0, n); /* Q = V_1(P0) = P0 */
79 :
80 : /* invariant: (P, Q) = (V_{k+1}(P0), V_k(P0)), start with k=1 */
81 4217 : size_e = mpz_sizeinbase (e, 2);
82 601928 : for (i = size_e - 1; i > 0;)
83 : {
84 597711 : if (ecm_tstbit (e, --i)) /* k -> 2k+1 */
85 : {
86 294226 : if (i) /* Q is not needed for last iteration */
87 : {
88 293941 : mpres_mul (Q, P, Q, n);
89 293941 : mpres_sub (Q, Q, P0, n);
90 : }
91 294226 : mpres_sqr (P, P, n);
92 294226 : mpres_sub_ui (P, P, 2, n);
93 : }
94 : else /* k -> 2k */
95 : {
96 303485 : mpres_mul (P, P, Q, n);
97 303485 : mpres_sub (P, P, P0, n);
98 303485 : if (i) /* Q is not needed for last iteration */
99 : {
100 299553 : mpres_sqr (Q, Q, n);
101 299553 : mpres_sub_ui (Q, Q, 2, n);
102 : }
103 : }
104 : }
105 :
106 4217 : mpres_set (P1, P, n);
107 4217 : mpz_add_ui (e, e, 1); /* recover original value of e */
108 : }
109 :
110 : /* Input: P0 is the initial point (sigma)
111 : n is the number to factor
112 : B1 is the stage 1 bound
113 : B1done: stage 1 was already done up to that limit
114 : go: if <> 1, group order to preload
115 : Output: a is the factor found, or the value at end of stage 1
116 : B1done is set to B1 if stage 1 completed normally,
117 : or to the largest prime processed if interrupted, but never
118 : to a smaller value than B1done was upon function entry.
119 : Return value: non-zero iff a factor was found.
120 : */
121 : static int
122 264 : pp1_stage1 (mpz_t f, mpres_t P0, mpmod_t n, double B1, double *B1done,
123 : mpz_t go, int (*stop_asap)(void), char *chkfilename)
124 : {
125 : double B0, p, q, r, last_chkpnt_p;
126 : mpz_t g;
127 : mpres_t P, Q;
128 : mpres_t R, S, T;
129 264 : int youpi = ECM_NO_FACTOR_FOUND;
130 : unsigned int max_size, size_n;
131 : long last_chkpnt_time;
132 : prime_info_t prime_info;
133 :
134 264 : mpz_init (g);
135 264 : mpres_init (P, n);
136 264 : mpres_init (Q, n);
137 264 : mpres_init (R, n);
138 264 : mpres_init (S, n);
139 264 : mpres_init (T, n);
140 :
141 264 : B0 = ceil (sqrt (B1));
142 :
143 264 : size_n = mpz_sizeinbase (n->orig_modulus, 2);
144 264 : max_size = L1 * size_n;
145 :
146 264 : if (mpz_cmp_ui (go, 1) > 0)
147 5 : pp1_mul (P0, P0, go, n, P, Q);
148 :
149 : /* suggestion from Peter Montgomery: start with exponent n^2-1,
150 : as factors of Lucas and Fibonacci number are either +/-1 (mod index),
151 : and so is n. Therefore, index will appear as a factor
152 : of n^2-1 and be included in stage 1.
153 : Do this only when n is composite, otherwise all tests with prime
154 : n factor of a Cunningham number will succeed in stage 1.
155 :
156 : As in P-1, for small overhead, use that trick only when lg(n) <= sqrt(B1).
157 : */
158 264 : if ((double) size_n <= B0 &&
159 144 : mpz_probab_prime_p (n->orig_modulus, PROBAB_PRIME_TESTS) == 0)
160 : {
161 26 : mpz_mul (g, n->orig_modulus, n->orig_modulus);
162 26 : mpz_sub_ui (g, g, 1);
163 26 : pp1_mul (P0, P0, g, n, P, Q);
164 : }
165 :
166 264 : mpz_set_ui (g, 1);
167 :
168 264 : last_chkpnt_p = 2.;
169 264 : last_chkpnt_time = cputime ();
170 : /* first loop through small primes <= sqrt(B1) */
171 264 : prime_info_init (prime_info);
172 95548 : for (p = 2.0; p <= B0; p = (double) getprime_mt (prime_info))
173 : {
174 299219 : for (q = 1, r = p; r <= B1; r *= p)
175 203935 : if (r > *B1done) q *= p;
176 95284 : mpz_mul_d (g, g, q, Q);
177 95284 : if (mpz_sizeinbase (g, 2) >= max_size)
178 : {
179 3953 : pp1_mul (P0, P0, g, n, P, Q);
180 3953 : mpz_set_ui (g, 1);
181 3953 : if (stop_asap != NULL && (*stop_asap) ())
182 : {
183 0 : interrupt:
184 0 : outputf (OUTPUT_NORMAL, "Interrupted at prime %.0f\n", p);
185 0 : if (p > *B1done)
186 0 : *B1done = p;
187 0 : goto clear_and_exit;
188 : }
189 : }
190 : }
191 :
192 264 : pp1_mul (P0, P0, g, n, P, Q);
193 :
194 : /* All primes sqrt(B1) < p <= B1 appear with exponent 1. All primes <= B1done
195 : are already included with exponent at least 1, so it's safe to skip
196 : ahead to B1done+1. */
197 :
198 2033517384 : while (p <= *B1done)
199 2033517120 : p = (double) getprime_mt (prime_info);
200 :
201 : /* then all primes > sqrt(B1) and taken with exponent 1 */
202 115425083 : for (; p <= B1; p = (double) getprime_mt (prime_info))
203 : {
204 115424819 : pp1_mul_prac (P0, (ecm_uint) p, n, P, Q, R, S, T);
205 :
206 115424819 : if (stop_asap != NULL && (*stop_asap) ())
207 0 : goto interrupt;
208 115811164 : if (chkfilename != NULL && p > last_chkpnt_p + 10000. &&
209 386345 : elltime (last_chkpnt_time, cputime ()) > CHKPNT_PERIOD)
210 : {
211 0 : writechkfile (chkfilename, ECM_PP1, p, n, NULL, P0, NULL, NULL);
212 0 : last_chkpnt_p = p;
213 0 : last_chkpnt_time = cputime ();
214 : }
215 : }
216 :
217 : /* If stage 1 finished normally, p is the smallest prime >B1 here.
218 : In that case, set to B1 */
219 264 : if (p > B1)
220 264 : p = B1;
221 :
222 264 : if (p > *B1done)
223 264 : *B1done = p;
224 :
225 264 : mpres_sub_ui (P, P0, 2, n);
226 264 : mpres_gcd (f, P, n);
227 264 : youpi = mpz_cmp_ui (f, 1);
228 :
229 264 : clear_and_exit:
230 264 : if (chkfilename != NULL)
231 5 : writechkfile (chkfilename, ECM_PP1, p, n, NULL, P0, NULL, NULL);
232 264 : prime_info_clear (prime_info); /* free the prime table */
233 264 : mpres_clear (Q, n);
234 264 : mpres_clear (R, n);
235 264 : mpres_clear (S, n);
236 264 : mpres_clear (T, n);
237 264 : mpz_clear (g);
238 264 : mpres_clear (P, n);
239 :
240 264 : return youpi;
241 : }
242 :
243 : /* checks if the factor p was found by P+1 or P-1 (when prime).
244 : a is the initial seed.
245 : */
246 : static void
247 214 : pp1_check_factor (mpz_t a, mpz_t p)
248 : {
249 214 : if (mpz_probab_prime_p (p, PROBAB_PRIME_TESTS))
250 : {
251 214 : mpz_mul (a, a, a);
252 214 : mpz_sub_ui (a, a, 4);
253 214 : if (mpz_jacobi (a, p) == 1)
254 4 : outputf (OUTPUT_NORMAL, "[factor found by P-1]\n");
255 : }
256 214 : }
257 :
258 :
259 : /******************************************************************************
260 : * *
261 : * Williams P+1 *
262 : * *
263 : ******************************************************************************/
264 :
265 : /* Input: p is the initial generator (x0), if 0 generate it at random.
266 : n is the number to factor (assumed to be odd)
267 : B1 is the stage 1 bound
268 : B2 is the stage 2 bound
269 : k is the number of blocks for stage 2
270 : verbose is the verbosity level
271 : Output: f is the factor found, p is the residue at end of stage 1
272 : Return value: non-zero iff a factor is found (1 for stage 1, 2 for stage 2)
273 : */
274 : int
275 269 : pp1 (mpz_t f, mpz_t p, mpz_t n, mpz_t go, double *B1done, double B1,
276 : mpz_t B2min_parm, mpz_t B2_parm, unsigned long k,
277 : int verbose, int repr, int use_ntt, FILE *os, FILE *es,
278 : char *chkfilename, char *TreeFilename, double maxmem,
279 : gmp_randstate_t rng, int (*stop_asap)(void))
280 : {
281 269 : int youpi = ECM_NO_FACTOR_FOUND;
282 : long st;
283 : mpres_t a;
284 : mpmod_t modulus;
285 : mpz_t B2min, B2; /* Local B2, B2min to avoid changing caller's values */
286 : faststage2_param_t faststage2_params;
287 269 : int twopass = 0;
288 : mpz_t p0;
289 :
290 : ASSERT (mpz_divisible_ui_p (n, 2) == 0);
291 :
292 269 : set_verbose (verbose);
293 269 : ECM_STDOUT = (os == NULL) ? stdout : os;
294 269 : ECM_STDERR = (es == NULL) ? stdout : es;
295 :
296 269 : st = cputime ();
297 :
298 269 : if (mpz_cmp_ui (p, 0) == 0)
299 50 : pp1_random_seed (p, n, rng);
300 :
301 269 : mpz_init_set (B2min, B2min_parm);
302 269 : mpz_init_set (B2, B2_parm);
303 :
304 : /* Set default B2. See ecm.c for comments */
305 269 : if (ECM_IS_DEFAULT_B2(B2))
306 40 : mpz_set_d (B2, pow (B1 * PP1FS2_COST, PM1FS2_DEFAULT_B2_EXPONENT));
307 :
308 : /* set B2min */
309 269 : if (mpz_sgn (B2min) < 0)
310 249 : mpz_set_d (B2min, B1);
311 :
312 : {
313 : long P;
314 269 : const unsigned long lmax = 1UL<<28; /* An upper bound */
315 : unsigned long lmax_NTT, lmax_noNTT;
316 :
317 269 : mpz_init (faststage2_params.m_1);
318 269 : faststage2_params.l = 0;
319 269 : faststage2_params.file_stem = TreeFilename;
320 :
321 : /* Find out what the longest transform length is we can do at all.
322 : If no maxmem is given, the non-NTT can theoretically do any length. */
323 :
324 269 : lmax_NTT = 0;
325 269 : if (use_ntt)
326 : {
327 215 : unsigned long t, t2 = 0;
328 : /* See what transform length that the NTT can handle (due to limited
329 : primes and limited memory) */
330 215 : t = mpzspm_max_len (n);
331 215 : lmax_NTT = MIN (lmax, t);
332 215 : if (maxmem != 0.)
333 : {
334 11 : t = pp1fs2_maxlen (double_to_size (maxmem), n, use_ntt, 0);
335 11 : t = MIN (t, lmax_NTT);
336 : /* Maybe the two pass variant lets us use a longer transform */
337 11 : t2 = pp1fs2_maxlen (double_to_size (maxmem), n, use_ntt, 1);
338 11 : t2 = MIN (t2, lmax_NTT);
339 11 : if (t2 > t)
340 : {
341 1 : t = t2;
342 1 : twopass = 1;
343 : }
344 11 : lmax_NTT = t;
345 : }
346 215 : outputf (OUTPUT_DEVVERBOSE, "NTT can handle lmax <= %lu\n", lmax_NTT);
347 : }
348 :
349 : /* See what transform length that the non-NTT code can handle */
350 269 : lmax_noNTT = lmax;
351 269 : if (maxmem != 0.)
352 : {
353 : unsigned long t;
354 21 : t = pp1fs2_maxlen (double_to_size (maxmem), n, 0, 0);
355 21 : lmax_noNTT = MIN (lmax_noNTT, t);
356 21 : outputf (OUTPUT_DEVVERBOSE, "non-NTT can handle lmax <= %lu\n",
357 : lmax_noNTT);
358 : }
359 :
360 269 : P = choose_P (B2min, B2, MAX(lmax_noNTT, lmax_NTT), k,
361 : &faststage2_params, B2min, B2, use_ntt, ECM_PP1);
362 269 : if (P == ECM_ERROR)
363 : {
364 5 : outputf (OUTPUT_ERROR,
365 : "Error: cannot choose suitable P value for your stage 2 "
366 : "parameters.\nTry a shorter B2min,B2 interval.\n");
367 5 : mpz_clear (faststage2_params.m_1);
368 5 : return ECM_ERROR;
369 : }
370 :
371 : /* See if the selected parameters let us use NTT or not */
372 264 : if (faststage2_params.l > lmax_NTT)
373 49 : use_ntt = 0;
374 :
375 264 : if (maxmem != 0.)
376 : {
377 : unsigned long MB;
378 : char *s;
379 21 : if (!use_ntt)
380 10 : s = "out";
381 11 : else if (twopass)
382 1 : s = " two pass";
383 : else
384 10 : s = " one pass";
385 :
386 21 : MB = pp1fs2_memory_use (faststage2_params.l, n, use_ntt, twopass)
387 : / 1048576;
388 21 : outputf (OUTPUT_VERBOSE, "Using lmax = %lu with%s NTT which takes "
389 : "about %luMB of memory\n", faststage2_params.l, s, MB);
390 : }
391 : }
392 :
393 : /* Print B1, B2, polynomial and x0 */
394 264 : print_B1_B2_poly (OUTPUT_NORMAL, ECM_PP1, B1, *B1done, B2min_parm, B2min,
395 : B2, 1, p, 0, 0, NULL, 0, 0);
396 :
397 : /* If we do a stage 2, print its parameters */
398 264 : if (mpz_cmp (B2, B2min) >= 0)
399 : {
400 : /* can't mix 64-bit types and mpz_t on win32 for some reason */
401 244 : outputf (OUTPUT_VERBOSE, "P = %" PRIu64 ", l = %" PRIu64 ", "
402 : "s_1 = %" PRIu64 ", k = s_2 = %" PRIu64,
403 : faststage2_params.P, faststage2_params.l,
404 : faststage2_params.s_1,faststage2_params.s_2);
405 244 : outputf (OUTPUT_VERBOSE, ", m_1 = %Zd\n",
406 : faststage2_params.m_1);
407 : }
408 :
409 264 : if (test_verbose (OUTPUT_VERBOSE))
410 : {
411 30 : if (mpz_sgn (B2min_parm) >= 0)
412 : {
413 0 : outputf (OUTPUT_VERBOSE,
414 : "Can't compute success probabilities for B1 <> B2min\n");
415 : }
416 : else
417 : {
418 30 : rhoinit (256, 10);
419 : /* If x0 is chosen randomly, the resulting group order will behave,
420 : on average, like for P-1, thus we use the same code as for P-1. */
421 30 : print_prob (B1, B2, 0, k, 1, go);
422 : }
423 : }
424 :
425 264 : mpmod_init (modulus, n, repr);
426 264 : mpres_init (a, modulus);
427 264 : mpres_set_z (a, p, modulus);
428 :
429 : /* since pp1_mul_prac takes an ecm_uint, we have to check
430 : that B1 <= ECM_UINT_MAX */
431 264 : if (B1 > (double) ECM_UINT_MAX)
432 : {
433 0 : outputf (OUTPUT_ERROR, "Error, maximal step1 bound for P+1 is %lu\n",
434 : ECM_UINT_MAX);
435 0 : youpi = ECM_ERROR;
436 0 : goto clear_and_exit;
437 : }
438 :
439 264 : if (B1 > *B1done || mpz_cmp_ui (go, 1) > 0)
440 264 : youpi = pp1_stage1 (f, a, modulus, B1, B1done, go, stop_asap,
441 : chkfilename);
442 :
443 264 : outputf (OUTPUT_NORMAL, "Step 1 took %ldms\n", elltime (st, cputime ()));
444 264 : if (test_verbose (OUTPUT_RESVERBOSE))
445 : {
446 : mpz_t t;
447 :
448 15 : mpz_init (t);
449 15 : mpres_get_z (t, a, modulus);
450 15 : outputf (OUTPUT_RESVERBOSE, "x=%Zd\n", t);
451 15 : mpz_clear (t);
452 : }
453 :
454 264 : mpz_init_set (p0, p);
455 :
456 : /* store in p the residue at end of stage 1 */
457 264 : mpres_get_z (p, a, modulus);
458 :
459 264 : if (stop_asap != NULL && (*stop_asap) ())
460 0 : goto clear_and_exit;
461 :
462 264 : if (youpi == ECM_NO_FACTOR_FOUND && mpz_cmp (B2, B2min) >= 0)
463 : {
464 234 : if (use_ntt)
465 187 : youpi = pp1fs2_ntt (f, a, modulus, &faststage2_params, twopass);
466 : else
467 47 : youpi = pp1fs2 (f, a, modulus, &faststage2_params);
468 : }
469 :
470 264 : if (youpi > 0 && test_verbose (OUTPUT_NORMAL))
471 214 : pp1_check_factor (p0, f); /* tell user if factor was found by P-1 */
472 :
473 264 : mpz_clear (p0);
474 :
475 264 : clear_and_exit:
476 264 : mpres_clear (a, modulus);
477 264 : mpmod_clear (modulus);
478 264 : mpz_clear (faststage2_params.m_1);
479 264 : mpz_clear (B2);
480 264 : mpz_clear (B2min);
481 :
482 264 : return youpi;
483 : }
|
