Annotation of OpenXM_contrib/pari-2.2/src/modules/mpqs.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: mpqs.c,v 1.11 2001/10/01 12:11:33 karim Exp $
2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /* Written by Thomas Papanikolaou and Xavier Roblot
17: *
18: * Implementation of the Self-Initializing Multi-Polynomial Quadratic Sieve
19: * based on code developed as part of the LiDIA project
20: * (http://www.informatik.tu-darmstadt.de/TI/LiDIA/Welcome.html)
21: */
22: #include "pari.h"
23:
24: #ifndef SEEK_SET
25: # define SEEK_SET 0
26: #endif
27:
28: #ifdef __CYGWIN32__
29: /* otherwise fseek() goes crazy due to silent \n <--> LF translations */
30: # define WRITE "wb"
31: # define READ "rb"
32: #else
33: # define WRITE "w"
34: # define READ "r"
35: #endif
36:
37: #ifdef MPQS_DEBUG_VERYVERBOSE
38: # ifndef MPQS_DEBUG_VERBOSE
39: # define MPQS_DEBUG_VERBOSE
40: # endif
41: #endif
42:
43: #ifdef MPQS_DEBUG_VERBOSE
44: # ifndef MPQS_DEBUG
45: # define MPQS_DEBUG
46: # endif
47: #endif
48:
49: #define MPQS_STRING_LENGTH (4 * 1024)
50: #define MPQS_CANDIDATE_ARRAY_SIZE 2000
51: #define MPQS_PRIMES_FOR_MULTIPLIER 5
52: /* `large primes' must be smaller than
53: * max(MPQS_LP_BOUND, largest_prime_in_FB) * MPQS_LP_FACTOR
54: */
55: #define MPQS_LP_BOUND 125000 /* works for 32 and 64bit */
56: #define MPQS_LP_FACTOR 80
57:
58: typedef ulong *mpqs_gauss_row;
59: typedef long *mpqs_gauss_indices;
60: typedef mpqs_gauss_row *mpqs_gauss_matrix;
61:
62: /* BEGIN: global variables to disappear as soon as possible */
63: /* names for temp. files, set in mpqs_get_filename */
64: static char *FREL_str;
65: static char *FNEW_str;
66: static char *LPREL_str;
67: static char *LPNEW_str;
68: static char *TMP_str;
69: static char *COMB_str;
70:
71: /* our own pointer to PARI's or to our own prime diffs table.
72: * NB the latter is never freed, unless we need to replace it with
73: * an even larger one. */
74: static byteptr mpqs_diffptr = NULL;
75: static long mpqs_prime_count = 0;
76:
77: /* END: global variables to disappear as soon as possible */
78:
79: /******************************/
80:
81: /**
82: ** determines a unique name for a file based on a short nickname
83: **/
84:
85: /* name is allocated on the stack */
86: static char *
87: mpqs_get_filename(char *s)
88: {
89: char *buf;
90:
91: s = pari_unique_filename(s);
92: buf = (char*) new_chunk((strlen(s) / sizeof(long)) + 2);
93: return strcpy(buf, s);
94: }
95:
96: /**
97: ** compares two `large prime' relations according to their
98: ** first element (the large prime itself).
99: **/
100:
101: static int
102: mpqs_relations_cmp(const void *a, const void *b)
103: {
104: char **sa = (char **) a;
105: char **sb = (char **) b;
106: long qa = strtol(*sa, NULL, 10);
107: long qb = strtol(*sb, NULL, 10);
108: /* atol() isn't entirely portable for the Full Relations case where the
109: strings of digits are too long to fit into a long --GN */
110: if (qa < qb) return -1;
111: else if (qa > qb) return 1;
112: else return strcmp(*sa, *sb);
113: }
114:
115: /**
116: ** Given a file "filename" containing full or `large prime' relations,
117: ** this function rearranges the file so that relations are sorted by
118: ** their first elements. Works also for sorting full relations.
119: ** Works in memory, discards duplicate lines, and overwrites the
120: ** original file.
121: **/
122:
123: #define bummer(fn) err(talker, "error whilst writing to file %s", fn)
124: #define min_bufspace 120 /* use new buffer when < min_bufspace left */
125: #define buflist_size 1024 /* size of list-of-buffers blocks */
126:
127: static long
128: mpqs_sort_lp_file(char *filename)
129: {
130: pariFILE *pTMP;
131: FILE *TMP;
132: char *old_s, **sort_table = (char**)avma, *buf, *cur_line;
133: static char **buflist_head = NULL;
134: char **buflist, **next_buflist;
135: long i, j, bufspace, length, count, av=avma;
136:
137: if (!buflist_head)
138: {
139: buflist_head = (char **) gpmalloc(buflist_size * sizeof(char *));
140: buflist = buflist_head;
141: *buflist++ = NULL; /* flag this as last and only buflist block */
142: }
143: else
144: buflist = buflist_head + 1;
145: /* this block will never be freed, and extra blocks may be allocated as
146: needed and linked ahead of it. NB it turns out that whilst extra
147: buflist blocks might have been needed when we were still sorting
148: entire FREL files (more than 1023 buffers, corresponding to about
149: 20000 lines of c.200 characters), they should never be touched now
150: that we only sort LPNEW and FNEW files, which are rather shorter.
151: But the code might as well stay around for paranoia, or for future
152: upgrades to handling even larger numbers (and factor bases and thus
153: relations files). It only costs one comparison per buffer allocation.
154: --GN */
155:
156: pTMP = pari_fopen(filename, READ);
157: TMP = pTMP->file;
158: /* get first buffer and read first line, if any, into it */
159: buf = (char *) gpmalloc(MPQS_STRING_LENGTH * sizeof(char));
160: cur_line = buf;
161: bufspace = MPQS_STRING_LENGTH;
162:
163: if (fgets(cur_line, bufspace, TMP) == NULL)
164: { /* file empty */
165: avma = av;
166: free(buf);
167: pari_fclose(pTMP);
168: return 0;
169: }
170: /* enter first buffer into buflist */
171: *buflist++ = buf; /* can't overflow the buflist block */
172: length = strlen(cur_line) + 1; /* count the \0 byte as well */
173: bufspace -= length;
174:
175: /* at start of loop, one line from the file is sitting in cur_line inside
176: buf, the next will go into cur_line + length, and there's room for
177: bufspace further characters in buf.
178: The loop reads another line if one exists, and if this overruns the
179: current buffer, it allocates a fresh one --GN */
180:
181: for (i=0, sort_table--; /* until end of file */; i++, sort_table--)
182: {
183: /* sort_table is allocated on the stack, 0x100 cells at a time. Hence the
184: * stack must be left alone in the rest of the loop to keep the array
185: * connected. In particular, buffers can't be new_chunk'ed --KB */
186: if ((i & 0xff) == 0) new_chunk(0x100);
187: *sort_table = cur_line;
188: cur_line += length;
189:
190: /* if very little room is left, allocate a fresh buffer before
191: attempting to read a line, and remember to free it if no
192: further line is forthcoming. This avoids some copying of
193: partial lines --GN */
194: if (bufspace < min_bufspace)
195: {
196: if (
197: #ifdef MPQS_DEBUG
198: 1 ||
199: #endif
200: DEBUGLEVEL >= 7)
201: fprintferr("MQPS: short of space -- another buffer for sorting\n");
202: buf = (char *) gpmalloc(MPQS_STRING_LENGTH * sizeof(char));
203: cur_line = buf;
204: bufspace = MPQS_STRING_LENGTH;
205: if (fgets(cur_line, bufspace, TMP) == NULL)
206: {
207: free(buf); break;
208: }
209: else
210: {
211: /* remember buffer for later deallocation */
212: if (buflist - buflist_head >= buflist_size)
213: {
214: /* need another buflist block */
215: next_buflist =
216: (char **) gpmalloc(buflist_size * sizeof(char *));
217: *next_buflist = (char *)buflist_head; /* link */
218: buflist_head = next_buflist;
219: buflist = buflist_head + 1;
220: }
221: *buflist++ = buf;
222: length = strlen(cur_line) + 1;
223: bufspace -= length;
224: continue;
225: }
226: }
227:
228: /* normal case: try fitting another line into the current buffer,
229: break if none exists */
230: if (fgets(cur_line, bufspace, TMP) == NULL) break;
231: length = strlen(cur_line) + 1;
232: bufspace -= length;
233:
234: /* check whether we got the entire line or only part of it */
235: if (bufspace == 0 && cur_line[length-2] != '\n')
236: {
237: long lg1;
238: if (
239: #ifdef MPQS_DEBUG
240: 1 ||
241: #endif
242: DEBUGLEVEL >= 7)
243: fprintferr("MQPS: line wrap -- another buffer for sorting\n");
244: buf = (char *) gpmalloc(MPQS_STRING_LENGTH * sizeof(char));
245: /* remember buffer for later deallocation */
246: if (buflist - buflist_head >= buflist_size)
247: {
248: /* need another buflist block */
249: next_buflist =
250: (char **) gpmalloc(buflist_size * sizeof(char *));
251: *next_buflist = (char *)buflist_head; /* link */
252: buflist_head = next_buflist;
253: buflist = buflist_head + 1;
254: }
255: *buflist++ = buf;
256:
257: /* copy what we've got to the new buffer */
258: (void)strcpy(buf, cur_line); /* cannot overflow */
259: cur_line = buf + length - 1; /* point at the \0 byte */
260: bufspace = MPQS_STRING_LENGTH - length + 1;
261: /* read remainder of line */
262: if (fgets(cur_line, bufspace, TMP) == NULL)
263: err(talker,"MPQS: relations file truncated?!\n");
264: lg1 = strlen(cur_line);
265: length += lg1; /* we already counted the \0 once */
266: bufspace -= (lg1 + 1); /* but here we must take it into account */
267: cur_line = buf; /* back up to the beginning of the line */
268: }
269: } /* for */
270:
271: pari_fclose(pTMP);
272:
273: #if 0 /* caught above, can no longer happen */
274: if (i == 0)
275: {
276: avma = av; return 0;
277: }
278: #endif
279:
280: /* sort the whole lot in place by swapping pointers */
281: qsort(sort_table, i, sizeof(char *), mpqs_relations_cmp);
282:
283: /* copy results back to the original file, skipping exact duplicates */
284: pTMP = pari_fopen(filename, WRITE); /* NOT safefopen */
285: TMP = pTMP->file;
286: old_s = sort_table[0];
287: if (fputs(sort_table[0], TMP) < 0)
288: bummer(filename);
289: count = 1;
290: for(j = 1; j < i; j++)
291: {
292: if (strcmp(old_s, sort_table[j]) != 0)
293: {
294: if (fputs(sort_table[j], TMP) < 0)
295: bummer(filename);
296: count++;
297: }
298: old_s = sort_table[j];
299: }
300: pari_fclose(pTMP);
301: if (
302: #ifdef MPQS_DEBUG
303: 1 ||
304: #endif
305: DEBUGLEVEL >= 6)
306: fprintferr("MPQS: done sorting one file.\n");
307:
308: /* deallocate buffers and any extraneous buflist blocks except the first */
309: while (*--buflist)
310: {
311: if (buflist != buflist_head) /* not a linkage pointer */
312: free((void *) *buflist); /* free a buffer */
313: else /* linkage pointer */
314: {
315: next_buflist = (char **)(*buflist);
316: free((void *)buflist_head); /* free a buflist block */
317: buflist_head = next_buflist;
318: buflist = buflist_head + buflist_size;
319: }
320: }
321: #if 0
322: for (j = 0; j < i; j++)
323: {
324: buf = (char *)(sort_table[j]) - 1;
325: if (*buf) /* check: deallocation flag or \0 ? */
326: free(buf);
327: }
328: #endif
329: avma = av; return count;
330: }
331:
332: /**
333: ** appends contents of file fp1 to file fp of name filename
334: ** (auxiliary routine for merge sort) and returns number of
335: ** lines copied.
336: ** fp assumed open for appending and fp1 for reading; fp will
337: ** be closed afterwards.
338: **/
339:
340: static long
341: mpqs_append_file(pariFILE *f, FILE *fp1)
342: {
343: FILE *fp = f->file;
344: char line[MPQS_STRING_LENGTH];
345: long c = 0;
346: while (fgets(line, MPQS_STRING_LENGTH, fp1) != NULL)
347: {
348: if (fputs(line, fp) < 0)
349: err(talker, "error whilst appending to file %s", f->name);
350: c++;
351: }
352: if (fflush(fp))
353: err(warner, "error whilst flushing file %s", f->name);
354: pari_fclose(f);
355: return c;
356: }
357:
358: /**
359: ** Does a merge-sort on the files LPREL and LPNEW;
360: ** assumes that LPREL and LPNEW are already sorted.
361: ** Creates/truncates the TMP file, writes result to it and closes it
362: ** (via a call to mpqs_append_file()).
363: ** Instead of LPREL, LPNEW we may also call this with FREL, FNEW.
364: ** In the latter case mode should be 1 (and we return the count of
365: ** all full relations), in the former 0 (and we return the count
366: ** of frels we expect to be able to combine out of the present lprels).
367: ** Moreover, if mode is 0, the combinable lprels are written out to a
368: ** separate file (COMB) which we also create and close (the caller
369: ** will find it exists iff our return value is nonzero), and we keep
370: ** only one occurrence of each `large prime' in TMP (i.e. in the future
371: ** LPREL file). --GN
372: **/
373:
374: #define swap_lines { \
375: line_tmp = line_new_old; \
376: line_new_old = line_new; \
377: line_new = line_tmp; \
378: }
379:
380: static long
381: mpqs_mergesort_lp_file0(FILE *LPREL, FILE *LPNEW, long mode)
382: {
383: /* TMP (renamed upon return) and COMB (deleted upon return) guaranteed not
384: * to exist yet --> safefopen for both temp files */
385: pariFILE *pTMP = pari_safefopen(TMP_str, WRITE);
386: FILE *TMP = pTMP->file;
387: pariFILE *pCOMB = NULL;
388: FILE *COMB = NULL; /* gcc -Wall */
389: char line1[MPQS_STRING_LENGTH], line2[MPQS_STRING_LENGTH];
390: char line[MPQS_STRING_LENGTH];
391: char *line_new = line1, *line_new_old = line2, *line_tmp;
392: long q_new, q_new_old = -1, q, i = 0, c = 0;
393: long comb_in_progress;
394:
395: /* if LPNEW is empty, copy LPREL to TMP. Could be done by a rename
396: if we didn't want to count the lines (again)... however, this
397: case will not normally happen */
398:
399: if (fgets(line_new, MPQS_STRING_LENGTH, LPNEW) == NULL)
400: {
401: i = mpqs_append_file(pTMP, LPREL);
402: return (mode ? i : 0);
403: }
404: /* we now have a line_new from LPNEW */
405:
406: /* if LPREL is empty, copy LPNEW to TMP... almost. */
407:
408: if (fgets(line, MPQS_STRING_LENGTH, LPREL) == NULL)
409: {
410: if (fputs(line_new, TMP) < 0)
411: bummer(TMP_str);
412: if (mode) /* full relations mode */
413: {
414: i = mpqs_append_file(pTMP, LPNEW);
415: return i + 1; /* COMB cannot be open here */
416: }
417:
418: /* LP mode: check for combinable relations */
419: q_new_old = atol(line_new);
420: /* we need to retain a copy of the old line just for a moment, because
421: we may yet have to write it to COMB, too. Do this simply by swapping
422: the two buffers */
423: swap_lines;
424: comb_in_progress = 0;
425: i = 0;
426:
427: while (fgets(line_new, MPQS_STRING_LENGTH, LPNEW) != NULL)
428: {
429: q_new = atol(line_new);
430: if (q_new_old == q_new)
431: {
432: /* found combinables, check whether we're already busy on this
433: particular `large prime' */
434: if (!comb_in_progress)
435: {
436: /* if not, write first line to COMB, creating and opening the
437: file first if it isn't open yet */
438: if (!pCOMB)
439: {
440: pCOMB = pari_safefopen(COMB_str, WRITE);
441: COMB = pCOMB->file;
442: }
443: if (fputs(line_new_old, COMB) < 0)
444: bummer(COMB_str);
445: comb_in_progress = 1;
446: }
447: /* in any case, write the current line, and count it */
448: if (fputs(line_new, COMB) < 0)
449: bummer(COMB_str);
450: i++;
451: }
452: else
453: { /* not combinable */
454: q_new_old = q_new;
455: comb_in_progress = 0;
456: /* and dump it to the TMP file */
457: if (fputs(line_new, TMP) < 0)
458: bummer(TMP_str);
459: /* and stash it away for a moment */
460: swap_lines;
461: comb_in_progress = 0;
462: }
463: } /* while */
464: if (pCOMB) pari_fclose(pCOMB);
465: pari_fclose(pTMP);
466: return i;
467: }
468:
469: /* normal case: both LPNEW and LPREL are not empty */
470: q_new = atol(line_new);
471: q = atol(line);
472:
473: for(;;) /* main merging loop */
474: {
475: comb_in_progress = 0;
476: i = 0;
477:
478: /* first the harder case: let LPNEW catch up with LPREL, and possibly
479: overtake it, checking for combinables coming from LPNEW alone */
480:
481: while (q > q_new)
482: {
483: if (mode || !comb_in_progress)
484: {
485: if (fputs(line_new, TMP) < 0)
486: bummer(TMP_str);
487: }
488: if (mode)
489: c++; /* in FREL mode, count lines written */
490: else if (!comb_in_progress)
491: {
492: q_new_old = q_new;
493: swap_lines;
494: }
495: if (fgets(line_new, MPQS_STRING_LENGTH, LPNEW) == NULL)
496: {
497: if (fputs(line, TMP) < 0)
498: bummer(TMP_str);
499: if (mode) c++; else c += i;
500: i = mpqs_append_file(pTMP, LPREL);
501: if (pCOMB) pari_fclose(pCOMB);
502: return (mode ? c + i : c);
503: }
504: q_new = atol(line_new);
505: if (mode) continue;
506: /* LP mode only: */
507: if (q_new_old != q_new)
508: { /* not combinable */
509: comb_in_progress = 0;
510: /* next loop iteration will deal with it, or the loop may end */
511: }
512: else
513: {
514: /* found combinables, check whether we're already busy on this
515: `large prime' */
516: if (!comb_in_progress)
517: {
518: if (!pCOMB)
519: {
520: pCOMB = pari_safefopen(COMB_str, WRITE);
521: COMB = pCOMB->file;
522: }
523: if (fputs(line_new_old, COMB) < 0)
524: bummer(COMB_str);
525: comb_in_progress = 1;
526: }
527: /* in any case, write the current line, and count it */
528: if (fputs(line_new, COMB) < 0)
529: bummer(COMB_str);
530: i++;
531: }
532: } /* while q > q_new */
533:
534: /* q <= q_new */
535:
536: if (!mode) c += i; /* accumulate count of combinables */
537: i = 0; /* and clear it */
538: comb_in_progress = 0; /* redundant */
539:
540: /* now let LPREL catch up with LPNEW, and possibly overtake it */
541:
542: while (q < q_new)
543: {
544: if (fputs(line, TMP) < 0)
545: bummer(TMP_str);
546: if (mode) c++;
547: if (fgets(line, MPQS_STRING_LENGTH, LPREL) == NULL)
548: {
549: if (fputs(line_new, TMP) < 0)
550: bummer(TMP_str);
551: i = mpqs_append_file(pTMP, LPNEW);
552: if (pCOMB) pari_fclose(pCOMB);
553: return (mode ? c + i + 1 : c);
554: }
555: else
556: q = atol(line);
557: }
558:
559: /* q >= q_new */
560:
561: /* Finally, it may happen that q == q_new, indicating combinables whose
562: `large prime' is already represented in LPREL, and appears now once
563: or more often in LPNEW. Thus in this sub-loop we advance LPNEW.
564: The `line' from LPREL is left alone, and will be written to TMP the
565: next time around the main for loop; we only write it to COMB here --
566: unless all we find is an exact duplicate of the line we already have,
567: that is. (There can be at most one such, and if so it will simply be
568: discarded, whether we are in LP or full relations mode.) */
569:
570: while (q == q_new)
571: {
572: if (strcmp(line_new, line) == 0)
573: {
574: /* duplicate -- move right ahead to the next LPNEW line */
575: ; /* do nothing here */
576: } /* end duplicate case */
577: else if (mode)
578: {
579: /* full relations mode: write line_new out first, keep line */
580: if (fputs(line_new, TMP) < 0)
581: bummer(TMP_str);
582: c++;
583: } /* end full relations case */
584: else
585: {
586: /* LP mode, and combinable relation */
587: if (!comb_in_progress)
588: {
589: if (!pCOMB)
590: {
591: pCOMB = pari_safefopen(COMB_str, WRITE);
592: COMB = pCOMB->file;
593: }
594: if (fputs(line, COMB) < 0)
595: bummer(COMB_str);
596: comb_in_progress = 1;
597: }
598: if (fputs(line_new, COMB) < 0)
599: bummer(COMB_str);
600: i++;
601: } /* end combinable LP case */
602: /* NB comb_in_progress is cleared by q_new becoming bigger than q,
603: and thus the current while loop terminating, the next time through
604: the main for loop */
605:
606: /* common ending: get another line_new, if any */
607: if (fgets(line_new, MPQS_STRING_LENGTH, LPNEW) == NULL)
608: {
609: if (fputs(line, TMP) < 0)
610: bummer(TMP_str);
611: if (mode) c++; else c += i;
612: i = mpqs_append_file(pTMP, LPREL);
613: if (pCOMB) pari_fclose(pCOMB);
614: return (mode ? c + i : c);
615: }
616: else
617: q_new = atol(line_new);
618: } /* while */
619:
620: if (!mode) c += i; /* accumulate count of combinables */
621: } /* for */
622: }
623:
624: static long
625: mpqs_mergesort_lp_file(char *REL_str, char *NEW_str, long mode)
626: {
627: pariFILE *pREL = pari_fopen(REL_str, READ);
628: pariFILE *pNEW = pari_fopen(NEW_str, READ);
629: long tp = mpqs_mergesort_lp_file0(pREL->file, pNEW->file, mode);
630:
631: pari_fclose(pREL);
632: pari_fclose(pNEW);
633: pari_unlink(REL_str);
634: if (rename(TMP_str,REL_str))
635: err(talker, "can\'t rename file %s to %s", TMP_str, REL_str);
636: if (
637: #ifdef MPQS_DEBUG
638: 1 ||
639: #endif
640: DEBUGLEVEL >= 6)
641: fprintferr("MPQS: renamed file %s to %s\n", TMP_str, REL_str);
642: return tp;
643: }
644:
645: /******************************/
646: /**
647: ** return next prime larger than p, using *primes_ptr on the
648: ** diffptr table first and pari's other wits after that
649: **/
650:
651: static byteptr
652: mpqs_iterate_primes(long *p, byteptr primes_ptr)
653: {
654: long prime = *p;
655: if (*primes_ptr)
656: prime += *primes_ptr++;
657: else
658: {
659: long av = avma;
660: prime = itos(nextprime(stoi(prime + 1)));
661: avma = av;
662: }
663: *p = prime;
664: return primes_ptr;
665: }
666:
667:
668: /**
669: ** return the multiplier k for N to use in MPQS
670: **/
671:
672: long
673: mpqs_find_k(GEN N, long tries)
674: {
675: static long cand_table[] = { 1,3,5,7,11,13,15,17,19,21,23,
676: 29,31,33,35,37,39,41,43,47,51,53,
677: 55,57,59,61,65,67,69,71,73,79,83,
678: 85,89,91,93,95,97,101,103,105,107,
679: 109,113,127,131,137,139};
680: long av = avma;
681: long best_k = 1, k, p, N_mod_4 = smodis(N, 4), x;
682: GEN kN;
683: double best_value = 1, value, dp;
684: long i, j;
685: byteptr primes_ptr;
686:
687: for (i = 0; i < MPQS_PRIMES_FOR_MULTIPLIER; i++)
688: {
689: k = cand_table[i];
690: x = (k * N_mod_4) % 4;
691: if (x == 1)
692: {
693: value = -0.7 * log2 ((double) k);
694: kN = mulis(N, k);
695: if (smodis(kN, 8) == 1)
696: value += 1.38629;
697:
698: j = 0; p = 0;
699: primes_ptr = diffptr;
700: while (j <= tries)
701: {
702: primes_ptr = mpqs_iterate_primes(&p, primes_ptr);
703: if (kross(smodis(kN, p), p) == 1)
704: {
705: j++;
706: dp = log2((double) p) / p;
707: value += (k % p == 0) ? dp : 2 * dp;
708: }
709: }
710:
711: if (value > best_value)
712: {
713: best_value = value;
714: best_k = k;
715: }
716: }
717: }
718: avma = av; return best_k;
719: }
720:
721: /******************************/
722: /**
723: ** guesstimate up to what size we're going to need precomputed
724: ** small primes
725: **/
726:
727: static long
728: mpqs_find_maxprime(long size)
729: {
730: double x;
731:
732: if (size < 16000)
733: return 176000;
734:
735: x = log((double)size);
736: x += log(x);
737: x -= 0.9427;
738: x *= size;
739:
740: return (long)x;
741: }
742:
743: /**
744: ** return the number of primes initialized in PARI
745: **/
746:
747: static long
748: mpqs_count_primes()
749: {
750: byteptr p = mpqs_diffptr;
751:
752: for ( ; *p; p++) /* empty */;
753: return (p - mpqs_diffptr);
754: }
755:
756: /**
757: ** create a factor base of size primes p_i with the
758: ** property that legendre(k*N, p_i) == 0, 1
759: ** FB[0] contains the size of the factor base; FB[FB[0]+1] is the
760: ** last prime p_i.
761: ** FB[1] contains -1 which allows to handle negative values of Q(x).
762: ** If a prime factor of N is found during the construction, it
763: ** is returned in f, otherwise f = 0 (such a factor will necessarily
764: ** fit into a C long).
765: **/
766:
767: static long*
768: mpqs_create_FB(long size, GEN kN, long k, long *f)
769: {
770: long i, kr, p = 0;
771: long *FB;
772: byteptr primes_ptr;
773:
774: FB = (long *) calloc(size + 3, sizeof(long));
775: if (!FB) err(memer);
776: FB[0] = size;
777: FB[1] = -1;
778: /* pick up PARI's current differences-of-small-primes array */
779: /* unless we already have our own. We must then reset it to NULL */
780: /* so we will not use a stale one the next time we are called */
781: if (!mpqs_diffptr) mpqs_diffptr = diffptr;
782:
783: if ((mpqs_prime_count ? mpqs_prime_count : mpqs_count_primes())
784: < 3 * size)
785: {
786: long newsize = 3 * mpqs_find_maxprime(size);
787: if (mpqs_diffptr != diffptr) free((void *) mpqs_diffptr);
788: if (DEBUGLEVEL >= 2)
789: {
790: fprintferr("MPQS: precomputing auxiliary primes up to %ld\n",
791: newsize);
792: }
793: mpqs_diffptr = initprimes(newsize);
794: /* count 'em once and remember the result */
795: mpqs_prime_count = mpqs_count_primes();
796: }
797:
798: if (
799: #ifdef MPQS_DEBUG
800: 1 ||
801: #endif
802: DEBUGLEVEL >= 7)
803: fprintferr("MPQS: FB [-1");
804: primes_ptr = mpqs_diffptr;
805: for (i = 2; i < size + 2; )
806: {
807: primes_ptr = mpqs_iterate_primes(&p, primes_ptr);
808:
809: /* testing p!=k was only correct because MPQS_PRIMES_FOR_MULTIPLIER
810: is so small that the composite values k > 1 will never be chosen.
811: Fixing this just in case someone increases that parameter one day...
812: --GN */
813: if (p > k || k%p)
814: {
815: kr = kross(smodis(kN, p), p);
816: if (kr != -1)
817: {
818: if (kr == 0)
819: {
820: if (
821: #ifdef MPQS_DEBUG
822: 1 ||
823: #endif
824: DEBUGLEVEL >= 7)
825: fprintferr(",%ld...] Wait a second --\n", p);
826: *f = p;
827: return FB;
828: }
829: if (
830: #ifdef MPQS_DEBUG
831: 1 ||
832: #endif
833: DEBUGLEVEL >= 7)
834: fprintferr(",%ld", p);
835: FB[i] = p, i++;
836: }
837: }
838: }
839:
840: if (
841: #ifdef MPQS_DEBUG
842: 1 ||
843: #endif
844: DEBUGLEVEL >= 7)
845: {
846: fprintferr("]\n");
847: flusherr();
848: }
849:
850: FB[i] = 0;
851: if (DEBUGLEVEL >= 6)
852: fprintferr("MPQS: last available index in FB is %ld\n", i - 1);
853:
854: *f = 0;
855: return FB;
856: }
857:
858: /******************************/
859:
860: /**
861: ** see below for the format of the parameter table.
862: ** These values are somewhat experimental and can be improved by
863: ** further extensive testing.
864: ** Tentative entries for 9...14 digits by GN, entry for 15 changed
865: ** Tentative entries for 81...100 digits by GN, entry for 74 changed
866: ** total_no_of_primes_for_A (penultimate entries in each row) increased
867: ** throughout, partly so as to offset the effect of changed bin_index
868: ** increment strategy --GN
869: ** Added entries to govern sort intervals, commented out unused numbers
870: ** of digits --GN
871: **/
872:
873: static double
874: mpqs_parameters[99][8] =
875: {
876: { /*9*/ 0.8, 900, 20, 3, 13, 2, 70, 8},
877: { /*10*/ 0.8, 900, 21, 3, 12, 2, 70, 8},
878: { /*11*/ 0.8, 920, 22, 3, 12, 2, 70, 6},
879: { /*12*/ 0.8, 960, 24, 3, 12, 2, 70, 6},
880: { /*13*/ 0.8, 1020, 26, 3, 12, 2, 70, 6},
881: { /*14*/ 0.8, 1100, 29, 3, 12, 2, 70, 6},
882: { /*15*/ 0.8, 1200, 32, 3, 12, 2, 60, 8},
883: { /*16*/ 0.8, 1400, 35, 3, 12, 2, 60, 8},
884: { /*17*/ 0.8, 3000, 40, 3, 12, 2, 60, 8},
885: { /*18*/ 0.8, 3000, 60, 3, 12, 2, 50, 10},
886: { /*19*/ 0.8, 3600, 80, 3, 13, 2, 50, 10},
887: { /*20*/ 0.8, 4000, 100, 3, 13, 2, 40, 10},
888: { /*21*/ 0.8, 4300, 100, 3, 13, 2, 40, 10},
889: { /*22*/ 0.8, 4500, 120, 3, 13, 3, 40, 10},
890: { /*23*/ 0.8, 4800, 140, 3, 14, 3, 30, 10},
891: { /*24*/ 0.8, 5000, 160, 3, 14, 4, 30, 10},
892: { /*25*/ 0.8, 5000, 180, 3, 14, 4, 30, 10},
893: { /*26*/ 0.9, 6000, 200, 5, 10, 4, 30, 10},
894: { /*27*/ 1.17, 6000, 220, 5, 10, 5, 30, 10},
895: { /*28*/ 1.17, 6500, 240, 5, 10, 5, 30, 10},
896: { /*29*/ 1.17, 6500, 260, 5, 10, 5, 30, 10},
897: { /*30*/ 1.36, 7000, 325, 5, 10, 5, 20, 10},
898: { /*31*/ 1.36, 7000, 355, 5, 10, 5, 20, 10},
899: { /*32*/ 1.36, 7500, 375, 5, 10, 5, 20, 10},
900: { /*33*/ 1.43, 7500, 400, 6, 11, 6, 20, 10},
901: { /*34*/ 1.43, 7500, 425, 6, 11, 6, 20, 10},
902: { /*35*/ 1.43, 7500, 550, 6, 11, 6, 20, 10},
903: { /*36*/ 1.43, 8000, 650, 6, 11, 6, 20, 10},
904: { /*37*/ 1.69, 9000, 750, 6, 11, 7, 20, 10},
905: { /*38*/ 1.69, 10000, 850, 6, 11, 7, 20, 10},
906: { /*39*/ 1.69, 11000, 950, 6, 11, 7, 20, 10},
907: { /*40*/ 1.69, 14000, 1000, 6, 11, 7, 20, 10},
908: { /*41*/ 1.69, 14000, 1150, 6, 11, 8, 10, 10},
909: { /*42*/ 1.69, 15000, 1300, 6, 11, 8, 10, 10},
910: { /*43*/ 1.69, 15000, 1600, 6, 11, 8, 10, 10},
911: { /*44*/ 1.69, 15000, 1900, 7, 12, 9, 10, 10},
912: { /*45*/ 1.69, 15000, 2200, 7, 12, 9, 10, 10},
913: { /*46*/ 1.69, 20000, 2500, 7, 12, 9, 10, 10},
914: { /*47*/ 1.69, 25000, 2500, 7, 12, 10, 10, 10},
915: { /*48*/ 1.69, 27500, 2700, 7, 12, 10, 10, 10},
916: { /*49*/ 1.69, 30000, 2800, 7, 12, 10, 10, 10},
917: { /*50*/ 1.75, 35000, 2900, 7, 12, 10, 10, 10},
918: { /*51*/ 1.75, 40000, 3000, 7, 12, 10, 10, 10},
919: { /*52*/ 1.85, 50000, 3200, 7, 12, 11, 10, 10},
920: { /*53*/ 1.85, 50000, 3500, 7, 12, 11, 10, 10},
921: { /*54*/ 1.95, 80000, 3800, 7, 13, 11, 10, 10},
922: { /*55*/ 1.95, 90000, 4100, 7, 13, 11, 10, 10},
923: { /*56*/ 1.95, 100000, 4400, 7, 13, 11, 10, 10},
924: { /*57*/ 1.95, 80000, 4700, 8, 14, 12, 10, 10},
925: { /*58*/ 1.95, 80000, 5000, 8, 14, 12, 10, 10},
926: { /*59*/ 2.15, 130000, 5500, 8, 14, 12, 10, 8},
927: { /*60*/ 2.15, 140000, 5800, 8, 14, 12, 10, 8},
928: { /*61*/ 2.15, 150000, 6100, 8, 14, 13, 10, 8},
929: { /*62*/ 2.15, 160000, 6400, 8, 14, 13, 10, 8},
930: { /*63*/ 2.35, 170000, 6700, 8, 14, 13, 10, 8},
931: { /*64*/ 2.35, 180000, 7000, 8, 14, 13, 10, 8},
932: { /*65*/ 2.35, 190000, 7300, 8, 14, 13, 10, 8},
933: { /*66*/ 2.35, 200000, 7600, 8, 14, 13, 10, 8},
934: { /*67*/ 2.4, 150000, 7900, 8, 14, 13, 10, 8},
935: { /*68*/ 2.4, 150000, 8200, 8, 15, 13, 10, 8},
936: { /*69*/ 2.4, 130000, 8600, 8, 15, 13, 8, 6},
937: { /*70*/ 2.45, 130000, 8800, 8, 15, 13, 8, 6},
938: { /*71*/ 2.45, 130000, 8800, 9, 16, 13, 8, 6},
939: { /*72*/ 2.4, 260000, 9400, 9, 16, 13, 5, 5},
940: { /*73*/ 2.4, 270000, 9700, 9, 16, 13, 5, 5},
941: { /*74*/ 2.4, 280000, 9800, 9, 16, 13, 5, 5},
942: { /*75*/ 2.6, 140000, 9000, 9, 17, 13, 5, 5},
943: { /*76*/ 2.6, 160000, 9400, 9, 17, 13, 5, 5},
944: { /*77*/ 2.6, 180000, 9600, 9, 17, 13, 5, 5},
945: { /*78*/ 2.6, 200000, 9800, 9, 17, 13, 5, 5},
946: { /*79*/ 2.65, 220000, 10000, 9, 18, 13, 5, 5},
947: { /*80*/ 2.65, 250000, 10500, 9, 18, 13, 5, 5},
948: /* entries below due to Thomas Denny */
949: { /*81*/ 3.1, 1500000, 16000, 9, 18, 15, 4, 4},
950: { /*82*/ 3.1, 1500000, 17000, 9, 18, 15, 4, 4},
951: { /*83*/ 3.1, 1500000, 18500, 9, 19, 16, 4, 4},
952: { /*84*/ 3.2, 1500000, 20000, 9, 19, 16, 4, 4},
953: { /*85*/ 3.2, 1500000, 21500, 9, 19, 16, 4, 4},
954: { /*86*/ 3.2, 1500000, 23000, 9, 19, 16, 4, 4},
955: { /*87*/ 3.2, 1500000, 25000, 9, 20, 17, 4, 4},
956: { /*88*/ 3.3, 1500000, 27000, 9, 20, 17, 4, 4},
957: { /*89*/ 3.3, 1500000, 28000, 9, 20, 17, 4, 4},
958: { /*90*/ 3.5, 2000000, 30500, 9, 20, 17, 2, 2},
959: { /*91*/ 3.6, 2000000, 32000, 9, 21, 18, 2, 2},
960: { /*92*/ 3.6, 2000000, 35000, 9, 21, 18, 2, 2},
961: { /*93*/ 3.7, 2000000, 37000, 9, 21, 18, 2, 2},
962: { /*94*/ 3.7, 2000000, 39500, 9, 22, 18, 2, 2},
963: { /*95*/ 3.7, 2500000, 41500, 9, 22, 18, 2, 2},
964: { /*96*/ 3.8, 2500000, 45000, 10, 23, 18, 2, 2},
965: { /*97*/ 3.8, 2500000, 47500, 10, 23, 18, 2, 2},
966: { /*98*/ 3.7, 3000000, 51000, 10, 24, 18, 2, 2},
967: { /*99*/ 3.8, 3000000, 53000, 10, 24, 18, 2, 2},
968: {/*100*/ 3.8, 3500000, 51000, 10, 25, 18, 2, 2},
969: {/*101*/ 3.8, 3500000, 54000, 10, 25, 18, 2, 2},
970: {/*102*/ 3.8, 3500000, 57000, 10, 26, 18, 2, 2},
971: {/*103*/ 3.9, 4000000, 61000, 10, 26, 18, 2, 2},
972: {/*104*/ 3.9, 4000000, 66000, 10, 27, 18, 2, 2},
973: {/*105*/ 3.9, 4000000, 70000, 10, 27, 18, 2, 2},
974: {/*106*/ 3.9, 4000000, 75000, 10, 28, 18, 2, 2},
975: {/*107*/ 3.9, 4000000, 80000, 10, 30, 18, 2, 2},
976: };
977:
978: /**
979: ** Return the appropriate parameters for the initialization of MPQS
980: **
981: ** the parameter table for mpqs has the following format:
982: **
983: ** [d -- decimal digits -- commented out]
984: ** e -- approximation accuracy
985: ** M -- size sieving interval
986: ** s -- size FB,
987: ** fA -- # prime factors in A
988: ** pA -- total # primes for determination of A,
989: ** [maxB -- how many Bs to use before choosing a new A -- determined by fA]
990: ** si -- initial sieving index
991: ** isoi -- initial sorting interval
992: ** fsoi -- followup sorting interval increment
993: **/
994:
995: static void
996: mpqs_read_parameters(ulong d, double *e, ulong *M, ulong *s, ulong *os,
997: ulong *fA, ulong *pA, long *maxB, ulong *si,
998: ulong *isoi, ulong *fsoi)
999: {
1000: long i = d;
1001: if (i < 9) i = 9;
1002: if (i > 107) i = 107;
1003: i -= 9;
1004:
1005: *e = mpqs_parameters[i][0];
1006: *M = (ulong)mpqs_parameters[i][1];
1007: *s = (ulong)mpqs_parameters[i][2];
1008: if (*s >= 200) *os = *s + 70;
1009: else *os = (ulong)(*s * 1.35);
1010: *fA = (ulong)mpqs_parameters[i][3];
1011: *pA = (ulong)mpqs_parameters[i][4];
1012: #if 0
1013: *maxB = (long) pow(2.0, mpqs_parameters[i][3] - 1.0);
1014: #else
1015: *maxB = 1l << (*fA - 1);
1016: #endif
1017: *si = (ulong)mpqs_parameters[i][5];
1018: *isoi = 10 * (ulong)mpqs_parameters[i][6]; /* percentages scaled by 10 */
1019: *fsoi = 10 * (ulong)mpqs_parameters[i][7];
1020: }
1021:
1022: /******************************/
1023:
1024: /**
1025: ** increment *x (!=0) to a larger value which has the same number of 1s
1026: ** in its binary representation. Cases which wrap around can be detected
1027: ** by the caller afterwards as long as we keep total_no_of_primes_for_A
1028: ** strictly less than BITS_IN_LONG.
1029: **
1030: ** Switch for the fast cases moved into here --GN
1031: **
1032: ** Changed switch to increment *x in all cases to the next larger number
1033: ** which (a) has the same count of 1 bits and (b) does not arise from the
1034: ** old value by moving a single 1 bit one position to the left (which was
1035: ** undesirable for the sieve). --GN based on discussion with TP
1036: **/
1037: static void
1038: mpqs_increment(ulong *x)
1039: {
1040: ulong r1_mask, r01_mask, slider=1;
1041:
1042: /* ultra-fast 32-way computed jump handles 22 out of 32 cases */
1043: switch (*x & 0x1f)
1044: {
1045: case 29:
1046: (*x)++; break; /* this does shift a single bit, but we'll
1047: postprocess this case */
1048: case 26:
1049: (*x) += 2; break; /* again */
1050: case 1: case 3: case 6: case 9: case 11:
1051: case 17: case 19: case 22: case 25: case 27:
1052: (*x) += 3; return;
1053: case 20:
1054: (*x) += 4; break; /* again */
1055: case 5: case 12: case 14: case 21:
1056: (*x) += 5; return;
1057: case 2: case 7: case 13: case 18: case 23:
1058: (*x) += 6; return;
1059: case 10:
1060: (*x) += 7; return;
1061: case 8:
1062: (*x) += 8; break; /* and again */
1063: case 4: case 15:
1064: (*x) += 12; return;
1065: default: /* 0, 16, 24, 28, 30, 31 */
1066: /* isolate rightmost 1 */
1067: r1_mask = ((*x ^ (*x - 1)) + 1) >> 1;
1068: /* isolate rightmost 1 which has a 0 to its left */
1069: r01_mask = ((*x ^ (*x + r1_mask)) + r1_mask) >> 2;
1070: /* simple cases. Both of these shift a single bit one position to the
1071: left, and will need postprocessing */
1072: if (r1_mask == r01_mask) { *x += r1_mask; break; }
1073: if (r1_mask == 1) { *x += r01_mask; break; }
1074: /* general case -- idea:
1075: add r01_mask, kill off as many 1 bits as possible to its right
1076: while at the same time filling in 1 bits from the LSB. */
1077: if (r1_mask == 2) { *x += (r01_mask>>1) + 1; return; }
1078: while (r01_mask > r1_mask && slider < r1_mask)
1079: {
1080: r01_mask >>= 1; slider <<= 1;
1081: }
1082: *x += r01_mask + slider - 1;
1083: return;
1084: }
1085: /* post-process all cases which couldn't be finalized above. If we get
1086: here, slider still has its original value. */
1087: r1_mask = ((*x ^ (*x - 1)) + 1) >> 1;
1088: r01_mask = ((*x ^ (*x + r1_mask)) + r1_mask) >> 2;
1089: if (r1_mask == r01_mask) { *x += r1_mask; return; }
1090: if (r1_mask == 1) { *x += r01_mask; return; }
1091: if (r1_mask == 2) { *x += (r01_mask>>1) + 1; return; }
1092: while (r01_mask > r1_mask && slider < r1_mask)
1093: {
1094: r01_mask >>= 1; slider <<= 1;
1095: }
1096: *x += r01_mask + slider - 1;
1097: return;
1098: }
1099:
1100: static long
1101: mpqs_invsmod(long b, long p)
1102: {
1103: long v1 = 0, v2 = 1, v3, r, oldp = p;
1104:
1105: while (b > 1)
1106: {
1107: v3 = v1 - (p / b) * v2; v1 = v2; v2 = v3;
1108: r = p % b; p = b; b = r;
1109: }
1110:
1111: if (v2 < 0) v2 += oldp;
1112: return v2;
1113: }
1114:
1115: /**
1116: ** compute coefficients of the sieving polynomial for self initializing
1117: ** variant. Coefficients A and B are returned and several tables are
1118: ** updated. See Thomas Sosnowski's diploma thesis.
1119: **/
1120:
1121: static void
1122: mpqs_self_init(GEN A, GEN B, GEN N, GEN kN, long *FB, long *sqrt_mod_p_kN,
1123: long *start_1, long *start_2, ulong no_of_primes_in_A,
1124: ulong total_no_of_primes_for_A, ulong M, long **inv_A_H_i,
1125: long *Q_prime, long *Q_prime_glob, GEN* H_i, long index_i,
1126: long start_index_FB_for_A, long *inv_A2, GEN inv_A4,
1127: ulong *bin_index, GEN *f)
1128: {
1129: long av;
1130: GEN p1, p2;
1131: ulong size_of_FB, nu_2;
1132: long i, j, p, M_mod_p, tmp, tmp1, tmp2;
1133:
1134: *f = NULL;
1135: if (index_i == 0)
1136: {
1137: /* compute first polynomial with new A */
1138: if (!*bin_index)
1139: {
1140: *bin_index = (1ul << no_of_primes_in_A) - 1;
1141: }
1142: else
1143: mpqs_increment(bin_index);
1144: if (*bin_index >= (1ul << total_no_of_primes_for_A)) /* wraparound */
1145: {
1146: if (DEBUGLEVEL >= 2)
1147: fprintferr("MPQS: bin_index wraparound\n");
1148: /* complain: caller should increase some parameters... */
1149: return;
1150: /* caller should repeat the above test; if we simply increase the
1151: total_no_of_primes_for_A, then bin_index now contains something
1152: suitable for carrying on -- we'll probably end up not using
1153: that particular index and skipping ahead to the next one, but
1154: this isn't a problem. OTOH total_no_of_primes_for_A should now
1155: be so large that this should never happen except for ridiculous
1156: input, like asking mpqs to factor 27, in which case it will act
1157: as a safeguard against an infinite loop --GN */
1158: }
1159:
1160: /* determine primes used for A in this iteration */
1161: for (i = 0, j = 0; (ulong)j < total_no_of_primes_for_A; j++)
1162: {
1163: if (*bin_index & (1 << j))
1164: {
1165: Q_prime[i] = Q_prime_glob[j];
1166: i++;
1167: }
1168: }
1169:
1170: /*
1171: * compute coefficient A using
1172: * A = prod (Q_prime[i], 0 < i < no_of_primes_in_A)
1173: */
1174:
1175: av = avma;
1176: p1 = stoi(Q_prime[0]);
1177: for (i = 1; (ulong)i < no_of_primes_in_A; i++)
1178: p1 = mulis(p1, Q_prime[i]);
1179: affii(p1, A);
1180: avma = av;
1181:
1182: /*
1183: * let bound := no_of_primes_in_A. Compute H_i, 0 <= i < bound.
1184: */
1185:
1186: for (i = 0; (ulong)i < no_of_primes_in_A; i++)
1187: {
1188: av = avma;
1189: p = Q_prime[i];
1190: p1 = divis(A, p);
1191: tmp = smodis(p1, p);
1192: p1 = mulis(p1, mpqs_invsmod(tmp, p));
1193: tmp = itos(mpsqrtmod(modis(kN, p), stoi(p)));
1194: p1 = mulis(p1, tmp);
1195: affii(modii(p1, A), H_i[i]);
1196: avma = av;
1197: }
1198:
1199: /*
1200: * Each B is of the form sum(v_i * H_i, 0 <= i < bound)
1201: * where each v_j is +1 or -1. This has to be done only once
1202: * (for index_i==0) for the coefficient A; if index_i > 0 then
1203: * there is a linear recursion for B
1204: */
1205:
1206: /* compute actual B coefficient, by taking all v_i == 1 */
1207: av = avma;
1208: p1 = H_i[0];
1209: for (i = 1; (ulong)i < no_of_primes_in_A; i++)
1210: p1 = addii(p1, H_i[i]);
1211: affii(p1, B);
1212: avma = av;
1213:
1214: /* ensure B = 1 mod 4 */
1215: if (mod2(B) == 0)
1216: {
1217: /* B += (A % 4) * A; */
1218: av = avma;
1219: p1 = addii(B, mulsi(mod4(A), A));
1220: affii(p1, B);
1221: avma = av;
1222: }
1223:
1224: /* compute inv_A2[i] = 1/(2*A) mod p_i */
1225: av = avma;
1226: p1 = shifti(A, 1);
1227: size_of_FB = FB[0] + 1;
1228: for (i = 2; (ulong)i <= size_of_FB; i++)
1229: inv_A2[i] = mpqs_invsmod(smodis(p1, FB[i]), FB[i]);
1230: avma = av;
1231:
1232: /* compute inv_A_H_i[i][j] = 1/A * H_i[i] mod p_j */
1233: for (i = 0; (ulong)i < no_of_primes_in_A; i++)
1234: {
1235: for (j = 2; (ulong)j <= size_of_FB; j++)
1236: {
1237: av = avma;
1238: p1 = mulis(H_i[i], inv_A2[j] << 1);
1239: inv_A_H_i[i][j] = smodis(p1, FB[j]);
1240: avma = av;
1241: }
1242: }
1243:
1244: /* compute the roots z1, z2, of the polynomial Q(x) mod p_j and
1245: * initialize start_1[i] with the first value p_i | Q(z1 + i p_j)
1246: * initialize start_2[i] with the first value p_i | Q(z2 + i p_j)
1247: */
1248: for (j = 2; (ulong)j <= size_of_FB; j++)
1249: {
1250: p = FB[j];
1251: M_mod_p = M % p;
1252:
1253: tmp = p - smodis(B, p); /* tmp = -B mod p */
1254:
1255: tmp1 = (tmp - sqrt_mod_p_kN[j]) % p;
1256: if (tmp1 < 0) tmp1 += p;
1257: tmp1 = mulssmod(tmp1, inv_A2[j], p);
1258: tmp1 = (tmp1 + M_mod_p) % p;
1259:
1260: tmp2 = (tmp + sqrt_mod_p_kN[j]) % p;
1261: if (tmp2 < 0) tmp2 += p;
1262: tmp2 = mulssmod(tmp2, inv_A2[j], p);
1263: tmp2 = (tmp2 + M_mod_p) % p;
1264:
1265: start_1[j] = (tmp1 < 0) ? tmp1 + p : tmp1;
1266: start_2[j] = (tmp2 < 0) ? tmp2 + p : tmp2;
1267: }
1268:
1269: /* determine 1/(4A) mod kN */
1270: av = avma;
1271: if (!invmod(shifti(A,2), kN, &p1)) /* --GN */
1272: {
1273: /* inv_A4 is a factor > 1 of kN. This cannot actually happen.
1274: Catch it anyway, it doesn't cost more than two comparisons
1275: against 0 / NULL */
1276: p1 = mppgcd(p1, N);
1277: *f = gerepileupto(av, p1);
1278: return; /* we certainly can't use the current poly */
1279: }
1280: affii(p1, inv_A4);
1281: avma = av;
1282: }
1283: else
1284: {
1285: /* no "real" computation -- use recursive formula
1286: * first: update of B, compute B[index_i], index_i > 0
1287: */
1288:
1289: /* compute nu_2 = nu_2(index_i) */
1290: nu_2 = 0;
1291: j = index_i;
1292: while ((j & 1) == 0)
1293: {
1294: nu_2++;
1295: j >>= 1;
1296: }
1297:
1298: av = avma;
1299: p1 = shifti(H_i[nu_2], 1);
1300:
1301: if ((((j+1)/2) & 1) == 1)
1302: {
1303: i = -1;
1304: p1 = subii(B, p1);
1305: }
1306: else
1307: {
1308: i = 1;
1309: p1 = addii(B, p1);
1310: }
1311: affii(p1, B);
1312: avma = av;
1313:
1314: /* determine new starting positions for sieving */
1315: size_of_FB = FB[0] + 1;
1316:
1317: if (i == -1)
1318: {
1319: for (j = 2; (ulong)j <= size_of_FB; j++)
1320: {
1321: p = FB[j];
1322: start_1[j] += inv_A_H_i[nu_2][j];
1323: if (start_1[j] >= p)
1324: start_1[j] -= p;
1325: start_2[j] += inv_A_H_i[nu_2][j];
1326: if (start_2[j] >= p)
1327: start_2[j] -= p;
1328: }
1329: }
1330: else
1331: {
1332: for (j = 2; (ulong)j <= size_of_FB; j++)
1333: {
1334: p = FB[j];
1335: start_1[j] -= inv_A_H_i[nu_2][j];
1336: if (start_1[j] < 0)
1337: start_1[j] += p;
1338: start_2[j] -= inv_A_H_i[nu_2][j];
1339: if (start_2[j] < 0)
1340: start_2[j] += p;
1341: }
1342: }
1343: }
1344:
1345: /* p=2 is a special case */
1346: if (FB[2] == 2)
1347: {
1348: start_1[2] = 0; start_2[2] = 1;
1349: }
1350:
1351:
1352: /* now compute zeros of polynomials that have only one zero mod p
1353: because p divides the coefficient A */
1354:
1355: /* compute coefficient -C */
1356: av = avma;
1357: p1 = subii(kN, sqri(B));
1358: p2 = divii(p1, shifti(A, 2));
1359:
1360: for (j = 1; (ulong)j <= total_no_of_primes_for_A; j++)
1361: if (*bin_index & (1 << (j-1)))
1362: {
1363: p = FB[start_index_FB_for_A + j];
1364: tmp = mpqs_invsmod(smodis(B, p), p);
1365: tmp2 = smodis (p2, p);
1366: tmp = mulssmod(tmp2, tmp, p);
1367: start_1[start_index_FB_for_A + j] =
1368: start_2[start_index_FB_for_A + j] = (tmp + M) % p;
1369: }
1370: avma = av;
1371:
1372: #ifdef MPQS_DEBUG
1373: {
1374: /*
1375: * check correctness of the roots z1, z2 mod p_i by evaluting
1376: * A * z1^2 + B * z1 + C mod p_i and checking against 0 mod p_i
1377: * C is written as (B*B-kN)/(4*A)
1378: */
1379: size_of_FB = FB[0] + 1;
1380: for (j = 3; j <= size_of_FB; j++)
1381: {
1382: av = avma;
1383: p = FB[j];
1384: M_mod_p = M % p;
1385:
1386: p1 = mulis(A, start_1[j] - M_mod_p);
1387: p1 = addii(p1, B);
1388: p1 = mulis(p1, start_1[j] - M_mod_p);
1389: p2 = divii(subii(sqri(B), kN), shifti(A, 2));
1390: p1 = addii(p1, p2);
1391: if (smodis(p1, p) != 0)
1392: {
1393: fprintferr("MPQS: (1) j = %ld p = %ld\n", j, p);
1394: fprintferr("MPQS: (1) index_i = %ld\n", index_i);
1395: fprintferr("MPQS: A = %Z\n", A);
1396: fprintferr("MPQS: B = %Z\n", B);
1397: fprintferr("MPQS: C = %Z\n", p1);
1398: fprintferr("MPQS: z1 = %ld\n", start_1[j] - M_mod_p);
1399: err(talker, "MPQS: self_init: found wrong polynomial in (1)");
1400: }
1401:
1402: p1 = mulis(A, start_2[j]-M_mod_p);
1403: p1 = addii(p1, B);
1404: p1 = mulis(p1, start_2[j]-M_mod_p);
1405: p2 = divii(subii(sqri(B), kN), shifti(A, 2));
1406: p1 = addii(p1, p2);
1407: if (smodis(p1, p) != 0)
1408: {
1409: fprintferr("MPQS: (2) j = %ld p = %ld\n", j, p);
1410: fprintferr("MPQS: (2) index_i = %ld\n", index_i);
1411: fprintferr("MPQS: A = %Z\n", A);
1412: fprintferr("MPQS: B = %Z\n", B);
1413: fprintferr("MPQS: C = %Z\n", p1);
1414: fprintferr("MPQS: z2 = %ld\n", start_2[j] - M_mod_p);
1415: err(talker, "MPQS: self_init: found wrong polynomial in (2)");
1416: }
1417: avma = av;
1418: }
1419: if (DEBUGLEVEL >= 4)
1420: fprintferr("MPQS: checking of roots of Q(x) was successful\n");
1421: }
1422: #endif
1423:
1424: }
1425:
1426: /******************************/
1427:
1428: /**
1429: ** Main sieving routine:
1430: **
1431: ** FB is a pointer to an array which holds the factor base
1432: ** log_FB is a pointer to an array which holds the approximations for
1433: ** the logarithms of the factor base elements
1434: ** start_1, start_2
1435: ** are arrays for starting points for different FB elements
1436: ** sieve_array
1437: ** points to a sieve array
1438: ** sieve_array_end
1439: ** points to the end of the sieve array
1440: ** M is the size of the sieving interval
1441: ** starting_sieving_index
1442: ** marks the first FB element which is used for sieving
1443: **
1444: **/
1445:
1446: #ifdef INLINE
1447: INLINE
1448: #endif
1449: void
1450: mpqs_sieve_p(unsigned char *begin, unsigned char *end,
1451: long p_times_4, long p, unsigned char log_p)
1452: {
1453: register unsigned char *e = end - p_times_4;
1454: while (e - begin >= 0) /* signed comparison --GN2000Jun02 */
1455: {
1456: (*begin) += log_p, begin += p;
1457: (*begin) += log_p, begin += p;
1458: (*begin) += log_p, begin += p;
1459: (*begin) += log_p, begin += p;
1460: }
1461: while (end - begin >= 0) /* again --GN2000Jun02 */
1462: (*begin) += log_p, begin += p;
1463: }
1464:
1465: static void
1466: mpqs_sieve(long *FB, unsigned char *log_FB, long *start_1, long *start_2,
1467: unsigned char *sieve_array, unsigned char *sieve_array_end,
1468: ulong M, ulong starting_sieving_index)
1469: {
1470: register long p_times_4, p, *ptr_FB, start1, start2, l;
1471: register unsigned char log_p;
1472:
1473: memset (sieve_array, 0, (M << 1) * sizeof (unsigned char));
1474:
1475: l = starting_sieving_index;
1476: ptr_FB = &FB[l];
1477:
1478: while ((p = *ptr_FB++) != 0)
1479: {
1480: p_times_4 = p * 4;
1481: log_p = log_FB[l];
1482: start1 = start_1[l];
1483: start2 = start_2[l];
1484:
1485: /* sieve with FB[l] from start_1[l] */
1486: /* if start1 != start2 sieve with FB[l] from start_2[l] */
1487: if (start1 == start2)
1488: mpqs_sieve_p(sieve_array + start1, sieve_array_end,
1489: p_times_4, p, log_p);
1490: else
1491: {
1492: mpqs_sieve_p(sieve_array + start1, sieve_array_end,
1493: p_times_4, p, log_p);
1494: mpqs_sieve_p(sieve_array + start2, sieve_array_end,
1495: p_times_4, p, log_p);
1496: }
1497:
1498: l++;
1499: }
1500: }
1501:
1502: /******************************/
1503:
1504: static long
1505: mpqs_eval_sieve(unsigned char *sieve_array, ulong M, long *candidates)
1506: {
1507: ulong M_2 = M << 1;
1508: int i = 0, count, q, r;
1509: unsigned char *ptr = sieve_array;
1510:
1511: r = M_2 % 4;
1512: q = M_2 - r;
1513: count = 0;
1514: for (i = 0; i < q; i += 4, ptr += 4)
1515: {
1516: if (*(ptr+0) > 128) { candidates[count] = i+0; count++; }
1517: if (*(ptr+1) > 128) { candidates[count] = i+1; count++; }
1518: if (*(ptr+2) > 128) { candidates[count] = i+2; count++; }
1519: if (*(ptr+3) > 128) { candidates[count] = i+3; count++; }
1520: }
1521: switch (r)
1522: {
1523: case 3:
1524: if (*(ptr+0) > 128) { candidates[count] = i+0; count++; }
1525: if (*(ptr+1) > 128) { candidates[count] = i+1; count++; }
1526: if (*(ptr+2) > 128) { candidates[count] = i+2; count++; }
1527: break;
1528: case 2:
1529: if (*(ptr+0) > 128) { candidates[count] = i+0; count++; }
1530: if (*(ptr+1) > 128) { candidates[count] = i+1; count++; }
1531: break;
1532: case 1:
1533: if (*(ptr+0) > 128) { candidates[count] = i+0; count++; }
1534: break;
1535: default:
1536: break;
1537: }
1538: candidates[count] = 0;
1539: return count;
1540: }
1541:
1542: /**
1543: ** Main relation routine:
1544: **
1545: ** FB is a pointer to an array which holds the factor base
1546: ** log_FB is a pointer to an array which holds the approximations for
1547: ** the logarithms of the factor base elements
1548: ** start_1, start_2
1549: ** are arrays for starting points for different FB elements
1550: ** M is the size of the sieving interval
1551: **
1552: **/
1553:
1554: #ifdef INLINE
1555: INLINE void
1556: mpqs_add_factor(char *s, ulong ei, ulong pi)
1557: {
1558: sprintf(s + strlen(s), " %lu %lu", ei, pi);
1559: }
1560: #else
1561: #define mpqs_add_factor(s,ei,pi) \
1562: sprintf(s + strlen(s), " %lu %lu", (ulong)ei, (ulong)pi)
1563: #endif
1564:
1565: /**
1566: ** divide and return the remainder. Leaves both the quotient and
1567: ** the remainder on the stack as GENs; the caller should clean up.
1568: **/
1569:
1570: #ifdef INLINE
1571: INLINE
1572: #else
1573: static
1574: #endif
1575: ulong
1576: mpqs_div_rem(GEN x, long y, GEN *q)
1577: {
1578: GEN r;
1579: *q = dvmdis(x, y, &r);
1580: if (signe(r)) return (ulong)(r[2]);
1581: return 0;
1582: }
1583:
1584: static long
1585: mpqs_eval_candidates(GEN A, GEN inv_A4, GEN B, GEN kN, long k,
1586: double sqrt_kN,
1587: long *FB, long *start_1, long *start_2,
1588: ulong M, long bin_index, long *candidates,
1589: long number_of_candidates, long lp_bound,
1590: long start_index_FB_for_A,
1591: FILE *FREL, FILE *LPREL)
1592: /* NB FREL, LPREL are actually FNEW, LPNEW when we get called */
1593: {
1594: GEN Qx, A_2x_plus_B, Y, Qx_div_p;
1595: double a, b, inv_2_a;
1596: long av, powers_of_2, p, tmp_p, remd_p, bi;
1597: long z1, z2, x, number_of_relations, x_minus_M;
1598: char *relations;
1599: ulong i, pi, ei, size_of_FB;
1600: int useless_cand;
1601: #ifdef MPQS_DEBUG_VERBOSE
1602: static char complaint[256], complaint0[256];
1603: #endif
1604:
1605: /* FIXME: this is probably right !!! */
1606:
1607: /*
1608: * compute the roots of Q(X); this will be used to find the sign of
1609: * Q(x) after reducing mod kN
1610: */
1611:
1612: a = gtodouble(A);
1613: inv_2_a = 1 / (2.0 * a);
1614: b = gtodouble(B);
1615: z1 = (long) ((-b - sqrt_kN) * inv_2_a);
1616: z2 = (long) ((-b + sqrt_kN) * inv_2_a);
1617:
1618: number_of_relations = 0;
1619: /* FIXME: this is definitely too loooooong ... */
1620: /* then take the size of the FB into account.
1621: We have in the worst case: a leading " 1 1", a trailing " 0\n" with
1622: the final NUL character, and in between up to size_of_FB pairs each
1623: consisting of an exponent, a subscript into FB, and two spaces.
1624: Moreover, subscripts into FB fit into 5 digits, and exponents fit
1625: into 3 digits with room to spare -- anything needing 3 or more digits
1626: for the subscript must come with an exponent of at most 2 digits.
1627: Moreover the product of the first 58 primes is larger than 10^110,
1628: so there cannot be more than 60 pairs in all, even if size_of_FB > 10^4.
1629: --GN */
1630: size_of_FB = FB[0]; /* one less than usually: don't count FB[1] */
1631: if (size_of_FB > 60) size_of_FB = 60;
1632: relations =
1633: (char *) gpmalloc((8 + size_of_FB * 9) * sizeof(char));
1634: #ifdef MPQS_DEBUG_VERBOSE
1635: complaint[0] = '\0';
1636: #endif
1637: #ifdef MPQS_DEBUG
1638: size_of_FB = FB[0] + 1; /* last useful subscript into FB */
1639: #endif
1640:
1641: i = 0;
1642: while (i < (ulong)number_of_candidates)
1643: {
1644: #ifdef MPQS_DEBUG_VERYVERBOSE /* this one really fills the screen... */
1645: fprintferr("%c", (char)('0' + i%10));
1646: #endif
1647: x = candidates[i++];
1648: relations[0] = '\0';
1649: av = avma;
1650:
1651: /* A_2x_plus_B = (A*(2x)+B), Qx = (A*(2x)+B)^2/(4*A) = Q(x) */
1652: x_minus_M = x - M;
1653: A_2x_plus_B = modii(addii(mulis(A, x_minus_M << 1), B), kN);
1654: Y = subii(kN, A_2x_plus_B);
1655: if (absi_cmp(A_2x_plus_B, Y) < 0) Y = A_2x_plus_B;
1656: /* absolute value of smallest absolute residue of A_2x_plus_B mod kN */
1657:
1658: Qx = modii(sqri(Y), kN);
1659:
1660: /* Most of the time, gcd(Qx, kN) will be either 1 or k. However,
1661: it may happen that Qx is a multiple of N, especially when N
1662: is small, and this will lead to havoc below -- so we have to be
1663: a little bit careful. Of course we cannot possibly afford to
1664: compute the gcd each time through this loop unless we are
1665: debugging... --GN */
1666: #ifdef MPQS_DEBUG
1667: {
1668: long av1 = avma, ks;
1669: GEN g = mppgcd(Qx, kN);
1670: /* if ((ks = kronecker(divii(Qx, g), divii(kN, g))) != 1) */
1671: if (is_pm1(g))
1672: {
1673: if ((ks = kronecker(Qx, kN)) != 1)
1674: {
1675: fprintferr("\nMPQS: 4*A*Q(x) = %Z\n", Qx);
1676: fprintferr("\tKronecker symbol %ld\n", ks);
1677: err(talker, "MPQS: 4*A*Q(x) is not a square (mod kN)");
1678: }
1679: }
1680: #ifdef MPQS_DEBUG_VERBOSE
1681: else if (cmpis(g,k) /* != 0 */ )
1682: {
1683: char *gs = GENtostr(g);
1684: sprintf(complaint, "\nMPQS: gcd(4*A*Q(x), kN) = %s\n", gs);
1685: free(gs);
1686: if (strcmp(complaint, complaint0) != 0)
1687: {
1688: fprintferr(complaint);
1689: strcpy(complaint0, complaint);
1690: }
1691: }
1692: #endif
1693: avma = av1;
1694: }
1695: #endif
1696:
1697: Qx = modii(mulii(Qx, inv_A4), kN);
1698:
1699: /* check the sign of Qx */
1700: if (z1 < x_minus_M && x_minus_M < z2)
1701: {
1702: Qx = subii(kN, Qx);
1703: /* i = 1, ei = 1, pi */
1704: mpqs_add_factor(relations, 1, 1);
1705: }
1706: if (!signe(Qx))
1707: {
1708: #ifdef MPQS_DEBUG_VERBOSE
1709: fprintferr("<+>");
1710: #endif
1711: avma = av;
1712: continue;
1713: }
1714:
1715: /* divide by powers of 2; note that we're really dealing with
1716: 4*A*Q(x), so we remember an extra factor 2^2 */
1717: powers_of_2 = vali(Qx);
1718: Qx = shifti(Qx, -powers_of_2);
1719: mpqs_add_factor(relations, powers_of_2 + 2, 2);
1720: if (!signe(Qx)) /* this shouldn't happen */
1721: {
1722: #ifdef MPQS_DEBUG_VERBOSE
1723: fprintferr("<*>");
1724: #endif
1725: avma = av;
1726: continue;
1727: }
1728:
1729: /* we handled the case p = 2 already */
1730: pi = 3;
1731: bi = bin_index;
1732: useless_cand = 0;
1733: #ifdef MPQS_DEBUG_VERBOSE
1734: fprintferr("a");
1735: #endif
1736:
1737: /* FB[3] .. FB[start_index_FB_for_A] do not divide A .
1738: * p = FB[start_index_FB_for_A+j+1] divides A (to the first power)
1739: * iff the 2^j bit in bin_index is set */
1740: while ((p = FB[pi]) != 0 && !useless_cand)
1741: {
1742: tmp_p = x % p;
1743:
1744: ei = 0;
1745: if (bi && pi > (ulong)start_index_FB_for_A)
1746: {
1747: ei = bi & 1; /* either 0 or 1 */
1748: bi >>= 1;
1749: }
1750:
1751: if (tmp_p == start_1[pi] || tmp_p == start_2[pi])
1752: {
1753: /* p divides Q(x) and possibly A */
1754: remd_p = mpqs_div_rem(Qx, p, &Qx_div_p);
1755: if (remd_p)
1756: {
1757: useless_cand = 1;
1758: break;
1759: }
1760: do
1761: {
1762: ei++;
1763: Qx = Qx_div_p;
1764: remd_p = mpqs_div_rem(Qx, p, &Qx_div_p);
1765: } while (remd_p == 0);
1766: }
1767: if (ei) /* p might divide A but not Q(x) */
1768: mpqs_add_factor(relations, ei, pi);
1769:
1770: pi++;
1771: }
1772:
1773: if (useless_cand)
1774: {
1775: #ifdef MPQS_DEBUG_VERBOSE
1776: fprintferr("\b<");
1777: #endif
1778: avma = av;
1779: continue;
1780: }
1781: #ifdef MPQS_DEBUG_VERBOSE
1782: fprintferr("\bb");
1783: #endif
1784:
1785: if (is_pm1(Qx))
1786: {
1787: char *Qxstring = GENtostr(Y);
1788: strcat(relations, " 0");
1789: fprintf(FREL, "%s :%s\n", Qxstring, relations);
1790: number_of_relations++;
1791:
1792: #ifdef MPQS_DEBUG
1793: {
1794: GEN Qx_2, prod_pi_ei, pi_ei;
1795: long lr = strlen(relations);
1796: char *s, *t = gpmalloc(lr+1);
1797: long pi, ei, av1 = avma;
1798:
1799: #ifdef MPQS_DEBUG_VERBOSE
1800: fprintferr("\b(");
1801: #endif
1802: Qx_2 = modii(sqri(Y), kN);
1803: prod_pi_ei = gun;
1804: strcpy(t, relations);
1805: s = strtok(relations, " \n");
1806: while (s != NULL)
1807: {
1808: ei = atol(s);
1809: if (ei == 0) break;
1810: s = strtok(NULL, " \n");
1811: pi = atol(s);
1812: pi_ei = powmodulo(stoi(FB[pi]), stoi(ei), kN);
1813: prod_pi_ei = modii(mulii(prod_pi_ei, pi_ei), kN);
1814: s = strtok(NULL, " \n");
1815: }
1816:
1817: if (!egalii(Qx_2, prod_pi_ei))
1818: {
1819: #ifdef MPQS_DEBUG_VERBOSE
1820: fprintferr("!)\n");
1821: #endif
1822: fprintferr("MPQS: %s :%s\n", Qxstring, t);
1823: fprintferr("\tQx_2 = %Z\n", Qx_2);
1824: fprintferr("\t rhs = %Z\n", prod_pi_ei);
1825: err(talker, "MPQS: wrong full relation found!!");
1826: }
1827: #ifdef MPQS_DEBUG_VERBOSE
1828: else
1829: fprintferr(":)");
1830: #endif
1831: avma = av1;
1832: free(t);
1833: }
1834: #endif
1835:
1836: free(Qxstring);
1837: }
1838: else if (cmpis(Qx, lp_bound) > 0)
1839: {
1840: /* TO BE DONE: check for double large prime */
1841: /* moved this else clause here where I think it belongs -- GN */
1842: #ifdef MPQS_DEBUG_VERBOSE
1843: fprintferr("\b.");
1844: #endif
1845: }
1846: else if (k==1 || cgcd(k, itos(Qx)) == 1)
1847: {
1848: /* if (mpqs_isprime(itos(Qx))) */
1849: char *Qxstring = GENtostr(Y);
1850: char *L1string = GENtostr(Qx);
1851: strcat(relations, " 0");
1852: fprintf(LPREL, "%s @ %s :%s\n", L1string, Qxstring, relations);
1853:
1854: #ifdef MPQS_DEBUG
1855: {
1856: GEN Qx_2, prod_pi_ei, pi_ei;
1857: long lr = strlen(relations);
1858: char *s, *t = gpmalloc(lr+1);
1859: long pi, ei, av1 = avma;
1860:
1861: #ifdef MPQS_DEBUG_VERBOSE
1862: fprintferr("\b(");
1863: #endif
1864: Qx_2 = modii(sqri(Y), kN);
1865: prod_pi_ei = gun;
1866: strcpy(t, relations);
1867: s = strtok(relations, " \n");
1868: while (s != NULL)
1869: {
1870: ei = atol(s);
1871: if (ei == 0) break;
1872: s = strtok(NULL, " \n");
1873: pi = atol(s);
1874: pi_ei = powmodulo(stoi(FB[pi]), stoi(ei), kN);
1875: prod_pi_ei = modii(mulii(prod_pi_ei, pi_ei), kN);
1876: s = strtok(NULL, " \n");
1877: }
1878: prod_pi_ei = modii(mulii(prod_pi_ei, Qx), kN);
1879:
1880: if (!egalii(Qx_2, prod_pi_ei))
1881: {
1882: #ifdef MPQS_DEBUG_VERBOSE
1883: fprintferr("!)\n");
1884: #endif
1885: fprintferr("MPQS: %s @ %s :%s\n",
1886: L1string, Qxstring, t);
1887: fprintferr("\tQx_2 = %Z\n", Qx_2);
1888: fprintferr("\t rhs = %Z\n", prod_pi_ei);
1889: err(talker, "MPQS: wrong large prime relation found!!");
1890: }
1891: #ifdef MPQS_DEBUG_VERBOSE
1892: else
1893: fprintferr(";)");
1894: #endif
1895: avma = av1;
1896: free(t);
1897: }
1898: #endif
1899:
1900: free(Qxstring);
1901: free(L1string);
1902: }
1903: #ifdef MPQS_DEBUG_VERBOSE
1904: else
1905: fprintferr("\b<k>");
1906: #endif
1907:
1908: avma = av;
1909: } /* while */
1910:
1911: #ifdef MPQS_DEBUG_VERBOSE
1912: fprintferr("\n");
1913: #endif
1914:
1915: free(relations);
1916: return number_of_relations;
1917: }
1918:
1919: /******************************/
1920:
1921: /**
1922: ** combines the large prime relations in COMB to full
1923: ** relations in FNEW. FNEW is assumed to be open for
1924: ** writing / appending.
1925: **/
1926:
1927: typedef struct
1928: {
1929: long q;
1930: char Y [MPQS_STRING_LENGTH];
1931: char ep [MPQS_STRING_LENGTH];
1932: } mpqs_lp_entry;
1933:
1934: static void
1935: mpqs_combine_exponents(long *ei, long ei_size, char *r1, char *r2)
1936: {
1937: char ej [MPQS_STRING_LENGTH], ek [MPQS_STRING_LENGTH];
1938: char *s;
1939: long e, p;
1940:
1941: memset(ei, 0, ei_size * sizeof(long));
1942: strcpy(ej, r1);
1943: strcpy(ek, r2);
1944:
1945: s = strtok(ej, " \n");
1946: while (s != NULL)
1947: {
1948: e = atol(s);
1949: if (e == 0)
1950: break;
1951: s = strtok(NULL, " \n");
1952: p = atol(s);
1953: ei[p] += e;
1954: s = strtok(NULL, " \n");
1955: }
1956:
1957: s = strtok(ek, " \n");
1958: while (s != NULL)
1959: {
1960: e = atol(s);
1961: if (e == 0)
1962: break;
1963: s = strtok(NULL, " \n");
1964: p = atol(s);
1965: ei[p] += e;
1966: s = strtok(NULL, " \n");
1967: }
1968:
1969: }
1970:
1971: static long
1972: mpqs_combine_large_primes(FILE *COMB, FILE *FNEW, long size_of_FB,
1973: GEN N, GEN kN, GEN *f
1974: #ifdef MPQS_DEBUG
1975: , long *FB
1976: #endif
1977: )
1978: {
1979: char buf [MPQS_STRING_LENGTH], *s1, *s2;
1980: char new_relation [MPQS_STRING_LENGTH];
1981: mpqs_lp_entry e[2]; /* we'll use the two alternatingly */
1982: long *ei, ei_size = size_of_FB + 2;
1983: long old_q;
1984: GEN inv_q, Y1, Y2, new_Y, new_Y1;
1985: long i, l, c = 0;
1986: long av = avma, av2;
1987:
1988: *f = NULL;
1989: if (fgets(buf, MPQS_STRING_LENGTH, COMB) == NULL)
1990: return 0; /* nothing there -- should not happen */
1991:
1992: /* put first lp relation in row 0 of e */
1993: s1 = buf; s2 = strchr(s1, ' '); *s2 = '\0';
1994: e[0].q = atol(s1);
1995: s1 = s2 + 3; s2 = strchr(s1, ' '); *s2 = '\0';
1996: strcpy(e[0].Y, s1);
1997: s1 = s2 + 3; s2 = strchr(s1, '\n'); *s2 = '\0';
1998: strcpy(e[0].ep, s1);
1999:
2000: i = 1; /* second relation will go into row 1 */
2001: old_q = e[0].q;
2002:
2003: while (!invmod(stoi(old_q), kN, &inv_q)) /* can happen --GN */
2004: {
2005: inv_q = mppgcd(inv_q, N);
2006: if (is_pm1(inv_q) || egalii(inv_q, N)) /* pity */
2007: {
2008: #ifdef MPQS_DEBUG
2009: fprintferr("MPQS: skipping relation with non-invertible q\n");
2010: #endif
2011: avma = av;
2012: if (fgets(buf, MPQS_STRING_LENGTH, COMB) == NULL)
2013: return 0;
2014: s1 = buf; s2 = strchr(s1, ' '); *s2 = '\0';
2015: e[0].q = atol(s1);
2016: s1 = s2 + 3; s2 = strchr(s1, ' '); *s2 = '\0';
2017: strcpy(e[0].Y, s1);
2018: s1 = s2 + 3; s2 = strchr(s1, '\n'); *s2 = '\0';
2019: strcpy(e[0].ep, s1);
2020: old_q = e[0].q;
2021: continue;
2022: }
2023: *f = gerepileupto(av, inv_q);
2024: return c;
2025: }
2026: Y1 = lisexpr(e[0].Y);
2027: av2 = avma; /* preserve inv_q and Y1 */
2028:
2029: ei = (long *) gpmalloc(ei_size * sizeof(long));
2030:
2031: while (fgets(buf, MPQS_STRING_LENGTH, COMB) != NULL)
2032: {
2033: s1 = buf; s2 = strchr(s1, ' '); *s2 = '\0';
2034: e[i].q = atol(s1);
2035: s1 = s2 + 3; s2 = strchr(s1, ' '); *s2 = '\0';
2036: strcpy(e[i].Y, s1);
2037: s1 = s2 + 3; s2 = strchr(s1, '\n'); *s2 = '\0';
2038: strcpy(e[i].ep, s1);
2039:
2040: if (e[i].q != old_q)
2041: {
2042: /* switch to combining a new bunch, swapping the rows */
2043: old_q = e[i].q;
2044: avma = av; /* discard old inv_q and Y1 */
2045: if (!invmod(stoi(old_q), kN, &inv_q)) /* can happen --GN */
2046: {
2047: inv_q = mppgcd(inv_q, N);
2048: if (is_pm1(inv_q) || egalii(inv_q, N)) /* pity */
2049: {
2050: #ifdef MPQS_DEBUG
2051: fprintferr("MPQS: skipping relation with non-invertible q\n");
2052: #endif
2053: avma = av;
2054: old_q = -1; /* sentinel */
2055: av2 = avma;
2056: continue; /* discard this combination */
2057: }
2058: *f = gerepileupto(av, inv_q);
2059: free (ei);
2060: return c;
2061: }
2062: Y1 = lisexpr(e[i].Y);
2063: i = 1 - i; /* subsequent relations go to other row */
2064: av2 = avma; /* preserve inv_q and Y1 */
2065: continue;
2066: }
2067: /* count and combine the two we've got, and continue in the same row */
2068: c++;
2069:
2070: mpqs_combine_exponents(ei, ei_size, e[1-i].ep, e[i].ep);
2071: if (DEBUGLEVEL >= 6)
2072: {
2073: fprintferr("MPQS: combining\n");
2074: fprintferr(" {%ld @ %s : %s}\n", old_q, e[1-i].Y, e[1-i].ep);
2075: fprintferr(" * {%ld @ %s : %s}\n", e[i].q, e[i].Y, e[i].ep);
2076: }
2077: Y2 = lisexpr(e[i].Y);
2078: new_Y = modii(mulii(mulii(Y1, Y2), inv_q), kN);
2079: new_Y1 = subii(kN, new_Y);
2080: if (absi_cmp(new_Y1, new_Y) < 0) new_Y = new_Y1;
2081: s1 = GENtostr(new_Y);
2082: strcpy(new_relation, s1);
2083: free(s1);
2084: strcat(new_relation, " :");
2085: if (ei[1] % 2)
2086: strcat(new_relation, " 1 1");
2087: for (l = 2; l < ei_size; l++)
2088: {
2089: if (ei[l])
2090: {
2091: sprintf(buf, " %ld %ld", ei[l], l);
2092: strcat(new_relation, buf);
2093: }
2094: }
2095: strcat(new_relation, " 0");
2096: if (DEBUGLEVEL >= 6) fprintferr(" == {%s}\n", new_relation);
2097: strcat(new_relation, "\n");
2098:
2099: #ifdef MPQS_DEBUG
2100: {
2101: char ejk [MPQS_STRING_LENGTH];
2102: GEN Qx_2, prod_pi_ei, pi_ei;
2103: char *s;
2104: long pi, exi, av1 = avma;
2105: Qx_2 = modii(sqri(new_Y), kN);
2106:
2107: strcpy(ejk, new_relation);
2108: s = strchr(ejk, ':') + 2;
2109: s = strtok(s, " \n");
2110:
2111: prod_pi_ei = gun;
2112: while (s != NULL)
2113: {
2114: exi = atol(s);
2115: if (exi == 0) break;
2116: s = strtok(NULL, " \n");
2117: pi = atol(s);
2118: pi_ei = powmodulo(stoi(FB[pi]), stoi(exi), kN);
2119: prod_pi_ei = modii(mulii(prod_pi_ei, pi_ei), kN);
2120: s = strtok(NULL, " \n");
2121: }
2122: avma = av1;
2123:
2124: if (!egalii(Qx_2, prod_pi_ei))
2125: {
2126: free(ei);
2127: err(talker,
2128: "MPQS: combined large prime relation is false");
2129: }
2130: }
2131: #endif
2132:
2133: if (fputs(new_relation, FNEW) < 0)
2134: {
2135: free(ei);
2136: bummer(FNEW_str);
2137: }
2138: avma = av2;
2139: } /* while */
2140:
2141: free(ei);
2142:
2143: if (DEBUGLEVEL >= 4)
2144: fprintferr("MPQS: combined %ld full relation%s\n",
2145: c, (c!=1 ? "s" : ""));
2146: return c;
2147: }
2148:
2149: /******************************/
2150:
2151: #ifdef LONG_IS_64BIT
2152:
2153: #define MPQS_GAUSS_BITS 64
2154: static unsigned long mpqs_mask_bit[] =
2155: {
2156: 0x8000000000000000UL, 0x4000000000000000UL,
2157: 0x2000000000000000UL, 0x1000000000000000UL,
2158: 0x0800000000000000UL, 0x0400000000000000UL,
2159: 0x0200000000000000UL, 0x0100000000000000UL,
2160: 0x0080000000000000UL, 0x0040000000000000UL,
2161: 0x0020000000000000UL, 0x0010000000000000UL,
2162: 0x0008000000000000UL, 0x0004000000000000UL,
2163: 0x0002000000000000UL, 0x0001000000000000UL,
2164: 0x0000800000000000UL, 0x0000400000000000UL,
2165: 0x0000200000000000UL, 0x0000100000000000UL,
2166: 0x0000080000000000UL, 0x0000040000000000UL,
2167: 0x0000020000000000UL, 0x0000010000000000UL,
2168: 0x0000008000000000UL, 0x0000004000000000UL,
2169: 0x0000002000000000UL, 0x0000001000000000UL,
2170: 0x0000000800000000UL, 0x0000000400000000UL,
2171: 0x0000000200000000UL, 0x0000000100000000UL,
2172: 0x0000000080000000UL, 0x0000000040000000UL,
2173: 0x0000000020000000UL, 0x0000000010000000UL,
2174: 0x0000000008000000UL, 0x0000000004000000UL,
2175: 0x0000000002000000UL, 0x0000000001000000UL,
2176: 0x0000000000800000UL, 0x0000000000400000UL,
2177: 0x0000000000200000UL, 0x0000000000100000UL,
2178: 0x0000000000080000UL, 0x0000000000040000UL,
2179: 0x0000000000020000UL, 0x0000000000010000UL,
2180: 0x0000000000008000UL, 0x0000000000004000UL,
2181: 0x0000000000002000UL, 0x0000000000001000UL,
2182: 0x0000000000000800UL, 0x0000000000000400UL,
2183: 0x0000000000000200UL, 0x0000000000000100UL,
2184: 0x0000000000000080UL, 0x0000000000000040UL,
2185: 0x0000000000000020UL, 0x0000000000000010UL,
2186: 0x0000000000000008UL, 0x0000000000000004UL,
2187: 0x0000000000000002UL, 0x0000000000000001UL
2188: };
2189:
2190: #else
2191:
2192: #define MPQS_GAUSS_BITS 32
2193: static unsigned long mpqs_mask_bit[] =
2194: {
2195: 0x80000000UL, 0x40000000UL, 0x20000000UL, 0x10000000UL,
2196: 0x08000000UL, 0x04000000UL, 0x02000000UL, 0x01000000UL,
2197: 0x00800000UL, 0x00400000UL, 0x00200000UL, 0x00100000UL,
2198: 0x00080000UL, 0x00040000UL, 0x00020000UL, 0x00010000UL,
2199: 0x00008000UL, 0x00004000UL, 0x00002000UL, 0x00001000UL,
2200: 0x00000800UL, 0x00000400UL, 0x00000200UL, 0x00000100UL,
2201: 0x00000080UL, 0x00000040UL, 0x00000020UL, 0x00000010UL,
2202: 0x00000008UL, 0x00000004UL, 0x00000002UL, 0x00000001UL
2203: };
2204:
2205: #endif
2206:
2207: static mpqs_gauss_matrix
2208: mpqs_gauss_create_matrix(long rows, long cols)
2209: {
2210: mpqs_gauss_matrix m;
2211: long i, j, words = cols / MPQS_GAUSS_BITS;
2212: if (cols % MPQS_GAUSS_BITS)
2213: words++;
2214: m = (mpqs_gauss_row *) gpmalloc(rows * sizeof(mpqs_gauss_row));
2215: for (i = 0; i < rows; i++)
2216: {
2217: m[i] = (ulong *) gpmalloc(words * sizeof(ulong));
2218: for (j = 0; j < words; j++)
2219: m[i][j] = 0UL;
2220: }
2221: return m;
2222: }
2223:
2224: static void
2225: mpqs_gauss_destroy_matrix(mpqs_gauss_matrix m, long rows, long cols)
2226: {
2227: long i;
2228: (void) cols;
2229: for (i = 0; i < rows; i++)
2230: free(m[i]);
2231: free(m);
2232: }
2233:
2234: static ulong
2235: mpqs_gauss_get_bit(mpqs_gauss_matrix m, long i, long j)
2236: {
2237: return m[i][j / MPQS_GAUSS_BITS] & mpqs_mask_bit[j % MPQS_GAUSS_BITS];
2238: }
2239:
2240: static void
2241: mpqs_gauss_set_bit(mpqs_gauss_matrix m, long i, long j)
2242: {
2243: m[i][j / MPQS_GAUSS_BITS] |= mpqs_mask_bit[j % MPQS_GAUSS_BITS];
2244: }
2245:
2246: static void
2247: mpqs_gauss_clear_bit(mpqs_gauss_matrix m, long i, long j)
2248: {
2249: m[i][j / MPQS_GAUSS_BITS] &= ~mpqs_mask_bit[j % MPQS_GAUSS_BITS];
2250: }
2251:
2252: /* output an mpqs_gauss_matrix in PARI format */
2253: /* TODO -- this should do a little intelligent wrapping... GN */
2254: static void
2255: mpqs_gauss_print_matrix(mpqs_gauss_matrix m, long rows, long cols)
2256: {
2257: long i, j;
2258: fprintferr("\n[");
2259: for (i = 0; i < rows; i++)
2260: {
2261: for (j = 0; j < cols - 1; j++)
2262: {
2263: if (mpqs_gauss_get_bit(m, i, j))
2264: fprintferr("1, ");
2265: else
2266: fprintferr("0, ");
2267: }
2268:
2269: if (mpqs_gauss_get_bit(m, i, j))
2270: fprintferr("1");
2271: else
2272: fprintferr("0");
2273: if (i != rows - 1)
2274: fprintferr("; ");
2275: }
2276: fprintferr("]\n");
2277: }
2278:
2279: /* x ^= y : row addition over F_2 */
2280: static void
2281: mpqs_gauss_add_rows(mpqs_gauss_row y, mpqs_gauss_row x, long k, long n)
2282: {
2283: long i, q, r;
2284: n = n - k;
2285: r = n % 8; q = n - r + k; i = 0 + k;
2286: for (; i < q; i += 8)
2287: {
2288: x[0+i] ^= y[0+i]; x[1+i] ^= y[1+i];
2289: x[2+i] ^= y[2+i]; x[3+i] ^= y[3+i];
2290: x[4+i] ^= y[4+i]; x[5+i] ^= y[5+i];
2291: x[6+i] ^= y[6+i]; x[7+i] ^= y[7+i];
2292: }
2293: switch (r)
2294: {
2295: case 7: x[i] ^= y[i], i++;
2296: case 6: x[i] ^= y[i], i++;
2297: case 5: x[i] ^= y[i], i++;
2298: case 4: x[i] ^= y[i], i++;
2299: case 3: x[i] ^= y[i], i++;
2300: case 2: x[i] ^= y[i], i++;
2301: case 1: x[i] ^= y[i], i++;
2302: default: ;
2303: }
2304: }
2305:
2306: /* create and read an mpqs_gauss_matrix from a relation file FREL (opened
2307: by caller). Also record the position of each relation in the file for
2308: later use
2309: rows = size_of_FB+1, cols = rel */
2310: static mpqs_gauss_matrix
2311: mpqs_gauss_read_matrix(FILE *FREL, long rows, long cols, long *fpos)
2312: {
2313: mpqs_gauss_matrix m;
2314: long i = 0, e, p;
2315: char buf[MPQS_STRING_LENGTH], *s;
2316:
2317: m = mpqs_gauss_create_matrix(rows, cols);
2318: if ((fpos[0] = ftell(FREL)) < 0)
2319: err(talker, "ftell error on full relations file");
2320: while (fgets(buf, MPQS_STRING_LENGTH, FREL) != NULL)
2321: {
2322: s = strchr(buf, ':') + 2;
2323: s = strtok(s, " \n");
2324: while (s != NULL)
2325: {
2326: e = atol(s);
2327: if (e == 0)
2328: break;
2329: s = strtok(NULL, " \n");
2330: p = atol(s);
2331: if (e & 1)
2332: mpqs_gauss_set_bit(m, p - 1, i);
2333: s = strtok(NULL, " \n");
2334: }
2335: i++;
2336: if (i < cols && (fpos[i] = ftell(FREL)) < 0)
2337: err(talker, "ftell error on full relations file");
2338: }
2339: if (i != cols)
2340: {
2341: fprintferr("MPQS: full relations file %s than expected",
2342: i > cols ? "longer" : "shorter");
2343: err(talker, "MPQS panicking");
2344: }
2345: return m;
2346: }
2347:
2348: /* compute the kernel of an mpqs_gauss_matrix over F_2
2349: uses the algorithm described by Knuth, Vol. 2 as modified by Henri Cohen
2350: and optimized by Christian Batut, Thomas Papanikolaou and Xavier Roblot.
2351: m = rows, n = cols */
2352: static mpqs_gauss_matrix
2353: mpqs_kernel(mpqs_gauss_matrix x, long m, long n, long *rank)
2354: {
2355: long k, i, j, t, r;
2356: mpqs_gauss_matrix ker_x;
2357: mpqs_gauss_indices c, d;
2358: long words = n / MPQS_GAUSS_BITS;
2359: if (n % MPQS_GAUSS_BITS)
2360: words++;
2361:
2362: c = (long *) gpmalloc(m * sizeof(long));
2363: for (k = 0; k < m; k++) c[k] = -1;
2364:
2365: d = (long *) gpmalloc(n * sizeof(long));
2366: r = 0;
2367:
2368: for (k = 0; k < n; k++)
2369: {
2370: j = 0;
2371: while (j < m && (c[j] >= 0 || mpqs_gauss_get_bit(x, j, k) == 0))
2372: j++;
2373:
2374: if (j == m)
2375: {
2376: d[k] = -1; /* no pivot found in column k */
2377: r++; /* column k is a vector of the kernel */
2378: }
2379: else
2380: {
2381: d[k] = j; /* the pivot in column k is at row j */
2382: c[j] = k; /* a pivot at position j was found in column k */
2383: for (t = 0; t < m; t++)
2384: {
2385: if (t != j)
2386: {
2387: if (mpqs_gauss_get_bit(x, t, k))
2388: mpqs_gauss_add_rows(x[j], x[t], k / MPQS_GAUSS_BITS, words);
2389: }
2390: }
2391: }
2392: }
2393:
2394: ker_x = mpqs_gauss_create_matrix(n, r);
2395:
2396: for (j = k = 0; j < r; j++, k++)
2397: {
2398: while (d[k] != -1) k++;
2399: for (i = 0; i < k; i++)
2400: {
2401: if (d[i] != -1)
2402: {
2403: if (mpqs_gauss_get_bit(x, d[i], k))
2404: mpqs_gauss_set_bit(ker_x, i, j);
2405: else
2406: mpqs_gauss_clear_bit(ker_x, i, j);
2407: }
2408: else
2409: mpqs_gauss_clear_bit(ker_x, i, j);
2410: }
2411: mpqs_gauss_set_bit(ker_x, k, j);
2412: }
2413:
2414: free(c);
2415: free(d);
2416:
2417: *rank = r;
2418: return ker_x;
2419: }
2420:
2421: /******************************/
2422:
2423: static GEN
2424: mpqs_add_relation(GEN Y_prod, GEN N_or_kN, long *ei, char *r1)
2425: {
2426: GEN Y, res;
2427: char relation [MPQS_STRING_LENGTH];
2428: char *s;
2429: long av = avma, e, p;
2430:
2431: strcpy(relation, r1);
2432: s = strchr(relation, ':') - 1;
2433: *s = '\0';
2434: Y = lisexpr(relation);
2435: res = modii(mulii(Y_prod, Y), N_or_kN);
2436:
2437: s = s + 3;
2438: s = strtok(s, " \n");
2439: while (s != NULL)
2440: {
2441: e = atol(s);
2442: if (e == 0)
2443: break;
2444: s = strtok(NULL, " \n");
2445: p = atol(s);
2446: ei[p] += e;
2447: s = strtok(NULL, " \n");
2448: }
2449: return gerepileupto(av,res);
2450: }
2451:
2452: static char*
2453: mpqs_get_relation(long pos, FILE *FREL)
2454: {
2455: static char buf [MPQS_STRING_LENGTH];
2456: if (fseek(FREL, pos, SEEK_SET))
2457: err(talker, "can\'t seek full relations file");
2458: if (fgets(buf, MPQS_STRING_LENGTH, FREL) == NULL)
2459: err(talker, "full relations file truncated?!");
2460: return buf;
2461: }
2462:
2463: /* the following two reside in src/basemath/ifactor1.c */
2464: long is_odd_power(GEN x, GEN *pt, long *mask);
2465: int miller(GEN n, long k);
2466:
2467: #define isprobableprime(n) (miller((n),17))
2468:
2469: static GEN
2470: mpqs_solve_linear_system(GEN kN, GEN N, long rel, long *FB, long size_of_FB)
2471: {
2472: pariFILE *pFREL;
2473: FILE *FREL;
2474: GEN X, Y_prod, N_or_kN, D1, base, res, new_res;
2475: long av=avma, av2, av3, lim, lim3, tetpil;
2476: long *fpos, *ei;
2477: long i, j, H_cols, H_rows;
2478: long res_last, res_next, res_size, res_max;
2479: mpqs_gauss_matrix m, ker_m;
2480: long done, flag, mask, rank;
2481:
2482: #ifdef MPQS_DEBUG
2483: N_or_kN = kN;
2484: #else
2485: N_or_kN = (DEBUGLEVEL >= 4 ? kN : N);
2486: #endif /* --GN */
2487: pFREL = pari_fopen(FREL_str, READ);
2488: FREL = pFREL->file;
2489: fpos = (long *) gpmalloc(rel * sizeof(long));
2490:
2491: m = mpqs_gauss_read_matrix(FREL, size_of_FB+1, rel, fpos);
2492: if (DEBUGLEVEL >= 7)
2493: {
2494: fprintferr("\\\\ MATRIX READ BY MPQS\nFREL=");
2495: mpqs_gauss_print_matrix(m, size_of_FB+1, rel);
2496: fprintferr("\n");
2497: }
2498:
2499: ker_m = mpqs_kernel(m, size_of_FB+1, rel, &rank);
2500: if (DEBUGLEVEL >= 4)
2501: {
2502: if (DEBUGLEVEL >= 7)
2503: {
2504: fprintferr("\\\\ KERNEL COMPUTED BY MPQS\nKERNEL=");
2505: mpqs_gauss_print_matrix(ker_m, rel, rank);
2506: fprintferr("\n");
2507: }
2508: /* put this here where the poor user has a chance of seeing it: */
2509: fprintferr("MPQS: Gauss done: kernel has rank %ld, taking gcds...\n",
2510: rank);
2511: }
2512:
2513: H_rows = rel;
2514: H_cols = rank;
2515:
2516: /* If the computed kernel turns out to be trivial, fail gracefully;
2517: main loop may try to find some more relations --GN */
2518: if (H_cols == 0)
2519: {
2520: if (DEBUGLEVEL >= 3)
2521: err(warner, "MPQS: no solutions found from linear system solver");
2522: /* ei has not yet been allocated */
2523: mpqs_gauss_destroy_matrix(m, size_of_FB+1, rel);
2524: mpqs_gauss_destroy_matrix(ker_m, rel, rank);
2525: free(fpos);
2526: pari_fclose(pFREL);
2527: avma = av;
2528: return NULL; /* no factors found */
2529: }
2530:
2531: /* If the rank is r, we can expect up to 2^r pairwise coprime factors,
2532: but it may well happen that a kernel basis vector contributes nothing
2533: new to the decomposition. We will allocate room for up to eight factors
2534: initially (certainly adequate when one or two basis vectors work),
2535: adjusting this down at the end to what we actually found, or up if
2536: we are very lucky and find more factors. In the upper half of our
2537: vector, we store information about which factors we know to be composite
2538: (zero) or believe to be composite ((long)NULL which is just 0) or suspect
2539: to be prime (un), or an exponent (deux or some t_INT) if it is a proper
2540: power */
2541: av2 = avma; lim = stack_lim(av2,1);
2542: if (rank > BITS_IN_LONG - 2)
2543: res_max = VERYBIGINT; /* this, unfortunately, is the common case */
2544: else
2545: res_max = 1L<<rank; /* max number of factors we can hope for */
2546: res_size = 8; /* no. of factors we can store in this res */
2547: res = cgetg(2*res_size+1, t_VEC);
2548: for (i=2*res_size; i; i--) res[i] = 0;
2549: res_next = res_last = 1;
2550:
2551: ei = (long *) gpmalloc((size_of_FB + 2) * sizeof(long));
2552:
2553: for (i = 0; i < H_cols; i++) /* loop over basis of kernel */
2554: {
2555: X = Y_prod = gun;
2556: memset(ei, 0, (size_of_FB + 2) * sizeof(long));
2557:
2558: av3 = avma; lim3 = stack_lim(av3,1);
2559: for (j = 0; j < H_rows; j++)
2560: {
2561: if (mpqs_gauss_get_bit(ker_m, j, i))
2562: Y_prod = mpqs_add_relation(Y_prod, N_or_kN, ei,
2563: mpqs_get_relation(fpos[j], FREL));
2564: if (low_stack(lim3, stack_lim(av3,1)))
2565: {
2566: if(DEBUGMEM>1) err(warnmem,"[1]: mpqs_solve_linear_system");
2567: Y_prod = gerepileupto(av3, Y_prod);
2568: }
2569: }
2570: Y_prod = gerepileupto(av3, Y_prod); /* unconditionally */
2571:
2572: av3 = avma; lim3 = stack_lim(av3,1);
2573: for (j = 2; j <= (size_of_FB + 1); j++)
2574: {
2575: if (ei[j] & 1)
2576: { /* this cannot happen --GN */
2577: mpqs_gauss_destroy_matrix(m, size_of_FB+1, rel);
2578: mpqs_gauss_destroy_matrix(ker_m, rel, rank);
2579: free(ei); free(fpos);
2580: fprintferr("MPQS: the combination of the relations is a nonsquare\n");
2581: err(bugparier, "factoring (MPQS)");
2582: }
2583: if (ei[j])
2584: {
2585: X = modii(mulii(X, powmodulo(stoi(FB[j]), stoi(ei[j]>>1), N_or_kN)),
2586: N_or_kN);
2587: if (low_stack(lim3, stack_lim(av3,1)))
2588: {
2589: if(DEBUGMEM>1) err(warnmem,"[2]: mpqs_solve_linear_system");
2590: X = gerepileupto(av3, X);
2591: }
2592: }
2593: }
2594: X = gerepileupto(av3, X); /* unconditionally */
2595: if (
2596: #ifdef MPQS_DEBUG
2597: 1 ||
2598: #endif
2599: DEBUGLEVEL >= 1)
2600: { /* this shouldn't happen either --GN */
2601: if (signe(modii(subii(sqri(X), sqri(Y_prod)), N_or_kN)))
2602: {
2603: fprintferr("MPQS: X^2 - Y^2 != 0 mod %s\n",
2604: (N_or_kN == kN ? "kN" : "N"));
2605: fprintferr("\tindex i = %ld\n", i);
2606: err(warner, "MPQS: wrong relation found after Gauss");
2607: }
2608: }
2609: /* if res_next < 3, we still haven't decomposed the original N, and
2610: will want both a gcd and its cofactor, overwriting N. Note that
2611: gcd(X+Y_prod,N)==1 forces gcd(X-Y_prod,N)==N, so we can skip X-Y_prod
2612: in such cases */
2613: done = 0; /* (re-)counts probably-prime factors
2614: (or powers whose bases we don't want to
2615: handle any further) */
2616: if (res_next < 3)
2617: {
2618: D1 = mppgcd(addii(X,Y_prod),N);
2619: if (!is_pm1(D1))
2620: {
2621: if ((flag = egalii(D1, N))) /* assignment */
2622: /* this one's useless, try the other one */
2623: D1 = mppgcd(subii(X,Y_prod),N);
2624: if (!flag || (!is_pm1(D1) && !egalii(D1,N)))
2625: { /* got something that works */
2626: if (DEBUGLEVEL >= 5)
2627: fprintferr("MPQS: splitting N after %ld kernel vector%s\n",
2628: i+1, (i? "s" : ""));
2629: (void)mpdivis(N, D1, N); /* divide N in place */
2630: res[1] = (long)N; /* we'll gcopy this anyway before exiting */
2631: res[2] = (long)D1;
2632: res_last = res_next = 3;
2633: if (res_max == 2) break; /* two out of two possible factors seen */
2634: mask = 7;
2635: /* (actually, 5th/7th powers aren't really worth the trouble. OTOH
2636: once we have the hooks for dealing with cubes, higher powers can
2637: be handled essentially for free) */
2638: if (isprobableprime(N)) { done++; res[res_size+1] = un; }
2639: else if (carrecomplet(N, &base))
2640: { /* squares could cost us a lot of time */
2641: affii(base, N); done++; res[res_size+1] = deux;
2642: if (DEBUGLEVEL >= 5)
2643: fprintferr("MPQS: decomposed a square\n");
2644: }
2645: else if ((flag = is_odd_power(N, &base, &mask))) /* assignment */
2646: {
2647: affii(base, N); done++; res[res_size+1] = (long)(stoi(flag));
2648: if (DEBUGLEVEL >= 5)
2649: fprintferr("MPQS: decomposed a %s\n",
2650: (flag == 3 ? "cube" :
2651: (flag == 5 ? "5th power" : "7th power")));
2652: }
2653: else res[res_size+1] = zero; /* known composite */
2654: mask = 7;
2655: if (isprobableprime(D1)) { done++; res[res_size+2] = un; }
2656: else if (carrecomplet(D1, &base))
2657: {
2658: res[2] = (long)base;
2659: done++; res[res_size+2] = deux;
2660: if (DEBUGLEVEL >= 5)
2661: fprintferr("MPQS: decomposed a square\n");
2662: }
2663: else if ((flag = is_odd_power(D1, &base, &mask))) /* assignment */
2664: {
2665: res[2] = (long)base;
2666: done++; res[res_size+2] = (long)(stoi(flag));
2667: if (DEBUGLEVEL >= 5)
2668: fprintferr("MPQS: decomposed a %s\n",
2669: (flag == 3 ? "cube" :
2670: (flag == 5 ? "5th power" : "7th power")));
2671: }
2672: else res[res_size+2] = zero; /* known composite */
2673: if (done == 2) break; /* both factors look prime or were powers */
2674: if (DEBUGLEVEL >= 5)
2675: fprintferr("MPQS: got two factors, looking for more...\n");
2676: /* and loop (after collecting garbage if necessary) */
2677: } /* else loop to next kernel basis vector */
2678: } /* D1!=1, nothing to be done when ==1 */
2679: }
2680: else /* we already have factors */
2681: {
2682: GEN X_plus_Y = addii(X, Y_prod);
2683: GEN X_minus_Y = NULL; /* will only be computed when needed */
2684: for (j=1; j < res_next; j++)
2685: { /* loop over known-composite factors */
2686: if (res[res_size+j] && res[res_size+j] != zero)
2687: {
2688: done++; continue; /* skip probable primes etc */
2689: }
2690: /* (actually, also skip square roots of squares etc. They are
2691: necessarily a lot smaller than the original N, and should
2692: be easy to deal with later) */
2693: av3 = avma;
2694: D1 = mppgcd(X_plus_Y, (GEN)(res[j]));
2695: if (is_pm1(D1)) continue; /* this one doesn't help us */
2696: if ((flag = egalii(D1, (GEN)(res[j]))))
2697: { /* bad one, try the other */
2698: avma = av3;
2699: if (!X_minus_Y) X_minus_Y = subii(X, Y_prod);
2700: av3 = avma;
2701: D1 = mppgcd(X_minus_Y, (GEN)(res[j]));
2702: }
2703: if (!flag || (!is_pm1(D1) && !egalii(D1, (GEN)(res[j]))))
2704: { /* got one which splits this factor */
2705: if (DEBUGLEVEL >= 5)
2706: fprintferr("MPQS: resplitting a factor after %ld kernel vectors\n",
2707: i+1); /* always plural */
2708: /* first make sure there's room for another factor */
2709: if (res_next > res_size)
2710: { /* need to reallocate (_very_ rare case) */
2711: long i1, new_size = 2*res_size;
2712: GEN new_res;
2713: if (new_size > res_max) new_size = res_max;
2714: new_res = cgetg(2*new_size+1, t_VEC);
2715: for (i1=2*new_size; i1>=res_next; i1--) new_res[i1] = 0;
2716: for (i1=1; i1<res_next; i1++)
2717: {
2718: new_res[i1] = res[i1]; /* factors */
2719: new_res[new_size+i1] = res[res_size+i1]; /* primality tags */
2720: }
2721: res = new_res; res_size = new_size; /* res_next unchanged */
2722: }
2723: /* now there is room; divide into existing factor and store the
2724: new gcd. Remove the known-composite flag from the old entry */
2725: (void)mpdivis((GEN)(res[j]), D1, (GEN)(res[j]));
2726: res[res_next] = (long)D1; res[res_size+j] = 0;
2727: if (++res_next > res_max) break; /* all possible factors seen */
2728: mask = 7; /* check for 3rd and 5th powers */
2729: if (isprobableprime((GEN)(res[j])))
2730: {
2731: done++; res[res_size+j] = un;
2732: }
2733: else if (carrecomplet((GEN)(res[j]), &base))
2734: {
2735: /* we no longer bother preserving the cloned N or what is left
2736: of it when we hit it here */
2737: res[j] = (long)base;
2738: done++; res[res_size+j] = deux;
2739: if (DEBUGLEVEL >= 5)
2740: fprintferr("MPQS: decomposed a square\n");
2741: }
2742: else if ((flag = is_odd_power((GEN)(res[j]), &base, &mask))) /* = */
2743: {
2744: res[j] = (long)base;
2745: done++; res[res_size+j] = (long)stoi(flag);
2746: if (DEBUGLEVEL >= 5)
2747: fprintferr("MPQS: decomposed a %s\n",
2748: (flag == 3 ? "cube" :
2749: (flag == 5 ? "5th power" : "7th power")));
2750: }
2751: else res[res_size+j] = zero; /* known composite */
2752: mask = 7;
2753: if (isprobableprime(D1))
2754: {
2755: done++; res[res_size+res_last] = un;
2756: }
2757: else if (carrecomplet(D1, &base))
2758: {
2759: res[res_last] = (long)base;
2760: done++; res[res_size+res_last] = deux;
2761: if (DEBUGLEVEL >= 5)
2762: fprintferr("MPQS: decomposed a square\n");
2763: }
2764: else if ((flag = is_odd_power(D1, &base, &mask))) /* assignment */
2765: {
2766: res[res_last] = (long)base;
2767: done++; res[res_size+res_last] = (long)stoi(flag);
2768: if (DEBUGLEVEL >= 5)
2769: fprintferr("MPQS: decomposed a %s\n",
2770: (flag == 3 ? "cube" :
2771: (flag == 5 ? "5th power" : "7th power")));
2772: }
2773: else res[res_size+res_last] = zero; /* known composite */
2774: } /* else it didn't help on this factor, try the next one... */
2775: else avma = av3; /* ...after destroying the gcds */
2776: } /* loop over known composite factors */
2777: if (res_next > res_last)
2778: {
2779: if (DEBUGLEVEL >= 5 && done < res_last)
2780: fprintferr("MPQS: got %ld factors%s\n",
2781: res_last,
2782: (done < res_last ? ", looking for more..." : ""));
2783: res_last = res_next;
2784: }
2785: /* we should break out of the outer loop when we have seen res_max
2786: factors, and also when all current factors are probable primes */
2787: if (res_next > res_max || done == res_next - 1) break;
2788: /* else continue loop over kernel basis */
2789: } /* end case of further splitting of existing factors */
2790: /* garbage collection */
2791: if (low_stack(lim, stack_lim(av2,1)))
2792: {
2793: long i1;
2794: if(DEBUGMEM>1) err(warnmem,"[3]: mpqs_solve_linear_system");
2795: tetpil=avma;
2796: /* gcopy would have a problem with our NULL pointers... */
2797: new_res = cgetg(lg(res), t_VEC);
2798: for (i1=2*res_size; i1>=res_next; i1--) new_res[i1] = 0;
2799: for (i1=1; i1<res_next; i1++)
2800: {
2801: new_res[i1] =
2802: (isonstack((GEN)(res[i1])) ? licopy((GEN)(res[i1])) : res[i1]);
2803: new_res[res_size+i1] = res[res_size+i1];
2804: }
2805: res = gerepile(av2, tetpil, new_res);
2806: /* discard X and Y_prod */
2807: }
2808: } /* for (loop over kernel basis) */
2809:
2810: if (res_next < 3) /* still no factors seen */
2811: {
2812: pari_fclose(pFREL);
2813: mpqs_gauss_destroy_matrix(m, size_of_FB+1, rel);
2814: mpqs_gauss_destroy_matrix(ker_m, rel, rank);
2815: free(ei); free(fpos);
2816: avma = av;
2817: return NULL; /* no factors found */
2818: }
2819: /* normal case: convert internal format to ifac format as described in
2820: src/basemath/ifactor1.c (triples of components: value, exponent, class
2821: [unknown, known composite, known prime]), clean up and return result */
2822: pari_fclose(pFREL);
2823: mpqs_gauss_destroy_matrix(m, size_of_FB+1, rel);
2824: mpqs_gauss_destroy_matrix(ker_m, rel, rank);
2825: free(ei); free(fpos);
2826: res_last = res_next - 1; /* number of distinct factors found */
2827: new_res = cgetg(3*res_last + 1, t_VEC);
2828: if (DEBUGLEVEL >= 6)
2829: fprintferr("MPQS: wrapping up vector of %ld factors\n", res_last);
2830: for (i=1,j=1; i <= res_last; i++)
2831: {
2832: new_res[j++] = /* value of factor */
2833: isonstack((GEN)(res[i])) ? licopy((GEN)(res[i])) : res[i];
2834: flag = res[res_size+i];
2835: new_res[j++] = /* exponent */
2836: flag ? /* flag was zero or un or ... */
2837: (flag == zero ? un :
2838: (isonstack((GEN)flag) ? licopy((GEN)flag) : flag)
2839: ) :
2840: un; /* flag was (long)NULL */
2841: new_res[j++] = /* class */
2842: flag == zero ? zero : /* known composite */
2843: (long)NULL; /* base of power or suspected prime --
2844: mark as `unknown' */
2845: if (DEBUGLEVEL >= 6)
2846: fprintferr("\tpackaging %ld: %Z ^%ld (%s)\n", i, res[i],
2847: itos((GEN)(new_res[j-2])),
2848: (flag == zero ? "comp." : "unknown"));
2849: }
2850: return gerepileupto(av, new_res);
2851: }
2852:
2853: /******************************/
2854:
2855: /* All percentages below are actually fixed-point quantities scaled by 10
2856: (value of 1 means 0.1%, 1000 means 100%). --GN */
2857:
2858: #define max_percentage 1500 /* give up when nothing found after ~1.5
2859: times the required number of relations has
2860: been computed (n might be a prime power
2861: or the parameters might be exceptionally
2862: unfortunate for it) --GN */
2863:
2864: /**
2865: ** Factors N using the self-initializing multipolynomial quadratic sieve
2866: ** (SIMPQS). Returns one of the two factors, or (usually) a vector of
2867: ** factors and exponents and information about which ones are definitely
2868: ** still composite, or NULL when something goes wrong or when we can't
2869: ** seem to make any headway.
2870: **/
2871:
2872: /* TODO: this function to be renamed mpqs_main() with several extra para-
2873: meters, with mpqs() as a wrapper for the standard case, so we can do
2874: partial runs distributed across several machines etc. (not necessarily
2875: from gp, but certainly from a small dedicated C program). --GN */
2876:
2877: GEN
2878: mpqs(GEN N)
2879: {
2880: GEN fact; /* will in the end hold our factor(s) */
2881: ulong decimal_digits; /* number of decimals of N */
2882: long tc; /* number of candidates found in one
2883: iteration */
2884: long t; /* number of full relations found in one
2885: iteration */
2886: long tp; /* number of full relations found in one
2887: iteration after combining large prime
2888: relations */
2889: long k; /* the multiplier */
2890: GEN kN; /* k * N */
2891: double sqrt_kN; /* double approximation of sqrt(k*N) */
2892: long *FB; /* the factor base */
2893: ulong size_of_FB; /* the size of the factor base */
2894: ulong FB_overshoot; /* the number of relations aimed for */
2895: double log_multiplier; /* maps logarithms of FB elements to
2896: unsigned chars */
2897: double tolerance; /* used in sieving with logarithms */
2898: unsigned char *log_FB; /* logarithms of FB elements mapped to
2899: unsigned chars */
2900: unsigned char *sieve_array; /* the sieve array */
2901: unsigned char *sieve_array_end; /* the end of the sieve array */
2902: ulong M; /* the sieve interval size [-M, M] */
2903: ulong starting_sieving_index; /* the first element in the FB used for
2904: sieving */
2905: long last_prime_in_FB; /* the biggest prime in FB */
2906: long *sqrt_mod_p_kN; /* array containing sqrt(K*N) mod p_i,
2907: where p_i are FB elements */
2908: GEN A, B, inv_A4; /* coefficients of polynomials in (SI) */
2909: long *Q_prime_glob, *Q_prime; /* used for self initialization (SI) */
2910: GEN *H_i; /* will be used to compute the consecutive
2911: B during (SI) */
2912: long start_index_FB_for_A; /* the first index in the FB used for
2913: computing A during (SI), minus 1 */
2914: ulong no_of_primes_in_A; /* number of primes dividing A
2915: during (SI) */
2916: ulong total_no_of_primes_for_A; /* number of primes used for finding A
2917: during (SI) */
2918: long max_no_of_B; /* how many B's to use before choosing
2919: a new A during (SI) */
2920: long index_i; /* denotes the current pair of coeffs
2921: (A, B_i). If i == 0, a new A is
2922: generated */
2923: ulong bin_index; /* used for choosing the prime factors for
2924: for the A coeffs in mpqs_self_init() */
2925: long i, p; /* counters */
2926:
2927: /* Let (z0, z1) be the roots of Q(x) = A x^2 + Bx +C mod p_i;
2928: then we know that
2929: Q(z1 + i k) == 0 mod p_i and Q(z2 + i k) == 0 mod p_i */
2930:
2931: long *start_1; /* stores the position where p_i
2932: divides Q(x) for the first time,
2933: using z1 */
2934: long *start_2; /* same, but using z2 */
2935: long *inv_A2; /* inv_A2[i] = 1/(2*A) mod p_i */
2936: long **inv_A_H_i; /* inv_A_H_i[i][j] = 1/A H_i[i] mod p_j */
2937: long *candidates; /* positions in the sieve_array which
2938: probably factor over the FB */
2939: long lp_bound; /* size limit for large primes */
2940: ulong sort_interval; /* these determine when to sort and merge */
2941: ulong followup_sort_interval; /* the temporary files (scaled percentages) */
2942: long total_full_relations = 0;
2943: long total_partial_relations = 0;
2944: pariFILE *pFREL;
2945: pariFILE *pFNEW;
2946: pariFILE *pLPREL;
2947: pariFILE *pLPNEW;
2948: pariFILE *pCOMB;
2949: FILE *FNEW;
2950: FILE *LPNEW;
2951: long percentage = 0; /* scaled by 10, see comment above */
2952: long total_candidates_number = 0;
2953: long useless_candidates = 0; /* another scaled percentage */
2954: long vain_iterations = 0;
2955: long good_iterations = 0;
2956: long iterations = 0;
2957: long av = avma;
2958:
2959: /* state flags for cleanup if we get interrupted --GN */
2960: static char all_clean = 1; /* set to 0 while mpqs() is busy */
2961:
2962: if (DEBUGLEVEL >= 4)
2963: {
2964: if (!all_clean) { /* TODO: clean up... */ }
2965: (void) timer2(); /* clear timer */
2966: fprintferr("MPQS: number to factor N = %Z\n", N);
2967: }
2968:
2969: all_clean = 0; /* activate clean-up trap */
2970: {
2971: char *s = GENtostr(N);
2972: decimal_digits = strlen(s);
2973: free(s);
2974: }
2975:
2976: if (decimal_digits >= 108) /* changed to match Denny's parameters */
2977: {
2978: err(warner, "MPQS: number too big to be factored with MPQS, giving up");
2979: avma = av;
2980: all_clean = 1;
2981: return NULL;
2982: }
2983:
2984: if (DEBUGLEVEL >= 4)
2985: fprintferr("MPQS: factoring number of %ld decimal digits\n",
2986: decimal_digits);
2987:
2988: if (decimal_digits >= 70)
2989: err(warner,
2990: "MPQS: the factorization of this number will take %s hours",
2991: decimal_digits >= 84 ? "many" : "several");
2992:
2993: k = mpqs_find_k(N, 5);
2994: if (DEBUGLEVEL >= 5)
2995: fprintferr("MPQS: found multiplier %ld for N\n", k);
2996:
2997: kN = mulis(N, k);
2998: {
2999: char *s = GENtostr(kN);
3000: decimal_digits = strlen(s);
3001: free(s);
3002: }
3003:
3004: if (DEBUGLEVEL >= 5)
3005: {
3006: fprintferr("MPQS: kN = %Z\n", kN);
3007: fprintferr("MPQS: kN has %ld decimal digits\n",
3008: decimal_digits);
3009: }
3010:
3011: mpqs_read_parameters(decimal_digits, &tolerance, &M, &size_of_FB,
3012: &FB_overshoot,
3013: &no_of_primes_in_A, &total_no_of_primes_for_A,
3014: &max_no_of_B, &starting_sieving_index,
3015: &sort_interval, &followup_sort_interval);
3016: {
3017: double bytesize = (size_of_FB + 1)/(8.*1048576.);
3018: bytesize *= FB_overshoot;
3019: if (bytesize > 32.)
3020: {
3021: err(warner,
3022: "MPQS: Gauss elimination will require more than 32MBy of memory");
3023: if (DEBUGLEVEL >= 1)
3024: fprintferr("\t(estimated memory needed: %4.1fMBy)\n", bytesize);
3025: }
3026: }
3027:
3028: if (DEBUGLEVEL >= 4)
3029: {
3030: fprintferr("MPQS: sieving interval = [%ld, %ld]\n",
3031: -(long)M, M);
3032: fprintferr("MPQS: size of factor base = %ld\n",
3033: size_of_FB);
3034: fprintferr("MPQS: striving for %ld relations\n",
3035: FB_overshoot);
3036: fprintferr("MPQS: first sorting at %ld%%, then every %3.1f%% / %3.1f%%\n",
3037: sort_interval/10, followup_sort_interval/10.,
3038: followup_sort_interval/20.);
3039: fprintferr("MPQS: initial sieving index = %ld\n",
3040: starting_sieving_index);
3041: }
3042: sieve_array =
3043: (unsigned char *) gpmalloc((M << 1) * sizeof(unsigned char));
3044: sieve_array_end = sieve_array + (M << 1) - 1;
3045:
3046: if (DEBUGLEVEL >= 5)
3047: fprintferr("MPQS: creating factor base FB of size = %ld\n",
3048: size_of_FB);
3049:
3050: FB = mpqs_create_FB(size_of_FB, kN, k, &p);
3051:
3052: /* the following this way round to avoid overflow: --GN */
3053: lp_bound = FB[FB[0]+1];
3054: if (DEBUGLEVEL >= 4)
3055: fprintferr("MPQS: largest prime in FB = %ld\n",
3056: lp_bound);
3057: if (lp_bound > MPQS_LP_BOUND)
3058: lp_bound = MPQS_LP_BOUND;
3059: lp_bound *= MPQS_LP_FACTOR;
3060: if (DEBUGLEVEL >= 4)
3061: fprintferr("MPQS: bound for `large primes\' = %ld\n",
3062: lp_bound);
3063:
3064: if (p != 0)
3065: {
3066: free(FB);
3067: free(sieve_array);
3068: if (DEBUGLEVEL >= 4)
3069: fprintferr("\nMPQS: found factor = %ld whilst creating factor base\n",
3070: p);
3071: if (mpqs_diffptr == diffptr) mpqs_diffptr = NULL;
3072: avma = av;
3073: /* TODO: free the FB etc!!! */
3074: all_clean = 1;
3075: return stoi(p);
3076: }
3077:
3078: if (DEBUGLEVEL >= 5)
3079: {
3080: fprintferr("MPQS: computing logarithm approximations for p_i in FB\n");
3081: fprintferr("MPQS: computing sqrt(k*N) mod p_i\n");
3082: }
3083: last_prime_in_FB = FB[FB[0] + 1];
3084: log_multiplier = (2 * (unsigned char) (0.5 * log2 (gtodouble(kN)) +
3085: log2 ((double) M) - tolerance *
3086: log2 ((double) last_prime_in_FB)));
3087: log_multiplier = 127.0 / log_multiplier;
3088:
3089: log_FB =
3090: (unsigned char *) gpmalloc((size_of_FB + 2) * sizeof(unsigned char));
3091: sqrt_mod_p_kN =
3092: (long *) gpmalloc((size_of_FB + 2) * sizeof(long));
3093:
3094: for (i = 2; (ulong)i < size_of_FB + 2; i++)
3095: {
3096: long av1 = avma;
3097: p = FB[i];
3098:
3099: /* compute the approximations of the logarithms of p_i */
3100: log_FB[i] =
3101: (unsigned char) (log_multiplier * log2 ((double) p) * 2);
3102:
3103: /* compute the x_i such that x_i^2 = (kN % p_i) mod p_i */
3104: sqrt_mod_p_kN[i] = itos(mpsqrtmod(modis(kN, p), stoi(p)));
3105: avma = av1;
3106: }
3107:
3108: if (DEBUGLEVEL >= 5)
3109: fprintferr("MPQS: allocating arrays for self-initialization\n");
3110:
3111: /* the size of coefficient A should be approximately */
3112: /* sqrt(kN)/M, so the size of the primes p dividing */
3113: /* A should be approximately (sqrt(kN/M))^(1/no_of_primes_in_A) */
3114:
3115: /* compute expected size of coefficients, which also defines the tolerance
3116: for the logarithmic sieve */
3117: sqrt_kN = sqrt(gtodouble(kN));
3118: tolerance = sqrt_kN / M;
3119: tolerance = pow(2.0, log2(tolerance) / no_of_primes_in_A);
3120: if (tolerance > FB[size_of_FB - 1]) /* ??? */
3121: err(talker, "MPQS: number of prime factors in A is too small");
3122:
3123: /* choose primes in FB to use for building the A coefficient */
3124: start_index_FB_for_A = 2;
3125: while (FB[start_index_FB_for_A] < tolerance)
3126: start_index_FB_for_A++;
3127:
3128: if (start_index_FB_for_A > 7)
3129: start_index_FB_for_A -= (no_of_primes_in_A >> 1);
3130:
3131: /* ensure that the primes to be used for A do occur in the FB */
3132: if (start_index_FB_for_A + total_no_of_primes_for_A >= size_of_FB)
3133: start_index_FB_for_A = size_of_FB - total_no_of_primes_for_A - 1;
3134:
3135: /* ensure that the multiplier k is smaller than the prime divisors of A */
3136: if (k >= FB[start_index_FB_for_A])
3137: {
3138: do
3139: start_index_FB_for_A++;
3140: while (k >= FB[start_index_FB_for_A]);
3141: if (start_index_FB_for_A + total_no_of_primes_for_A >= size_of_FB)
3142: err(talker,
3143: "MPQS: number of primes for A is too large, or FB too small");
3144: /* this shouldn't happen */
3145: }
3146:
3147: /* now a selection of total_no_of_primes_for_A consecutive primes
3148: * p = FB[start_index_FB_for_A+1], ... is chosen from the factor base */
3149:
3150: /* (no need to copy them around as long as they are consecutive entries
3151: in FB --GN) */
3152: Q_prime =
3153: (long *) gpmalloc(no_of_primes_in_A * sizeof(long));
3154: #if 1
3155: Q_prime_glob = &(FB[start_index_FB_for_A + 1]);
3156: #else
3157: Q_prime_glob =
3158: (long *) gpmalloc(total_no_of_primes_for_A * sizeof(long));
3159: for (i = 0; i < total_no_of_primes_for_A; i++)
3160: Q_prime_glob[i] = FB[start_index_FB_for_A + i + 1];
3161: #endif
3162:
3163: /* used for choosing the correct A coeffs in mpqs_self_init() */
3164: bin_index = 0;
3165:
3166: /* will be used to compute the consecutive B during
3167: self-initialization */
3168: H_i = (GEN *) gpmalloc(no_of_primes_in_A * sizeof(GEN));
3169: for (i = 0; (ulong)i < no_of_primes_in_A; i++)
3170: H_i[i] = cgeti(2 * total_no_of_primes_for_A);
3171:
3172: inv_A_H_i =
3173: (long **) gpmalloc(total_no_of_primes_for_A * sizeof(long *));
3174: for (i = 0; (ulong)i < total_no_of_primes_for_A; i++)
3175: inv_A_H_i[i] = (long *) gpmalloc((size_of_FB + 2) * sizeof(long));
3176:
3177: start_1 = (long *) gpmalloc((size_of_FB + 2) * sizeof(long));
3178: start_2 = (long *) gpmalloc((size_of_FB + 2) * sizeof(long));
3179:
3180: /* the next will hold 1/(2*A) mod p_i, p_i in FB */
3181: inv_A2 = (long *) gpmalloc((size_of_FB + 2) * sizeof(long));
3182:
3183: candidates = (long *) gpmalloc(MPQS_CANDIDATE_ARRAY_SIZE * sizeof(long));
3184:
3185: if (DEBUGLEVEL >= 5)
3186: {
3187: fprintferr("MPQS: index range of primes for A: [%ld, %ld]\n",
3188: start_index_FB_for_A + 1,
3189: start_index_FB_for_A + total_no_of_primes_for_A);
3190: fprintferr("MPQS: coefficients A will be built from %ld primes each\n",
3191: no_of_primes_in_A);
3192: }
3193:
3194: /* CLEAN UP define */
3195:
3196: #define CLEANUP() \
3197: free(FB); \
3198: free(sieve_array); \
3199: free(log_FB); \
3200: free(sqrt_mod_p_kN); \
3201: free(Q_prime); \
3202: /* free(Q_prime_glob); */ \
3203: free(H_i); \
3204: for (i = 0; (ulong)i < total_no_of_primes_for_A; i++) \
3205: free(inv_A_H_i[i]); \
3206: free(inv_A_H_i); \
3207: free(start_1); \
3208: free(start_2); \
3209: free(inv_A2); \
3210: free(candidates);
3211:
3212: /*
3213: * main loop which
3214: * - computes polynomials and their zeros (self initialization)
3215: * - does the sieving
3216: * - tests candidates of the sieve array
3217: */
3218:
3219: /* denotes the current pair of coeffs (A, B_i). If i == 0 a new A is
3220: generated */
3221: index_i = -1;
3222:
3223: /* the coefficients of the polynomial currently used for sieving */
3224: /* A will be at most as long as the product of
3225: no_of_primes_in_A word-sized primes (was total_..., GN) */
3226: /* B needs at most total_no_of_primes_for_A + 1 words [CHECK] --GN */
3227: /* Take the double is just for safety ... */
3228: A = cgeti(2 * no_of_primes_in_A);
3229: B = cgeti(2 * total_no_of_primes_for_A + 1);
3230:
3231: /* 1/4A mod kN */
3232: inv_A4 = cgeti(lg(kN));
3233:
3234: if (DEBUGLEVEL >= 5)
3235: fprintferr("MPQS: starting main loop\n");
3236: /* compute names for the temp files we'll need */
3237: FREL_str = mpqs_get_filename("FREL");
3238: FNEW_str = mpqs_get_filename("FNEW");
3239: LPREL_str = mpqs_get_filename("LPREL");
3240: LPNEW_str = mpqs_get_filename("LPNEW");
3241: COMB_str = mpqs_get_filename("COMB");
3242: TMP_str = mpqs_get_filename("LPTMP");
3243:
3244: /* just to create the temp files */
3245: pFREL = pari_safefopen(FREL_str, WRITE); pari_fclose(pFREL);
3246: pLPREL = pari_safefopen(LPREL_str, WRITE); pari_fclose(pLPREL);
3247:
3248: pFNEW = pari_safefopen(FNEW_str, WRITE); FNEW = pFNEW->file;
3249: pLPNEW= pari_safefopen(LPNEW_str, WRITE); LPNEW = pLPNEW->file;
3250: for(;;)
3251: {
3252: /* at start of loop, FNEW and LPNEW are open for writing */
3253: iterations++;
3254: /* when all of the B's have already been used, choose new A ;
3255: this is indicated by setting index_i to 0 */
3256: if (index_i == (max_no_of_B - 1))
3257: index_i = 0;
3258: else
3259: index_i++;
3260:
3261: /* self initialization: compute polynomial and its zeros */
3262:
3263: mpqs_self_init(A, B, N, kN, FB, sqrt_mod_p_kN, start_1, start_2,
3264: no_of_primes_in_A, total_no_of_primes_for_A,
3265: M, inv_A_H_i, Q_prime, Q_prime_glob, H_i,
3266: index_i, start_index_FB_for_A, inv_A2, inv_A4,
3267: &bin_index, &fact);
3268: if (bin_index >= (1ul << total_no_of_primes_for_A)) /* wraparound */
3269: {
3270: /* TODO: increase some parameters. For the moment, simply give up */
3271: CLEANUP();
3272: pari_fclose(pFNEW);
3273: pari_fclose(pLPNEW);
3274: /* FREL, LPREL are closed at this point */
3275: pari_unlink(FREL_str);
3276: pari_unlink(FNEW_str);
3277: pari_unlink(LPREL_str);
3278: pari_unlink(LPNEW_str);
3279: if (mpqs_diffptr == diffptr) mpqs_diffptr = NULL;
3280: all_clean = 1;
3281: avma=av; return NULL;
3282: }
3283: if (fact) /* A4 not coprime to kN, fact divides N */
3284: { /* (cannot actually happen) */
3285: if (is_pm1(fact) || egalii(fact,N))
3286: continue; /* cannot use the current polynomial */
3287: CLEANUP();
3288: pari_fclose(pFNEW);
3289: pari_fclose(pLPNEW);
3290: pari_unlink(FREL_str);
3291: pari_unlink(FNEW_str);
3292: pari_unlink(LPREL_str);
3293: pari_unlink(LPNEW_str);
3294: if (mpqs_diffptr == diffptr) mpqs_diffptr = NULL;
3295: all_clean = 1;
3296: if (DEBUGLEVEL >= 4)
3297: {
3298: fprintferr("MPQS: whilst trying to invert A4 mod kN,\n");
3299: fprintferr("\tfound factor = %Z\n", fact);
3300: flusherr();
3301: }
3302: return gerepileupto(av, fact);
3303: }
3304:
3305: if (DEBUGLEVEL >= 6)
3306: {
3307: if (!index_i)
3308: fprintferr("MPQS: chose prime pattern 0x%lX for A\n",
3309: bin_index);
3310: if (signe(B) < 0)
3311: {
3312: setsigne(B,1);
3313: fprintferr("MPQS: chose Q_%ld(x) = %Z x^2 - %Z x + C\n",
3314: index_i, A, B);
3315: setsigne(B,-1);
3316: }
3317: else
3318: {
3319: fprintferr("MPQS: chose Q_%ld(x) = %Z x^2 + %Z x + C\n",
3320: index_i, A, B);
3321: }
3322: }
3323:
3324: mpqs_sieve(FB, log_FB, start_1, start_2, sieve_array,
3325: sieve_array_end, M, starting_sieving_index);
3326:
3327: tc = mpqs_eval_sieve(sieve_array, M, candidates);
3328: total_candidates_number += tc;
3329: if (DEBUGLEVEL >= 6)
3330: fprintferr("MPQS: found %lu candidate%s\n",
3331: tc, (tc==1? "" : "s"));
3332:
3333: if (tc)
3334: {
3335: good_iterations++;
3336: t = mpqs_eval_candidates(A, inv_A4, B, kN, k, sqrt_kN, FB,
3337: start_1, start_2,
3338: M, bin_index, candidates, tc, lp_bound,
3339: start_index_FB_for_A, FNEW, LPNEW);
3340: }
3341: else
3342: t = 0;
3343:
3344: total_full_relations += t;
3345: percentage =
3346: (long)((1000.0 * total_full_relations) / FB_overshoot);
3347: useless_candidates =
3348: (long)((1000.0 * (total_candidates_number - total_full_relations))
3349: / (total_candidates_number ? total_candidates_number : 1));
3350:
3351: if ((ulong)percentage < sort_interval)
3352: continue; /* most main loops end here... */
3353:
3354: /* Extra processing when we have completed a sort interval: */
3355: if (DEBUGLEVEL >= 3)
3356: {
3357: if (DEBUGLEVEL >= 4)
3358: fprintferr("\nMPQS: passing the %3.1f%% checkpoint, time = %ld ms\n",
3359: sort_interval/10., timer2());
3360: else
3361: fprintferr("\nMPQS: passing the %3.1f%% checkpoint\n",
3362: sort_interval/10.);
3363: flusherr();
3364: }
3365: /* sort LPNEW and merge it into LPREL, diverting combinables into COMB */
3366: pari_fclose(pLPNEW);
3367: mpqs_sort_lp_file(LPNEW_str);
3368: tp = mpqs_mergesort_lp_file(LPREL_str, LPNEW_str, 0);
3369: pLPNEW = pari_fopen(LPNEW_str, WRITE); /* NOT safefopen */
3370: LPNEW = pLPNEW->file;
3371:
3372: /* combine whatever there is to be combined */
3373: if (tp > 0)
3374: {
3375: /* build full relations out of large prime relations */
3376: pCOMB = pari_fopen(COMB_str, READ);
3377: tp = mpqs_combine_large_primes(pCOMB->file, FNEW, size_of_FB, N, kN, &fact
3378: #ifdef MPQS_DEBUG
3379: , FB
3380: #endif
3381: );
3382: pari_fclose(pCOMB);
3383: pari_unlink(COMB_str);
3384: /* now FREL, LPREL are closed and FNEW, LPNEW are still open */
3385: if (fact) /* factor found during combining */
3386: {
3387: /* clean up */
3388: CLEANUP();
3389: pari_fclose(pLPNEW);
3390: pari_fclose(pFNEW);
3391: pari_unlink(FREL_str);
3392: pari_unlink(FNEW_str);
3393: pari_unlink(LPREL_str);
3394: pari_unlink(LPNEW_str);
3395: if (DEBUGLEVEL >= 4)
3396: {
3397: fprintferr("\nMPQS: split N whilst combining, time = %ld ms\n",
3398: timer2());
3399: fprintferr("MPQS: found factor = %Z\n", fact);
3400: }
3401: if (mpqs_diffptr == diffptr) mpqs_diffptr = NULL;
3402: all_clean = 1;
3403: return gerepileupto(av, fact);
3404: }
3405: }
3406:
3407: /* sort FNEW and merge it into FREL */
3408: pari_fclose(pFNEW);
3409: mpqs_sort_lp_file(FNEW_str);
3410: total_full_relations = mpqs_mergesort_lp_file(FREL_str, FNEW_str, 1);
3411: /* this being the definitive count (combinables combined, and
3412: duplicates removed) */
3413: /* FNEW stays closed until we know whether we need to reopen it for
3414: another iteration */
3415:
3416: total_partial_relations += tp;
3417:
3418: /* note that due to the removal of duplicates, percentage may actually
3419: decrease at this point. This may look funny in the diagnostics but
3420: is nothing to worry about, since we _are_ making progress neverthe-
3421: less. --GN */
3422: percentage =
3423: (long)((1000.0 * total_full_relations) / FB_overshoot);
3424: vain_iterations =
3425: (long)((1000.0 * (iterations - good_iterations))
3426: / iterations);
3427:
3428: /* TODO: the following could be improved (extrapolate from past
3429: experience how many combined full relations we can expect in
3430: addition to those we're finding directly) --GN */
3431: if (percentage >= 840)
3432: {
3433: if (percentage >= 980)
3434: sort_interval = percentage + 10;
3435: else
3436: sort_interval = percentage + (followup_sort_interval >> 1);
3437: if (sort_interval >= 1000 && percentage < 1000)
3438: sort_interval = 1000;
3439: }
3440: else
3441: {
3442: sort_interval = percentage + followup_sort_interval;
3443: }
3444:
3445: if (DEBUGLEVEL >= 4)
3446: {
3447: fprintferr("MPQS: done sorting%s, time = %ld ms\n",
3448: tp > 0 ? " and combining" : "",
3449: timer2());
3450: fprintferr("MPQS: found %3.1f%% of the required relations\n",
3451: percentage/10.);
3452: if (DEBUGLEVEL >= 5)
3453: {
3454: /* total_full_relations are always plural */
3455: fprintferr("MPQS: found %ld full relations\n",
3456: total_full_relations);
3457: fprintferr("MPQS: (%ld of these from partial relations)\n",
3458: total_partial_relations);
3459: fprintferr("MPQS: %4.1f%% useless candidates\n",
3460: useless_candidates/10.);
3461: fprintferr("MPQS: %4.1f%% of the iterations yielded no candidates\n",
3462: vain_iterations/10.);
3463: fprintferr("MPQS: next checkpoint at %3.1f%%\n",
3464: sort_interval/10.);
3465: }
3466: }
3467:
3468: if (percentage < 1000)
3469: {
3470: pFNEW = pari_fopen(FNEW_str, WRITE); /* NOT safefopen */
3471: FNEW = pFNEW->file;
3472: /* at this point, LPNEW and FNEW are again open for writing */
3473: continue; /* main loop */
3474: }
3475:
3476: /* percentage >= 1000, which implies total_full_relations > size_of_FB:
3477: try finishing it off */
3478:
3479: /* solve the system over F_2 */
3480: if (DEBUGLEVEL >= 4)
3481: fprintferr("\nMPQS: starting Gauss over F_2 on %ld distinct relations\n",
3482: total_full_relations);
3483: fact = mpqs_solve_linear_system(kN, N,
3484: total_full_relations,
3485: FB, size_of_FB);
3486:
3487: /* if solution found, clean up and return */
3488: if (fact)
3489: {
3490: /* clean up */
3491: CLEANUP();
3492: pari_fclose(pLPNEW);
3493: pari_unlink(FREL_str);
3494: pari_unlink(FNEW_str);
3495: pari_unlink(LPREL_str);
3496: pari_unlink(LPNEW_str);
3497: if (DEBUGLEVEL >= 4)
3498: {
3499: fprintferr("\nMPQS: time in Gauss and gcds = %ld ms\n",
3500: timer2());
3501: if (typ(fact) == t_INT)
3502: fprintferr("MPQS: found factor = %Z\n", fact);
3503: else
3504: {
3505: long j, nf=(lg(fact)-1)/3;
3506: if (nf==2)
3507: fprintferr("MPQS: found factors = %Z\n\tand %Z\n",
3508: fact[4], fact[1]);
3509: else
3510: {
3511: fprintferr("MPQS: found %ld factors =\n", nf);
3512: for (j=nf; j; j--)
3513: fprintferr("\t%Z%s\n", fact[3*j-2], (j>1? "," : ""));
3514: }
3515: }
3516: }
3517: if (mpqs_diffptr == diffptr) mpqs_diffptr = NULL;
3518: all_clean = 1;
3519: /* nuisance: fact may not be safe for a gcopy(), and thus
3520: gerepilecopy(av, fact) will usually segfault on
3521: one of the NULL markers. However, it is already a nice
3522: connected object, and it resides already the top of the
3523: stack, so... :^) --GN */
3524: return gerepileupto(av, fact);
3525: }
3526: else
3527: {
3528: if (DEBUGLEVEL >= 4)
3529: {
3530: fprintferr("\nMPQS: time in Gauss and gcds = %ld ms\n",
3531: timer2());
3532: fprintferr("MPQS: no factors found.\n");
3533: if (percentage <= max_percentage)
3534: fprintferr("\nMPQS: restarting sieving ...\n");
3535: else
3536: fprintferr("\nMPQS: giving up.\n");
3537: }
3538: if (percentage > max_percentage)
3539: {
3540: /* clean up */
3541: CLEANUP();
3542: pari_fclose(pLPNEW);
3543: pari_unlink(FREL_str);
3544: pari_unlink(FNEW_str);
3545: pari_unlink(LPREL_str);
3546: pari_unlink(LPNEW_str);
3547: if (mpqs_diffptr == diffptr) mpqs_diffptr = NULL;
3548: all_clean = 1;
3549: avma=av; return NULL;
3550: }
3551: pFNEW = pari_fopen(FNEW_str, WRITE); /* NOT safefopen */
3552: FNEW = pFNEW->file;
3553: }
3554: } /* main loop */
3555: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>