[BACK]Return to mpqs.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari-2.2 / src / modules

Annotation of OpenXM_contrib/pari-2.2/src/modules/mpqs.c, Revision 1.1.1.1

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

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>