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: mpqs.cc,v 2.17 2006/03/06 12:08:36 lidiaadm Exp $
12 //
13 // Author : Volker Mueller (VM), based on an implementation
14 // of Thomas Sosnowski
15 // Changes : See CVS log
16 //
17 // CAVEAT: This MPQS implementation is superseded by
18 // single_factor<bigint>::mpqs(). In the near future, it will
19 // delegate all calls to single_factor<bigint>.
20 // rational_factorization will eventually be removed.
21 //
22 //==============================================================================================
23
24 // Description : see the diploma thesis of Thomas Sosnowski
25
26 #ifdef HAVE_CONFIG_H
27 #include "config.h"
28 #endif
29 #include "LiDIA/rational_factorization.h"
30 #include "LiDIA/random_generator.h"
31 #include "LiDIA/udigit.h"
32 #include "LiDIA/lanczos.h"
33 #include "LiDIA/lidia_signal.h"
34 #include "LiDIA/mpqs_timing.h"
35 #include "LiDIA/arith.inl"
36
37 #include <cmath>
38 #include <cstdlib>
39 #include <cstdio>
40 #include <fstream>
41
42
43 #ifdef LIDIA_NAMESPACE
44 namespace LiDIA
45 {
46 #endif
47
48
49
50 #ifdef DEBUG // for testing correctness of fulls, partials
51 static bigint global_kN; // and combined relations
52 static int *global_FB;
53 #endif
54
55 void concat (char *, char *, char *); // functions for file handling.
56 void merge_unique (char *, char *, char *); // in future these functions
57 char *get_name (const char *); // should be collected in an own
58 void sort_file_lp (FILE *, char *); // class for LP variations
59 void sort_file_line (char *, char *);
60
61 const unsigned int STLEN = 6000;
62 const unsigned int CAND_NUMBER = 700;
63 typedef unsigned char SIEBTYP;
64
65
66 //******************************************************
67 // our own signal handler
68 //******************************************************
69
70 #define lidia_error_handler_n(f, m)\
71 { char *s; \
72 s = get_name("RELATIONS"); std::remove(s); \
73 s = get_name("LANCZOS_FORMAT"); std::remove(s); \
74 s = get_name("NEW_PARTIALS"); std::remove(s); \
75 s = get_name("NEW_FULLS"); std::remove(s); \
76 s = get_name("LP_RELATIONS"); std::remove(s); \
77 s = get_name("TEMPORARY"); std::remove(s); \
78 s = get_name("SORT_LINE0"); std::remove(s); \
79 s = get_name("SORT_LINE1"); std::remove(s); \
80 s = get_name("SORT_LP0"); std::remove(s); \
81 s = get_name("SORT_LP1"); std::remove(s); \
82 lidia_error_handler(f, m); }
83
84
85
LIDIA_SIGNAL_FUNCTION(stop_rf_mpqs)86 LIDIA_SIGNAL_FUNCTION (stop_rf_mpqs)
87 {
88 lidia_error_handler_n ("rational_factorization",
89 "mpqs::the factorization was interrupted");
90 }
91
92
93 // the parameter table for mpqs has the following format:
94 // dd, approximation accuracy, size sieving interval, size FB,
95 // # prime factors in A, total # primes for determination of A,
96 // starting sieving index
97 //
98 // these values are somewhat experimental and can be improved by
99 // heavy testing.
100
101 const float rational_factorization::qs_params[66][7] = {
102 {15, 0.8, 1200, 35, 3, 10, 2},
103 {16, 0.8, 1400, 35, 3, 10, 2},
104 {17, 0.8, 3000, 40, 3, 10, 2},
105 {18, 0.8, 3000, 60, 3, 10, 2},
106 {19, 0.8, 3600, 80, 3, 10, 2},
107 {20, 0.8, 4000, 100, 3, 10, 2},
108 {21, 0.8, 4250, 100, 3, 10, 2},
109 {22, 0.8, 4500, 120, 3, 10, 3},
110 {23, 0.8, 4750, 140, 3, 10, 3},
111 {24, 0.8, 5000, 160, 3, 12, 4},
112 {25, 0.8, 5000, 180, 3, 12, 4},
113 {26, 0.9, 6000, 200, 5, 8, 4},
114 {27, 1.17, 6000, 220, 5, 8, 5},
115 {28, 1.17, 6500, 240, 5, 8, 5},
116 {29, 1.17, 6500, 260, 5, 8, 5},
117 {30, 1.36, 7000, 325, 5, 8, 5},
118 {31, 1.36, 7000, 355, 5, 8, 5},
119 {32, 1.36, 7500, 375, 5, 8, 5},
120 {33, 1.43, 7500, 400, 6, 8, 6},
121 {34, 1.43, 7500, 425, 6, 8, 6},
122 {35, 1.43, 7500, 550, 6, 10, 6},
123 {36, 1.43, 8000, 650, 6, 10, 6},
124 {37, 1.69, 9000, 750, 6, 10, 7},
125 {38, 1.69, 10000, 850, 6, 10, 7},
126 {39, 1.69, 11000, 950, 6, 10, 7},
127 {40, 1.69, 14000, 1000, 6, 10, 7},
128 {41, 1.69, 14000, 1150, 6, 10, 8},
129 {42, 1.69, 15000, 1300, 6, 10, 8},
130 {43, 1.69, 15000, 1600, 6, 10, 8},
131 {44, 1.69, 15000, 1900, 7, 10, 9},
132 {45, 1.69, 15000, 2200, 7, 10, 9},
133 {46, 1.69, 20000, 2500, 7, 10, 9},
134 {47, 1.69, 25000, 2500, 7, 10, 10},
135 {48, 1.69, 27500, 2700, 7, 10, 10},
136 {49, 1.69, 30000, 2800, 7, 10, 10},
137 {50, 1.75, 35000, 2900, 7, 11, 10},
138 {51, 1.75, 40000, 3000, 7, 11, 10},
139 {52, 1.85, 50000, 3200, 7, 11, 11},
140 {53, 1.85, 50000, 3500, 7, 11, 11},
141 {54, 1.95, 80000, 3800, 7, 12, 11},
142 {55, 1.95, 90000, 4100, 7, 12, 11},
143 {56, 1.95, 100000, 4400, 7, 12, 11},
144 {57, 1.95, 80000, 4700, 8, 12, 12},
145 {58, 1.95, 80000, 5000, 8, 12, 12},
146 {59, 2.15, 130000, 5500, 8, 12, 12},
147 {60, 2.15, 140000, 5800, 8, 12, 12},
148 {61, 2.15, 150000, 6100, 8, 13, 13},
149 {62, 2.15, 160000, 6400, 8, 13, 13},
150 {63, 2.35, 170000, 6700, 8, 13, 13},
151 {64, 2.35, 180000, 7000, 8, 13, 13},
152 {65, 2.35, 190000, 7300, 8, 13, 13},
153 {66, 2.35, 200000, 7600, 8, 13, 13},
154 {67, 2.4, 150000, 7900, 8, 13, 13},
155 {68, 2.4, 150000, 8200, 8, 14, 13},
156 {69, 2.4, 130000, 8600, 8, 14, 13},
157 {70, 2.45, 130000, 8800, 8, 14, 13},
158 {71, 2.45, 130000, 8800, 9, 14, 13},
159 {72, 2.4, 260000, 9400, 9, 14, 13},
160 {73, 2.4, 270000, 9700, 9, 14, 13},
161 {74, 2.4, 280000, 9000, 9, 14, 13},
162 {75, 2.6, 140000, 9000, 9, 14, 13},
163 {76, 2.6, 160000, 9400, 9, 14, 13},
164 {77, 2.6, 180000, 9600, 9, 14, 13},
165 {78, 2.6, 200000, 9800, 9, 14, 13},
166 {79, 2.65, 220000, 10000, 9, 14, 13},
167 {80, 2.65, 250000, 10500, 9, 14, 13}
168 };
169
170
171
172 //************************************************************
173 // compute the number of occurances of char c in the string s
174 //************************************************************
175
number_c(char * s,char c)176 inline unsigned int number_c (char *s, char c)
177 {
178 unsigned int l = 0;
179
180 while (*s++)
181 if (*s == c)
182 l++;
183
184 return (l);
185 }
186
187
188
189 //************************************************************
190 // transform relations from internal format into format needed
191 // in lanczos algorithm.
192 //************************************************************
193
194
transform_relations()195 static void transform_relations ()
196 {
197 char *s1;
198 long val;
199
200 s1 = get_name ("RELATIONS"); // file with relations in internal format
201 FILE *fpin = fopen (s1, "r");
202
203 s1 = get_name ("LANCZOS_FORMAT"); // file with relations in lanczos format
204 FILE *fpout = fopen (s1, "w");
205
206 char *p, *pp, zeile[STLEN], string[STLEN];
207 int i = 0, j, number_entries;
208
209 if (!fpin || !fpout)
210 lidia_error_handler_n ("rational_factorization",
211 "mpqs::transform_relations::can't open files");
212
213 while (fgets (zeile, STLEN, fpin))
214 {
215 if ((pp = (char *) strchr (zeile, ':')))
216 {
217 pp += 2;
218 number_entries = number_c (pp, ' ') + 2; // count # entries per row
219 string[0] = '\0';
220
221 p = (char *) strtok (pp, " \n");
222
223 while (p != NULL)
224 {
225 j = std::strtol (p,NULL,10); // handle only odd exponents
226 if (j & 1)
227 {
228 if(strlen(string) > STLEN - 3)
229 {
230 lidia_error_handler_n("rational_factorization",
231 "mpqs::transform_relations::insufficent string length");
232 }
233 strcat(string, " 1");
234 p = (char *) strtok (NULL, " \n");
235 val = std::strtol (p,NULL,10);
236 if(val == LONG_MAX)
237 {
238 lidia_error_handler_n ("rational_factorization",
239 "mpqs::transform_relations::error in file content");
240
241 }
242 else
243 {
244 char valbuffer[STLEN];
245 int vallength = snprintf(valbuffer, STLEN, " %ld", val);
246 if(STLEN <= strlen(string) + vallength)
247 {
248 lidia_error_handler_n("rational_factorization",
249 "mpqs::transform_relations::insufficent string length");
250 }
251 strcat(string, valbuffer);
252 }
253 p = (char *) strtok (NULL, " \n");
254 }
255 else
256 {
257 number_entries -= 2;
258 p = (char *) strtok (NULL, " \n");
259 p = (char *) strtok (NULL, " \n");
260 }
261 }
262 fprintf (fpout, "%d %d 0%s 0\n", ++i, number_entries, string);
263 }
264 }
265 fflush(fpout);
266 fclose (fpin);
267 fclose (fpout);
268 }
269
270
271
272
273
274 //********************************************************************
275 // initialize mpqs parameters from precomputed table
276 //********************************************************************
277
278 void rational_factorization::
qs_read_par(unsigned int stellen,double & T,unsigned int & M,unsigned int & size_FB,unsigned int & P_ONCE,unsigned int & POLY,unsigned int & P_TOTAL,unsigned int & smallstart)279 qs_read_par (unsigned int stellen, double &T, unsigned int &M,
280 unsigned int &size_FB, unsigned int &P_ONCE,
281 unsigned int &POLY, unsigned int &P_TOTAL,
282 unsigned int &smallstart)
283 {
284
285 unsigned int i;
286
287 if (stellen < 15)
288 i = 0;
289 else
290 i = stellen - 15;
291
292 T = qs_params[i][1];
293 M = static_cast < unsigned int >(qs_params[i][2]); // length of sieve array is 2*M
294
295 size_FB = static_cast < unsigned int >(qs_params[i][3]); // size of factor base
296
297 P_ONCE = static_cast < unsigned int >(qs_params[i][4]); // number of prime factors of
298
299 // coefficient A
300
301 P_TOTAL = static_cast < unsigned int >(qs_params[i][5]); // total number of primes,
302
303 // the P_ONCEs are taken from
304
305 POLY = static_cast < unsigned int >(std::pow (2.0, static_cast < double >(P_ONCE - 1))); // number of coefficients B that
306
307 // are used for
308 smallstart = static_cast < unsigned int >(qs_params[i][6]); // begin sieving
309
310 // with FB[smallstart]
311
312 }
313
314
315
316 //************************************************************
317 // assume that s is a relation of the form "... : FB_i exp_i ..."
318 // and exp is a vector of exponents.
319 // the function scans s and adds (sign == 1) or subtracts
320 // (sign == -1) s to the array exp
321 //
322 // -->used in combination of LP relations
323 //************************************************************
324
add_sub_exp(char * s,int * exp,int sign)325 inline void add_sub_exp (char *s, int *exp, int sign)
326 {
327 char *pp, *p, hilf[STLEN];
328 long e;
329 long val;
330
331 strcpy (hilf, s);
332
333 if ((pp = strchr (hilf, ':'))) {
334 pp += 2;
335
336 p = strtok (pp, " \n");
337 while (p != NULL) {
338 e = std::strtol (p,NULL,10);
339 if (e==LONG_MAX) {
340 lidia_error_handler_n ("rational_factorization",
341 "mpqs::add_sub_exp::error in input file (invalid exponent)");
342 }
343 else if(!e) {
344 break;
345 }
346 else
347 {
348 p = strtok (NULL, " \n");
349 val = std::strtol (p,NULL,10);
350 if (val == LONG_MAX || val < 0) {
351 lidia_error_handler_n ("rational_factorization",
352 "mpqs::add_sub_exp::error in input file (invalid prime)");
353 }
354 else {
355 exp[val] += sign * e;
356 }
357 p = strtok (NULL, " \n");
358 }
359 }
360 }
361 }
362
363
364
365 //******************************************************************
366 // count full relations which can be combined from partial relations
367 //******************************************************************
368
count_partials()369 static unsigned long count_partials ()
370 {
371 char *s1, *s2, *s3;
372 FILE *fp;
373 char line[STLEN];
374
375 s1 = get_name ("NEW_PARTIALS"); // new found LP relations
376 s2 = get_name ("TEMPORARY"); // temporary file
377 s3 = get_name ("LP_RELATIONS"); // holds already processed LP relations
378
379 if ((fp = fopen (s1, "r")) != NULL) // if there are new LP relations
380 {
381 sort_file_lp (fp, s2); // sort file with partial relations with key LP
382 fclose (fp);
383 std::remove (s1);
384 merge_unique (s3, s2, s1); // merge file to already processed LP
385 std::remove (s3);
386 rename (s1, s3);
387 }
388
389 fp = fopen (s3, "r");
390
391 if (fp == NULL) // no file --> no LP relations can be combined
392 return 0;
393
394 // now the counting phase starts
395
396 long counter = 0, large, large_old = -1, i = 1;
397
398 while (fgets (line, STLEN, fp)) {
399 // idea: count all LP relations with same LP
400 large = std::strtol (line,NULL,10);
401
402 if (large == large_old)
403 {
404 counter += i;
405 i++;
406 }
407 else
408 {
409 large_old = large;
410 i = 1;
411 }
412 }
413 fclose (fp);
414 return (counter);
415 }
416
417
418
419 //**************************************************************
420 // combine LP relations to full relations
421 //**************************************************************
422
zusam(unsigned int size_exp,int * FB,const bigint & kN)423 static void zusam (unsigned int size_exp, int * FB, const bigint & kN)
424 {
425 char *s1, *s2, *s3;
426
427 s1 = get_name ("LP_RELATIONS"); // file with LP relations
428 s2 = get_name ("TEMPORARY"); // temporary file
429 sort_file_line (s1, s2);
430
431 FILE *fi = fopen (s2, "r");
432
433 if (fi == NULL)
434 return;
435
436 FILE *fneu;
437 char line[50][STLEN], faktor[STLEN], Hstring[STLEN], Hstring2[STLEN];
438 bigint h, h_old, h_inv, neu;
439 unsigned int large=0, large_old=0, i;
440 int *exp = NULL;
441 long pos, ii, j, ss = 0;
442 bigint H1, H2, H3;
443
444 s2 = get_name ("NEW_FULLS"); // file for new full relations
445 fneu = fopen (s2, "w");
446 if (fneu == NULL)
447 lidia_error_handler_n ("rational_factorization",
448 "zusam::can't open output file");
449
450 exp = new int[size_exp];
451 memset(exp,0,size_exp * sizeof(int));
452 memset(line,0,sizeof(line));
453
454 while (fgets (line[ss & 1], STLEN, fi))
455 {
456 sscanf (line[ss & 1], "%d @ %*s ", &large);
457 ss++;
458
459 if (large == large_old) // two LP relations with same large prime
460 { // -->combination is done
461 pos = 2;
462 while (large == large_old && fgets (line[pos], STLEN, fi))
463 {
464 if (pos >= 50)
465 {
466 lidia_error_handler_n ("rational_factorization",
467 "zusam::too many entries with same large");
468 break;
469 }
470 sscanf (line[pos], "%d @ %*s ", &large);
471 pos++;
472 }
473 pos--;
474
475 for (ii = 0; ii < pos - 1; ii++)
476 {
477 Hstring2[0] = '\0';
478 sscanf (line[ii], "%d @ %s ", &large_old, Hstring2);
479
480 for (j = ii + 1; j < pos; j++)
481 {
482 Hstring[0] = '\0';
483 sscanf (line[j], "%d @ %s ", &large, Hstring);
484 string_to_bigint(Hstring, H1);
485 string_to_bigint(Hstring2, H2);
486 H1 = (H1 * H1) % kN;
487 H2 = (H2 * H2) % kN;
488
489 memset (exp, 0, size_exp * sizeof (int));
490 add_sub_exp (line[j], exp, 1);
491 add_sub_exp (line[ii], exp, -1);
492 memset (faktor, 0, STLEN);
493
494 if (exp[1] != 0)
495 {
496 sprintf (faktor, "%s 1 1", faktor);
497 H2.negate();
498 }
499
500 for (i = 2; i < size_exp; i++)
501 {
502 if (exp[i] != 0)
503 {
504 sprintf (faktor, "%s %d %d", faktor, exp[i], i);
505
506 if (exp[i] > 0)
507 {
508 power (H3, bigint (FB[i]), exp[i]);
509 H2 = (H2 * H3) % kN;
510 }
511 else
512 {
513 power (H3, bigint (FB[i]), -exp[i]);
514 H1 = (H1 * H3) % kN;
515 }
516 }
517
518 if (faktor[0] == (char) NULL)
519 continue;
520 }
521
522 if (H1 != H2 && H1 != kN - H2)
523 fprintf (fneu, "%s %s : %s 0\n", Hstring, Hstring2,
524 faktor);
525
526 #ifdef DEBUG
527 if ((H1 - H2) % global_kN != 0)
528 std::cout << "\n ERROR found: " << std::flush;
529 #endif
530 }
531 }
532 strcpy (line[(ss + 1) & 1], line[pos]);
533 sscanf (line[pos], "%d @ %*s ", &large_old);
534 }
535 else
536 large_old = large;
537 }
538
539 fclose (fi);
540 fclose (fneu);
541 delete[]exp;
542
543 s3 = get_name ("TEMPORARY");
544 std::remove (s3);
545 s1 = get_name ("RELATIONS");
546
547 concat (s1, s2, s3); // concat fulls and newly combined LP relations
548 std::remove (s1);
549 std::remove (s2);
550 sort_file_line (s3, s1); // automatically delete multiple relations
551 }
552
553
554
555 //**********************************************************
556 // determine approximate running time for mpqs on i digit
557 // number, if print = true, then print result on stdout
558 // seems to be very out-of-date !!
559 //**********************************************************
560
zeitqs(unsigned int dec_size,bool print)561 double rational_factorization::zeitqs (unsigned int dec_size, bool print)
562 {
563 double sec;
564 int day, st, min;
565 double timings_celeron[] = {3.2, 3.4, 4.2, 4.6, 5.5, 8.0, 9.3, 11.5, 14.1,
566 17.4, 17.8, 23.9, 29.5, 40.9, 61.7, 109.5, 142.6, 199.6, 264.2};
567
568 // average timings for integers >= 39, <= 57 digits
569 // on Celeron
570
571 if (dec_size <= 38) // "const time" for small numbers
572 sec = 2 * mpqs_machine_factor; // mpqs_machine_factor is computed externally
573 else
574 {
575 if (dec_size <= 57)
576 sec = timings_celeron[dec_size - 39] * mpqs_machine_factor;
577 else
578 sec = std::pow (2.0, LiDIA::log2(timings_celeron[57-39])
579 + (dec_size - 57)*0.53) * mpqs_machine_factor;
580 }
581
582 if (print == false)
583 return sec;
584
585 st = static_cast < int >(sec / 3600.0);
586 min = static_cast < int >(sec / 60.0);
587
588 if (st < 1)
589 {
590 if (min < 1)
591 std::cout << static_cast < int >(sec) << " seconds \n";
592 else
593 std::cout << min << " minutes " <<
594 static_cast < int > (sec - min * 60.0) <<" seconds\n";
595 }
596 else
597 {
598 day = static_cast < int >(st / 24.0);
599 min -= 60 * st;
600
601 if (day < 1)
602 std::cout << st << " hours " << min << " minutes \n";
603 else
604 {
605 st -= 24 * day;
606 std::cout << day << " days " << st <<" hours ";
607 std::cout << min << " minutes";
608 }
609 }
610 return (sec);
611 }
612
613
614
615 //*******************************************************************
616 // main sieving routine:
617 //
618 // FB is a pointer to an array which holds the factor basis
619 // LOGP is a pointer to an array which holds the approximations for
620 // the logarithms of the factor basis elements
621 // START1, START2 are arrays for starting points for different FB elements
622 // sieb points to a sieve array
623 // ende points to the end of the sieve array
624 // M is the size of the sieving interval
625 // CANDIDATE is an array which is filled with candidates which might split
626 // smallstart marks the first FB element which is used for sieving
627 //
628 //******************************************************************
629
630 // FIXME: not portable!
631 #ifdef WORDS_BIGENDIAN
632 // big endian
633 # define INT32_OCTET_0_0x80 0x80000000U
634 # define INT32_OCTET_1_0x80 0x00800000U
635 # define INT32_OCTET_2_0x80 0x00008000U
636 # define INT32_OCTET_3_0x80 0x00000080U
637 #else
638 // little endian
639 # define INT32_OCTET_0_0x80 0x00000080U
640 # define INT32_OCTET_1_0x80 0x00008000U
641 # define INT32_OCTET_2_0x80 0x00800000U
642 # define INT32_OCTET_3_0x80 0x80000000U
643 #endif
644
645 static void
qs_sieve_interval(int * FB,SIEBTYP * LOGP,int * START1,int * START2,SIEBTYP * sieb,SIEBTYP * ende,unsigned int M,int * CANDIDATE,unsigned int smallstart)646 qs_sieve_interval (int *FB, SIEBTYP * LOGP, int *START1,
647 int *START2, SIEBTYP * sieb, SIEBTYP * ende,
648 unsigned int M, int *CANDIDATE, unsigned int smallstart)
649 {
650 register int p, l, *fbp, *lsieb = (int *) sieb; // FIXME: SIEBTYP and int are incompatible types
651 register SIEBTYP logp;
652 register SIEBTYP *begin;
653
654 register int x, counter = 0, M_2 = M << 1;
655 register int oldstart1;
656
657 memset (sieb, 0, (M_2) * sizeof (SIEBTYP));
658
659 fbp = &FB[smallstart];
660 l = smallstart;
661
662 while ((p = *fbp++) != 0) {
663 logp = LOGP[l];
664 begin = sieb + START1[l];
665 oldstart1 = START1[l];
666
667 for (;;) {
668 // sieving with FB[l] from START1[l]
669 if (begin <= ende) {
670 (*begin) += logp;
671 begin += p;
672 } else
673 break;
674 }
675
676 if (oldstart1 != START2[l]) {
677 begin = sieb + START2[l];
678 for (;;) {
679 if (begin <= ende) {
680 (*begin) += logp;
681 begin += p;
682 } else
683 break;
684 }
685 }
686 l++;
687 }
688
689 l = 0;
690 while (l < M_2) {
691 if ((p = (*lsieb)) & (0x80808080)) {
692 // check whether at least one of
693 // four consequent sieve entries
694 // is a candidate
695
696 x = l;
697
698 if (p & (INT32_OCTET_0_0x80 | INT32_OCTET_1_0x80)) {
699 if (p & (INT32_OCTET_0_0x80)) {
700 CANDIDATE[counter++] = x;
701 if (p & (INT32_OCTET_1_0x80))
702 CANDIDATE[counter++] = x + 1;
703 } else
704 CANDIDATE[counter++] = x + 1;
705
706 if (p & (INT32_OCTET_2_0x80))
707 CANDIDATE[counter++] = x + 2;
708
709 if (p & (INT32_OCTET_3_0x80))
710 CANDIDATE[counter++] = x + 3;
711 } else {
712 if (p & (INT32_OCTET_2_0x80))
713 CANDIDATE[counter++] = x + 2;
714 else
715 CANDIDATE[counter++] = x + 3;
716 }
717 }
718 lsieb++;
719 l += 4;
720 }
721 CANDIDATE[counter] = 0;
722 }
723
724
725
726 //******************************************************************
727 // compute optimal multiplier in the set cand of candidates and return
728 // this value
729 //******************************************************************
730
731 int rational_factorization::
compute_multiplier(const bigint & N,int bis,ecm_primes & prim)732 compute_multiplier (const bigint & N, int bis, ecm_primes & prim)
733 {
734 int cand[5] = {1, 3, 5, 7, 11};
735
736 bigint kN;
737 register unsigned long plauf;
738 register int p, j, i, k = 1, nmod4;
739 double wert, bestwert = 1, plus;
740
741 nmod4 = static_cast < int >(N.least_significant_digit ()) & 0x3;
742
743 for (j = 0; j <= 4; j++)
744 {
745 if ((((p = cand[j]) * nmod4) & 0x00000003) != 1)
746 continue;
747
748 wert = -0.7 * LiDIA::log2 (static_cast < double >(p));
749
750 multiply (kN, N, p);
751
752 if ((kN.least_significant_digit () & 0x7) == 1)
753 wert += 1.38629;
754
755 plauf = prim.getprimes ();
756 i = 0;
757
758 while (i <= bis)
759 {
760 if (legendre
761 (static_cast <
762 int >(remainder (kN, static_cast < long >(plauf))),
763 static_cast < int >(plauf)) == 1)
764 {
765 i++;
766 plus =
767 LiDIA::log2 (static_cast < double >(plauf)) / static_cast <
768 double >(plauf);
769 if ((p % plauf) == 0)
770 wert += plus;
771 else
772 wert += 2 * plus;
773 }
774
775 plauf = prim.getprimes ();
776 }
777
778 if (wert > bestwert) {
779 bestwert = wert;
780 k = p;
781 }
782 prim.resetprimes (1);
783 }
784 return (k);
785 }
786
787
788
789 //******************************************************************
790 // create the factor basis
791 //******************************************************************
792
793
794 int rational_factorization::
create_FB(unsigned int size,const bigint & kN,int ** FB,ecm_primes & prim)795 create_FB (unsigned int size, const bigint & kN, int **FB,
796 ecm_primes & prim)
797 {
798 register unsigned int osize, p;
799 register int *fbb;
800
801 if (!(*FB = new int[size + 3]))
802 lidia_error_handler_n ("rational_factorization",
803 "mpqs::create_FB::can't allocate memory");
804
805 fbb = FB[0];
806
807 *fbb++ = static_cast < int >(size);
808
809 *fbb++ = -1;
810
811 osize = 0;
812 prim.resetprimes (1);
813 p = 2;
814
815 while (osize < size)
816 {
817 if (legendre
818 (static_cast < int >(remainder (kN, static_cast < long >(p))),
819 p) != -1)
820 {
821 // p in factor base!
822 *fbb++ = p;
823 osize++;
824 }
825 p = prim.getprimes ();
826 }
827
828 *fbb = 0;
829 return (0);
830 }
831
832
833
834 //**************************************************************
835 // determine the number of ones in binary representation of bi
836 //**************************************************************
837
count_ones(unsigned int bi)838 inline int count_ones (unsigned int bi)
839 {
840 int s = 0;
841
842 while (bi > 0)
843 {
844 if (bi & 1)
845 s++;
846 bi >>= 1;
847 }
848 return s;
849 }
850
851
852
853 //*****************************************************************
854 // compute coefficients of sieving polynomial for self initializing
855 // variant. Coefficients A and B are returned abd several tables are
856 // updated -->see Thomas Sosnowskis diploma thesis.
857 //*****************************************************************
858
859 static void
compute_coeff(bigint & A,bigint & B,const bigint & kN,int * FB,int * SQRTkN,int * START1,int * START2,int P_ONCE,int P_TOTAL,int SIEBSTART,int ** vorb,int * Q_prime,int * Q_prime_glob,bigint * BG,lidia_size_t index_i,int start_fb,int * a_inv,bigint & A4_inverse,unsigned int & bin_index)860 compute_coeff (bigint & A, bigint & B, const bigint & kN, int *FB,
861 int *SQRTkN, int *START1, int *START2,
862 int P_ONCE, int P_TOTAL, int SIEBSTART,
863 int **vorb, int *Q_prime, int *Q_prime_glob,
864 bigint * BG, lidia_size_t index_i,
865 int start_fb, int *a_inv, bigint & A4_inverse,
866 unsigned int &bin_index)
867 {
868 register int p, size_FB;
869 int SIEBS, tmp, tmp1, tmp2;
870 lidia_size_t j, nu_2, i;
871 bigint reserve, TMP;
872
873 if (index_i == 0)
874 {
875 bin_index++;
876
877 while (count_ones (bin_index) != P_ONCE)
878 bin_index++;
879
880 i = 0;
881 for (j = 0; j < P_TOTAL; j++) // determine primes used for A
882 { // in this iteration
883 if (bin_index & (1 << j)) {
884 Q_prime[i] = Q_prime_glob[j];
885 i++;
886 }
887 }
888
889 A.assign (Q_prime[0]); // compute coefficient A
890 for (i = 1; i < P_ONCE; i++)
891 multiply (A, A, Q_prime[i]);
892
893 shift_left (A4_inverse, A, 2);
894
895 // compute BG[0] to BG[P_ONCE-1]
896 // each B is of the form
897 // v_0*BG[0]+v_1*BG[1]+...+v_(P_ONCE-1)*BG[(P_ONCE-1)],
898 // where each v_j is +1 or -1
899
900 // this has to be done only once (for index_i==0) for the
901 // coefficient A; if index_i > 0 then there is a linear
902 // recursion for B
903
904 for (i = 0; i < P_ONCE; i++) {
905 p = Q_prime[i];
906 divide (reserve, A, p);
907 tmp =
908 static_cast < int >(remainder (reserve, static_cast < long >(p)));
909 if (tmp < 0)
910 tmp += p;
911 multiply (reserve, reserve, invert (tmp, p));
912 tmp =
913 ressol (static_cast <
914 int >(remainder (kN, static_cast < long >(p))), p);
915 if (tmp < 0)
916 tmp += p;
917 multiply (reserve, reserve, tmp);
918 remainder (BG[i], reserve, A);
919 }
920
921 B.assign (BG[0]); // compute actual B coefficient
922 for (i = 1; i < P_ONCE; i++)
923 add (B, B, BG[i]);
924
925 if (B.is_even ()) // assure B = 1 mod 4
926 add (B, (A.least_significant_digit () & 3) * A, B);
927
928 size_FB = FB[0] + 1; // a_inv[i] = 1/(2*A) mod p_i
929
930 for (i = 2; i <= size_FB; i++)
931 a_inv[i] =
932 invert (static_cast <
933 int >(remainder (A << 1, static_cast < long >(FB[i]))),
934 FB[i]);
935
936 for (i = 0; i < P_ONCE; i++) {
937 // vorb[i][j] = 1/A * B[i] mod p_j
938 for (j = 2; j <= size_FB; j++) {
939 p = FB[j];
940 multiply (reserve, BG[i], a_inv[j] << 1);
941 if ((tmp =
942 static_cast <
943 int >(remainder (reserve, static_cast < long >(p)))) <0)
944 tmp += p;
945
946 vorb[i][j] = tmp;
947 }
948 }
949
950 for (j = 2; j <= size_FB; j++) {
951 p = FB[j];
952 SIEBS = SIEBSTART % p;
953
954 tmp = static_cast < int >(remainder (-B, p));
955
956 tmp1 = (tmp - SQRTkN[j]) % p;
957 if (tmp1 < 0)
958 tmp1 += p;
959 tmp = (tmp + SQRTkN[j]) % p;
960 if (tmp < 0)
961 tmp += p;
962
963 tmp2 =
964 static_cast <
965 int
966 >(multiply_mod
967 (static_cast < udigit > (tmp), udigit (a_inv[j]),
968 static_cast < udigit > (p)));
969
970 tmp2 = (tmp2 + SIEBS) % p;
971 if (tmp2 < 0)
972 START1[j] = tmp2 + p;
973 else
974 START1[j] = tmp2;
975
976 tmp2 =
977 static_cast <
978 int
979 >(multiply_mod
980 (static_cast < udigit > (tmp1), udigit (a_inv[j]),
981 static_cast < udigit > (p)));
982 tmp2 = (tmp2 + SIEBS) % p;
983 if (tmp2 < 0)
984 START2[j] = tmp2 + p;
985 else
986 START2[j] = tmp2;
987 }
988
989 TMP = xgcd_left (A4_inverse, A4_inverse, kN); // determine 1/(4A) mod kN
990 }
991
992 else // no "real" computation -- use recursive formula
993 { // first: update of B, compute B[index_i], index_i > 0
994
995 nu_2 = 0; // nu_2 = nu_2(index_i)
996 j = index_i;
997 while ((j & 1) == 0) {
998 nu_2++;
999 j >>= 1;
1000 }
1001
1002 shift_left (TMP, BG[nu_2], 1);
1003
1004 if ((((j + 1) / 2) & 1) == 1) {
1005 i = -1;
1006 subtract (B, B, TMP);
1007 } else {
1008 i = 1;
1009 add (B, B, TMP);
1010 }
1011
1012 size_FB = FB[0] + 1; // determine new starting positions
1013 // for sieving
1014
1015 if (i == -1) {
1016 for (j = 2; j <= size_FB; j++) {
1017 p = FB[j];
1018 START1[j] += vorb[nu_2][j];
1019 if (START1[j] >= p)
1020 START1[j] -= p;
1021 START2[j] += vorb[nu_2][j];
1022 if (START2[j] >= p)
1023 START2[j] -= p;
1024 }
1025 } else {
1026 for (j = 2; j <= size_FB; j++) {
1027 p = FB[j];
1028 START1[j] -= vorb[nu_2][j];
1029 if (START1[j] < 0)
1030 START1[j] += p;
1031 START2[j] -= vorb[nu_2][j];
1032 if (START2[j] < 0)
1033 START2[j] += p;
1034 }
1035 }
1036 }
1037
1038 if (FB[2] == 2) // note special situation for p = 2
1039 {
1040 START1[2] = 1;
1041 START2[2] = 1;
1042 }
1043
1044 // now compute zeros of polynomials that have only one zero mod p
1045 // because p divides coefficient A
1046
1047 square (reserve, B); // compute coefficient -C
1048 subtract (reserve, kN, reserve);
1049 shift_left (TMP, A, 2);
1050 divide (TMP, reserve, TMP);
1051
1052 for (j = 1; j <= P_TOTAL; j++)
1053 if (bin_index & (1 << (j - 1))) {
1054 p = FB[start_fb + j];
1055 tmp =
1056 invert (static_cast <
1057 int >(remainder (B, static_cast < long >(p))), p);
1058 if (tmp < 0)
1059 tmp += p;
1060 tmp2 =
1061 static_cast < int >(remainder (TMP, static_cast < long >(p)));
1062 if (tmp2 < 0)
1063 tmp2 += p;
1064
1065 tmp = static_cast <int>(multiply_mod
1066 (static_cast < udigit > (tmp2),
1067 static_cast < udigit > (tmp),
1068 static_cast < udigit > (p)));
1069 START1[start_fb + j] = START2[start_fb + j] =
1070 (tmp + SIEBSTART) % p;
1071 }
1072
1073 /*
1074 #ifdef DEBUG // check correctness of roots mod p
1075 if(FB[2] == 2)
1076 j = 3;
1077 else
1078 j = 2;
1079 for (; j <= FB[0] + 1; j++)
1080 {
1081 p = FB[j];
1082 SIEBS = SIEBSTART % p;
1083 if (((A * (START1[j] - SIEBS) + B) * (START1[j] - SIEBS) +
1084 (B * B - kN) / (4 * A)) % p != 0)
1085 lidia_error_handler ("rational_factorization", "mpqs::compute_coeff::"
1086 "found wrong polynomial in (1)");
1087
1088 if (((A * (START2[j] - SIEBS) + B) * (START2[j] - SIEBS) +
1089 (B * B - kN) / (4 * A)) % p != 0)
1090 lidia_error_handler ("rational_factorization","mpqs::compute_coeff::"
1091 "found wrong polynomial in (2)");
1092 }
1093 #endif
1094 */
1095 }
1096
1097
1098
1099 //-------------------------------------------------------------------
1100 // insert ul into string *p, add blank and return pointer to first
1101 // char after blank --> much faster than sprintf
1102
insert_at(char * p,unsigned long n)1103 inline char *insert_at (char *p, unsigned long n)
1104 {
1105 register int c, i, j, e;
1106
1107 i = 0;
1108 do {
1109 p[i++] = (char) (n % 10 + '0');
1110 } while ((n /= 10) > 0);
1111 e = i;
1112
1113 if (e > 1)
1114 for (i = 0, j = e - 1; i < j; i++, j--) {
1115 c = *(p + i);
1116 *(p + i) = *(p + j);
1117 *(p + j) = c;
1118 }
1119 p[e] = ' ';
1120 return (p + e + 1);
1121 }
1122
1123
1124
1125 //***********************************************************************
1126 // testing routine which filters correct full and LP relations out of
1127 // candidates found in the sieving step.
1128 //***********************************************************************
1129
1130 static int
teste(const bigint & kN,const bigint & A,int * FB,int * START1,int * START2,char * faktor,unsigned int M,double d_wurz,int * Q_prime,const bigint & B,unsigned int start_fb,unsigned int P_ONCE,unsigned int P_TOTAL,int * CANDIDATE,unsigned int smallstart,const bigint & A4_inverse,FILE * fpfull,FILE * fppart)1131 teste (const bigint & kN, const bigint & A, int *FB, int *START1,
1132 int *START2, char *faktor, unsigned int M, double d_wurz,
1133 int *Q_prime, const bigint & B, unsigned int start_fb,
1134 unsigned int P_ONCE, unsigned int P_TOTAL, int *CANDIDATE,
1135 unsigned int smallstart, const bigint & A4_inverse,
1136 FILE * fpfull, FILE * fppart)
1137 {
1138 int small_value = 0;
1139 unsigned int fak_i, ii, M_2, counter, upper_bound;
1140 bigint H, Qx, TMP;
1141 int x = 0, p, vorber, divides = 0, N1, N2;
1142 long rest;
1143 int geteilt = 0, ready, rest_i;
1144 double a, b;
1145 char Hstring[STLEN];
1146 char *faktorp;
1147
1148 // compute the roots of the polynomial
1149 b = dbl (B);
1150 a = dbl (A);
1151 a *= 2;
1152 N1 = static_cast < int >((-b - d_wurz) / a);
1153 N2 = static_cast < int >((-b + d_wurz) / a);
1154
1155 M_2 = M << 1;
1156 upper_bound = start_fb + P_TOTAL + 1;
1157 counter = 0;
1158
1159 while ((x = CANDIDATE[counter++]) != 0) {
1160 // while there are candidates to test
1161 x -= M;
1162 multiply (Qx, A, x << 1);
1163 add (Qx, Qx, B);
1164 div_rem (TMP, H, Qx, kN); // needed for writing relation to file
1165 if (H.is_negative ())
1166 H.negate ();
1167 square (Qx, H);
1168 div_rem (TMP, Qx, Qx, kN);
1169 multiply (Qx, Qx, A4_inverse);
1170 div_rem (TMP, Qx, Qx, kN);
1171
1172 faktor[0] = ' ';
1173 faktorp = faktor + 1;
1174
1175 if (Qx.is_negative () && (x <= N1 || x >= N2))
1176 add (Qx, kN, Qx);
1177 else if (Qx.is_positive () && (N1 < x && x < N2)) {
1178 subtract (Qx, kN, Qx);
1179 faktorp = insert_at (faktorp, 1);
1180 faktorp = insert_at (faktorp, 1);
1181 }
1182
1183 if (Qx.is_negative ()) {
1184 faktorp = insert_at (faktorp, 1);
1185 faktorp = insert_at (faktorp, 1);
1186 Qx.absolute_value ();
1187 }
1188
1189 while (Qx.is_even ()) {
1190 divides++;
1191 Qx.divide_by_2 ();
1192 }
1193
1194 if (divides) {
1195 faktorp = insert_at (faktorp, divides + 2); // '+2' because of 4*A
1196 faktorp = insert_at (faktorp, 2);
1197 divides = 0;
1198 } else {
1199 faktorp = insert_at (faktorp, 2);
1200 faktorp = insert_at (faktorp, 2);
1201 }
1202
1203 fak_i = 2;
1204 divides = 0;
1205 ready = 0;
1206
1207 while (fak_i++ < smallstart - 1) {
1208 p = FB[fak_i];
1209
1210 vorber = (M + x) % p;
1211
1212 if ((vorber == START1[fak_i]) || (vorber == START2[fak_i])) {
1213 do {
1214 div_rem (TMP, rest, Qx, static_cast < long >(p));
1215
1216 if (rest == 0) {
1217 divides++;
1218 Qx.assign (TMP);
1219 }
1220 } while (rest == 0);
1221
1222 faktorp = insert_at (faktorp, divides);
1223 faktorp = insert_at (faktorp, fak_i);
1224 divides = 0;
1225 }
1226 }
1227
1228 while ((p = FB[fak_i]) != 0) {
1229 vorber = (M + x) % p;
1230
1231 if ((fak_i <= upper_bound) && (fak_i > start_fb)) {
1232 ii = ready;
1233 while (ii < P_ONCE) {
1234 if (p == Q_prime[ii]) {
1235 geteilt = 1;
1236 ready = ii + 1;
1237 break;
1238 }
1239 ii++;
1240 }
1241 }
1242
1243 if (geteilt == 1) {
1244 // p divides coefficient A ...
1245 if (vorber == START1[fak_i]) {
1246 // ... and divides Qx
1247 do {
1248 div_rem (TMP, rest, Qx, static_cast < long >(p));
1249
1250 if (rest == 0) {
1251 divides++;
1252 Qx.assign (TMP);
1253 }
1254 } while (rest == 0);
1255 }
1256 faktorp = insert_at (faktorp, divides + 1);
1257 faktorp = insert_at (faktorp, fak_i);
1258
1259 geteilt = 0;
1260 divides = 0;
1261 } else {
1262 if ((vorber == START1[fak_i]) || (vorber == START2[fak_i])) {
1263 do {
1264 div_rem (TMP, rest, Qx, static_cast < long >(p));
1265
1266 if (rest == 0) {
1267 divides++;
1268 Qx.assign (TMP);
1269 }
1270 } while (rest == 0);
1271
1272 faktorp = insert_at (faktorp, divides);
1273 faktorp = insert_at (faktorp, fak_i);
1274 divides = 0;
1275 }
1276 }
1277 fak_i++;
1278 }
1279
1280
1281 if (Qx.is_one ()) {
1282 // full relation found
1283 small_value++;
1284 *(faktorp - 1) = '\0';
1285
1286 bigint_to_string (H, Hstring);
1287 fprintf (fpfull, "%s 1 : %s 0\n", Hstring, faktor);
1288
1289 /*
1290 #ifdef DEBUG // test correctness of full relation
1291 char *p;
1292 bigint h1, h2;
1293 int e;
1294
1295 h1.assign (1);
1296 p = (char *) strtok (faktor, " \n");
1297 while (p != NULL) {
1298 e = std::strtol (p,NULL,10);
1299 if (!e)
1300 break;
1301 p = (char *) strtok (NULL, " \n");
1302 power (h2, bigint (FB[std::strtol (p,NULL,10)]), e);
1303 multiply (h1, h1, h2);
1304 p = (char *) strtok (NULL, " \n");
1305 }
1306
1307 if ((H * H - h1) == 0)
1308 cout<<"\nERROR: FULL RELATION IN Z !!";
1309
1310
1311 if ((H * H - h1) % kN != 0)
1312 lidia_error_handler ("rational_factorization",
1313 "mpqs::teste::found wrong full relation");
1314 #endif
1315 */
1316 }
1317 else
1318 if ((Qx.intify (rest_i)) == 0) // LP relation found
1319 if (rest_i < 10000000) {
1320 *(faktorp - 1) = '\0';
1321 bigint_to_string (H, Hstring);
1322 fprintf (fppart, "%9d @ %s : %s 0\n", rest_i, Hstring, faktor);
1323
1324 /*
1325 #ifdef DEBUG // test correctness of partial relation
1326 char *p;
1327 bigint h1, h2;
1328 int e;
1329
1330 h1.assign (Qx);
1331 p = (char *) strtok (faktor, " \n");
1332 while (p != NULL) {
1333 e = std::strtol (p,NULL,10);
1334 if (!e)
1335 break;
1336 p = (char *) strtok (NULL, " \n");
1337 power (h2, bigint (FB[std::strtol (p,NULL,10)]), e);
1338 multiply (h1, h1, h2);
1339 p = (char *) strtok (NULL, " \n");
1340 }
1341
1342 if ((H * H - h1) == 0)
1343 cout<<"\nERROR: PARTIAL RELATION IN Z !!";
1344
1345
1346 if ((H * H - h1) % kN != 0)
1347 lidia_error_handler ("rational_factoriztion",
1348 "mpqs::teste::found wrong partial relation");
1349 #endif
1350 */
1351 }
1352 }
1353 return (small_value);
1354 }
1355
1356
1357
1358 //*****************************************************************
1359 // after LP relations have been combined, start linear system solver
1360 // and use solutions of linear system to determine factors of kN
1361 //*****************************************************************
1362
1363 bool rational_factorization::
qs_build_factors(const bigint & N,const bigint & kN,unsigned int index,int * FB)1364 qs_build_factors (const bigint & N,
1365 const bigint & kN, unsigned int index, int *FB)
1366 {
1367 rf_single_factor fact;
1368 FILE *faktorenmatrix;
1369 bigint X_quad, Y_quad, hilf, hilf2;
1370 int counter = 0, k, end_FB = FB[0] + 2, e;
1371 bool erg1, erg2;
1372 bool found = false;
1373 size_t i, zeilenindex = 0;
1374 char zeile[STLEN], Hstring[STLEN], Hstring2[STLEN], rest[STLEN];
1375 long *expo;
1376 long exp_i;
1377 char *p, *pp, *s1;
1378
1379 if (!(expo = new long[end_FB + 1]))
1380 lidia_error_handler_n ("rational_factorization",
1381 "mpqs::qs_build_factors::can't allocate memory");
1382
1383 memset (expo, 0, (end_FB + 1) * sizeof (long));
1384 transform_relations();
1385
1386 if (info)
1387 std::cout << "\n\nLinear system built " << std::flush;
1388
1389 // here starts the linear system solver
1390
1391 s1 = get_name ("LANCZOS_FORMAT");
1392 lanczos_sparse_matrix solve_matrix(s1);
1393
1394 preprocess pre;
1395 postprocess post;
1396
1397 std::auto_ptr<index_list> correction_list = pre.process(solve_matrix);
1398
1399 lanczos lan(solve_matrix);
1400
1401 do {
1402 lan.solve();
1403
1404 #ifdef DEBUG
1405 if (lan.get_result_rank() <= 0)
1406 lidia_error_handler("rational_factorization",
1407 "mpqs::Lanczos found no solution!");
1408 #endif
1409
1410 }
1411 while(lan.get_result_rank() <= 0);
1412
1413 std::auto_ptr<lanczos_vector_block> solution = post.process(lan.get_result(),
1414 *correction_list);
1415
1416 lanczos_vector_block::result_vector_type lanczos_result = solution->result();
1417
1418 #ifdef DEBUG
1419 std::cout << "\nVerifying Lanczos solution ... " << std::flush;
1420
1421 lanczos_sparse_matrix *matrix;
1422 lanczos_vector_block *vector;
1423
1424 s1 = get_name ("LANCZOS_FORMAT");
1425 matrix = new lanczos_sparse_matrix (s1);
1426 if (!matrix)
1427 lidia_error_handler("rational_factorization",
1428 "mpqs::Error in Input Lanczos Matrix");
1429
1430 vector = new lanczos_vector_block (matrix->number_of_columns ());
1431 vector->read(lanczos_result);
1432
1433 unsigned long i2, j, h, g;
1434 lanczos_vector_block *result_vec;
1435
1436 result_vec = new lanczos_vector_block (matrix->number_of_rows ());
1437 result_vec->clear ();
1438 i2 = 0;
1439
1440 while (i2 < matrix->number_of_columns ())
1441 {
1442 j = 0;
1443 h = matrix->get_vector (i2).get_number_of_entries ();
1444 while (j < h) {
1445 g = matrix->get_vector (i2).get_entry (j);
1446 result_vec->put_row (g, vector->get_row (i2) ^
1447 (result_vec->get_row (g)));
1448 j++;
1449 }
1450 i2++;
1451 }
1452 if (!result_vec->is_zero ())
1453 lidia_error_handler("rational_factorization",
1454 "mpqs::Solution of Lanczos is wrong !!");
1455 else
1456 std::cout << "\nSolution of Lanczos is correct \n " << std::flush;
1457
1458 delete (matrix);
1459 delete (vector);
1460 delete (result_vec);
1461 #endif
1462
1463 if (info)
1464 std::cout << "\nand solved with Block-Lanczos Algorithm\n" << std::flush;
1465
1466 s1 = get_name ("RELATIONS");
1467 if (!(faktorenmatrix = fopen (s1, "r")))
1468 lidia_error_handler_n ("rational_factorization",
1469 "mpqs::qs_build_factors::can't open relation"
1470 " file");
1471
1472 //--------------------------------------------------------------
1473 // now the combination of relations with these solutions starts
1474
1475 int no_solution = 0;
1476
1477 while (!found) // use solutions to find X, Y
1478 // with X^2 = Y^2 mod N
1479 {
1480 size_t lanczos_number_of_entries;
1481
1482 X_quad.assign_one ();
1483 Y_quad.assign_one ();
1484
1485 counter = 1;
1486 fgets (zeile, STLEN, faktorenmatrix); // skip .. relations
1487
1488 if (lanczos_result[no_solution].empty())
1489 {
1490 break;
1491 }
1492 else
1493 lanczos_number_of_entries = lanczos_result[no_solution][0];
1494
1495 for (i = 1; i <= lanczos_number_of_entries; i++)
1496 {
1497 zeilenindex = lanczos_result[no_solution][i];
1498
1499 while (counter <= zeilenindex)
1500 {
1501 fgets (zeile, STLEN, faktorenmatrix);
1502 counter++;
1503 }
1504
1505 sscanf (zeile, "%s %s : %s", Hstring, Hstring2, rest);
1506
1507 string_to_bigint (Hstring, hilf);
1508 string_to_bigint (Hstring2, hilf2);
1509
1510 if (!hilf2.is_one ())
1511 {
1512 multiply (X_quad, X_quad, hilf2);
1513 remainder (X_quad, X_quad, kN);
1514 }
1515
1516 multiply (Y_quad, Y_quad, hilf);
1517 remainder (Y_quad, Y_quad, kN);
1518
1519 // read single exponents and sum them up in array for later use
1520
1521 if ((pp = (char *) strchr (zeile, ':'))) {
1522 pp += 2;
1523 p = (char *) strtok (pp, " \n");
1524 while (p != NULL) {
1525 e = std::strtol (p,NULL,10);
1526 if (!e)
1527 break;
1528 p = (char *) strtok (NULL, " \n");
1529 expo[std::strtol (p,NULL,10)] += e;
1530 p = (char *) strtok (NULL, " \n");
1531 }
1532 }
1533 } // solution read
1534 no_solution ++;
1535
1536
1537 for (i = 2; i <= end_FB; i++)
1538 {
1539 exp_i = expo[i];
1540
1541 if (exp_i > 0)
1542 {
1543 exp_i >>= 1;
1544
1545 while (exp_i != 0)
1546 {
1547 exp_i--;
1548 multiply (X_quad, X_quad, FB[i]);
1549 }
1550 if (X_quad.abs_compare (kN) >= 0)
1551 remainder (X_quad, X_quad, kN);
1552 }
1553 else
1554 if (exp_i < 0)
1555 {
1556 exp_i = -((-exp_i) >> 1);
1557
1558 while (exp_i != 0)
1559 {
1560 exp_i++;
1561 multiply (Y_quad, Y_quad, FB[i]);
1562 }
1563 if (Y_quad.abs_compare (kN) >= 0)
1564 remainder (Y_quad, Y_quad, kN);
1565 }
1566 }
1567
1568 add (hilf, X_quad, Y_quad);
1569 subtract (Y_quad, Y_quad, X_quad);
1570
1571 X_quad = gcd (hilf, N);
1572 Y_quad = gcd (Y_quad, N);
1573
1574 if (X_quad.is_one () || Y_quad.is_one ()) {
1575 // no proper factor found
1576 fseek (faktorenmatrix, 0, SEEK_SET);
1577 memset (expo, 0, (end_FB + 1) * sizeof (long));
1578 continue;
1579 }
1580 else {
1581 if((fact.single_exponent =
1582 is_power(fact.single_base, X_quad)) == 0) {
1583 fact.single_base.assign(X_quad);
1584 fact.single_exponent = 1;
1585 }
1586 fact.factor_state = is_prime(fact.single_base) ?
1587 rf_single_factor::prime : rf_single_factor::not_prime;
1588 factors[index] = fact;
1589
1590 if((fact.single_exponent =
1591 is_power(fact.single_base, Y_quad)) == 0) {
1592 fact.single_base.assign(Y_quad);
1593 fact.single_exponent = 1;
1594 }
1595 fact.factor_state = is_prime(fact.single_base) ?
1596 rf_single_factor::prime : rf_single_factor::not_prime;
1597 factors[no_of_comp()] = fact;
1598
1599 found = true;
1600 fseek (faktorenmatrix, 0, 0);
1601 memset (expo, 0, (end_FB + 1) * sizeof (long));
1602 }
1603 }
1604 fclose (faktorenmatrix);
1605
1606 delete[] expo;
1607
1608 return (!found);
1609 }
1610
1611
1612
1613 //*****************************************************************
1614 // main routine which factors one component of a rational_factorization
1615 // with MPQS, calls all relevant functions in a loop until one factor was
1616 // found.
1617 //******************************************************************
1618
mpqs_impl(lidia_size_t index,ecm_primes & prim)1619 rational_factorization & rational_factorization::mpqs_impl(lidia_size_t index,
1620 ecm_primes & prim) {
1621 bigint N (factors[index].single_base);
1622 bigint A, B, kN, A4_inverse;
1623 bigint *BG;
1624 long tmp, i, k, p;
1625 unsigned int stellenzahl;
1626 unsigned long counter_treff = 0;
1627 FILE *fpfull, *fppart;
1628
1629 unsigned int M, size_FB, smalls = 0;
1630 int small_value;
1631 int *SQRTkN, *FB, *START1, *a_inv;
1632 int *CANDIDATE, *Q_prime, *Q_prime_glob;
1633 int *START2;
1634 int start_fb;
1635 unsigned int P_TOTAL, POLY, P_ONCE, smallstart, prozent, vergleich = 10;
1636 SIEBTYP *sieb = NULL, *ende = NULL, *LOGP = NULL;
1637 double d_wurz, T, LOGMUL;
1638 int **vorb;
1639 lidia_size_t index_i;
1640 lidia_size_t added_relations;
1641 char *s1, *faktor;
1642 unsigned long last_cnt=0;
1643
1644 i = is_power (A, N); // perfect power handling
1645 if (i > 0) {
1646 rf_single_factor fact;
1647
1648 if (is_prime (A, 5))
1649 {
1650 fact.single_base = A;
1651 fact.single_exponent =
1652 factors[index].single_exponent * static_cast < int >(i);
1653 fact.factor_state = rf_single_factor::prime;
1654 factors[index] = fact;
1655 }
1656 else
1657 {
1658 rational_factorization A_fact;
1659
1660 A_fact.assign (A);
1661 A_fact.factor ();
1662
1663 int j, kk = A_fact.no_of_comp ();
1664
1665 for (j = 0; j < kk; j++)
1666 A_fact.factors[j].single_exponent *= static_cast < int >(i);
1667
1668 divide (*this, *this, rational_factorization (N));
1669 multiply (*this, *this, A_fact);
1670 }
1671 return *this;
1672 }
1673
1674
1675 if (info) {
1676 std::cout << "\n\nQS -- quadratic sieve\n == \n\n";
1677 std::
1678 cout << "number to factor: " << N << " (" << decimal_length (N) <<
1679 ")\n\n";
1680 std::cout.flush ();
1681 }
1682
1683 k = compute_multiplier (N, 5, prim);
1684 multiply (kN, N, k);
1685
1686 if (info)
1687 std::cout << "Multiplier: " << k << "\n" << std::flush;
1688 stellenzahl = decimal_length (kN);
1689
1690 if (info) {
1691 std::cout << "\nestimated running time (user-time): ";
1692 zeitqs (stellenzahl, true);
1693 std::cout.flush ();
1694 }
1695
1696 if ((stellenzahl > 80))
1697 {
1698 lidia_warning_handler ("rational_factorization",
1699 "mpqs::Input Number too big to be factored"
1700 " on one machine in reasonable time");
1701 stellenzahl = 80;
1702 }
1703
1704 qs_read_par (stellenzahl, T, M, size_FB, P_ONCE, POLY, P_TOTAL,
1705 smallstart);
1706
1707 if (info) {
1708 std::cout << "\nSieve Interval: [-" << M << " , " << M << "]";
1709 std::cout << "\n# Factor Basis: " << size_FB << "\n";
1710 std::cout.flush ();
1711 }
1712
1713 added_relations = static_cast < lidia_size_t >(size_FB * 0.05);
1714 if (added_relations < 50)
1715 added_relations = 50;
1716
1717 create_FB (size_FB, kN, &FB, prim);
1718
1719 #ifdef DEBUG
1720 global_kN = kN;
1721 global_FB = FB;
1722 #endif
1723
1724 BG = new bigint[P_ONCE + 1];
1725 sieb = new SIEBTYP[M << 1];
1726 faktor = new char[STLEN];
1727
1728 LOGP = new SIEBTYP[size_FB + 2];
1729 SQRTkN = new int[size_FB + 2];
1730 START1 = new int[size_FB + 2];
1731 START2 = new int[size_FB + 2];
1732 a_inv = new int[size_FB + 2];
1733 CANDIDATE = new int[CAND_NUMBER];
1734 Q_prime_glob = new int[P_TOTAL];
1735 Q_prime = new int[P_ONCE];
1736
1737 if (sieb == NULL || faktor == NULL || LOGP == NULL || SQRTkN == NULL ||
1738 START1 == NULL || START2 == NULL || a_inv == NULL || BG == NULL ||
1739 CANDIDATE == NULL || Q_prime_glob == NULL || Q_prime == NULL)
1740 lidia_error_handler_n ("rational_factorization",
1741 "mpqs::can't allocate memory");
1742
1743 vorb = new int *[P_TOTAL];
1744
1745 for (i = 0; i < static_cast < int >(P_TOTAL); i++)
1746 {
1747 if (!(vorb[i] = new int[size_FB + 2]))
1748 lidia_error_handler_n ("rational_factorization",
1749 "mpqs::allocate vorb");
1750 }
1751
1752 ende = sieb + (M << 1) - 1;
1753
1754 // determine approximations for log(p) for FB elements p and sqrt(kN)
1755 // mod p (used for fast computation of sieving start points
1756
1757 LOGMUL = (2 * static_cast < SIEBTYP > (0.5 * LiDIA::log2 (dbl (kN)) +
1758 LiDIA::log2 (static_cast <
1759 double >(M)) -T *
1760 LiDIA::log2 (static_cast <
1761 double >(FB[FB[0] + 1]))));
1762 LOGMUL = 127.0 / LOGMUL;
1763
1764 tmp = size_FB + 2;
1765 for (i = 2; i < tmp; i++) {
1766 p = FB[i];
1767 LOGP[i] =
1768 static_cast < SIEBTYP >
1769 (LOGMUL * LiDIA::log2 (static_cast < double >(p)) * 2);
1770
1771 if ((SQRTkN[i] =
1772 ressol (static_cast <
1773 int >(remainder (kN, static_cast < long >(p))),
1774 static_cast < int >(p))) < 0)
1775 SQRTkN[i] += static_cast < int >(p); // compute sqrt(kN) modulo different moduli p
1776 }
1777
1778 d_wurz = std::sqrt (dbl (kN));
1779
1780 // the size of coefficient A should be approximately
1781 // sqrt(kN)/M, so the size of the primes p dividing
1782 // A should be approximately (sqrt(kN/M))^(1/P_ONCE)
1783
1784 T = d_wurz / M;
1785 T = std::pow (2.0, LiDIA::log2 (T) / P_ONCE);
1786
1787 if (T > FB[size_FB - 1])
1788 lidia_error_handler_n ("rational_factorization",
1789 "mpqs: P_ONCE too small");
1790
1791 i = 2;
1792 while (FB[i] < T)
1793 i++;
1794 start_fb = static_cast < int >(i); // P_TOTAL consecutive primes p[start_fb], ...,
1795
1796 // are chosen from the factor basis
1797
1798 if (i > 7)
1799 start_fb -= (P_ONCE >> 1);
1800
1801 if (k >= FB[start_fb]) // Multiplier must not occur in factor basis
1802 while (k >= FB[start_fb])
1803 start_fb++;
1804
1805
1806 for (i = 0; i < static_cast < int >(P_TOTAL); i++) {
1807 Q_prime_glob[i] = FB[start_fb + i + 1]; // collect prime numbers which
1808 // will build the A coefficents
1809 }
1810
1811 index_i = -1;
1812
1813 lidia_signal sig1 (LIDIA_SIGTERM, stop_rf_mpqs), sig2 (LIDIA_SIGINT,
1814 stop_rf_mpqs);
1815 lidia_signal sig3 (LIDIA_SIGHUP, stop_rf_mpqs), sig4 (LIDIA_SIGSEGV,
1816 stop_rf_mpqs);
1817
1818 s1 = get_name ("LP_RELATIONS");
1819 fpfull = fopen (s1, "w");
1820 fprintf (fpfull, "\n");
1821 fclose (fpfull);
1822
1823 s1 = get_name ("RELATIONS");
1824 if (!(fpfull = fopen (s1, "a")))
1825 lidia_error_handler_n ("rational_factorization",
1826 "mpqs::can't open file RELATIONS");
1827
1828 fprintf (fpfull, "\n");
1829
1830 s1 = get_name ("NEW_PARTIALS");
1831 if (!(fppart = fopen (s1, "a")))
1832 lidia_error_handler_n ("rational_factorization",
1833 "mpqs::can't open file NEW_PARTIALS");
1834
1835 unsigned int bin_index = (1 << P_ONCE) - 1; // variable used for
1836 //choosing the
1837 // correct A coeffs in compute_coeff
1838
1839 while (true) // central loop of
1840 // - computing polynomials and zeros
1841 // - sieving
1842 // - testing candidates of the sieve array
1843 {
1844 if (index_i == static_cast < int >(POLY) - 1) // when all of the B's have already
1845 // been used, choose new A
1846 index_i = 0;
1847 else
1848 index_i++;
1849
1850
1851 compute_coeff (A, B, kN, FB, SQRTkN, START1, START2, P_ONCE, P_TOTAL,
1852 M, vorb, Q_prime, Q_prime_glob, BG, index_i,
1853 start_fb, a_inv, A4_inverse, bin_index);
1854
1855 qs_sieve_interval (FB, LOGP, START1, START2, sieb, ende, M,
1856 CANDIDATE, smallstart);
1857
1858 small_value = teste (kN, A, FB, START1, START2, faktor, M, d_wurz, Q_prime,
1859 B, start_fb, P_ONCE, P_TOTAL, CANDIDATE,
1860 smallstart, A4_inverse, fpfull, fppart);
1861
1862 smalls += small_value;
1863
1864 prozent =
1865 static_cast <
1866 unsigned int >((static_cast < double >(smalls + counter_treff)
1867 / static_cast < double >(size_FB)) *100);
1868
1869 if (prozent >= vergleich)
1870 {
1871 fflush (fppart);
1872 fclose (fppart);
1873
1874 counter_treff = count_partials ();
1875
1876 prozent =
1877 static_cast <
1878 unsigned int >((static_cast < double >(smalls + counter_treff)
1879 / static_cast < double >(size_FB)) *100);
1880
1881 if (info && (last_cnt != (smalls + counter_treff)) )
1882 {
1883 std::cout << smalls +
1884 counter_treff << " (" << prozent << "%) relations ";
1885 std::cout << "found. (" << smalls << " by fulls, "
1886 << counter_treff;
1887 std::cout << " by partials)\n" << std::flush;
1888 }
1889
1890 last_cnt = smalls+counter_treff;
1891
1892 while (vergleich <= prozent)
1893 {
1894 if (vergleich >= 90 || stellenzahl >= 75)
1895 vergleich += 5;
1896 else
1897 vergleich += 10;
1898 }
1899
1900 s1 = get_name ("NEW_PARTIALS");
1901 if (!(fppart = fopen (s1, "a")))
1902 lidia_error_handler_n ("rational_factorization",
1903 "mpqs::can't open file NEW_PARTIALS");
1904 }
1905
1906 if ((static_cast < long >(smalls) + counter_treff) >
1907 (static_cast < long >(size_FB) + added_relations) )
1908 {
1909 fflush (fppart);
1910 fflush (fpfull);
1911 fclose (fppart);
1912 fclose (fpfull);
1913
1914 // need to include the partials that we have already.
1915 counter_treff = count_partials ();
1916
1917 if (info && (last_cnt != (smalls + counter_treff)) )
1918 {
1919 prozent = static_cast < unsigned int >
1920 ((static_cast < double >(smalls + counter_treff)
1921 / static_cast < double >(size_FB)) *100);
1922
1923 std::cout << smalls +
1924 counter_treff << " (" << prozent << "%) relations ";
1925 std::cout << "found. (" << smalls << " by fulls, "
1926 << counter_treff;
1927 std::cout << " by partials)\n" << std::flush;
1928 }
1929 last_cnt = smalls+counter_treff;
1930
1931 zusam (size_FB + 2, FB, kN);
1932
1933 if (qs_build_factors (N, kN, index, FB) == false) {
1934 break;
1935 }
1936 else {
1937 added_relations = static_cast <long> (static_cast < double >
1938 (smalls + counter_treff) * 1.05) - size_FB;
1939 if (info) {
1940 std::cout << "\nNo non-trivial congruence found -->";
1941 std::cout << " Restart sieving ...\n" << std::flush;
1942 }
1943
1944 s1 = get_name ("RELATIONS");
1945 if (!(fpfull = fopen (s1, "a")))
1946 lidia_error_handler_n ("rational_factorization",
1947 "mpqs::can't open file RELATIONS");
1948
1949 s1 = get_name ("NEW_PARTIALS");
1950 if (!(fppart = fopen (s1, "a")))
1951 lidia_error_handler_n ("rational_factorization",
1952 "mpqs::can't open file "
1953 "NEW_PARTIALS");
1954 }
1955 }
1956 }
1957
1958 s1 = get_name ("RELATIONS");
1959 std::remove (s1);
1960 s1 = get_name ("LANCZOS_FORMAT");
1961 std::remove (s1);
1962 s1 = get_name ("NEW_PARTIALS");
1963 std::remove (s1);
1964 s1 = get_name ("NEW_FULLS");
1965 std::remove (s1);
1966 s1 = get_name ("LP_RELATIONS");
1967 std::remove (s1);
1968 s1 = get_name ("TEMPORARY");
1969 std::remove (s1);
1970 s1 = get_name ("SORT_LINE0");
1971 std::remove (s1);
1972 s1 = get_name ("SORT_LINE1");
1973 std::remove (s1);
1974 s1 = get_name ("SORT_LP0");
1975 std::remove (s1);
1976 s1 = get_name ("SORT_LP1");
1977 std::remove (s1);
1978
1979 delete[]LOGP;
1980 delete[]SQRTkN;
1981 delete[]START1;
1982 delete[]START2;
1983 delete[]sieb;
1984 delete[]CANDIDATE;
1985 delete[]a_inv;
1986 delete[]FB;
1987 delete[]faktor;
1988 delete[]Q_prime_glob;
1989 delete[]Q_prime;
1990 delete[]BG;
1991
1992 for (i = 0; i < static_cast < int >(P_TOTAL); i++)
1993 delete[]vorb[i];
1994 delete[]vorb;
1995
1996 return *this;
1997 }
1998
1999
2000
2001 //*********************************************************
2002 // mpqs function for component index of a factorization
2003 //*********************************************************
2004
2005
2006 rational_factorization & rational_factorization::
mpqs_comp(lidia_size_t index)2007 mpqs_comp (lidia_size_t index)
2008 {
2009 if ((index < 0) || (index > no_of_comp ()))
2010 lidia_error_handler_n ("rational_factorization",
2011 "mpqs::index out of range");
2012
2013 if (factors[index].factor_state == rf_single_factor::prime ||
2014 is_prime (factors[index].single_base, 8))
2015 {
2016 if (info)
2017 std::cout << "\nprime number : " << factors[index].
2018 single_base << "\n";
2019 factors[index].factor_state = rf_single_factor::prime;
2020 return *this;
2021 }
2022 long B = 200000;
2023
2024 if (decimal_length (factors[index].single_base) > 65)
2025 B = 500000;
2026
2027 ecm_primes prim (2, B, 100000);
2028
2029 trialdiv (index, 1, 180000, prim);
2030 prim.resetprimes (1);
2031
2032 if (is_prime_factor (index)) // fact. found
2033 return *this; // with TD
2034
2035 mpqs (index, prim);
2036 compose ();
2037 return *this;
2038 }
2039
2040
2041
mpqs(lidia_size_t index,ecm_primes & prim)2042 rational_factorization & rational_factorization::mpqs(lidia_size_t index,
2043 ecm_primes & prim) {
2044 // Compute complete factorization
2045 rf_single_factor& factor_ref = factors[index];
2046 rational_factorization f(factor_ref.single_base, 1);
2047
2048 bool failed = false;
2049 while(!(failed || f.is_prime_factorization())) {
2050 failed = true;
2051 lidia_size_t f_len = f.no_of_comp();
2052 for(lidia_size_t i = 0; i < f_len; ++i) {
2053 if(f.factors[i].factor_state == rf_single_factor::dont_know) {
2054 f.factors[i].factor_state =
2055 is_prime(f.factors[i].single_base) ?
2056 rf_single_factor::prime : rf_single_factor::not_prime;
2057 }
2058 if(f.factors[i].factor_state == rf_single_factor::not_prime) {
2059 prim.resetprimes(1);
2060 f.mpqs_impl(i, prim);
2061 f.refine();
2062 failed = false;
2063 break;
2064 }
2065 }
2066 }
2067
2068 lidia_size_t last_f_index = f.no_of_comp() - 1;
2069 for(lidia_size_t i = 0; i <= last_f_index; ++i) {
2070 f.factors[i].single_exponent *= factor_ref.single_exponent;
2071 }
2072 factor_ref = f.factors[last_f_index];
2073 f.factors.remove_from(last_f_index);
2074 factors.concat(factors, f.factors);
2075
2076 return *this;
2077 }
2078
2079 //********************************************************************
2080 // mpqs function for rational integers
2081 //********************************************************************
2082
mpqs(const bigint & N)2083 rational_factorization mpqs (const bigint & N)
2084 {
2085 rational_factorization f (N);
2086
2087 f.verbose (0);
2088 f.mpqs_comp (0);
2089 return f;
2090 }
2091
2092
2093
2094 #ifdef LIDIA_NAMESPACE
2095 } // end of namespace LiDIA
2096 #endif
2097