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