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