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