1 //==============================================================================================
2 //
3 // This file is part of LiDIA --- a library for computational number theory
4 //
5 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved.
6 //
7 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
8 //
9 //----------------------------------------------------------------------------------------------
10 //
11 // $Id$
12 //
13 // Author : Andreas Mueller (AM), Volker Mueller (VM)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18 // Description : see the diploma thesis of the first author
19
20
21 #ifdef HAVE_CONFIG_H
22 # include "config.h"
23 #endif
24 #include "LiDIA/rational_factorization.h"
25 #include "LiDIA/timer.h"
26 #include "LiDIA/random_generator.h"
27
28
29
30 #ifdef LIDIA_NAMESPACE
31 namespace LiDIA {
32 #endif
33
34
35
36 const unsigned int rational_factorization::ecm_params[29][3] = {
37 {151, 1009, 3},
38 {211, 2003, 4},
39 {283, 3203, 5},
40 {349, 4297, 5},
41 {411, 8861, 6},
42 {659, 17981, 8},
43 {997, 31489, 10},
44 {1439, 54617, 13},
45 {2111, 89501, 17},
46 {3067, 144961, 21},
47 {4391, 245719, 27},
48 {6287, 382919, 33},
49 {8839, 577721, 42},
50 {12241, 845951, 54},
51 {16339, 1285633, 68},
52 {21977, 1788863, 88},
53 {29860, 2508269, 111},
54 {41011, 3550163, 137},
55 {53214, 5139913, 172},
56 {71537, 6996393, 215},
57 {97499, 9650917, 261},
58 {125653, 13767227, 325},
59 {168273, 18372001, 398},
60 {223559, 24108067, 489},
61 {280891, 33036147, 612},
62 {376513, 43533099, 736},
63 {489249, 54671451, 909},
64 {631461, 77351879, 1088},
65 {830049, 97901685, 1315}
66 };
67
68
69
70 unsigned int rational_factorization::
ecm_read_max(int stell)71 ecm_read_max(int stell) // returns the limit for prime numbers
72 { // to execute the given strategy
73 return (ecm_params[stell-6][1]);
74 }
75
76
77
78 int rational_factorization::
ecm_job_planing(int strat[30],int buf[30][5])79 ecm_job_planing(int strat[30], int buf[30][5])
80 { // read the parameters for ECM for
81 // a given strategy
82 int i = 0, j = strat[0];
83 unsigned int job = 1;
84
85 if (info) {
86 std::cout << "Strategy for ECM : \n\n";
87 std::cout << " # | | prime limits | #curves \n";
88 std::cout << " digits | curves | step 1 | cont | total \n";
89 std::cout << "-------------------------------------------------------\n";
90 std::cout.flush();
91 }
92
93 while (j) {
94 j -= 6;
95
96 buf[i][0] = ecm_params[j][0];
97 buf[i][1] = buf[i][0];
98 buf[i][2] = ecm_params[j][1];
99 buf[i][3] = ecm_params[j][2];
100 buf[i][4] = buf[i][3];
101
102 if (info) printf(" %2d | %4d -", strat[i], job);
103
104 job += buf[i][3];
105
106 if (info)
107 printf(" %4d |%9d | %8d | %4d\n", job-1, buf[i][0], buf[i][2], buf[i][3]);
108
109 j = strat[++i];
110 }
111
112 if (info) {
113 fflush(stdout);
114 std::cout << "========================================================\n";
115 std::cout.flush();
116 }
117
118 return(i);
119 }
120
121
122
mecfind(ec_curve & me,ec_point_M & mp,bigint & factor,unsigned long B,unsigned long C,ecm_primes & prim)123 void mecfind(ec_curve & me, ec_point_M & mp, bigint & factor, unsigned long B, unsigned long C, ecm_primes & prim)
124 {
125 // first step of ECM:
126 // compute mp <-- K * mp
127 // K is the product of all primes powers up to B
128 unsigned long p;
129
130 bigmod a_2, m1, inv_4(4L);
131
132 long k, count = 0;
133 unsigned long sqrt_C = static_cast<long>(std::sqrt (static_cast<double>(C)));
134
135 inv_4.invert(0); // precompute the constant a_2 = (me.a+2)/4
136 add(m1, me.A(), 2L); // for the given curve me
137 multiply(a_2, inv_4, m1);
138
139 for (k = 0; k < 31; k++)
140 multiply_by_2(mp, mp, a_2); // mp = 2^30 * mp
141
142 // successive computation of mp = p^k * mp
143 // for all primes p = 3, ..., B and
144 // k integer with p^k <= C && p^(k+1) > C
145 while ((p = prim.getprimes()) <= sqrt_C) {
146 // loop: exponent k > 1
147 k = static_cast<long>(std::log(static_cast<double>(C))/std::log(static_cast<double>(p)));
148 multiply(mp, mp, a_2, p, power(p, k));
149 }
150
151 do {
152 // loop: exponent k = 1
153 multiply(mp, mp, a_2, p, p);
154
155 if (count++ > 1000) {
156 m1.assign(mp.Z());
157 factor = m1.invert(1);
158 if (!factor.is_one())
159 return;
160
161 count = 0;
162 multiply(m1, m1, mp.X());
163 mp.assign(m1, bigmod(1L));
164 }
165 } while ((p = prim.getprimes()) < B);
166
167 m1.assign(mp.Z());
168 factor = m1.invert(1);
169 if (!factor.is_one())
170 return;
171
172 multiply(m1, m1, mp.X());
173 mp.assign(m1, bigmod(1L));
174 }
175
176
177
cont(ec_point_W & wp,bigint & factor,unsigned long B,unsigned long D,ecm_primes & prim)178 void cont(ec_point_W & wp, bigint & factor, unsigned long B, unsigned long D, ecm_primes & prim)
179 {
180 register long w, u, v, count;
181
182 unsigned long p;
183
184 bigmod * x, h, c(1L);
185
186 ec_point_W wq;
187
188 // compute the number of babysteps w
189 // and an initial value v for then giantsteps
190 u = static_cast<long>(std::ceil(std::sqrt(static_cast<double>(D))));
191
192 w = (u / 30) * 30;
193 if (w + 30 - u < u - w)
194 w += 30;
195 D = w * w;
196 v = static_cast<long>(B) / w;
197
198 x = new bigmod[w+2];
199
200 // precompute babysteps:
201 // k * wp for k = 1, .., w
202
203 for (count = 1; count <= w; count += 2) {
204 add(wq, wp, wq, factor);
205 if (!factor.is_one()) {
206 delete [] x;
207 return;
208 }
209
210 if (gcd(w, count) == 1)
211 x[count] = wq.X();
212
213 add(wq, wp, wq, factor);
214 if (!factor.is_one()) {
215 delete [] x;
216 return;
217 }
218 }
219
220 // computed: wq = w * wp
221 // initialization for giantsteps:
222 // wp = v * wq = v * w * wp
223
224 multiply(wp, wq, v, factor);
225 if (!factor.is_one()) {
226 delete [] x;
227 return;
228 }
229
230 count = 1;
231 p = prim.getprimes();
232 // giantsteps
233 while (p < D && p > 1) {
234 if ((u = v*w - p) < 1) {
235 v++;
236 add(wp, wp, wq, factor);
237 if (!factor.is_one()) {
238 delete [] x;
239 return;
240 }
241
242 u = v*w - p;
243 if (u < 0)
244 u += w;
245 }
246
247 subtract(h, wp.X(), x[u]);
248
249 multiply(c, c, h);
250
251 if (count++ >= 500) {
252 count = 1;
253 factor = c.invert(1);
254 if (!factor.is_one()) {
255 delete [] x;
256 return;
257 }
258 }
259 p = prim.getprimes();
260 }
261
262 factor = c.invert(1);
263 delete [] x;
264 }
265
266
267
268 void rational_factorization::
ecm(lidia_size_t index,int ecm_restlength,int jobs_avail,int job_buffer[30][5],ecm_primes & prim)269 ecm(lidia_size_t index, int ecm_restlength, int jobs_avail, int job_buffer[30][5], ecm_primes & prim)
270 {
271 bigint N (factors[index].single_base);
272 int exp = factors[index].single_exponent;
273 lidia_size_t n = no_of_comp();
274
275 unsigned int B, C, D;
276 int uc, current_length, count = 0, m, job_nr = 0, job_grp = 0, qs_2;
277
278 char *n_string = new char[(N.bit_length() / 3) + 10];
279
280 random_generator rg;
281
282 current_length = bigint_to_string(N, n_string);
283 uc = ((current_length-30) / 3) << 1;
284
285 if (current_length <= ecm_restlength && uc <= 0) {
286 if (info) {
287 std::cout << "\nUsing MPQS for further factorization ....\n";
288 std::cout.flush();
289 }
290
291 delete [] n_string;
292 return;
293 }
294
295 if ((factors[index].factor_state >> 2) > 0) {
296 int params_for_computed = ecm_params[(factors[index].factor_state >> 2)-6][0];
297
298 while (job_buffer[job_grp][0] < params_for_computed)
299 job_grp ++;
300
301 jobs_avail -= job_grp;
302 }
303
304 B = job_buffer[job_grp][0]; // first job
305 C = job_buffer[job_grp][1];
306 D = job_buffer[job_grp][2];
307
308 bigint old_modul = bigmod::modulus();
309 bigmod::set_modulus(N);
310
311 rf_single_factor fact;
312 bigint factor, Q, R;
313 factor.assign_one();
314
315 ec_curve ell_curve;
316 ec_point_M mp;
317 ec_point_W wp;
318
319 timer ti;
320 qs_2 = static_cast<int>(zeitqs(current_length)) >> 1;
321
322 ti.start_timer();
323
324 while (!N.is_one() && jobs_avail > 0) {
325 while (jobs_avail > 0) {
326 prim.resetprimes(1);
327
328 if (job_buffer[job_grp][3] < 1) {
329 if (--jobs_avail <= 0)
330 break;
331
332 job_grp++;
333 B = job_buffer[job_grp][0]; // next job
334 C = job_buffer[job_grp][1];
335 D = job_buffer[job_grp][2];
336 }
337
338 job_buffer[job_grp][3]--;
339
340 if (info) {
341 std::cout << "\n" << ++job_nr << ".curve ";
342 std::cout.flush();
343 }
344
345 rg >> m;
346
347 mecgen16(ell_curve, mp, factor, m);
348 if (!factor.is_one()) break;
349
350 mecfind(ell_curve, mp, factor, B, C, prim);
351 if (!factor.is_one()) break;
352
353 factor = trans(ell_curve, wp, mp);
354 if (!factor.is_one()) break;
355
356 cont(wp, factor, B, D, prim);
357 if (!factor.is_one()) break;
358
359 ti.stop_timer();
360 ti.cont_timer();
361
362 if ((current_length <= ecm_restlength && --uc <= 0) ||
363 ((ti.user_time()+ti.sys_time())/100 > qs_2)) {
364 fact.factor_state = rf_single_factor::not_prime;
365
366 fact.single_base = N;
367 fact.single_exponent = exp;
368 factors[index] = fact;
369
370 delete [] n_string;
371 if (!old_modul.is_zero()) bigmod::set_modulus(old_modul);
372 return;
373 }
374 }
375
376 if (jobs_avail > 0) {
377 if (!(factor.is_one() || factor.is_zero() || factor == N)) {
378 div_rem(Q, R, N, factor);
379
380 while (R.is_zero()) {
381 count++;
382 N.assign(Q);
383 div_rem(Q, R, N, factor);
384 }
385
386 if (count) {
387 if (info) {
388 if (count > 1)
389 std::cout << "\nfactor " << factor << " ^ " << count;
390 else
391 std::cout << "\nfactor " << factor;
392 std::cout.flush();
393 }
394
395 if (is_prime(factor, 8))
396 fact.factor_state = rf_single_factor::prime;
397 else
398 fact.factor_state = rf_single_factor::not_prime;
399
400 fact.single_base = factor;
401 fact.single_exponent = count * exp;
402
403 count = 0; // ts
404
405 if (N.is_one()) {
406 factors[index] = fact;
407 break;
408 }
409 else
410 factors[n++] = fact;
411
412 if (is_prime(N, 8)) {
413 fact.factor_state = rf_single_factor::prime;
414
415 fact.single_base = N;
416 fact.single_exponent = exp;
417 factors[index] = fact;
418
419 if (info) {
420 std::cout << "\nprime number " << N << "\n";
421 std::cout.flush();
422 }
423
424 N.assign_one();
425 break;
426 }
427
428 current_length = bigint_to_string(N, n_string);
429 uc = ((current_length-30) / 3) << 1;
430
431 qs_2 = static_cast<int>(zeitqs(current_length)) >> 1;
432 if ((current_length <= ecm_restlength && uc <= 0) ||
433 ((ti.user_time()+ti.sys_time())/100 > qs_2)) {
434 fact.factor_state = rf_single_factor::not_prime;
435 fact.single_base = N;
436 fact.single_exponent = exp;
437 factors[index] = fact;
438
439 delete [] n_string;
440 if (!old_modul.is_zero()) bigmod::set_modulus(old_modul);
441 return;
442 }
443
444 bigmod::set_modulus(N);
445
446 }
447 }
448 }
449 }
450
451
452 if (!N.is_one()) {
453 int i = 0;
454
455 unsigned int table = job_buffer[job_grp][0];
456
457 while (ecm_params[i][0] != table)
458 i++;
459
460 i += 6;
461
462 fact.single_base = N;
463 fact.single_exponent = exp;
464 fact.factor_state = (rf_single_factor::not_prime) | (i << 2);
465 factors[index] = fact;
466 }
467
468 delete [] n_string;
469 if (!old_modul.is_zero()) bigmod::set_modulus(old_modul);
470 }
471
472
473
474 void rational_factorization::
ecm(lidia_size_t index,int jobs_avail,int job_buffer[30][5],ecm_primes & prim)475 ecm(lidia_size_t index, int jobs_avail, int job_buffer[30][5], ecm_primes & prim)
476 {
477 bigint N (factors[index].single_base);
478 int exp = factors[index].single_exponent;
479 lidia_size_t n = no_of_comp();
480
481 unsigned int B, C, D;
482 int count = 0, m, job_nr = 0, job_grp = 0;
483
484 random_generator rg;
485
486 if ((factors[index].factor_state >> 2) > 0) {
487 int params_for_computed = ecm_params[(factors[index].factor_state >> 2)-6][0];
488
489 while (job_buffer[job_grp][0] < params_for_computed)
490 job_grp ++;
491
492 jobs_avail -= job_grp;
493 }
494
495 B = job_buffer[job_grp][0]; // first job
496 C = job_buffer[job_grp][1];
497 D = job_buffer[job_grp][2];
498
499 bigint old_modul = bigmod::modulus();
500 bigmod::set_modulus(N);
501
502 rf_single_factor fact;
503 bigint factor, Q, R;
504 factor.assign_one();
505
506 ec_curve ell_curve;
507 ec_point_M mp;
508 ec_point_W wp;
509
510 while (!N.is_one() && jobs_avail > 0) {
511 while (jobs_avail > 0) {
512 prim.resetprimes(1);
513
514 if (info) {
515 std::cout << "\n" << ++job_nr << ".curve ";
516 std::cout.flush();
517 }
518
519 if (job_buffer[job_grp][3] <= 0) {
520 if (--jobs_avail <= 0)
521 break;
522
523 job_grp++;
524 B = job_buffer[job_grp][0]; // next job
525 C = job_buffer[job_grp][1];
526 D = job_buffer[job_grp][2];
527 }
528 else
529 job_buffer[job_grp][3]--;
530
531 rg >> m;
532
533 mecgen16(ell_curve, mp, factor, m);
534 if (!factor.is_one()) break;
535
536 mecfind(ell_curve, mp, factor, B, C, prim);
537 if (!factor.is_one()) break;
538
539 factor = trans(ell_curve, wp, mp);
540 if (!factor.is_one()) break;
541
542 cont(wp, factor, B, D, prim);
543 if (!factor.is_one()) break;
544 }
545
546 if (jobs_avail > 0) {
547 if (!(factor.is_one() || factor.is_zero() || factor == N)) {
548 div_rem(Q, R, N, factor);
549
550 while (R.is_zero()) {
551 count++;
552 N.assign(Q);
553 div_rem(Q, R, N, factor);
554 }
555
556 if (count) {
557 if (info) {
558 if (count > 1)
559 std::cout << "\nfactor " << factor << " ^ " << count;
560 else
561 std::cout << "\nfactor " << factor;
562 std::cout.flush();
563 }
564
565 if (is_prime(factor, 8))
566 fact.factor_state = rf_single_factor::prime;
567 else
568 fact.factor_state = rf_single_factor::not_prime;
569
570 fact.single_base = factor;
571 fact.single_exponent = count * exp;
572
573 count = 0; // ts
574
575 if (N.is_one()) {
576 factors[index] = fact;
577 break;
578 }
579 else
580 factors[n++] = fact;
581
582 if (is_prime(N, 8)) {
583 fact.factor_state = rf_single_factor::prime;
584
585 fact.single_base = N;
586 fact.single_exponent = exp;
587 factors[index] = fact;
588
589 if (info) {
590 std::cout << "\nprime number " << N << "\n";
591 std::cout.flush();
592 }
593
594 N.assign(1);
595 break;
596 }
597
598 count = 0;
599 bigmod::set_modulus(N);
600 }
601 }
602 }
603 }
604
605 if (!N.is_one()) {
606 int i = 0;
607
608 unsigned int table = job_buffer[job_grp][0];
609
610 while (ecm_params[i][0] != table)
611 i++;
612
613 i += 6;
614
615 fact.single_base = N;
616 fact.single_exponent = exp;
617 fact.factor_state = (rf_single_factor::not_prime) | (i << 2);
618 factors[index] = fact;
619 }
620
621 if (!old_modul.is_zero()) bigmod::set_modulus(old_modul);
622 }
623
624
625
626 void rational_factorization::
ecm_comp(lidia_size_t index,int upper_bound,int lower_bound,int step)627 ecm_comp(lidia_size_t index, int upper_bound, int lower_bound, int step)
628 {
629 if ((index< 0) || (index >= no_of_comp())) {
630 lidia_error_handler("rational_factorization", "ecm_comp::index out of range");
631 return;
632 }
633
634 if (rf_single_factor::state(factors[index].factor_state) == rf_single_factor::prime) {
635 if (info)
636 std::cout << "index " << index << " : prime number " << factors[index].single_base << "\n";
637 return;
638 }
639
640 if (rf_single_factor::state(factors[index].factor_state) == rf_single_factor::dont_know) {
641 if (is_prime(factors[index].single_base, 8)) {
642 if (info)
643 std::cout << "index " << index << " : prime number " << factors[index].single_base << "\n";
644
645 factors[index].factor_state = rf_single_factor::prime;
646 return;
647 }
648 else
649 factors[index].factor_state = rf_single_factor::not_prime;
650 }
651
652 if (step <= 0) {
653 lidia_warning_handler("rational_factorization", "ecm_comp::step <= 0");
654 step = 3;
655 }
656
657 if (upper_bound > 34) {
658 lidia_warning_handler("rational_factorization", "ecm_comp::upper_bound > 34");
659 upper_bound = 34;
660 }
661
662 if (lower_bound < 6) {
663 lidia_warning_handler("rational_factorization", "ecm_comp::lower_bound < 6");
664 lower_bound = 6;
665 }
666
667 if (lower_bound > upper_bound) {
668 lidia_warning_handler("rational_factorization", "ecm_comp::lower_bound > upper_bound");
669 upper_bound = lower_bound;
670 }
671
672 int k, D, jobs_avail, job_buffer[30][5], strategy[30];
673 int n = ((factors[index].single_base).bit_length() / 3) + 10;
674
675 if (upper_bound == 34) {
676 char *n_string;
677 n_string = new char[n];
678 upper_bound = (bigint_to_string(factors[index].single_base, n_string) >> 1) + 1;
679 delete [] n_string;
680 if (upper_bound > 34) upper_bound = 34;
681 }
682
683 strategy[0] = lower_bound;
684 n = lower_bound + step;
685 k = 1;
686
687 while (n < upper_bound) {
688 strategy[k++] = n;
689 n += step;
690 }
691
692 if (lower_bound < upper_bound)
693 strategy[k++] = upper_bound;
694
695 strategy[k] = 0;
696
697 D = ecm_read_max(upper_bound);
698
699 ecm_primes prim(1, static_cast<unsigned int>(D+200), 200000);
700
701 jobs_avail = ecm_job_planing(strategy, job_buffer);
702
703 ecm(index, jobs_avail, job_buffer, prim);
704
705 compose();
706 }
707
708
709
710 void rational_factorization::
ecm(int upper_bound,int lower_bound,int step)711 ecm(int upper_bound, int lower_bound, int step)
712 {
713 if (is_prime_factorization())
714 return;
715 lidia_size_t index, len = no_of_comp();
716
717 int D, k, jobs_avail, job_buffer[30][5], strategy[30];
718 int n = ((factors[len-1].single_base).bit_length() / 3) + 10;
719
720 if (step <= 0) {
721 lidia_warning_handler("rational_factorization", "ecm::step < 0");
722 step = 3;
723 }
724
725 if (upper_bound > 34) {
726 lidia_warning_handler("rational_factorization", "ecm::upper_bound > 34");
727 upper_bound = 34;
728 }
729
730 if (lower_bound < 6) {
731 lidia_warning_handler("rational_factorization", "ecm::lower_bound < 6");
732 lower_bound = 6;
733 }
734
735 if (lower_bound > upper_bound) {
736 lidia_warning_handler("rational_factorization", "ecm::lower_bound > upper_bound");
737 upper_bound = lower_bound;
738 }
739
740 if (upper_bound == 34) {
741 char *n_string;
742 n_string = new char[n];
743 upper_bound = (bigint_to_string(factors[len-1].single_base, n_string) >> 1) + 1;
744 delete [] n_string;
745 if (upper_bound > 34) upper_bound = 34;
746 if (upper_bound < 6) upper_bound = 6;
747 }
748
749 D = ecm_read_max(upper_bound);
750
751 ecm_primes prim(1, D+200, 200000);
752
753 strategy[0] = lower_bound; k = 1;
754 n = lower_bound + step;
755
756 while (n < upper_bound) {
757 strategy[k++] = n;
758 n += step;
759 }
760
761 if (lower_bound < upper_bound) {
762 strategy[k] = upper_bound;
763 strategy[k+1] = 0;
764 }
765 else strategy[k] = 0;
766
767 if (info) {
768 std::cout << "\n ECM with same strategy for each number";
769 std::cout << "\n --------------------------------------\n\n";
770 }
771
772 jobs_avail = ecm_job_planing(strategy, job_buffer);
773
774 for (index = 0; index < len; index++) {
775 if (info) {
776 std::cout << "\nindex " << index;
777 std::cout << "\n--------\n";
778 }
779
780 if (rf_single_factor::state(factors[index].factor_state) == rf_single_factor::prime) {
781 if (info)
782 std::cout << "index " << index << " : prime number " << factors[index].single_base << "\n";
783 continue;
784 }
785
786 if (rf_single_factor::state(factors[index].factor_state) == rf_single_factor::dont_know) {
787 if (is_prime(factors[index].single_base, 8)) {
788 if (info)
789 std::cout << "index " << index << " : prime number " << factors[index].single_base << "\n";
790
791 factors[index].factor_state = rf_single_factor::prime;
792
793 continue;
794 }
795 else
796 factors[index].factor_state = rf_single_factor::not_prime;
797 }
798
799 ecm(index, jobs_avail, job_buffer, prim);
800
801 for (D = 0; D < jobs_avail; D++)
802 job_buffer[D][3] = job_buffer[D][4];
803 }
804
805 compose();
806 }
807
808
809
ecm(const bigint & N,int upper_bound,int lower_bound,int step)810 rational_factorization ecm (const bigint &N, int upper_bound, int lower_bound, int step)
811 {
812 rational_factorization f(N);
813
814 f.ecm(upper_bound, lower_bound, step);
815
816 return f;
817 }
818
819
820
821 #ifdef LIDIA_NAMESPACE
822 } // end of namespace LiDIA
823 #endif
824