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