Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - modules - mpqs.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16746-c2cb716) Lines: 920 1136 81.0 %
Date: 2014-08-31 Functions: 38 38 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 447 706 63.3 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000  The PARI group.
       2                 :            : 
       3                 :            : This file is part of the PARI/GP package.
       4                 :            : 
       5                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       6                 :            : terms of the GNU General Public License as published by the Free Software
       7                 :            : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8                 :            : ANY WARRANTY WHATSOEVER.
       9                 :            : 
      10                 :            : Check the License for details. You should have received a copy of it, along
      11                 :            : with the package; see the file 'COPYING'. If not, write to the Free Software
      12                 :            : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13                 :            : 
      14                 :            : /* Written by Thomas Papanikolaou and Xavier Roblot
      15                 :            :  *
      16                 :            :  * Implementation of the Self-Initializing Multi-Polynomial Quadratic Sieve
      17                 :            :  * based on code developed as part of the LiDIA project
      18                 :            :  * (http://www.informatik.tu-darmstadt.de/TI/LiDIA/)
      19                 :            :  *
      20                 :            :  * Extensively modified by The PARI group.
      21                 :            :  */
      22                 :            : /* Notation commonly used in this file, and sketch of algorithm:
      23                 :            :  *
      24                 :            :  * Given an odd integer N > 1 to be factored, we throw in a small odd and
      25                 :            :  * squarefree multiplier k so as to make kN congruent 1 mod 4 and to have
      26                 :            :  * many small primes over which X^2 - kN splits.  We compute a factor base
      27                 :            :  * FB of such primes, and then essentially look for values x0 such that
      28                 :            :  * Q0(x0) = x0^2 - kN can be decomposed over this factor base, up to a
      29                 :            :  * possible factor dividing k and a possible "large prime".  Relations
      30                 :            :  * involving the latter can be combined into full relations (working mod N
      31                 :            :  * now) which don't, and full relations, by Gaussian elimination over the
      32                 :            :  * 2-element field for the exponent vectors, will finally lead us to an
      33                 :            :  * expression X^2 - Y^2 divisible by N and hopefully to a nontrivial
      34                 :            :  * splitting when we compute gcd(X +- Y, N).  Note that this can never
      35                 :            :  * split prime powers.  (Any odd prime dividing the X - Y factor, say, will
      36                 :            :  * divide it to the same power as it divides N.)
      37                 :            :  *
      38                 :            :  * Candidates x0 are found by sieving along arithmetic progressions modulo
      39                 :            :  * the small primes in the factor base, and evaluation of candidates picks
      40                 :            :  * out those x0 where many of these APs happen to coincide, resulting in a
      41                 :            :  * highly divisible Q0(x0).
      42                 :            :  *
      43                 :            :  * The Multi-Polynomial version improves this by choosing a modest subset of
      44                 :            :  * factor base primes and forcing these to divide Q0(x).  Let A be the product
      45                 :            :  * of the chosen primes, and write Q(x) = Q0(2Ax + B) = (2Ax + B)^2 - kN =
      46                 :            :  * 4A(Ax^2 + Bx + C), where B has been suitably chosen.  For each A, there
      47                 :            :  * are 2^omega_A possible values for B when A is the product of omega_A
      48                 :            :  * distinct primes, but we'll use only half of these, since the other half
      49                 :            :  * is more easily covered by exploiting the symmetry x <-> -x of the original
      50                 :            :  * quadratic.  The "Self-Initializating" bit refers to the fact that switching
      51                 :            :  * from one B to the next can be done very fast, whereas switching to the
      52                 :            :  * next A involves some recomputation from scratch.  (C is never needed
      53                 :            :  * explicitly except for debug diagnostics.)  Thus we can very quickly run
      54                 :            :  * through an entire "cohort" of polynomials sharing the same A.
      55                 :            :  *
      56                 :            :  * The sieve now ranges over values x0 such that |x0| < M  (we use x = x0 + M
      57                 :            :  * as the non-negative array subscript).  The coefficients A are chosen so
      58                 :            :  * that A*M is approximately sqrt(kN).  Then |B| will be bounded by about
      59                 :            :  * (j+4)*A, and |C| = -C will be about (M/4)*sqrt(kN), so Q(x0)/(4A) will
      60                 :            :  * take values roughly between -|C| and 3|C|.
      61                 :            :  *
      62                 :            :  * There are numerous refinements to this basic outline  (e.g. it is more
      63                 :            :  * efficient to _not_ use the very smallest primes in the FB for sieving,
      64                 :            :  * incorporating them only after selecting candidates).  The substition of
      65                 :            :  * 2Ax+B into X^2 - kN, with odd B, forces 2 to occur;  when kN is 1 mod 8,
      66                 :            :  * it will always occur at least to the 3rd power;  when kN = 5 mod 8, it
      67                 :            :  * will always occur exactly to the 2nd power.  We never sieve on 2 and we
      68                 :            :  * always pull out the power of 2 directly  (which is easy, of course).
      69                 :            :  * The prime factor(s) of k will show up whenever 2Ax + B has a factor in
      70                 :            :  * common with k;  we don't sieve on these either but can easily recognize
      71                 :            :  * them in a candidate.
      72                 :            :  */
      73                 :            : 
      74                 :            : #include "pari.h"
      75                 :            : #include "paripriv.h"
      76                 :            : 
      77                 :            : static char *
      78                 :    5845638 : paristrtok_r(char *str, const char *delim, char **saveptr)
      79                 :            : {
      80                 :            :   char *res;
      81         [ +  + ]:    5845638 :   if (!str) str = *saveptr;
      82                 :    5845638 :   str += strspn(str, delim);
      83         [ -  + ]:    5845638 :   if (!*str) return NULL;
      84                 :    5845638 :   res = str;
      85                 :    5845638 :   str += strcspn(str, delim);
      86         [ +  + ]:    5845638 :   if (*str) *str++ = 0;
      87                 :    5845638 :   *saveptr = str;
      88                 :    5845638 :   return res;
      89                 :            : }
      90                 :            : 
      91                 :            : /** DEBUG **/
      92                 :            : /* #define MPQS_DEBUG_VERBOSE 1 */
      93                 :            : /* histograms are pretty, but don't help performance after all (see below) */
      94                 :            : /* #define MPQS_USE_HISTOGRAMS */
      95                 :            : 
      96                 :            : #include "mpqs.h"
      97                 :            : 
      98                 :            : /*********************************************************************/
      99                 :            : /**                                                                 **/
     100                 :            : /**                         INITIAL SIZING                          **/
     101                 :            : /**                                                                 **/
     102                 :            : /*********************************************************************/
     103                 :            : /* number of decimal digits of argument - for parameter choosing and for
     104                 :            :  * diagnostics */
     105                 :            : static long
     106                 :        770 : decimal_len(GEN N)
     107                 :            : {
     108                 :        770 :   pari_sp av = avma;
     109                 :        770 :   long d = strlen(itostr(N));
     110                 :        770 :   avma = av; return d;
     111                 :            : }
     112                 :            : 
     113                 :            : /* To be called after choosing k and putting kN into the handle:
     114                 :            :  * Pick up the requested parameter set for the given size of kN in decimal
     115                 :            :  * digits and fill in various fields in the handle.  Return 0 when kN is
     116                 :            :  * too large, 1 when we're ok. */
     117                 :            : static int
     118                 :        385 : mpqs_set_parameters(mpqs_handle_t *h)
     119                 :            : {
     120                 :            :   long i;
     121                 :            :   const mpqs_parameterset_t *P;
     122                 :            :   double mb;
     123                 :            : 
     124                 :        385 :   h->digit_size_kN = decimal_len(h->kN);
     125         [ -  + ]:        385 :   if (h->digit_size_kN <= 9)
     126                 :          0 :     i = 0;
     127         [ -  + ]:        385 :   else if (h->digit_size_kN > MPQS_MAX_DIGIT_SIZE_KN)
     128                 :          0 :     return 0;
     129                 :            :   else
     130                 :        385 :     i = h->digit_size_kN - 9;
     131                 :            : 
     132                 :            :   /* (cf PARI bug#235) the following has always been, and will remain,
     133                 :            :    * a moving target... increased thresholds from 64, 80 to 79, 86
     134                 :            :    * respectively --GN20050601.  Note that the new values correspond to
     135                 :            :    * kN having >= 86 or >= 95 decimal digits, respectively.  Note also
     136                 :            :    * that the current sizing parameters for 90 or more digits are based
     137                 :            :    * on 100% theory and 0% practice. */
     138         [ -  + ]:        385 :   if (i >= 79)
     139         [ #  # ]:          0 :     pari_warn(warner, "MPQS: factoring this number will take %s hours:\nN = %Ps",
     140                 :            :         i >= 86 ? "many": "several", h->N);
     141                 :            : 
     142         [ -  + ]:        385 :   if (DEBUGLEVEL >= 5)
     143                 :            :   {
     144                 :          0 :     err_printf("MPQS: kN = %Ps\n", h->kN);
     145                 :          0 :     err_printf("MPQS: kN has %ld decimal digits\n", h->digit_size_kN);
     146                 :            :   }
     147                 :            : 
     148                 :        385 :   P = &(mpqs_parameters[i]);
     149                 :        385 :   h->tolerance        = P->tolerance;
     150                 :        385 :   h->lp_scale         = P->lp_scale;
     151                 :            :   /* make room for prime factors of k if any: */
     152                 :        385 :   h->size_of_FB       = P->size_of_FB + h->_k->omega_k;
     153                 :            :   /* for the purpose of Gauss elimination etc., prime factors of k behave
     154                 :            :    * like real FB primes, so take them into account when setting the goal: */
     155                 :        770 :   h->target_no_rels   = (h->size_of_FB >= 200 ?
     156         [ +  + ]:        385 :                          h->size_of_FB + 70 :
     157                 :        185 :                          (mpqs_int32_t)(h->size_of_FB * 1.35));
     158                 :        385 :   h->M                = P->M;
     159                 :        385 :   h->omega_A          = P->omega_A;
     160                 :        385 :   h->no_B             = 1UL << (P->omega_A - 1);
     161                 :        385 :   h->pmin_index1      = P->pmin_index1;
     162                 :            :   /* certain subscripts into h->FB should also be offset by omega_k: */
     163                 :        385 :   h->index0_FB        = 3 + h->_k->omega_k;
     164                 :            :   /* following are converted from % to parts per thousand: */
     165                 :        385 :   h->first_sort_point = 10 * P->first_sort_point;
     166                 :        385 :   h->sort_pt_interval = 10 * P->sort_pt_interval;
     167                 :            : 
     168                 :        385 :   mb = (h->size_of_FB + 1)/(8.*1048576.) * h->target_no_rels;
     169         [ -  + ]:        385 :   if (mb > 128.)
     170                 :            :   {
     171                 :          0 :     pari_warn(warner,
     172                 :            :         "MPQS: Gauss elimination will require more than\n\t128MBy of memory");
     173         [ #  # ]:          0 :     if (DEBUGLEVEL >= 1)
     174                 :          0 :       err_printf("\t(estimated memory needed: %4.1fMBy)\n", mb);
     175                 :            :   }
     176                 :            : 
     177                 :        385 :   return 1;
     178                 :            : }
     179                 :            : 
     180                 :            : /*********************************************************************/
     181                 :            : /**                                                                 **/
     182                 :            : /**                       OBJECT HOUSEKEEPING                       **/
     183                 :            : /**                                                                 **/
     184                 :            : /*********************************************************************/
     185                 :            : 
     186                 :            : /* The sub-constructors for the pieces of the handle will be called in the
     187                 :            :  * same order as their appearance here, and the later ones in part rely on
     188                 :            :  * the earlier ones having filled in some fields.
     189                 :            :  * There's a single destructor to handle all cleanup at the end  (except
     190                 :            :  * for mpqs() itself resetting avma). */
     191                 :            : 
     192                 :            : /* main handle constructor */
     193                 :            : static mpqs_handle_t *
     194                 :        385 : mpqs_handle_ctor(GEN N)
     195                 :            : {
     196                 :        385 :   mpqs_handle_t *h = (mpqs_handle_t *) pari_calloc(sizeof(mpqs_handle_t));
     197                 :        385 :   h->N = N;
     198                 :            : #ifdef MPQS_DEBUG_VERBOSE
     199                 :            :   err_printf("MPQS DEBUG: created handle @0x%p\n", (void *)h);
     200                 :            : #endif
     201                 :        385 :   return h;
     202                 :            : }
     203                 :            : 
     204                 :            : /* factor base constructor. Really a home-grown memalign(3c) underneath.
     205                 :            :  * We don't want FB entries to straddle L1 cache line boundaries, and
     206                 :            :  * malloc(3c) only guarantees alignment adequate for all primitive data
     207                 :            :  * types of the platform ABI - typically to 8 or 16 byte boundaries.
     208                 :            :  * Also allocate the inv_A_H array.
     209                 :            :  * The FB array pointer is returned for convenience */
     210                 :            : static mpqs_FB_entry_t *
     211                 :        385 : mpqs_FB_ctor(mpqs_handle_t *h)
     212                 :            : {
     213                 :            :   /* leave room for slots 0, 1, and sentinel slot at the end of the array */
     214                 :        385 :   long size_FB_chunk = (h->size_of_FB + 3) * sizeof(mpqs_FB_entry_t);
     215                 :            :   /* like FB, except this one does not have a sentinel slot at the end */
     216                 :        385 :   long size_IAH_chunk = (h->size_of_FB + 2) * sizeof(mpqs_inv_A_H_t);
     217                 :        385 :   char *fbp = (char*)pari_malloc(size_FB_chunk + 64);
     218                 :        385 :   char *iahp = (char*)pari_malloc(size_IAH_chunk + 64);
     219                 :            :   long fbl, iahl;
     220                 :            : 
     221                 :        385 :   h->FB_chunk = (void *)fbp;
     222                 :        385 :   h->invAH_chunk = (void *)iahp;
     223                 :            :   /* round up to next higher 64-bytes-aligned address */
     224                 :        385 :   fbl = (((long)fbp) + 64) & ~0x3FL;
     225                 :            :   /* and put the actual array there */
     226                 :        385 :   h->FB = (mpqs_FB_entry_t *)fbl;
     227                 :            : 
     228                 :        385 :   iahl = (((long)iahp) + 64) & ~0x3FL;
     229                 :        385 :   h->inv_A_H = (mpqs_inv_A_H_t *)iahl;
     230                 :        385 :   return (mpqs_FB_entry_t *)fbl;
     231                 :            : }
     232                 :            : 
     233                 :            : /* sieve array constructor;  also allocates the candidates array, the
     234                 :            :  * histograms, and temporary storage for relations under construction */
     235                 :            : static void
     236                 :        385 : mpqs_sieve_array_ctor(mpqs_handle_t *h)
     237                 :            : {
     238                 :        385 :   long size = (h->M << 1) + 1;
     239                 :        385 :   mpqs_int32_t size_of_FB = h->size_of_FB;
     240                 :            : 
     241                 :        385 :   h->sieve_array = (unsigned char *) pari_malloc(size * sizeof(unsigned char));
     242                 :        385 :   h->sieve_array_end = h->sieve_array + size - 2;
     243                 :        385 :   h->sieve_array_end[1] = 255; /* sentinel */
     244                 :        385 :   h->candidates = (long *)pari_malloc(MPQS_CANDIDATE_ARRAY_SIZE * sizeof(long));
     245                 :            : 
     246                 :            :   /* Room needed for string representation of a relation - worst case:
     247                 :            :    * + leading " 1 1"
     248                 :            :    * + trailing " 0\n" with final NUL character
     249                 :            :    * + in between up to size_of_FB pairs each consisting of an exponent, a
     250                 :            :    *   subscript into FB, and two spaces.
     251                 :            :    * Subscripts into FB fit into 5 digits, and exponents fit into 3 digits
     252                 :            :    * with room to spare -- anything needing 3 or more digits for the
     253                 :            :    * subscript must come with an exponent of at most 2 digits. Moreover the
     254                 :            :    * product of the first 58 primes is larger than 10^110  (and the righthand
     255                 :            :    * sides of proto-relations are much smaller than kN: on the order of
     256                 :            :    * M*sqrt(kN)),  so there cannot be more than 60 pairs in all, even if
     257                 :            :    * size_of_FB > 10^5. --GN */
     258                 :            : 
     259                 :            :   /* whereas mpqs_self_init() uses size_of_FB+1, we just use the size as
     260                 :            :    * it is, not counting FB[1], to start off the following estimate */
     261         [ +  + ]:        385 :   if (size_of_FB > 60) size_of_FB = 60;
     262                 :        385 :   h->relations = (char *) pari_malloc((8 + size_of_FB * 9) * sizeof(char));
     263                 :            :   /* and for tracking which primes occur in the current relation: */
     264                 :        385 :   h->relaprimes = (long *) pari_malloc((size_of_FB << 1) * sizeof(long));
     265                 :            : 
     266                 :            : #ifdef MPQS_USE_HISTOGRAMS
     267                 :            :   /* histograms to be used only when kN isn't very small */
     268                 :            :   if (h->size_of_FB > MPQS_MIN_SIZE_FB_FOR_HISTO) {
     269                 :            :     h->do_histograms = 1;
     270                 :            :     h->histo_full = (long *) pari_calloc(128 * sizeof(long));
     271                 :            :     h->histo_lprl = (long *) pari_calloc(128 * sizeof(long));
     272                 :            :     h->histo_drop = (long *) pari_calloc(128 * sizeof(long));
     273                 :            :   }
     274                 :            : #endif
     275                 :        385 : }
     276                 :            : 
     277                 :            : /* mpqs() calls the following (after recording avma) to allocate GENs for
     278                 :            :  * the current polynomial and self-initialization scratchpad data on the
     279                 :            :  * PARI stack.  This space is released by mpqs() itself at the end. */
     280                 :            : static void
     281                 :        385 : mpqs_poly_ctor(mpqs_handle_t *h)
     282                 :            : {
     283                 :            :   mpqs_int32_t i;
     284                 :        385 :   long size_per = h->omega_A * sizeof(mpqs_per_A_prime_t);
     285                 :            : 
     286                 :        385 :   h->per_A_pr = (mpqs_per_A_prime_t *) pari_calloc(size_per);
     287                 :            :   /* Sizing:  A is the product of omega_A primes, each well below word
     288                 :            :    * size.
     289                 :            :    * |B| is bounded by (omega_A + 4) * A, so can have at most one word
     290                 :            :    * more, and that's generous.
     291                 :            :    * |C| is less than A*M^2, so can take at most two words more than A.
     292                 :            :    * The array H holds residues modulo A, so the same size as used for A
     293                 :            :    * is sufficient. */
     294                 :        385 :   h->A = cgeti(h->omega_A + 2);
     295                 :        385 :   h->B = cgeti(h->omega_A + 3);
     296                 :            : #ifdef MPQS_DEBUG
     297                 :            :   h->C = cgeti(h->omega_A + 4);
     298                 :            : #endif
     299         [ +  + ]:       1860 :   for (i = 0; i < h->omega_A; i++)
     300                 :       1475 :     h->per_A_pr[i]._H = cgeti(h->omega_A + 2);
     301                 :            :   /* the handle starts out all zero, so in particular bin_index and index_i
     302                 :            :    * are initially 0.
     303                 :            :    * TODO: index_j currently initialized in mqps() but this is going to
     304                 :            :    * change. */
     305                 :        385 : }
     306                 :            : 
     307                 :            : /* main handle destructor, also cleans up all other allocated pieces
     308                 :            :  * (except for stuff created on the PARI stack which the caller should
     309                 :            :  * deal with by resetting avma) */
     310                 :            : static void
     311                 :        385 : mpqs_handle_dtor(mpqs_handle_t *h)
     312                 :            : {
     313                 :            : #define myfree(x) if(x) pari_free((void*)x)
     314         [ +  - ]:        385 :   myfree((h->per_A_pr));
     315         [ +  - ]:        385 :   myfree((h->relaprimes));
     316         [ +  - ]:        385 :   myfree(h->relations);
     317                 :            : 
     318                 :            : #ifdef MPQS_USE_HISTOGRAMS
     319                 :            :   myfree((h->histo_drop));
     320                 :            :   myfree((h->histo_lprl));
     321                 :            :   myfree((h->histo_full));
     322                 :            : #endif
     323                 :            : 
     324         [ +  - ]:        385 :   myfree((h->candidates));
     325         [ +  - ]:        385 :   myfree((h->sieve_array));
     326         [ +  - ]:        385 :   myfree((h->invAH_chunk));
     327         [ +  - ]:        385 :   myfree((h->FB_chunk));
     328         [ +  - ]:        385 :   myfree(h);
     329                 :        385 : }
     330                 :            : 
     331                 :            : /* TODO: relationsdb handle */
     332                 :            : 
     333                 :            : /*********************************************************************/
     334                 :            : /**                                                                 **/
     335                 :            : /**                        FACTOR BASE SETUP                        **/
     336                 :            : /**                                                                 **/
     337                 :            : /*********************************************************************/
     338                 :            : /* fill in the best-guess multiplier k for N. We force kN = 1 mod 4.
     339                 :            :  * Caller should proceed to fill in kN */
     340                 :            : static ulong
     341                 :        385 : mpqs_find_k(mpqs_handle_t *h)
     342                 :            : {
     343                 :        385 :   const pari_sp av = avma;
     344                 :        385 :   const long N_mod_8 = mod8(h->N), N_mod_4 = N_mod_8 & 3;
     345                 :            :   forprime_t S;
     346                 :            :   struct {
     347                 :            :     const mpqs_multiplier_t *_k;
     348                 :            :     long np; /* number of primes in factorbase so far for this k */
     349                 :            :     double value; /* the larger, the better */
     350                 :            :   } cache[MPQS_POSSIBLE_MULTIPLIERS];
     351                 :        385 :   ulong p, i, j, nbk = MPQS_POSSIBLE_MULTIPLIERS;
     352                 :            : 
     353         [ +  - ]:       3705 :   for (i = j = 0; i < sizeof(cand_multipliers)/sizeof(mpqs_multiplier_t); i++)
     354                 :            :   {
     355                 :       3705 :     const mpqs_multiplier_t *cand_k = &cand_multipliers[i];
     356                 :       3705 :     long k = cand_k->k;
     357                 :            :     double v;
     358         [ +  + ]:       3705 :     if ((k & 3) != N_mod_4) continue; /* want kN = 1 (mod 4) */
     359                 :       1925 :     v = -0.35 * log2((double)k);
     360         [ +  + ]:       1925 :     if ((k & 7) == N_mod_8) v += LOG2; /* kN = 1 (mod 8) */
     361                 :       1925 :     cache[j].np = 0;
     362                 :       1925 :     cache[j]._k = cand_k;
     363                 :       1925 :     cache[j].value = v;
     364         [ +  + ]:       1925 :     if (++j == nbk) break; /* enough */
     365                 :            :   }
     366                 :        385 :   nbk = j;
     367                 :        385 :   u_forprime_init(&S, 2, ULONG_MAX);
     368         [ +  - ]:       7135 :   while ( (p = u_forprime_next(&S)) )
     369                 :            :   {
     370                 :       7135 :     ulong Np = umodiu(h->N, p);
     371                 :       7135 :     long kroNp, seen = 0;
     372         [ -  + ]:       7135 :     if (!Np) return p;
     373                 :       7135 :     kroNp = krouu(Np, p);
     374         [ +  + ]:      42810 :     for (i = 0; i < nbk; i++)
     375                 :            :     {
     376         [ +  + ]:      35675 :       if (cache[i].np > MPQS_MULTIPLIER_SEARCH_DEPTH) continue;
     377                 :      24020 :       seen++;
     378         [ +  + ]:      24020 :       if (krouu(cache[i]._k->k % p, p) == kroNp) /* kronecker(k*N, p)=1 */
     379                 :            :       {
     380                 :      11550 :         cache[i].value += log2((double) p)/p;
     381                 :      11550 :         cache[i].np++;
     382                 :            :       }
     383                 :            :     }
     384         [ +  + ]:       7135 :     if (!seen) break; /* we're gone through SEARCH_DEPTH primes for all k */
     385                 :            :   }
     386         [ -  + ]:        385 :   if (!p) pari_err_OVERFLOW("mpqs_find_k [ran out of primes]");
     387                 :            :   {
     388                 :        385 :     long best_i = 0;
     389                 :        385 :     double v = cache[0].value;
     390         [ +  + ]:       1925 :     for (i = 1; i < nbk; i++)
     391         [ +  + ]:       1540 :       if (cache[i].value > v) { best_i = i; v = cache[i].value; }
     392                 :        385 :     h->_k = cache[best_i]._k; avma = av; return 0;
     393                 :            :   }
     394                 :            : }
     395                 :            : 
     396                 :            : /******************************/
     397                 :            : /* Create a factor base of 'size' primes p_i such that legendre(k*N, p_i) != -1
     398                 :            :  * We could have shifted subscripts down from their historical arrangement,
     399                 :            :  * but this seems too risky for the tiny potential gain in memory economy.
     400                 :            :  * The real constraint is that the subscripts of anything which later shows
     401                 :            :  * up at the Gauss stage must be nonnegative, because the exponent vectors
     402                 :            :  * there use the same subscripts to refer to the same FB entries.  Thus in
     403                 :            :  * particular, the entry representing -1 could be put into FB[0], but could
     404                 :            :  * not be moved to FB[-1] (although mpqs_FB_ctor() could be easily adapted
     405                 :            :  * to support negative subscripts).-- The historically grown layout is:
     406                 :            :  * FB[0] is unused.
     407                 :            :  * FB[1] is not explicitly used but stands for -1.
     408                 :            :  * FB[2] contains 2 (always).
     409                 :            :  * Before we are called, the size_of_FB field in the handle will already have
     410                 :            :  * been adjusted by _k->omega_k, so there's room for the primes dividing k,
     411                 :            :  * which when present will occupy FB[3] and following.
     412                 :            :  * The "real" odd FB primes begin at FB[h->index0_FB].
     413                 :            :  * FB[size_of_FB+1] is the last prime p_i.
     414                 :            :  * FB[size_of_FB+2] is a sentinel to simplify some of our loops.
     415                 :            :  * Thus we allocate size_of_FB+3 slots for FB.
     416                 :            :  *
     417                 :            :  * If a prime factor of N is found during the construction, it is returned
     418                 :            :  * in f, otherwise f = 0. */
     419                 :            : 
     420                 :            : /* returns the FB array pointer for convenience */
     421                 :            : static mpqs_FB_entry_t *
     422                 :        385 : mpqs_create_FB(mpqs_handle_t *h, ulong *f)
     423                 :            : {
     424                 :        385 :   const pari_sp av = avma;
     425                 :        385 :   mpqs_int32_t size = h->size_of_FB;
     426                 :            :   long i;
     427                 :        385 :   mpqs_uint32_t k = h->_k->k;
     428                 :            :   mpqs_FB_entry_t *FB;
     429                 :            :   forprime_t S;
     430                 :            : 
     431                 :        385 :   FB = mpqs_FB_ctor(h);
     432                 :            : 
     433                 :        385 :   FB[2].fbe_p = 2;
     434                 :            :   /* the fbe_logval and the fbe_sqrt_kN for 2 are never used */
     435                 :        385 :   FB[2].fbe_flags = MPQS_FBE_CLEAR;
     436                 :        385 :   (void)u_forprime_init(&S, 3, ULONG_MAX);
     437                 :            : 
     438                 :            :   /* the first loop executes h->_k->omega_k = 0, 1, or 2 times */
     439         [ +  + ]:        565 :   for (i = 3; i < h->index0_FB; i++)
     440                 :            :   {
     441                 :        180 :     mpqs_uint32_t kp = (ulong)h->_k->kp[i-3];
     442         [ -  + ]:        180 :     if (MPQS_DEBUGLEVEL >= 7) err_printf(",<%lu>", (ulong)kp);
     443                 :        180 :     FB[i].fbe_p = kp;
     444                 :            :     /* we *could* flag divisors of k here, but so far I see no need,
     445                 :            :      * and no flags bit has been assigned for the purpose */
     446                 :        180 :     FB[i].fbe_flags = MPQS_FBE_CLEAR;
     447                 :        180 :     FB[i].fbe_flogp = (float) log2((double) kp);
     448                 :        180 :     FB[i].fbe_sqrt_kN = 0;
     449                 :            :   }
     450                 :            :   /* now i == h->index0_FB */
     451         [ +  + ]:     202520 :   while (i < size + 2)
     452                 :            :   {
     453                 :     202135 :     ulong p = u_forprime_next(&S);
     454 [ +  + ][ +  + ]:     202135 :     if (p > k || k % p)
     455                 :            :     {
     456                 :     201955 :       ulong kN_mod_p = umodiu(h->kN, p);
     457                 :     201955 :       long kr = krouu(kN_mod_p, p);
     458         [ +  + ]:     201955 :       if (kr != -1)
     459                 :            :       {
     460         [ -  + ]:     101740 :         if (kr == 0) { *f = p; return FB; }
     461                 :     101740 :         FB[i].fbe_p = (mpqs_uint32_t) p;
     462                 :     101740 :         FB[i].fbe_flags = MPQS_FBE_CLEAR;
     463                 :            :         /* dyadic logarithm of p; single precision suffices */
     464                 :     101740 :         FB[i].fbe_flogp = (float) log2((double)p);
     465                 :            :         /* cannot yet fill in fbe_logval because the scaling multiplier
     466                 :            :          * depends on the largest prime in FB, as yet unknown */
     467                 :            : 
     468                 :            :         /* x such that x^2 = kN (mod p_i) */
     469                 :     101740 :         FB[i++].fbe_sqrt_kN = (mpqs_uint32_t)Fl_sqrt(kN_mod_p, p);
     470                 :            :       }
     471                 :            :     }
     472                 :            :   }
     473                 :        385 :   avma = av;
     474         [ -  + ]:        385 :   if (MPQS_DEBUGLEVEL >= 7)
     475                 :            :   {
     476                 :          0 :     err_printf("MPQS: FB [-1,2");
     477         [ #  # ]:          0 :     for (i = 3; i < h->index0_FB; i++) err_printf(",<%lu>", FB[i].fbe_p);
     478         [ #  # ]:          0 :     for (; i < size + 2; i++) err_printf(",%lu", FB[i].fbe_p);
     479                 :          0 :     err_printf("]\n");
     480                 :            :   }
     481                 :            : 
     482                 :        385 :   FB[i].fbe_p = 0;              /* sentinel */
     483                 :        385 :   h->largest_FB_p = FB[i-1].fbe_p; /* at subscript size_of_FB + 1 */
     484                 :            : 
     485                 :            :   /* locate the smallest prime that will be used for sieving */
     486         [ +  - ]:        905 :   for (i = h->index0_FB; FB[i].fbe_p != 0; i++)
     487         [ +  + ]:        905 :     if (FB[i].fbe_p >= h->pmin_index1) break;
     488                 :        385 :   h->index1_FB = i;
     489                 :            :   /* with our parameters this will never fall of the end of the FB */
     490                 :        385 :   *f = 0; return FB;
     491                 :            : }
     492                 :            : 
     493                 :            : /*********************************************************************/
     494                 :            : /**                                                                 **/
     495                 :            : /**                      MISC HELPER FUNCTIONS                      **/
     496                 :            : /**                                                                 **/
     497                 :            : /*********************************************************************/
     498                 :            : 
     499                 :            : /* Effect of the following:  multiplying the base-2 logarithm of some
     500                 :            :  * quantity by log_multiplier will rescale something of size
     501                 :            :  *    log2 ( sqrt(kN) * M / (largest_FB_prime)^tolerance )
     502                 :            :  * to 232.  Note that sqrt(kN) * M is just A*M^2, the value our polynomials
     503                 :            :  * take at the outer edges of the sieve interval.  The scale here leaves
     504                 :            :  * a little wiggle room for accumulated rounding errors from the approximate
     505                 :            :  * byte-sized scaled logarithms for the factor base primes which we add up
     506                 :            :  * in the sieving phase.-- The threshold is then chosen so that a point in
     507                 :            :  * the sieve has to reach a result which, under the same scaling, represents
     508                 :            :  *    log2 ( sqrt(kN) * M / (largest_FB_prime)^tolerance )
     509                 :            :  * in order to be accepted as a candidate. */
     510                 :            : /* The old formula was...
     511                 :            :  *   log_multiplier =
     512                 :            :  *      127.0 / (0.5 * log2 (handle->dkN)
     513                 :            :  *               + log2((double)M)
     514                 :            :  *               - tolerance * log2((double)handle->largest_FB_p)
     515                 :            :  *               );
     516                 :            :  * and we used to use this with a constant threshold of 128. */
     517                 :            : 
     518                 :            : /* NOTE: We used to divide log_multiplier by an extra factor 2, and in
     519                 :            :  * compensation we were multiplying by 2 when the fbe_logp fields were being
     520                 :            :  * filled in, making all those bytes even.  Tradeoff: the extra bit of
     521                 :            :  * precision is helpful, but interferes with a possible sieving optimization
     522                 :            :  * (artifically shift right the logp's of primes in A, and just run over both
     523                 :            :  * arithmetical progressions  (which coincide in this case)  instead of
     524                 :            :  * skipping the second one, to avoid the conditional branch in the
     525                 :            :  * mpqs_sieve() loops).  We could still do this, but might lose a little bit
     526                 :            :  * accuracy for those primes.  Probably no big deal. */
     527                 :            : static void
     528                 :        385 : mpqs_set_sieve_threshold(mpqs_handle_t *h)
     529                 :            : {
     530                 :        385 :   mpqs_FB_entry_t *FB = h->FB;
     531                 :            :   long i;
     532                 :            :   double log_maxval;
     533                 :            :   double log_multiplier;
     534                 :            : 
     535                 :        385 :   h->l2sqrtkN = 0.5 * log2(h->dkN);
     536                 :        385 :   h->l2M = log2((double)h->M);
     537                 :        385 :   log_maxval = h->l2sqrtkN + h->l2M - MPQS_A_FUDGE;
     538                 :        385 :   log_multiplier = 232.0 / log_maxval;
     539                 :        385 :   h->sieve_threshold =
     540                 :        385 :     (unsigned char) (log_multiplier *
     541                 :            :                      (log_maxval
     542                 :        385 :                       - h->tolerance * log2((double)h->largest_FB_p)
     543                 :            :                       )
     544                 :            :                      ) + 1;
     545                 :            :   /* That "+ 1" really helps - we may want to tune towards somewhat smaller
     546                 :            :    * tolerances  (or introduce self-tuning one day)... */
     547                 :            : 
     548                 :            :   /* If this turns out to be <128, scream loudly.
     549                 :            :    * That means that the FB or the tolerance or both are way too
     550                 :            :    * large for the size of kN.  (Normally, the threshold should end
     551                 :            :    * up in the 150...170 range.) */
     552         [ -  + ]:        385 :   if (h->sieve_threshold < 128) {
     553                 :          0 :     h->sieve_threshold = 128;
     554                 :          0 :     pari_warn(warner,
     555                 :            :         "MPQS: sizing out of tune, FB size or tolerance\n\ttoo large");
     556                 :            :   }
     557                 :            : 
     558                 :            :   /* Now fill in the byte-sized approximate scaled logarithms of p_i */
     559         [ -  + ]:        385 :   if (DEBUGLEVEL >= 5)
     560                 :            :   {
     561                 :          0 :     err_printf("MPQS: computing logarithm approximations for p_i in FB\n");
     562                 :            :   }
     563         [ +  + ]:     102125 :   for (i = h->index0_FB; i < h->size_of_FB + 2; i++)
     564                 :            :   {
     565                 :     203480 :     FB[i].fbe_logval =
     566                 :     101740 :       (unsigned char) (log_multiplier * FB[i].fbe_flogp);
     567                 :            :   }
     568                 :        385 : }
     569                 :            : 
     570                 :            : /* Given the partially populated handle, find the optimum place in the FB
     571                 :            :  * to pick prime factors for A from.  The lowest admissible subscript is
     572                 :            :  * index0_FB, but unless kN is very small, we stay away a bit from that.
     573                 :            :  * The highest admissible is size_of_FB + 1, where the largest FB prime
     574                 :            :  * resides.  The ideal corner is about (sqrt(kN)/M) ^ (1/omega_A),
     575                 :            :  * so that A will end up of size comparable to sqrt(kN)/M;  experimentally
     576                 :            :  * it seems desirable to stay slightly below this.  Moreover, the selection
     577                 :            :  * of the individual primes happens to err on the large side, for which we
     578                 :            :  * compensate a bit, using the (small positive) quantity MPQS_A_FUDGE.
     579                 :            :  * We rely on a few auxiliary fields in the handle to be already set by
     580                 :            :  * mqps_set_sieve_threshold() before we are called.
     581                 :            :  * Return 1 on success, and 0 otherwise. */
     582                 :            : static int
     583                 :        385 : mpqs_locate_A_range(mpqs_handle_t *h)
     584                 :            : {
     585                 :            :   /* i will be counted up to the desirable index2_FB + 1, and omega_A is never
     586                 :            :    * less than 3, and we want
     587                 :            :    *   index2_FB - (omega_A - 1) + 1 >= index0_FB + omega_A - 3,
     588                 :            :    * so: */
     589                 :        385 :   long i = h->index0_FB + 2*(h->omega_A) - 4;
     590                 :            :   double l2_target_pA;
     591                 :        385 :   mpqs_FB_entry_t *FB = h->FB;
     592                 :            : 
     593                 :        385 :   h->l2_target_A = (h->l2sqrtkN - h->l2M - MPQS_A_FUDGE);
     594                 :        385 :   l2_target_pA = h->l2_target_A / h->omega_A;
     595                 :            : 
     596                 :            :   /* find the sweet spot, normally shouldn't take long */
     597 [ +  - ][ +  + ]:      10435 :   while ((FB[i].fbe_p != 0) && (FB[i].fbe_flogp <= l2_target_pA)) i++;
     598                 :            : 
     599                 :            : #ifdef MPQS_DEBUG_LOCATE_A_RANGE
     600                 :            :   err_printf("MPQS DEBUG: omega_A=%ld, index0=%ld, i=%ld\n",
     601                 :            :              (long) h->omega_A, (long) h->index0_FB, i);
     602                 :            : #endif
     603                 :            : 
     604                 :            :   /* check whether this hasn't walked off the top end... */
     605                 :            :   /* The following should actually NEVER happen. */
     606         [ -  + ]:        385 :   if (i > h->size_of_FB - 3)
     607                 :            :   { /* this isn't going to work at all. */
     608                 :          0 :     pari_warn(warner,
     609                 :            :         "MPQS: sizing out of tune, FB too small or\n\tway too few primes in A");
     610                 :          0 :     return 0;
     611                 :            :   }
     612                 :        385 :   h->index2_FB = i - 1;
     613                 :            : #ifdef MPQS_DEBUG_LOCATE_A_RANGE
     614                 :            :   err_printf("MPQS DEBUG: index2_FB = %ld\n", i - 1);
     615                 :            : #endif
     616                 :            :   /* GN20050723
     617                 :            :    * assert: index0_FB + (omega_A - 3) [the lowest FB subscript eligible to
     618                 :            :    * be used in picking primes for A]  plus  (omega_A - 2)  does not exceed
     619                 :            :    * index2_FB  [the subscript from which the choice of primes for A starts,
     620                 :            :    * putting omega_A - 1 of them at or below index2_FB, and the last and
     621                 :            :    * largest one above, cf. mpqs_si_choose_primes() below].
     622                 :            :    * Moreover, index2_FB indicates the last prime below the ideal size, unless
     623                 :            :    * (when kN is very small) the ideal size was too small to use. */
     624                 :            : 
     625                 :        385 :   return 1;
     626                 :            : }
     627                 :            : 
     628                 :            : 
     629                 :            : /*********************************************************************/
     630                 :            : /**                                                                 **/
     631                 :            : /**                HISTOGRAMS AND THRESHOLD FEEDBACK                **/
     632                 :            : /**                                                                 **/
     633                 :            : /*********************************************************************/
     634                 :            : 
     635                 :            : #ifdef MPQS_USE_HISTOGRAMS
     636                 :            : 
     637                 :            : /* The histogram-related code is left in this file, but all under the
     638                 :            :  * above #ifdef, and disabled by default.  I'm finding that:
     639                 :            :  * - merely keeping the numbers updated in mpqs_eval_cand() below  (and
     640                 :            :  *   keeping the "negligible" 1.5 or 3KBys' worth of extra arrays in use)
     641                 :            :  *   causes us to run quite noticeably slower: 8-10% for a 73-digit number,
     642                 :            :  * - mpqs_eval_cand() has already become so much faster than it used to be
     643                 :            :  *   that raising the threshold to get rid of many low-valued unpromising
     644                 :            :  *   candidates does not save any significant time, and even losing a pretty
     645                 :            :  *   small number of additional LP relations actually harms us by lowering
     646                 :            :  *   the efficiency of LP relation combining.
     647                 :            :  * (The first point might be due merely to code bloat and less effective
     648                 :            :  * compiler optimizations - I'm not sure about that.)
     649                 :            :  * Just for getting a visual impression of how the sieve is performing,
     650                 :            :  * however, this is nice to have available.  Just turn on the #define at
     651                 :            :  * the top of the file and recompile with it. --GN2005-02-03
     652                 :            :  */
     653                 :            : 
     654                 :            : /* histogram evaluation happens very infrequently if at all.  So we'll do
     655                 :            :  * all the adding and putting-into-relation-with here, while mpqs_eval_cand()
     656                 :            :  * merely bumps one cell at a time and doesn't keep running totals. */
     657                 :            : 
     658                 :            : static void
     659                 :            : mpqs_print_histo(mpqs_handle_t *h)
     660                 :            : {
     661                 :            :   long i, tot = 0;
     662                 :            : 
     663                 :            :   if (!h->do_histograms) return;
     664                 :            : 
     665                 :            :   err_printf("\nMPQS: values from sieve vs. distribution of evaluated candidates:\n");
     666                 :            :   err_printf("   val  ___full __lprel ___none ___total\n");
     667                 :            :   for (i = 127; i >= 0; i--)
     668                 :            :   {
     669                 :            :     long rowtot = h->histo_full[i] + h->histo_lprl[i] + h->histo_drop[i];
     670                 :            :     tot += rowtot;
     671                 :            :     if ((rowtot > 0) || (i == h->sieve_threshold))
     672                 :            :       err_printf("%s[%3d] %7ld %7ld %7ld %8ld\n",
     673                 :            :                  i + 128 == h->sieve_threshold ? "^-" : "  ", i + 128,
     674                 :            :                  h->histo_full[i], h->histo_lprl[i], h->histo_drop[i],
     675                 :            :                  rowtot);
     676                 :            :   }
     677                 :            :   err_printf("        (total evaluated candidates: %ld)\n", tot);
     678                 :            : }
     679                 :            : 
     680                 :            : /* evaluation/feedback heuristics:
     681                 :            :  * First of all, refuse to draw any conclusions unless and until there's
     682                 :            :  * enough material to be statistically significant.
     683                 :            :  * Second, after sifting through the histo arrays, the new threshold is
     684                 :            :  * set to the minimum of the following three quantities:
     685                 :            :  * - the position where going down from the top value the histo_full
     686                 :            :  *   totals first exceed 100 - MPQS_HISTO_FREL_QUANTILE percent of all
     687                 :            :  *   full relations found so far by direct sieving.  (I.e. if the quantile
     688                 :            :  *   is 4, we want to keep all rows which account for 96% of all frels
     689                 :            :  *   obtained from the sieve.  Note that once we increase the threshold,
     690                 :            :  *   further counts will be biased against smaller values;  but we normally
     691                 :            :  *   don't expect to do many adjustments.)
     692                 :            :  * - the position where, going down from the top towards smaller values,
     693                 :            :  *   the cumulative number of useless candidates in histo_drop first exceeds
     694                 :            :  *   MPQS_HISTO_DROP_LIMIT times the number of useful ones.  I.e. when that
     695                 :            :  *   limit is 2, we're aiming for at least about 1/3 of all candidates coming
     696                 :            :  *   from the sieve to result in usable relations.
     697                 :            :  * - one less than the position where the histo_lprl count first falls below
     698                 :            :  *   MPQS_HISTO_LPREL_BASEFLOW times the number of useless candidates.  This
     699                 :            :  *   one will be capable of lowering the current threshold  (but never below
     700                 :            :  *   128).
     701                 :            :  * FIXME: For Double Large Prime mode, this will need to be seriously reworked.
     702                 :            :  */
     703                 :            : 
     704                 :            : /* This function returns 1 when it actually did something and decided on a
     705                 :            :  * good threshold  (although possibly the same as it was before),  -1 when
     706                 :            :  * there was nothing to do and never will be ("don't call us again"), 0
     707                 :            :  * when the caller should retry somewhat later.  Note that mpqs() already
     708                 :            :  * knows the total number of candidates generated so far  (from the return
     709                 :            :  * values of mpqs_eval_sieve()),  and won't call us too early;  but we also
     710                 :            :  * insist on minimal standards for the column sums.  Conversely when we ever
     711                 :            :  * lower the threshold, we ask for a re-evaluation later on.
     712                 :            :  * NB With the present accounting, once the threshold has been raised, it
     713                 :            :  * won't ever be lowered again, since the increasing counts above it will
     714                 :            :  * totally swamp the few earlier measurements below which can no longer
     715                 :            :  * grow.  So we might chop off those accumulating loops at the current sieve
     716                 :            :  * threshold. */
     717                 :            : static int
     718                 :            : mpqs_eval_histograms(mpqs_handle_t *h)
     719                 :            : {
     720                 :            :   long tot_full = 0, tot_lprl = 0, tot_drop = 0, total = 0;
     721                 :            :   long target_full, i;
     722                 :            :   int th_full, th_base, th_drop;
     723                 :            :   int th = h->sieve_threshold - 128;
     724                 :            : 
     725                 :            :   if (!h->do_histograms) return -1;
     726                 :            : 
     727                 :            :   /* first compute column sums */
     728                 :            :   for (i = 127; i >= 0; i--)
     729                 :            :   {
     730                 :            :     tot_full += h->histo_full[i];
     731                 :            :     tot_lprl += h->histo_lprl[i];
     732                 :            :     tot_drop += h->histo_drop[i];
     733                 :            :   }
     734                 :            :   total = tot_full + tot_lprl + tot_drop;
     735                 :            :   if ((total < MPQS_MIN_CANDS_FOR_HISTO) ||
     736                 :            :       (tot_full < MPQS_MIN_FRELS_FOR_HISTO))
     737                 :            :     return 0;                   /* too early to call the race */
     738                 :            : 
     739                 :            :   th_full = th_drop = th_base = -1;
     740                 :            :   /* find the full relations quantile point */
     741                 :            :   target_full = tot_full - (tot_full * MPQS_HISTO_FREL_QUANTILE) / 100.;
     742                 :            : 
     743                 :            :   tot_full = 0;
     744                 :            :   for (i = 127; i >= th; i--)
     745                 :            :   {
     746                 :            :     if ((tot_full += h->histo_full[i]) >= target_full)
     747                 :            :     {
     748                 :            :       th_full = i; break;
     749                 :            :     }
     750                 :            :   }
     751                 :            : 
     752                 :            :   /* find the "lp relations baseflow" point */
     753                 :            :   for (i = 127; i >= th; i--)
     754                 :            :   {
     755                 :            :     if (h->histo_lprl[i] + 1 <
     756                 :            :         MPQS_HISTO_LPREL_BASEFLOW * h->histo_drop[i])
     757                 :            :     {
     758                 :            :       th_base = i; break;
     759                 :            :     }
     760                 :            :   }
     761                 :            : 
     762                 :            :   /* find the wastefulness point */
     763                 :            :   tot_lprl = 0; tot_drop = 0;
     764                 :            :   for (i = 127; i >= th; i--)
     765                 :            :   {
     766                 :            :     tot_lprl += h->histo_full[i] + h->histo_lprl[i];
     767                 :            :     tot_drop += h->histo_drop[i];
     768                 :            :     if (tot_drop >
     769                 :            :         MPQS_HISTO_DROP_LIMIT * (tot_lprl + 1))
     770                 :            :     {
     771                 :            :       th_drop = i; break;
     772                 :            :     }
     773                 :            :   }
     774                 :            :   /* if these loops found nothing, then th_(full|base|drop) will still be -1.
     775                 :            :    * We won't tighten the sieve, but th_base would tell us we should loosen
     776                 :            :    * it  (reluctantly). */
     777                 :            : 
     778                 :            :   if (MPQS_DEBUGLEVEL >= 5)
     779                 :            :   {
     780                 :            :     mpqs_print_histo(h);
     781                 :            :     if (th_full >= 0)
     782                 :            :       err_printf("MPQS: threshold estimate for full rels: %d\n",
     783                 :            :                  th_full + 128);
     784                 :            :     if (th_drop >= 0)
     785                 :            :       err_printf("MPQS: threshold estimate for useful candidates: %d\n",
     786                 :            :                  th_drop + 128);
     787                 :            :   }
     788                 :            : 
     789                 :            :   /* any reason to open up the sieve?  wait until a minimal number of lprels
     790                 :            :    * at the threshold has been seen before going down... */
     791                 :            :   if ((th > 0) && (th_base <= th) &&
     792                 :            :       (h->histo_lprl[th] > (MPQS_MIN_FRELS_FOR_HISTO * 3.5)) )
     793                 :            :   {
     794                 :            :     h->sieve_threshold = th + 127;
     795                 :            :     if (MPQS_DEBUGLEVEL >= 4)
     796                 :            :       err_printf("MPQS: loosening sieve tolerance, new threshold %d\n",
     797                 :            :                  h->sieve_threshold);
     798                 :            :     return 0;                   /* this should be re-examined after a bit */
     799                 :            :   }
     800                 :            :   /* otherwise, any reason to tighten it? */
     801                 :            :   th = (th_full < th_drop ? th_full : th_drop) + 128;
     802                 :            :   if (th > h->sieve_threshold)
     803                 :            :   {
     804                 :            :     h->sieve_threshold = th;
     805                 :            :     if (MPQS_DEBUGLEVEL >= 4)
     806                 :            :       err_printf("MPQS: tightening sieve tolerance, new threshold %d\n",
     807                 :            :                  h->sieve_threshold);
     808                 :            :   }
     809                 :            :   /* maybe also loosen it if th_drop persistently stays below th... */
     810                 :            :   return 1;                     /* wait a good while before rechecking */
     811                 :            : }
     812                 :            : #endif
     813                 :            : 
     814                 :            : /*********************************************************************/
     815                 :            : /**                                                                 **/
     816                 :            : /**           RELATIONS AS STRINGS AND RELATIONS DATABASE           **/
     817                 :            : /**                                                                 **/
     818                 :            : /*********************************************************************/
     819                 :            : 
     820                 :            : /* determines a unique name for a file based on a short nickname
     821                 :            :  * name is allocated on the stack */
     822                 :            : static char *
     823                 :       2310 : mpqs_get_filename(char *dir, const char *s)
     824                 :            : {
     825                 :       2310 :   char *buf = stack_malloc(strlen(dir) + strlen(s) + 2);
     826                 :            : #if defined(__EMX__) || defined(WINCE)
     827                 :            :   sprintf(buf, "%s\\%s", dir,s);
     828                 :            : #else
     829                 :       2310 :   sprintf(buf, "%s/%s", dir,s);
     830                 :            : #endif
     831                 :       2310 :   return buf;
     832                 :            : }
     833                 :            : 
     834                 :            : /* compares two `large prime' relations according to their first element
     835                 :            :  * (the large prime itself). */
     836                 :            : static int
     837                 :     714645 : mpqs_relations_cmp(const void *a, const void *b)
     838                 :            : {
     839                 :     714645 :   char **sa = (char**) a;
     840                 :     714645 :   char **sb = (char**) b;
     841                 :     714645 :   long qa = strtol(*sa, NULL, 10);
     842                 :     714645 :   long qb = strtol(*sb, NULL, 10);
     843                 :            :   /* atol() isn't entirely portable for the Full Relations case where the
     844                 :            :      strings of digits are too long to fit into a long --GN */
     845         [ +  + ]:     714645 :   if (qa < qb) return -1;
     846         [ +  + ]:     479546 :   else if (qa > qb) return 1;
     847                 :     714645 :   else return strcmp(*sa, *sb);
     848                 :            : }
     849                 :            : 
     850                 :            : static void
     851                 :    1116084 : pari_fputs(char *s, pariFILE *f)
     852                 :            : {
     853         [ -  + ]:    1116084 :   if (fputs(s, f->file) < 0) pari_err_FILE("output file [fputs]", f->name);
     854                 :    1116084 : }
     855                 :            : #define min_bufspace 120UL /* use new buffer when < min_bufspace left */
     856                 :            : #define buflist_size 1024  /* size of list-of-buffers blocks */
     857                 :            : 
     858                 :            : /* Given a file "filename" containing full or `large prime' relations,
     859                 :            :  * rearrange the file so that relations are sorted by their first elements.
     860                 :            :  * Works also for sorting full relations. Works in memory, discards duplicate
     861                 :            :  * lines, and overwrites the original file. */
     862                 :            : static long
     863                 :       6466 : mpqs_sort_lp_file(char *filename)
     864                 :            : {
     865                 :            :   pariFILE *pTMP;
     866                 :            :   FILE *TMP;
     867                 :            :   char *old_s, *buf, *cur_line;
     868                 :            :   char **sort_table, **buflist, **next_buflist, **buflist_head;
     869                 :            :   long i, j, count;
     870                 :            :   size_t length, bufspace;
     871                 :       6466 :   pari_sp av=avma;
     872                 :            : 
     873                 :       6466 :   buflist_head = (char**) stack_malloc(buflist_size * sizeof(char*));
     874                 :       6466 :   buflist = buflist_head;
     875                 :       6466 :   *buflist++ = NULL; /* flag this as last and only buflist block */
     876                 :            :   /* extra blocks may be allocated as needed and linked ahead of
     877                 :            :    * buflist_head.  NB: whilst extra buflist blocks might have been
     878                 :            :    * needed when we were still sorting entire FREL files (more than 1023
     879                 :            :    * buffers, corresponding to about 20000 lines of ~200 characters), they
     880                 :            :    * should never be touched now that we only sort LPNEW and FNEW files, which
     881                 :            :    * are rather shorter. But the code might as well stay around for future
     882                 :            :    * upgrades to handling even larger numbers (and factor bases and thus
     883                 :            :    * relations files).  It costs one comparison per buffer allocation. --GN */
     884                 :            : 
     885                 :       6466 :   pTMP = pari_fopen_or_fail(filename, READ);
     886                 :       6466 :   TMP = pTMP->file;
     887                 :            :   /* get first buffer and read first line, if any, into it */
     888                 :       6466 :   buf = (char*) pari_malloc(MPQS_STRING_LENGTH * sizeof(char));
     889                 :       6466 :   cur_line = buf;
     890                 :       6466 :   bufspace = MPQS_STRING_LENGTH;
     891                 :            : 
     892         [ +  + ]:       6466 :   if (fgets(cur_line, bufspace, TMP) == NULL)
     893                 :            :   { /* file empty */
     894                 :       2926 :     pari_free(buf); pari_fclose(pTMP);
     895                 :       2926 :     avma = av; return 0;
     896                 :            :   }
     897                 :            :   /* enter first buffer into buflist */
     898                 :       3540 :   *buflist++ = buf; /* can't overflow the buflist block */
     899                 :       3540 :   length = strlen(cur_line) + 1; /* count the \0 byte as well */
     900                 :       3540 :   bufspace -= length;
     901                 :            : 
     902                 :       3540 :   sort_table = (char**)avma;
     903                 :            :   /* at start of loop, one line from the file is sitting in cur_line inside buf,
     904                 :            :    * the next will go into cur_line + length, and there's room for bufspace
     905                 :            :    * further characters in buf. The loop reads another line if one exists, and
     906                 :            :    * if this overruns the current buffer, it allocates a fresh one --GN */
     907                 :       3540 :   for (i=0, sort_table--; /* until end of file */; i++, sort_table--)
     908                 :            :   { /* sort_table is allocated on the stack, 0x100 cells at a time. Hence the
     909                 :            :      * stack must be left alone in the rest of the loop to keep the array
     910                 :            :      * connected. In particular, buffers can't be new_chunk'ed */
     911         [ +  + ]:     158640 :     if ((i & 0xff) == 0) (void)new_chunk(0x100);
     912                 :     158640 :     *sort_table = cur_line;
     913                 :     158640 :     cur_line += length;
     914                 :            : 
     915                 :            :     /* if little room is left, allocate a fresh buffer before attempting to
     916                 :            :      * read a line, and remember to free it if no further line is forthcoming.
     917                 :            :      * This avoids some copying of partial lines --GN */
     918         [ +  + ]:     158640 :     if (bufspace < min_bufspace)
     919                 :            :     {
     920         [ -  + ]:       1960 :       if (MPQS_DEBUGLEVEL >= 7)
     921                 :          0 :         err_printf("MQPS: short of space -- another buffer for sorting\n");
     922                 :       1960 :       buf = (char*) pari_malloc(MPQS_STRING_LENGTH * sizeof(char));
     923                 :       1960 :       cur_line = buf;
     924                 :       1960 :       bufspace = MPQS_STRING_LENGTH;
     925         [ +  + ]:       1960 :       if (fgets(cur_line, bufspace, TMP) == NULL) { pari_free(buf); break; }
     926                 :            : 
     927                 :            :       /* remember buffer for later deallocation */
     928         [ -  + ]:       1885 :       if (buflist - buflist_head >= buflist_size)
     929                 :            :       { /* need another buflist block */
     930                 :          0 :         next_buflist = (char**) pari_malloc(buflist_size * sizeof(char*));
     931                 :          0 :         *next_buflist = (char*)buflist_head; /* link */
     932                 :          0 :         buflist_head = next_buflist;
     933                 :          0 :         buflist = buflist_head + 1;
     934                 :            :       }
     935                 :       1885 :       *buflist++ = buf;
     936                 :       1885 :       length = strlen(cur_line) + 1;
     937                 :       1885 :       bufspace -= length; continue;
     938                 :            :     }
     939                 :            : 
     940                 :            :     /* normal case:  try fitting another line into the current buffer */
     941         [ +  + ]:     156680 :     if (fgets(cur_line, bufspace, TMP) == NULL) break; /* none exists */
     942                 :     153215 :     length = strlen(cur_line) + 1;
     943                 :     153215 :     bufspace -= length;
     944                 :            : 
     945                 :            :     /* check whether we got the entire line or only part of it */
     946 [ +  + ][ +  - ]:     153215 :     if (bufspace == 0 && cur_line[length-2] != '\n')
     947                 :            :     {
     948                 :            :       size_t lg1;
     949         [ -  + ]:         30 :       if (MPQS_DEBUGLEVEL >= 7)
     950                 :          0 :         err_printf("MQPS: line wrap -- another buffer for sorting\n");
     951                 :         30 :       buf = (char*) pari_malloc(MPQS_STRING_LENGTH * sizeof(char));
     952                 :            :       /* remember buffer for later deallocation */
     953         [ -  + ]:         30 :       if (buflist - buflist_head >= buflist_size)
     954                 :            :       { /* need another buflist block */
     955                 :          0 :         next_buflist = (char**)pari_malloc(buflist_size * sizeof(char*));
     956                 :          0 :         *next_buflist = (char*)buflist_head; /* link */
     957                 :          0 :         buflist_head = next_buflist;
     958                 :          0 :         buflist = buflist_head + 1;
     959                 :            :       }
     960                 :         30 :       *buflist++ = buf;
     961                 :            : 
     962                 :            :       /* copy what we've got to the new buffer */
     963                 :         30 :       (void)strcpy(buf, cur_line); /* cannot overflow */
     964                 :         30 :       cur_line = buf + length - 1; /* point at the \0 byte */
     965                 :         30 :       bufspace = MPQS_STRING_LENGTH - length + 1;
     966                 :            :       /* read remainder of line */
     967         [ -  + ]:         30 :       if (fgets(cur_line, bufspace, TMP) == NULL)
     968                 :          0 :         pari_err_FILE("TMP file [fgets]", pTMP->name);
     969                 :         30 :       lg1 = strlen(cur_line);
     970                 :         30 :       length += lg1; /* we already counted the \0 once */
     971                 :         30 :       bufspace -= (lg1 + 1); /* but here we must take it into account */
     972                 :         30 :       cur_line = buf; /* back up to the beginning of the line */
     973                 :            :     }
     974                 :     155100 :   } /* for */
     975                 :            : 
     976                 :       3540 :   pari_fclose(pTMP);
     977                 :            : 
     978                 :            :   /* sort the whole lot in place by swapping pointers */
     979                 :       3540 :   qsort(sort_table, i, sizeof(char*), mpqs_relations_cmp);
     980                 :            : 
     981                 :            :   /* copy results back to the original file, skipping exact duplicates */
     982                 :       3540 :   pTMP = pari_fopen_or_fail(filename, WRITE);
     983                 :       3540 :   old_s = sort_table[0];
     984                 :       3540 :   pari_fputs(sort_table[0], pTMP);
     985                 :       3540 :   count = 1;
     986         [ +  + ]:     155100 :   for(j = 1; j < i; j++)
     987                 :            :   {
     988         [ +  + ]:     151560 :     if (strcmp(old_s, sort_table[j]))
     989                 :            :     {
     990                 :     151510 :       pari_fputs(sort_table[j], pTMP);
     991                 :     151510 :       count++;
     992                 :            :     }
     993                 :     151560 :     old_s = sort_table[j];
     994                 :            :   }
     995                 :       3540 :   pari_fclose(pTMP);
     996         [ -  + ]:       3540 :   if (MPQS_DEBUGLEVEL >= 6) err_printf("MPQS: done sorting one file.\n");
     997                 :            : 
     998                 :            :   /* deallocate buffers and any extraneous buflist blocks except the first */
     999         [ +  + ]:       8995 :   while (*--buflist)
    1000                 :            :   {
    1001         [ +  - ]:       5455 :     if (buflist != buflist_head) /* not a linkage pointer */
    1002                 :       5455 :       pari_free((void*) *buflist);   /* free a buffer */
    1003                 :            :     else
    1004                 :            :     { /* linkage pointer */
    1005                 :          0 :       next_buflist = (char**)(*buflist);
    1006                 :          0 :       pari_free((void*)buflist_head); /* free a buflist block */
    1007                 :          0 :       buflist_head = next_buflist;
    1008                 :          0 :       buflist = buflist_head + buflist_size;
    1009                 :            :     }
    1010                 :            :   }
    1011                 :       6466 :   avma = av; return count;
    1012                 :            : }
    1013                 :            : 
    1014                 :            : /* appends contents of file fp1 to f (auxiliary routine for merge sort) and
    1015                 :            :  * returns number of lines copied. Close f afterwards */
    1016                 :            : static long
    1017                 :       6436 : mpqs_append_file(pariFILE *f, FILE *fp1)
    1018                 :            : {
    1019                 :       6436 :   FILE *fp = f->file;
    1020                 :            :   char line[MPQS_STRING_LENGTH];
    1021                 :       6436 :   long c = 0;
    1022         [ +  + ]:     290859 :   while (fgets(line, MPQS_STRING_LENGTH, fp1)) { pari_fputs(line, f); c++; }
    1023         [ -  + ]:       6436 :   if (fflush(fp)) pari_warn(warner, "error whilst flushing file %s", f->name);
    1024                 :       6436 :   pari_fclose(f); return c;
    1025                 :            : }
    1026                 :            : 
    1027                 :            : /* Merge-sort on the files LPREL and LPNEW; assumes that LPREL and LPNEW are
    1028                 :            :  * already sorted. Creates/truncates the TMP file, writes result to it and
    1029                 :            :  * closes it (via mpqs_append_file()). Instead of LPREL, LPNEW we may also call
    1030                 :            :  * this with FREL, FNEW. In the latter case pCOMB should be NULL (and we
    1031                 :            :  * return the count of all full relations), in the former non-NULL (and we
    1032                 :            :  * return the count of frels we expect to be able to combine out of the
    1033                 :            :  * present lprels). If pCOMB is non-NULL, the combinable lprels are written
    1034                 :            :  * out to this separate file.
    1035                 :            :  * We keep only one occurrence of each `large prime' in TMP (i.e. in the
    1036                 :            :  * future LPREL file). --GN */
    1037                 :            : 
    1038                 :            : #define swap_lines() { char *line_tmp;\
    1039                 :            :   line_tmp = line_new_old; \
    1040                 :            :   line_new_old = line_new; \
    1041                 :            :   line_new = line_tmp; }
    1042                 :            : 
    1043                 :            : static long
    1044                 :       6466 : mpqs_mergesort_lp_file0(FILE *LPREL, FILE *LPNEW, pariFILE *pCOMB,
    1045                 :            :                         pariFILE *pTMP)
    1046                 :            : {
    1047                 :            :   char line1[MPQS_STRING_LENGTH], line2[MPQS_STRING_LENGTH];
    1048                 :            :   char line[MPQS_STRING_LENGTH];
    1049                 :       6466 :   char *line_new = line1, *line_new_old = line2;
    1050                 :       6466 :   long q_new, q_new_old = -1, q, i = 0, c = 0;
    1051                 :            :   long comb_in_progress;
    1052                 :            : 
    1053         [ +  + ]:       6466 :   if ( !fgets(line_new, MPQS_STRING_LENGTH, LPNEW) )
    1054                 :            :   { /* LPNEW is empty: copy LPREL to TMP. Could be done by a rename if we
    1055                 :            :      * didn't want to count the lines (again)... however, this case will not
    1056                 :            :      * normally happen */
    1057                 :       2926 :     i = mpqs_append_file(pTMP, LPREL);
    1058         [ -  + ]:       2926 :     return pCOMB ? 0 : i;
    1059                 :            :   }
    1060                 :            :   /* we now have a line_new from LPNEW */
    1061                 :            : 
    1062         [ +  + ]:       3540 :   if (!fgets(line, MPQS_STRING_LENGTH, LPREL))
    1063                 :            :   { /* LPREL is empty: copy LPNEW to TMP... almost. */
    1064                 :        415 :     pari_fputs(line_new, pTMP);
    1065         [ +  + ]:        415 :     if (!pCOMB)
    1066                 :            :     { /* full relations mode */
    1067                 :        385 :       i = mpqs_append_file(pTMP, LPNEW);
    1068                 :        385 :       return i + 1;
    1069                 :            :     }
    1070                 :            : 
    1071                 :            :     /* LP mode:  check for combinable relations */
    1072                 :         30 :     q_new_old = atol(line_new);
    1073                 :            :     /* we need to retain a copy of the old line just for a moment, because we
    1074                 :            :      * may yet have to write it to pCOMB. Do this by swapping the two buffers */
    1075                 :         30 :     swap_lines();
    1076                 :         30 :     comb_in_progress = 0;
    1077                 :         30 :     i = 0;
    1078                 :            : 
    1079         [ +  + ]:       3870 :     while (fgets(line_new, MPQS_STRING_LENGTH, LPNEW))
    1080                 :            :     {
    1081                 :       3840 :       q_new = atol(line_new);
    1082         [ +  + ]:       3840 :       if (q_new_old == q_new)
    1083                 :            :       { /* found combinables, check whether we're already busy on this
    1084                 :            :            particular `large prime' */
    1085         [ +  - ]:         65 :         if (!comb_in_progress)
    1086                 :            :         { /* if not, write first line to pCOMB, creating and opening the
    1087                 :            :            * file first if it isn't open yet */
    1088                 :         65 :           pari_fputs(line_new_old, pCOMB);
    1089                 :         65 :           comb_in_progress = 1;
    1090                 :            :         }
    1091                 :            :         /* in any case, write the current line, and count it */
    1092                 :         65 :         pari_fputs(line_new, pCOMB);
    1093                 :         65 :         i++;
    1094                 :            :       }
    1095                 :            :       else
    1096                 :            :       { /* not combinable */
    1097                 :       3775 :         q_new_old = q_new;
    1098                 :       3775 :         comb_in_progress = 0;
    1099                 :            :         /* and dump it to the TMP file */
    1100                 :       3775 :         pari_fputs(line_new, pTMP);
    1101                 :            :         /* and stash it away for a moment */
    1102                 :       3775 :         swap_lines();
    1103                 :            :       }
    1104                 :            :     } /* while */
    1105                 :         30 :     pari_fclose(pTMP); return i;
    1106                 :            :   }
    1107                 :            : 
    1108                 :            :   /* normal case: both LPNEW and LPREL are not empty */
    1109                 :       3125 :   q_new = atol(line_new);
    1110                 :       3125 :   q = atol(line);
    1111                 :            : 
    1112                 :            :   for(;;)
    1113                 :            :   { /* main merging loop */
    1114                 :      71990 :     i = comb_in_progress = 0;
    1115                 :            : 
    1116                 :            :     /* first the harder case:  let LPNEW catch up with LPREL, and possibly
    1117                 :            :        overtake it, checking for combinables coming from LPNEW alone */
    1118         [ +  + ]:     151014 :     while (q > q_new)
    1119                 :            :     {
    1120 [ +  + ][ +  + ]:      81015 :       if (!pCOMB || !comb_in_progress) pari_fputs(line_new, pTMP);
    1121         [ +  + ]:      81015 :       if (!pCOMB) c++; /* in FREL mode, count lines written */
    1122         [ +  + ]:      22508 :       else if (!comb_in_progress)
    1123                 :            :       {
    1124                 :      22343 :         q_new_old = q_new;
    1125                 :      22343 :         swap_lines();
    1126                 :            :       }
    1127         [ +  + ]:      81015 :       if (!fgets(line_new, MPQS_STRING_LENGTH, LPNEW))
    1128                 :            :       {
    1129                 :       1991 :         pari_fputs(line, pTMP);
    1130         [ +  + ]:       1991 :         if (!pCOMB) c++; else c += i;
    1131                 :       1991 :         i = mpqs_append_file(pTMP, LPREL);
    1132         [ +  + ]:       1991 :         return pCOMB? c: c + i;
    1133                 :            :       }
    1134                 :      79024 :       q_new = atol(line_new);
    1135         [ +  + ]:      79024 :       if (!pCOMB) continue;
    1136                 :            : 
    1137                 :            :       /* LP mode only: */
    1138         [ +  + ]:      22286 :       if (q_new_old != q_new) /* not combinable */
    1139                 :      22121 :         comb_in_progress = 0; /* next loop will deal with it, or loop may end */
    1140                 :            :       else
    1141                 :            :       { /* found combinables, check whether we're already busy on this
    1142                 :            :            `large prime' */
    1143         [ +  + ]:        165 :         if (!comb_in_progress)
    1144                 :            :         {
    1145                 :        155 :           pari_fputs(line_new_old, pCOMB);
    1146                 :        155 :           comb_in_progress = 1;
    1147                 :            :         }
    1148                 :            :         /* in any case, write the current line, and count it */
    1149                 :        165 :         pari_fputs(line_new, pCOMB);
    1150                 :        165 :         i++;
    1151                 :            :       }
    1152                 :            :     } /* while q > q_new */
    1153                 :            : 
    1154                 :            :     /* q <= q_new */
    1155                 :            : 
    1156         [ +  + ]:      69999 :     if (pCOMB) c += i;   /* accumulate count of combinables */
    1157                 :      69999 :     i = 0;               /* and clear it */
    1158                 :      69999 :     comb_in_progress = 0;/* redundant */
    1159                 :            : 
    1160                 :            :     /* now let LPREL catch up with LPNEW, and possibly overtake it */
    1161         [ +  + ]:     614916 :     while (q < q_new)
    1162                 :            :     {
    1163                 :     545283 :       pari_fputs(line, pTMP);
    1164         [ +  + ]:     545283 :       if (!pCOMB) c++;
    1165         [ +  + ]:     545283 :       if (!fgets(line, MPQS_STRING_LENGTH, LPREL))
    1166                 :            :       {
    1167                 :        366 :         pari_fputs(line_new, pTMP);
    1168                 :        366 :         i = mpqs_append_file(pTMP, LPNEW);
    1169         [ +  + ]:        366 :         return pCOMB? c: c + i + 1;
    1170                 :            :       }
    1171                 :            :       else
    1172                 :     544917 :         q = atol(line);
    1173                 :            :     }
    1174                 :            : 
    1175                 :            :     /* q >= q_new */
    1176                 :            : 
    1177                 :            :     /* Finally, it may happen that q == q_new, indicating combinables whose
    1178                 :            :      * `large prime' is already in LPREL, and appears now once or more often in
    1179                 :            :      * LPNEW. Thus in this sub-loop we advance LPNEW. The `line' from LPREL is
    1180                 :            :      * left alone, and will be written to TMP the next time around the main for
    1181                 :            :      * loop; we only write it to pCOMB here -- unless all we find is an exact
    1182                 :            :      * duplicate of the line we already have, that is. (There can be at most
    1183                 :            :      * one such, and if so it is simply discarded.) */
    1184         [ +  + ]:     107108 :     while (q == q_new)
    1185                 :            :     {
    1186         [ +  + ]:      38243 :       if (!strcmp(line_new, line))
    1187                 :            :       { /* duplicate -- move right ahead to the next LPNEW line */
    1188                 :            :         ;/* do nothing here */
    1189                 :            :       }
    1190         [ +  + ]:      37648 :       else if (!pCOMB)
    1191                 :            :       { /* full relations mode: write line_new out first, keep line */
    1192                 :      35203 :         pari_fputs(line_new, pTMP);
    1193                 :      35203 :         c++;
    1194                 :            :       }
    1195                 :            :       else
    1196                 :            :       { /* LP mode, and combinable relation */
    1197         [ +  + ]:       2445 :         if (!comb_in_progress)
    1198                 :            :         {
    1199                 :       2390 :           pari_fputs(line, pCOMB);
    1200                 :       2390 :           comb_in_progress = 1;
    1201                 :            :         }
    1202                 :       2445 :         pari_fputs(line_new, pCOMB);
    1203                 :       2445 :         i++;
    1204                 :            :       }
    1205                 :            :       /* NB comb_in_progress is cleared by q_new becoming bigger than q, thus
    1206                 :            :        * the current while loop terminating, the next time through the main for
    1207                 :            :        * loop */
    1208                 :            : 
    1209                 :            :       /* common ending: get another line_new, if any */
    1210         [ +  + ]:      38243 :       if (!fgets(line_new, MPQS_STRING_LENGTH, LPNEW))
    1211                 :            :       {
    1212                 :        768 :         pari_fputs(line, pTMP);
    1213         [ +  - ]:        768 :         if (!pCOMB) c++; else c += i;
    1214                 :        768 :         i = mpqs_append_file(pTMP, LPREL);
    1215         [ +  - ]:        768 :         return pCOMB? c: c + i;
    1216                 :            :       }
    1217                 :            :       else
    1218                 :      37475 :         q_new = atol(line_new);
    1219                 :            :     } /* while */
    1220                 :            : 
    1221         [ +  + ]:      68865 :     if (pCOMB) c += i; /* accumulate count of combinables */
    1222                 :      75331 :   }
    1223                 :            : }
    1224                 :            : 
    1225                 :            : static long
    1226                 :       6466 : mpqs_mergesort_lp_file(char *REL_str, char *NEW_str, char *TMP_str, pariFILE *pCOMB)
    1227                 :            : {
    1228                 :       6466 :   pariFILE *pREL = pari_fopen_or_fail(REL_str, READ);
    1229                 :       6466 :   pariFILE *pNEW = pari_fopen_or_fail(NEW_str, READ);
    1230                 :       6466 :   pariFILE *pTMP = pari_fopen_or_fail(TMP_str, WRITE);
    1231                 :            :   long tp;
    1232                 :            : 
    1233                 :       6466 :   tp = mpqs_mergesort_lp_file0(pREL->file, pNEW->file, pCOMB, pTMP);
    1234                 :       6466 :   pari_fclose(pREL);
    1235                 :       6466 :   pari_fclose(pNEW);
    1236                 :       6466 :   pari_unlink(REL_str);
    1237         [ -  + ]:       6466 :   if (rename(TMP_str,REL_str)) pari_err_FILE("output file [rename]", REL_str);
    1238         [ -  + ]:       6466 :   if (MPQS_DEBUGLEVEL >= 6)
    1239                 :          0 :     err_printf("MPQS: renamed file %s to %s\n", TMP_str, REL_str);
    1240                 :       6466 :   return tp;
    1241                 :            : }
    1242                 :            : 
    1243                 :            : 
    1244                 :            : /*********************************************************************/
    1245                 :            : /**                                                                 **/
    1246                 :            : /**                       SELF-INITIALIZATION                       **/
    1247                 :            : /**                                                                 **/
    1248                 :            : /*********************************************************************/
    1249                 :            : 
    1250                 :            : #ifdef MPQS_DEBUG
    1251                 :            : /* Debug-only helper routine: check correctness of the root z mod p_i
    1252                 :            :  * by evaluting A * z^2 + B * z + C mod p_i  (which should be 0).
    1253                 :            :  * C is written as (B*B-kN)/(4*A) */
    1254                 :            : static void
    1255                 :            : check_root(mpqs_handle_t *h, long p, long start)
    1256                 :            : {
    1257                 :            :   long z = start - ((long)(h->M) % p);
    1258                 :            :   if (smodis(addii(h->C, mului(z, addii(h->B, mului(z, h->A)))), p))
    1259                 :            :   {
    1260                 :            :     err_printf("MPQS: p = %ld\n", p);
    1261                 :            :     err_printf("MPQS: A = %Ps\n", h->A);
    1262                 :            :     err_printf("MPQS: B = %Ps\n", h->B);
    1263                 :            :     err_printf("MPQS: C = %Ps\n", h->C);
    1264                 :            :     err_printf("MPQS: z = %ld\n", z);
    1265                 :            :     pari_err_BUG("MPQS: self_init: found wrong polynomial");
    1266                 :            :   }
    1267                 :            : }
    1268                 :            : #endif
    1269                 :            : 
    1270                 :            : /* There are four parts to self-initialization, which are exercised at
    1271                 :            :  * somewhat different times:
    1272                 :            :  * - choosing a new coefficient A  (selecting the prime factors to go into it,
    1273                 :            :  *   and doing the required bookkeeping in the FB entries, including clearing
    1274                 :            :  *   out the flags from the previous cohort), together with:
    1275                 :            :  * - doing the actual computations associated with a new A
    1276                 :            :  * - choosing a new B keeping the same A (very much simpler and quicker)
    1277                 :            :  * - and a small common bit that needs to happen in both cases.
    1278                 :            :  * As to the first item, the new scheme works as follows:
    1279                 :            :  * We pick omega_A - 1 prime factors for A below the index2_FB point which
    1280                 :            :  * marks their ideal size, and one prime above this point, choosing the
    1281                 :            :  * latter so as to get log2(A) as close as possible to l2_target_A.
    1282                 :            :  * The lower prime factors are chosen using bit patterns of constant weight,
    1283                 :            :  * gradually moving away from index2_FB towards smaller FB subscripts.
    1284                 :            :  * If this bumps into index0_FB  (might happen for very small input),  we
    1285                 :            :  * back up by increasing index2_FB by two, and from then on choosing only
    1286                 :            :  * bit patterns with either or both of their bottom bits set, so at least
    1287                 :            :  * one of the omega_A - 1 smaller prime factor will be beyond the original
    1288                 :            :  * index2_FB point.  In this way we avoid re-using A's which had already
    1289                 :            :  * been done.
    1290                 :            :  * (The choice of the upper "flyer" prime is of course constrained by the
    1291                 :            :  * size of the FB, which normally should never be anywhere close to becoming
    1292                 :            :  * a problem.  In unfortunate cases, e.g. for very small kN, we might have
    1293                 :            :  * to live with a rather non-optimal choice.
    1294                 :            :  * Then again, MPQS as such is surprisingly robust.  One day, I had got the
    1295                 :            :  * order of entries in mpqs_parameterset_t mixed up, and was running on a
    1296                 :            :  * smallish N with a huuuuge factor base and extremely tiny sieve interval,
    1297                 :            :  * and it still managed to factor it fairly quickly...)
    1298                 :            :  *
    1299                 :            :  * Mathematically, there isn't much more to this than the usual formula for
    1300                 :            :  * solving a quadratic  (over the field of p elements for each prime p in
    1301                 :            :  * the FB which doesn't divide A),  solving a linear equation for each of
    1302                 :            :  * the primes which do divide A, and precomputing differences between roots
    1303                 :            :  * mod p so we can adjust the roots quickly when we change B.
    1304                 :            :  * See Thomas Sosnowski's diploma thesis.
    1305                 :            :  */
    1306                 :            : 
    1307                 :            : /* Helper function:
    1308                 :            :  * Increment *x (!=0) to a larger value which has the same number of 1s in its
    1309                 :            :  * binary representation.  Wraparound can be detected by the caller as long as
    1310                 :            :  * we keep total_no_of_primes_for_A strictly less than BITS_IN_LONG.
    1311                 :            :  *
    1312                 :            :  * Changed switch to increment *x in all cases to the next larger number
    1313                 :            :  * which (a) has the same count of 1 bits and (b) does not arise from the
    1314                 :            :  * old value by moving a single 1 bit one position to the left  (which was
    1315                 :            :  * undesirable for the sieve). --GN based on discussion with TP */
    1316                 :            : INLINE void
    1317                 :       2163 : mpqs_increment(mpqs_uint32_t *x)
    1318                 :            : {
    1319                 :       2163 :   mpqs_uint32_t r1_mask, r01_mask, slider=1UL;
    1320                 :            : 
    1321                 :            :   /* 32-way computed jump handles 22 out of 32 cases */
    1322   [ +  -  +  +  :       2163 :   switch (*x & 0x1F)
          +  +  -  -  +  
                      + ]
    1323                 :            :   {
    1324                 :            :   case 29:
    1325                 :         15 :     (*x)++; break; /* shifts a single bit, but we postprocess this case */
    1326                 :            :   case 26:
    1327                 :          0 :     (*x) += 2; break; /* again */
    1328                 :            :   case 1: case 3: case 6: case 9: case 11:
    1329                 :            :   case 17: case 19: case 22: case 25: case 27:
    1330                 :       1262 :     (*x) += 3; return;
    1331                 :            :   case 20:
    1332                 :         15 :     (*x) += 4; break; /* again */
    1333                 :            :   case 5: case 12: case 14: case 21:
    1334                 :         53 :     (*x) += 5; return;
    1335                 :            :   case 2: case 7: case 13: case 18: case 23:
    1336                 :        555 :     (*x) += 6; return;
    1337                 :            :   case 10:
    1338                 :          0 :     (*x) += 7; return;
    1339                 :            :   case 8:
    1340                 :          0 :     (*x) += 8; break; /* and again */
    1341                 :            :   case 4: case 15:
    1342                 :         65 :     (*x) += 12; return;
    1343                 :            :   default: /* 0, 16, 24, 28, 30, 31 */
    1344                 :            :     /* isolate rightmost 1 */
    1345                 :        198 :     r1_mask = ((*x ^ (*x - 1)) + 1) >> 1;
    1346                 :            :     /* isolate rightmost 1 which has a 0 to its left */
    1347                 :        198 :     r01_mask = ((*x ^ (*x + r1_mask)) + r1_mask) >> 2;
    1348                 :            :     /* simple cases.  Both of these shift a single bit one position to the
    1349                 :            :        left, and will need postprocessing */
    1350         [ -  + ]:        198 :     if (r1_mask == r01_mask) { *x += r1_mask; break; }
    1351         [ +  + ]:        198 :     if (r1_mask == 1) { *x += r01_mask; break; }
    1352                 :            :     /* general case -- idea: add r01_mask, kill off as many 1 bits as possible
    1353                 :            :      * to its right while at the same time filling in 1 bits from the LSB. */
    1354         [ +  + ]:        183 :     if (r1_mask == 2) { *x += (r01_mask>>1) + 1; return; }
    1355 [ +  + ][ +  + ]:        369 :     while (r01_mask > r1_mask && slider < r1_mask)
    1356                 :            :     {
    1357                 :        246 :       r01_mask >>= 1; slider <<= 1;
    1358                 :            :     }
    1359                 :        123 :     *x += r01_mask + slider - 1;
    1360                 :        123 :     return;
    1361                 :            :   }
    1362                 :            :   /* post-process all cases which couldn't be finalized above.  If we get
    1363                 :            :      here, slider still has its original value. */
    1364                 :         45 :   r1_mask = ((*x ^ (*x - 1)) + 1) >> 1;
    1365                 :         45 :   r01_mask = ((*x ^ (*x + r1_mask)) + r1_mask) >> 2;
    1366         [ -  + ]:         45 :   if (r1_mask == r01_mask) { *x += r1_mask; return; }
    1367         [ +  + ]:         45 :   if (r1_mask == 1) { *x += r01_mask; return; }
    1368         [ +  + ]:         30 :   if (r1_mask == 2) { *x += (r01_mask>>1) + 1; return; }
    1369 [ +  + ][ +  - ]:         30 :   while (r01_mask > r1_mask && slider < r1_mask)
    1370                 :            :   {
    1371                 :         15 :     r01_mask >>= 1; slider <<= 1;
    1372                 :            :   }
    1373                 :         15 :   *x += r01_mask + slider - 1;
    1374                 :       2163 :   return;
    1375                 :            : }
    1376                 :            : 
    1377                 :            : /* self-init (1): advancing the bit pattern, and choice of primes for A.
    1378                 :            :  * The first time this is called, it finds h->bin_index == 0, which tells us
    1379                 :            :  * to initialize things from scratch.  On later occasions, we need to begin
    1380                 :            :  * by clearing the MPQS_FBE_DIVIDES_A bit in the fbe_flags of the former
    1381                 :            :  * prime factors of A.  We use, of course, the per_A_pr array for finding
    1382                 :            :  * them.  Upon successful return, that array will have been filled in, and
    1383                 :            :  * the flag bits will have been turned on again in the right places.
    1384                 :            :  * We return 1 (true) when we could set things up successfully, and 0 when
    1385                 :            :  * we found we'd be using more bits to the left in bin_index than we have
    1386                 :            :  * matching primes for in the FB.  In the latter case, bin_index will be
    1387                 :            :  * zeroed out, index2_FB will be incremented by 2, index2_moved will be
    1388                 :            :  * turned on, and the caller, after checking that index2_FB has not become
    1389                 :            :  * too large, should just call us again, which then is guaranteed to succeed:
    1390                 :            :  * we'll start again with a right-justified sequence of 1 bits in bin_index,
    1391                 :            :  * now interpreted as selecting primes relative to the new index2_FB. */
    1392                 :            : #ifndef MPQS_DEBUG_SI_CHOOSE_PRIMES
    1393                 :            : #  define MPQS_DEBUG_SI_CHOOSE_PRIMES 0
    1394                 :            : #endif
    1395                 :            : INLINE int
    1396                 :       2548 : mpqs_si_choose_primes(mpqs_handle_t *h)
    1397                 :            : {
    1398                 :       2548 :   mpqs_FB_entry_t *FB = h->FB;
    1399                 :       2548 :   mpqs_per_A_prime_t *per_A_pr = h->per_A_pr;
    1400                 :       2548 :   double l2_last_p = h->l2_target_A;
    1401                 :       2548 :   mpqs_int32_t omega_A = h->omega_A;
    1402                 :            :   int i, j, v2, prev_last_p_idx;
    1403                 :       2548 :   int room = h->index2_FB - h->index0_FB - omega_A + 4;
    1404                 :            :   /* GN 20050723:  I.e., index2_FB minus (index0_FB + omega_A - 3) plus 1
    1405                 :            :    * The notion of room here (cf mpqs_locate_A_range() above) is the number
    1406                 :            :    * of primes at or below index2_FB which are eligible for A.
    1407                 :            :    * At the very least, we need omega_A - 1 of them, and it is guaranteed
    1408                 :            :    * by mpqs_locate_A_range() that at least this many are available when we
    1409                 :            :    * first get here.  The logic here ensures that the lowest FB slot used
    1410                 :            :    * for A is never less than index0_FB + omega_A - 3.  In other words, when
    1411                 :            :    * omega_A == 3 (very small kN), we allow ourselves to reach all the way
    1412                 :            :    * down to index0_FB;  otherwise, we keep away from it by at least one
    1413                 :            :    * position.  For omega_A >= 4 this avoids situations where the selection
    1414                 :            :    * of the smaller primes here has advanced to a lot of very small ones, and
    1415                 :            :    * the single last larger one has soared away to bump into the top end of
    1416                 :            :    * the FB. */
    1417                 :            :   mpqs_uint32_t room_mask;
    1418                 :            :   mpqs_int32_t p;
    1419                 :            :   ulong bits;
    1420                 :            : 
    1421                 :            :   /* XXX also clear the index_j field here? */
    1422         [ +  + ]:       2548 :   if (h->bin_index == 0)
    1423                 :            :   {
    1424                 :            :     /* first time here, or after increasing index2_FB, initialize to a pattern
    1425                 :            :      * of omega_A - 1 consecutive right-justified 1 bits.
    1426                 :            :      * Caller will have ensured that there are enough primes for this in the
    1427                 :            :      * FB below index2_FB. */
    1428                 :        385 :     h->bin_index = (1UL << (omega_A - 1)) - 1;
    1429                 :        385 :     prev_last_p_idx = 0;
    1430                 :            :   }
    1431                 :            :   else
    1432                 :            :   {
    1433                 :            :     /* clear out the old flags */
    1434         [ +  + ]:      11189 :     for (i = 0; i < omega_A; i++)
    1435                 :       9026 :       MPQS_FLG(i) &= ~MPQS_FBE_DIVIDES_A;
    1436                 :       2163 :     prev_last_p_idx = MPQS_I(omega_A-1);
    1437                 :            : 
    1438                 :            :     /* find out how much maneuvering room we have before we're up against
    1439                 :            :      * the index0_FB wall */
    1440         [ +  + ]:       2163 :     if (room > 30) room = 30;
    1441                 :       2163 :     room_mask = ~((1UL << room) - 1);
    1442                 :            : 
    1443                 :            :     /* bump bin_index to the next acceptable value.  If index2_moved is off,
    1444                 :            :      * call mpqs_increment() just once;  otherwise, repeat until there's
    1445                 :            :      * something in the least significant 2 bits - this to ensure that we
    1446                 :            :      * never re-use an A which we'd used before increasing index2_FB - but
    1447                 :            :      * also stop if something shows up in the forbidden bits on the left
    1448                 :            :      * where we'd run out of bits or out of subscripts  (i.e. walk beyond
    1449                 :            :      * index0_FB + omega_A - 3). */
    1450                 :       2163 :     mpqs_increment(&h->bin_index);
    1451         [ -  + ]:       2163 :     if (h->index2_moved)
    1452                 :            :     {
    1453         [ #  # ]:          0 :       while ((h->bin_index & (room_mask | 0x3)) == 0)
    1454                 :          0 :         mpqs_increment(&h->bin_index);
    1455                 :            :     }
    1456                 :            :     /* ok so did we fall off the edge on the left? */
    1457         [ -  + ]:       2163 :     if ((h->bin_index & room_mask) != 0)
    1458                 :            :     {
    1459                 :            :       /* Yes.  Turn on the index2_moved flag in the handle */
    1460                 :          0 :       h->index2_FB += 2;        /* caller to check this isn't too large!!! */
    1461                 :          0 :       h->index2_moved = 1;
    1462                 :          0 :       h->bin_index = 0;
    1463         [ #  # ]:          0 :       if (MPQS_DEBUG_SI_CHOOSE_PRIMES || (MPQS_DEBUGLEVEL >= 5))
    1464                 :          0 :         err_printf("MPQS: wrapping, more primes for A now chosen near FB[%ld] = %ld\n",
    1465                 :          0 :                    (long)h->index2_FB,
    1466                 :          0 :                    (long)FB[h->index2_FB].fbe_p);
    1467                 :          0 :       return 0;                 /* back off - caller should retry */
    1468                 :            :     }
    1469                 :            :   }
    1470                 :            :   /* assert: we aren't occupying any of the room_mask bits now, and if
    1471                 :            :    * index2_moved had already been on, at least one of the two LSBs is on */
    1472                 :       2548 :   bits = h->bin_index;
    1473         [ -  + ]:       2548 :   if (MPQS_DEBUG_SI_CHOOSE_PRIMES || (MPQS_DEBUGLEVEL >= 6))
    1474                 :          0 :     err_printf("MPQS: new bit pattern for primes for A: 0x%lX\n", bits);
    1475                 :            : 
    1476                 :            :   /* map bits to FB subscripts, counting downward with bit 0 corresponding
    1477                 :            :    * to index2_FB, and accumulate logarithms against l2_last_p */
    1478                 :       2548 :   j = h->index2_FB;
    1479                 :       2548 :   v2 = vals((long)bits);
    1480         [ +  + ]:       2548 :   if (v2) { j -= v2; bits >>= v2; }
    1481         [ +  - ]:       7953 :   for (i = omega_A - 2; i >= 0; i--)
    1482                 :            :   {
    1483                 :       7953 :     MPQS_I(i) = j;
    1484                 :       7953 :     l2_last_p -= MPQS_LP(i);
    1485                 :       7953 :     MPQS_FLG(i) |= MPQS_FBE_DIVIDES_A;
    1486                 :       7953 :     bits &= ~1UL;
    1487         [ +  + ]:       7953 :     if (!bits) break;           /* that was the i=0 iteration */
    1488                 :       5405 :     v2 = vals((long)bits);
    1489                 :       5405 :     j -= v2;
    1490                 :       5405 :     bits >>= v2;
    1491                 :            :   }
    1492                 :            :   /* Choose the larger prime.  Note we keep index2_FB <= size_of_FB - 3 */
    1493         [ +  - ]:      35201 :   for (j = h->index2_FB + 1; (p = FB[j].fbe_p) != 0; j++)
    1494                 :            :   {
    1495         [ +  + ]:      35201 :     if (FB[j].fbe_flogp > l2_last_p) break;
    1496                 :            :   }
    1497                 :            :   /* GN 20050724: The following trick avoids generating a relatively large
    1498                 :            :    * proportion of duplicate relations when the last prime happens to fall
    1499                 :            :    * into an area where there are large gaps from one FB prime to the next,
    1500                 :            :    * and would otherwise often be repeated  (so that successive A's would
    1501                 :            :    * wind up too similar to each other).  While this trick isn't perfect,
    1502                 :            :    * it seems to get rid of a major part of the potential duplication. */
    1503 [ +  - ][ +  + ]:       2548 :   if ((p != 0) && (j == prev_last_p_idx))
    1504                 :            :   {
    1505                 :        363 :     j++; p = FB[j].fbe_p;
    1506                 :            :   }
    1507                 :       5096 :   MPQS_I(omega_A - 1) = (p == 0 ? /* did we fall off the end of the FB? */
    1508         [ -  + ]:       2548 :                          h->size_of_FB + 1 : /* then improvise */
    1509                 :            :                          j);
    1510                 :       2548 :   MPQS_FLG(omega_A - 1) |= MPQS_FBE_DIVIDES_A;
    1511                 :            : 
    1512         [ -  + ]:       2548 :   if (MPQS_DEBUG_SI_CHOOSE_PRIMES || (MPQS_DEBUGLEVEL >= 6))
    1513                 :            :   {
    1514                 :          0 :     err_printf("MPQS: chose primes for A");
    1515         [ #  # ]:          0 :     for (i = 0; i < omega_A; i++)
    1516                 :            :     {
    1517         [ #  # ]:          0 :       err_printf(" FB[%ld]=%ld%s",
    1518                 :          0 :                  (long) MPQS_I(i),
    1519                 :          0 :                  (long) MPQS_AP(i),
    1520                 :          0 :                  i < omega_A - 1 ? "," : "\n");
    1521                 :            :     }
    1522                 :            :   }
    1523                 :       2548 :   return 1;
    1524                 :            : }
    1525                 :            : 
    1526                 :            : 
    1527                 :            : 
    1528                 :            : /* compute coefficients of the sieving polynomial for self initializing
    1529                 :            :  * variant. Coefficients A and B are returned and several tables are
    1530                 :            :  * updated.   */
    1531                 :            : /* A and B are kept on the PARI stack in preallocated GENs.  So is C when
    1532                 :            :  * we're compiled for debugging. */
    1533                 :            : static void
    1534                 :      25375 : mpqs_self_init(mpqs_handle_t *h)
    1535                 :            : {
    1536                 :      25375 :   const ulong size_of_FB = h->size_of_FB + 1;
    1537                 :      25375 :   mpqs_FB_entry_t *FB = h->FB;
    1538                 :      25375 :   mpqs_inv_A_H_t *inv_A_H = h->inv_A_H;
    1539                 :      25375 :   const pari_sp av = avma;
    1540                 :            :   GEN p1, p2;
    1541                 :      25375 :   GEN A = h->A;
    1542                 :      25375 :   GEN B = h->B;
    1543                 :      25375 :   mpqs_per_A_prime_t *per_A_pr = h->per_A_pr;
    1544                 :            :   long i, j;
    1545                 :            :   long inv_A2;
    1546                 :            :   ulong p;
    1547                 :            : 
    1548                 :            : #ifdef MPQS_DEBUG_AVMA
    1549                 :            :   err_printf("MPQS DEBUG: enter self init, avma = 0x%lX\n", (ulong)avma);
    1550                 :            : #endif
    1551                 :            : 
    1552                 :            :   /* when all of the B's have already been used, choose new A ;
    1553                 :            :      this is indicated by setting index_j to 0 */
    1554         [ +  + ]:      25375 :   if (++h->index_j == (mpqs_uint32_t)h->no_B)
    1555                 :            :   {
    1556                 :       2163 :     h->index_j = 0;
    1557                 :       2163 :     h->index_i++; /* count finished A's */
    1558                 :            :   }
    1559                 :            : 
    1560         [ +  + ]:      25375 :   if (h->index_j == 0)
    1561                 :            :   /* "mpqs_self_init_A()" */
    1562                 :            :   { /* compute first polynomial with new A */
    1563         [ -  + ]:       2548 :     if (!mpqs_si_choose_primes(h))
    1564                 :            :     {
    1565                 :            :       /* We ran out of room towards small primes, and index2_FB was raised.
    1566                 :            :        * Check that we're still ok in that direction before re-trying the
    1567                 :            :        * operation, which then is guaranteed to succeed. The invariant we
    1568                 :            :        * maintain towards the top end is that h->size_of_FB - h->index2_FB >= 3,
    1569                 :            :        * but note that our size_of_FB is one larger. */
    1570                 :            : 
    1571                 :            :       /* "throw exception up to caller." ( bin_index set to impossible value
    1572                 :            :        * 0 by mpqs_si_choose_primes() */
    1573         [ #  # ]:      25375 :       if (size_of_FB - h->index2_FB < 4) return;
    1574                 :          0 :       (void) mpqs_si_choose_primes(h);
    1575                 :            :     }
    1576                 :            :     /* assert: bin_index and per_A_pr are now populated with consistent
    1577                 :            :      * values */
    1578                 :            : 
    1579                 :            :     /* compute A = product of omega_A primes given by bin_index */
    1580                 :       2548 :     p1 = NULL;
    1581         [ +  + ]:      13049 :     for (i = 0; i < h->omega_A; i++)
    1582                 :            :     {
    1583                 :      10501 :       p = (ulong) MPQS_AP(i);
    1584         [ +  + ]:      10501 :       p1 = p1 ? muliu(p1, p): utoipos(p);
    1585                 :            :     }
    1586                 :       2548 :     affii(p1, A); avma = av;
    1587                 :            : 
    1588                 :            :     /* Compute H[i], 0 <= i < omega_A.  Also compute the initial
    1589                 :            :      * B = sum(v_i*H[i]), by taking all v_i = +1 */
    1590                 :            :     /* TODO: following needs to be changed later for segmented FB and sieve
    1591                 :            :      * interval, where we'll want to precompute several B's. */
    1592                 :       2548 :     p2 = NULL;
    1593         [ +  + ]:      13049 :     for (i = 0; i < h->omega_A; i++)
    1594                 :            :     {
    1595                 :      10501 :       p = (ulong) MPQS_AP(i);
    1596                 :      10501 :       p1 = divis(A, (long)p);
    1597                 :      10501 :       p1 = muliu(p1, Fl_inv(umodiu(p1, p), p));
    1598                 :      10501 :       p1 = muliu(p1, MPQS_SQRT(i));
    1599                 :      10501 :       affii(remii(p1, A), MPQS_H(i));
    1600         [ +  + ]:      10501 :       p2 = p2 ? addii(p2, MPQS_H(i)) : MPQS_H(i);
    1601                 :            :     }
    1602                 :       2548 :     affii(p2, B);
    1603                 :       2548 :     avma = av;                  /* must happen outside the loop */
    1604                 :            : 
    1605                 :            :     /* ensure B = 1 mod 4 */
    1606         [ +  + ]:       2548 :     if (mod2(B) == 0)
    1607                 :       1296 :       affii(addii(B, mului(mod4(A), A)), B); /* B += (A % 4) * A; */
    1608                 :            : 
    1609                 :       2548 :     p1 = shifti(A, 1);
    1610                 :            :     /* compute the roots z1, z2, of the polynomial Q(x) mod p_j and
    1611                 :            :      * initialize start1[i] with the first value p_i | Q(z1 + i p_j)
    1612                 :            :      * initialize start2[i] with the first value p_i | Q(z2 + i p_j)
    1613                 :            :      * The following loop "miraculously" does The Right Thing for the
    1614                 :            :      * primes dividing k (where sqrt_kN is 0 mod p).  Primes dividing A
    1615                 :            :      * are skipped here, and are handled further down in the common part
    1616                 :            :      * of SI. */
    1617         [ +  + ]:     955807 :     for (j = 3; (ulong)j <= size_of_FB; j++)
    1618                 :            :     {
    1619                 :            :       ulong mb, tmp1, tmp2, m;
    1620         [ +  + ]:     953259 :       if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
    1621                 :     942758 :       p = (ulong)FB[j].fbe_p; m = h->M % p;
    1622                 :     942758 :       inv_A2 = Fl_inv(umodiu(p1, p), p); /* = 1/(2*A) mod p_j */
    1623         [ +  + ]:     942758 :       mb = umodiu(B, p); if (mb) mb = p - mb;
    1624                 :            :       /* mb = -B mod p */
    1625                 :     942758 :       tmp1 = Fl_sub(mb, FB[j].fbe_sqrt_kN, p);
    1626                 :     942758 :       tmp1 = Fl_mul(tmp1, inv_A2, p);
    1627                 :     942758 :       FB[j].fbe_start1 = (mpqs_int32_t)Fl_add(tmp1, m, p);
    1628                 :            : 
    1629                 :     942758 :       tmp2 = Fl_add(mb, FB[j].fbe_sqrt_kN, p);
    1630                 :     942758 :       tmp2 = Fl_mul(tmp2, inv_A2, p);
    1631                 :     942758 :       FB[j].fbe_start2 = (mpqs_int32_t)Fl_add(tmp2, m, p);
    1632         [ +  + ]:    4593726 :       for (i = 0; i < h->omega_A - 1; i++)
    1633                 :            :       {
    1634         [ +  + ]:    3650968 :         ulong h = umodiu(MPQS_H(i), p) << 1; if (h > p) h -= p;
    1635                 :    3650968 :         MPQS_INV_A_H(i,j) = Fl_mul(h, inv_A2, p); /* 1/A * H[i] mod p_j */
    1636                 :            :       }
    1637                 :            :     }
    1638                 :            :   }
    1639                 :            :   else
    1640                 :            :   /* "mpqs_self_init_B()" */
    1641                 :            :   { /* no "real" computation -- use recursive formula */
    1642                 :            :     /* The following exploits that B is the sum of omega_A terms +-H[i].
    1643                 :            :      * Each time we switch to a new B, we choose a new pattern of signs;
    1644                 :            :      * the precomputation of the inv_A_H array allows us to change the
    1645                 :            :      * two arithmetic progressions equally fast.  The choice of sign
    1646                 :            :      * patterns does *not* follow the bit pattern of the ordinal number
    1647                 :            :      * of B in the current cohort;  rather, we use a Gray code, changing
    1648                 :            :      * only one sign each time.  When the i-th rightmost bit of the new
    1649                 :            :      * ordinal number index_j of B is 1, the sign of H[i] is changed;
    1650                 :            :      * the next bit to the left tells us whether we should be adding or
    1651                 :            :      * subtracting the difference term.  We never need to change the sign
    1652                 :            :      * of H[omega_A-1] (the topmost one), because that would just give us
    1653                 :            :      * the same sieve items Q(x) again with the opposite sign of x.  This
    1654                 :            :      * is why we only precomputed inv_A_H up to i = omega_A - 2. */
    1655                 :            : 
    1656                 :      22827 :     ulong v2 = 0;               /* 2-valuation of h->index_j */
    1657                 :            : 
    1658                 :      22827 :     j = h->index_j;
    1659                 :            :     /* could use vals() here, but we need to right shift the bit pattern
    1660                 :            :      * anyway in order to find out which inv_A_H entries must be added to or
    1661                 :            :      * subtracted from the modular roots */
    1662         [ +  + ]:      38260 :     while ((j & 1) == 0) { v2++; j >>= 1; }
    1663                 :            : 
    1664                 :            :     /* v2 = v_2(index_j), determine new starting positions for sieving */
    1665                 :      22827 :     p1 = shifti(MPQS_H(v2), 1);
    1666         [ +  + ]:      22827 :     if (j & 2)
    1667                 :            :     { /* j = 3 mod 4 */
    1668         [ +  + ]:    7174276 :       for (j = 3; (ulong)j <= size_of_FB; j++)
    1669                 :            :       {
    1670         [ +  + ]:    7164226 :         if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
    1671                 :    7115266 :         p = (ulong)FB[j].fbe_p;
    1672                 :    7115266 :         FB[j].fbe_start1 = Fl_sub(FB[j].fbe_start1, MPQS_INV_A_H(v2,j), p);
    1673                 :    7115266 :         FB[j].fbe_start2 = Fl_sub(FB[j].fbe_start2, MPQS_INV_A_H(v2,j), p);
    1674                 :            :       }
    1675                 :      10050 :       p1 = addii(B, p1);
    1676                 :            :     }
    1677                 :            :     else
    1678                 :            :     { /* j = 1 mod 4 */
    1679         [ +  + ]:    8212148 :       for (j = 3; (ulong)j <= size_of_FB; j++)
    1680                 :            :       {
    1681         [ +  + ]:    8199371 :         if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
    1682                 :    8139128 :         p = (ulong)FB[j].fbe_p;
    1683                 :    8139128 :         FB[j].fbe_start1 = Fl_add(FB[j].fbe_start1, MPQS_INV_A_H(v2,j), p);
    1684                 :    8139128 :         FB[j].fbe_start2 = Fl_add(FB[j].fbe_start2, MPQS_INV_A_H(v2,j), p);
    1685                 :            :       }
    1686                 :      12777 :       p1 = subii(B, p1);
    1687                 :            :     }
    1688                 :      22827 :     affii(p1, B);
    1689                 :            :   }
    1690                 :      25375 :   avma = av;
    1691                 :            : 
    1692                 :            :   /* p=2 is a special case.  start1[2], start2[2] are never looked at,
    1693                 :            :    * so don't bother setting them. */
    1694                 :            : 
    1695                 :            :   /* "mpqs_self_init_common()" */
    1696                 :            : 
    1697                 :            :   /* now compute zeros of polynomials that have only one zero mod p
    1698                 :            :      because p divides the coefficient A */
    1699                 :            : 
    1700                 :            :   /* compute coefficient -C */
    1701                 :      25375 :   p1 = diviiexact(subii(h->kN, sqri(B)), shifti(A, 2));
    1702                 :            : 
    1703         [ +  + ]:     145079 :   for (i = 0; i < h->omega_A; i++)
    1704                 :            :   {
    1705                 :            :     ulong tmp, s;
    1706                 :     119704 :     p = (ulong) MPQS_AP(i);
    1707                 :     119704 :     tmp = Fl_div(umodiu(p1, p), umodiu(B, p), p); s = (tmp + h->M) % p;
    1708                 :     119704 :     FB[MPQS_I(i)].fbe_start1 = (mpqs_int32_t)s;
    1709                 :     119704 :     FB[MPQS_I(i)].fbe_start2 = (mpqs_int32_t)s;
    1710                 :            :   }
    1711                 :            : 
    1712         [ -  + ]:      25375 :   if (MPQS_DEBUGLEVEL >= 6)
    1713                 :            :   {
    1714                 :            :     /* must happen before resetting avma, because of the absi() */
    1715         [ #  # ]:          0 :     err_printf("MPQS: chose Q_%ld(x) = %Ps x^2 %c %Ps x + C\n",
    1716                 :          0 :                (long) h->index_j, h->A,
    1717                 :          0 :                signe(h->B) < 0? '-': '+', absi(h->B));
    1718                 :            :   }
    1719                 :            : 
    1720                 :      25375 :   avma = av;
    1721                 :            : 
    1722                 :            : #ifdef MPQS_DEBUG
    1723                 :            :   /* stash C into the handle.  Since check_root() is the only thing which
    1724                 :            :    * uses it, and only for debugging, C doesn't even exist as a field in
    1725                 :            :    * the handle unless we're built with MPQS_DEBUG. */
    1726                 :            :   affii(negi(p1), h->C);
    1727                 :            :   for (j = 3; j <= size_of_FB; j++)
    1728                 :            :   {
    1729                 :            :     check_root(h, FB[j].fbe_p, FB[j].fbe_start1);
    1730                 :            :     check_root(h, FB[j].fbe_p, FB[j].fbe_start2); avma = av;
    1731                 :            :   }
    1732                 :            :   if (DEBUGLEVEL >= 6)
    1733                 :            :     PRINT_IF_VERBOSE("MPQS: checking of roots of Q(x) was successful\n");
    1734                 :            : #endif
    1735                 :            : 
    1736                 :            : #ifdef MPQS_DEBUG_AVMA
    1737                 :            :   err_printf("MPQS DEBUG: leave self init, avma = 0x%lX\n", (ulong)avma);
    1738                 :            : #endif
    1739                 :            : }
    1740                 :            : 
    1741                 :            : /*********************************************************************/
    1742                 :            : /**                                                                 **/
    1743                 :            : /**                           THE SIEVE                             **/
    1744                 :            : /**                                                                 **/
    1745                 :            : /*********************************************************************/
    1746                 :            : 
    1747                 :            : /* Main sieving routine:
    1748                 :            :  * p4 = 4*p (used for loop unrolling)
    1749                 :            :  * log_p: approximation for log(p)
    1750                 :            :  * begin: points to a sieve array
    1751                 :            :  * end: points to the end of the sieve array
    1752                 :            :  * starting_sieving_index: marks the first FB element used for sieving */
    1753                 :            : INLINE void
    1754                 :   32363840 : mpqs_sieve_p(unsigned char *begin, unsigned char *end,
    1755                 :            :              long p4, long p, unsigned char log_p)
    1756                 :            : {
    1757                 :   32363840 :   register unsigned char *e = end - p4;
    1758                 :            :   /* Loop unrolled some time ago. It might be better to let the compiler worry
    1759                 :            :    * about *this* kind of optimization, based on its knowledge of whatever
    1760                 :            :    * useful tricks the machine instruction set architecture is offering
    1761                 :            :    * ("speculative loads" being the buzzword). --GN */
    1762         [ +  + ]:  149507332 :   while (e - begin >= 0) /* signed comparison */
    1763                 :            :   {
    1764                 :  117143492 :     (*begin) += log_p, begin += p;
    1765                 :  117143492 :     (*begin) += log_p, begin += p;
    1766                 :  117143492 :     (*begin) += log_p, begin += p;
    1767                 :  117143492 :     (*begin) += log_p, begin += p;
    1768                 :            :   }
    1769         [ +  + ]:  107932297 :   while (end - begin >= 0) /* again */
    1770                 :   75568457 :     (*begin) += log_p, begin += p;
    1771                 :   32363840 : }
    1772                 :            : 
    1773                 :            : static void
    1774                 :      25375 : mpqs_sieve(mpqs_handle_t *h)
    1775                 :            : {
    1776                 :      25375 :   long p, l = h->index1_FB;
    1777                 :            :   mpqs_FB_entry_t *ptr_FB;
    1778                 :      25375 :   unsigned char *sieve_array = h->sieve_array;
    1779                 :      25375 :   unsigned char *sieve_array_end = h->sieve_array_end;
    1780                 :            : 
    1781         [ +  + ]:   16267147 :   for (ptr_FB = &(h->FB[l]); (p = ptr_FB->fbe_p) != 0; ptr_FB++, l++)
    1782                 :            :   {
    1783                 :   16241772 :     unsigned char log_p = ptr_FB->fbe_logval;
    1784                 :   16241772 :     long start1 = ptr_FB->fbe_start1;
    1785                 :   16241772 :     long start2 = ptr_FB->fbe_start2;
    1786                 :            : 
    1787                 :            :     /* sieve with FB[l] from start_1[l] */
    1788                 :            :     /* if start1 != start2 sieve with FB[l] from start_2[l] */
    1789                 :            :     /* Maybe it is more efficient not to have a conditional branch in
    1790                 :            :      * the present loop body, and instead to right-shift log_p one bit
    1791                 :            :      * based on a flag bit telling us that we're on a one-root prime?
    1792                 :            :      * And instead roll the two invocations of mpqs_sieve_p into one. */
    1793                 :   16241772 :     mpqs_sieve_p(sieve_array + start1, sieve_array_end, p << 2, p, log_p);
    1794         [ +  + ]:   16241772 :     if (start1 != start2)
    1795                 :   16122068 :       mpqs_sieve_p(sieve_array + start2, sieve_array_end, p << 2, p, log_p);
    1796                 :            :   }
    1797                 :      25375 : }
    1798                 :            : 
    1799                 :            : /******************************/
    1800                 :            : 
    1801                 :            : /* Could make shameless use of the fact that M is divisible by 4, but
    1802                 :            :  * let the compiler worry about loop unrolling.  Indeed I wonder whether
    1803                 :            :  * modern compilers woudln't have an easier time optimizing this if it
    1804                 :            :  * were written as array accesses.  Doing so. */
    1805                 :            : static long
    1806                 :      25375 : mpqs_eval_sieve(mpqs_handle_t *h)
    1807                 :            : {
    1808                 :      25375 :   long x = 0, count = 0, M_2 = h->M << 1;
    1809                 :            :   /* TODO: replace the following by an auto-adjusting threshold driven
    1810                 :            :    * by histogram yield measurements */
    1811                 :      25375 :   unsigned char th = h->sieve_threshold;
    1812                 :      25375 :   unsigned char *sieve_array = h->sieve_array;
    1813                 :      25375 :   long *candidates = h->candidates;
    1814                 :            : 
    1815                 :            :   /* The following variation on the original is marginally faster with a
    1816                 :            :    * good optimizing compiler.  Exploiting the sentinel, we don't need to
    1817                 :            :    * check for x < M_2 in the inner while loop - this more than makes up
    1818                 :            :    * for the "lack" of explicit unrolling.  Optimizations like loop
    1819                 :            :    * unrolling are best left to the compiler anyway... */
    1820         [ +  - ]:     183168 :   while (count < MPQS_CANDIDATE_ARRAY_SIZE - 1)
    1821                 :            :   {
    1822         [ +  + ]:  459565575 :     while (sieve_array[x] < th) x++;
    1823         [ +  + ]:     183168 :     if (x >= M_2) break;
    1824                 :     157793 :     candidates[count++] = x++;
    1825                 :            :   }
    1826                 :            : 
    1827                 :      25375 :   candidates[count] = 0; return count;
    1828                 :            : }
    1829                 :            : 
    1830                 :            : /*********************************************************************/
    1831                 :            : /**                                                                 **/
    1832                 :            : /**                     CONSTRUCTING RELATIONS                      **/
    1833                 :            : /**                                                                 **/
    1834                 :            : /*********************************************************************/
    1835                 :            : 
    1836                 :            : /* Main relation routine */
    1837                 :            : static void
    1838                 :    2090978 : mpqs_add_factor(char **last, ulong ei, ulong pi) {
    1839                 :    2090978 :   sprintf(*last, " %lu %lu", ei, pi);
    1840                 :    2090978 :   *last += strlen(*last);
    1841                 :    2090978 : }
    1842                 :            : 
    1843                 :            : /* concatenate " 0" */
    1844                 :            : static void
    1845                 :     155965 : mpqs_add_0(char **last) {
    1846                 :     155965 :   char *s = *last;
    1847                 :     155965 :   *s++ = ' ';
    1848                 :     155965 :   *s++ = '0';
    1849                 :     155965 :   *s++ = 0; *last = s;
    1850                 :     155965 : }
    1851                 :            : 
    1852                 :            : #ifdef MPQS_DEBUG
    1853                 :            : static GEN
    1854                 :            : mpqs_factorback(mpqs_handle_t *h, char *relations)
    1855                 :            : {
    1856                 :            :   char *s, *t = stack_strdup(relations), *tok;
    1857                 :            :   GEN N = h->N, prod = gen_1;
    1858                 :            :   long i;
    1859                 :            :   mpqs_FB_entry_t *FB = h->FB;
    1860                 :            : 
    1861                 :            :   s = paristrtok_r(t, " \n", &tok);
    1862                 :            :   while (s != NULL)
    1863                 :            :   {
    1864                 :            :     long e = atol(s); if (!e) break;
    1865                 :            :     s = paristrtok_r(NULL, " \n", &tok);
    1866                 :            :     i = atol(s);
    1867                 :            :     /* special case -1 */
    1868                 :            :     if (i == 1) { prod = Fp_neg(prod,N); s = paristrtok_r(NULL, " \n", &tok); continue; }
    1869                 :            :     prod = Fp_mul(prod, Fp_powu(utoipos(FB[i].fbe_p), e, N), N);
    1870                 :            :     s = paristrtok_r(NULL, " \n", &tok);
    1871                 :            :   }
    1872                 :            :   return prod;
    1873                 :            : }
    1874                 :            : #endif
    1875                 :            : 
    1876                 :            : /* NB FREL, LPREL are actually FNEW, LPNEW when we get called */
    1877                 :            : static long
    1878                 :      25185 : mpqs_eval_cand(mpqs_handle_t *h, long number_of_cand,
    1879                 :            :                FILE *FREL, FILE *LPREL)
    1880                 :            : {
    1881                 :            :   pari_sp av;
    1882                 :      25185 :   long number_of_relations = 0;
    1883                 :      25185 :   char *relations = h->relations;
    1884                 :      25185 :   long *relaprimes = h->relaprimes;
    1885                 :            :   ulong i, pi;
    1886                 :      25185 :   mpqs_FB_entry_t *FB = h->FB;
    1887                 :      25185 :   GEN A = h->A;
    1888                 :      25185 :   GEN B = h->B;                 /* we don't need coefficient C here */
    1889                 :            :   int pii;
    1890                 :      25185 :   long *candidates = h->candidates;
    1891                 :            : 
    1892                 :      25185 :   av = avma;
    1893                 :            : #ifdef MPQS_DEBUG_AVMA
    1894                 :            :   err_printf("MPQS DEBUG: enter eval cand, avma = 0x%lX\n", (ulong)avma);
    1895                 :            : #endif
    1896         [ +  + ]:     182978 :   for (i = 0; i < (ulong)number_of_cand; i++, avma = av)
    1897                 :            :   {
    1898                 :            :     GEN Qx, Qx_part, A_2x_plus_B, Y;
    1899                 :            :     long powers_of_2, p;
    1900                 :     157793 :     long x = candidates[i];
    1901                 :     157793 :     long x_minus_M = x - h->M;
    1902                 :     157793 :     char *relations_end = relations;
    1903                 :     157793 :     int relaprpos = 0;
    1904                 :            : 
    1905                 :            : #ifdef MPQS_DEBUG_AVMA
    1906                 :            :     err_printf("MPQS DEBUG: eval loop 1, avma = 0x%lX\n", (ulong)avma);
    1907                 :            : #endif
    1908                 :            : 
    1909                 :     157793 :     *relations_end = 0;
    1910                 :            : #ifdef MPQS_DEBUG_VERYVERBOSE
    1911                 :            :     err_printf("%c", (char)('0' + i%10));
    1912                 :            : #endif
    1913                 :            : 
    1914                 :            :     /* A_2x_plus_B = (A*(2x)+B), Qx = (A*(2x)+B)^2/(4*A) = Q(x) */
    1915                 :     157793 :     A_2x_plus_B = addii(mulis(A, x_minus_M << 1), B);
    1916                 :     157793 :     Y = absi(A_2x_plus_B);
    1917                 :            : 
    1918                 :     157793 :     Qx = subii(sqri(A_2x_plus_B), h->kN);
    1919                 :            : 
    1920                 :            : #ifdef MPQS_DEBUG_AVMA
    1921                 :            :     err_printf("MPQS DEBUG: eval loop 2, avma = 0x%lX\n", (ulong)avma);
    1922                 :            : #endif
    1923                 :            :     /* When N is relatively small, it may happen that Qx is outright
    1924                 :            :      * divisible by N at this point.  In any case, when no extensive prior
    1925                 :            :      * trial division / Rho / ECM had been attempted, gcd(Qx,N) may turn
    1926                 :            :      * out to be a nontrivial factor of N  (larger than what the FB contains
    1927                 :            :      * or we'd have found it already, but possibly smaller than the large-
    1928                 :            :      * prime bound).  This is too rare to check for here in the inner loop,
    1929                 :            :      * but it will be caught if such an LP relation is ever combined with
    1930                 :            :      * another. */
    1931                 :            : 
    1932                 :            :     /* Qx cannot possibly vanish here */
    1933         [ -  + ]:     157793 :     if (!signe(Qx)) { PRINT_IF_VERBOSE("<+>"); continue; }
    1934         [ +  + ]:     157793 :     else if (signe(Qx) < 0) {
    1935                 :      86544 :       setabssign(Qx);
    1936                 :      86544 :       mpqs_add_factor(&relations_end, 1, 1); /* i = 1, ei = 1, pi */
    1937                 :            :     }
    1938                 :            : 
    1939                 :            :     /* divide by powers of 2;  we're really dealing with 4*A*Q(x), so we
    1940                 :            :      * always have at least 2^2 here, and at least 2^3 when kN is 1 mod 4 */
    1941                 :     157793 :     powers_of_2 = vali(Qx);
    1942                 :     157793 :     Qx = shifti(Qx, -powers_of_2);
    1943                 :     157793 :     mpqs_add_factor(&relations_end, powers_of_2, 2);
    1944                 :            : 
    1945                 :            :     /* That has dealt with a possible -1 and the power of 2.  First pass
    1946                 :            :      * over odd primes in FB: pick up all possible divisors of Qx including
    1947                 :            :      * those sitting in k or in A, and remember them in relaprimes. Do not
    1948                 :            :      * yet worry about possible repeated factors, these will be found in the
    1949                 :            :      * second pass. */
    1950                 :     157793 :     Qx_part = A;
    1951                 :            : 
    1952                 :            :     /* The first pass recognizes divisors of A by their corresponding flags
    1953                 :            :      * bit in the FB entry.  (Divisors of k require no special treatment at
    1954                 :            :      * this stage.)  We construct a preliminary table of FB subscripts and
    1955                 :            :      * "exponents" of the FB primes which divide Qx.  (We store subscripts
    1956                 :            :      * rather than the primes themselves because the string representation
    1957                 :            :      * of a relation is in terms of the subscripts.)
    1958                 :            :      * We must distinguish three cases so we can do the right thing in the
    1959                 :            :      * 2nd pass: prime not in A which divides Qx, prime in A which does not
    1960                 :            :      * divide Qx/A, prime in A which does divide Qx/A. The first and third
    1961                 :            :      * kinds need checking for repeated factors, the second kind doesn't. The
    1962                 :            :      * first and second kinds contribute 1 to the exponent in the relation,
    1963                 :            :      * the 3rd kind contributes 2. We store 1,0,2 respectively in these three
    1964                 :            :      * cases.
    1965                 :            :      * Factors in common with k are much simpler - if they occur, they occur
    1966                 :            :      * exactly to the first power, and this makes no difference in the first
    1967                 :            :      * pass - here they behave just like every normal odd factor base prime.
    1968                 :            :      */
    1969                 :            : 
    1970         [ +  + ]:   96870585 :     for (pi = 3; (p = FB[pi].fbe_p); pi++)
    1971                 :            :     {
    1972                 :   96712792 :       long tmp_p = x % p;
    1973                 :   96712792 :       ulong ei = 0;
    1974                 :            : 
    1975                 :            :       /* Here we use that MPQS_FBE_DIVIDES_A equals 1. */
    1976                 :   96712792 :       ei = FB[pi].fbe_flags & MPQS_FBE_DIVIDES_A;
    1977                 :            : 
    1978 [ +  + ][ +  + ]:   96712792 :       if (tmp_p == FB[pi].fbe_start1 || tmp_p == FB[pi].fbe_start2)
    1979                 :            :       { /* p divides Q(x)/A (and possibly A), 1st or 3rd case */
    1980                 :    1134321 :         relaprimes[relaprpos++] = pi;
    1981                 :    1134321 :         relaprimes[relaprpos++] = 1 + ei;
    1982                 :    1134321 :         Qx_part = muliu(Qx_part, p);
    1983                 :            :       }
    1984         [ +  + ]:   95578471 :       else if (ei)
    1985                 :            :       { /* p divides A but does not divide Q(x)/A, 2nd case */
    1986                 :     712320 :         relaprimes[relaprpos++] = pi;
    1987                 :     712320 :         relaprimes[relaprpos++] = 0;
    1988                 :            :       }
    1989                 :            :     }
    1990                 :            : 
    1991                 :            :     /* We have now accumulated the known factors of Qx except for possible
    1992                 :            :      * repeated factors and for possible large primes.  Divide off what we
    1993                 :            :      * have.  (This is faster than dividing off A and each prime separately.)
    1994                 :            :      */
    1995                 :     157793 :     Qx = diviiexact(Qx, Qx_part);
    1996                 :            :     /* (ToDo: MPQS_DEBUG sanity check...) */
    1997                 :            : 
    1998                 :            : #ifdef MPQS_DEBUG_AVMA
    1999                 :            :     err_printf("MPQS DEBUG: eval loop 3, avma = 0x%lX\n", (ulong)avma);
    2000                 :            : #endif
    2001                 :            : 
    2002                 :            :     /* second pass - deal with any repeated factors, and write out the string
    2003                 :            :      * representation of the tentative relation. At this point, the only
    2004                 :            :      * primes which can occur again in the adjusted Qx are those in relaprimes
    2005                 :            :      * which are followed by 1 or 2.  We must pick up those followed by a 0,
    2006                 :            :      * too, though. */
    2007                 :            :     PRINT_IF_VERBOSE("a");
    2008         [ +  + ]:    2004434 :     for (pii = 0; pii < relaprpos; pii+=2)
    2009                 :            :     {
    2010                 :            :       long remd_p;
    2011                 :    1846641 :       ulong ei = relaprimes[pii+1];
    2012                 :            :       GEN Qx_div_p;
    2013                 :    1846641 :       pi = relaprimes[pii];
    2014                 :            : 
    2015                 :            :       /* Here, prime factors of k go their separate way.  We could have
    2016                 :            :        * introduced another FB entry flag for marking them, but it is easier
    2017                 :            :        * to identify them just by their position before index0_FB. */
    2018         [ +  + ]:    1846641 :       if ((mpqs_int32_t)pi < h->index0_FB) {
    2019                 :            : #ifdef MPQS_DEBUG
    2020                 :            :         PRINT_IF_VERBOSE("\bk!");
    2021                 :            : #endif
    2022                 :      25587 :         mpqs_add_factor(&relations_end, 1, pi);
    2023                 :      25587 :         continue;
    2024                 :            :       }
    2025         [ +  + ]:    1821054 :       if (ei == 0) /* p divides A and that was it */
    2026                 :            :       {
    2027                 :     712320 :         mpqs_add_factor(&relations_end, 1, pi);
    2028                 :     712320 :         continue;
    2029                 :            :       }
    2030                 :    1108734 :       p = FB[pi].fbe_p;
    2031                 :            : #ifdef MPQS_DEBUG_CANDIDATE_EVALUATION
    2032                 :            :       err_printf("MPQS DEBUG: Qx=%Ps p=%ld\n", Qx, (long)p);
    2033                 :            : #endif
    2034                 :            :       /* otherwise p might still divide the current adjusted Qx. Try it... */
    2035                 :            :       /* NOTE: break out of loop when remaining Qx is 1 ?  Or rather, suppress
    2036                 :            :        * the trial divisions, since we still need to write our string.
    2037                 :            :        * Actually instead of testing for 1, test whether Qx is smaller than p;
    2038                 :            :        * cf Karim's mail from 20050124.  If it is, without being 1, then it has
    2039                 :            :        * a common factor with k.  But those factors are soon going to have
    2040                 :            :        * disappeared before we get here.  However, inserting
    2041                 :            :        * an explicit if (!is_pm1(Qx)) here did not help any. */
    2042                 :    1108734 :       Qx_div_p = divis_rem(Qx, p, &remd_p);
    2043         [ +  + ]:    1162336 :       while (remd_p == 0) {
    2044                 :      53602 :         ei++; Qx = Qx_div_p;
    2045                 :      53602 :         Qx_div_p = divis_rem(Qx, p, &remd_p);
    2046                 :            :       }
    2047                 :    1846641 :       mpqs_add_factor(&relations_end, ei, pi);
    2048                 :            :     }
    2049                 :            : 
    2050                 :            : #ifdef MPQS_DEBUG_AVMA
    2051                 :            :     err_printf("MPQS DEBUG: eval loop 4, avma = 0x%lX\n", (ulong)avma);
    2052                 :            : #endif
    2053                 :            :     PRINT_IF_VERBOSE("\bb");
    2054         [ +  + ]:     157793 :     if (is_pm1(Qx))
    2055                 :            :     {
    2056                 :     126655 :       mpqs_add_0(&relations_end);
    2057                 :     126655 :       fprintf(FREL, "%s :%s\n", itostr(Y), relations);
    2058                 :     126655 :       number_of_relations++;
    2059                 :            : #ifdef MPQS_USE_HISTOGRAMS
    2060                 :            :       /* bump full relations counter at candidate's value */
    2061                 :            :       if (h->do_histograms) h->histo_full[sa[x]-128]++;
    2062                 :            : #endif
    2063                 :            : 
    2064                 :            : #ifdef MPQS_DEBUG
    2065                 :            :       {
    2066                 :            :         pari_sp av1 = avma;
    2067                 :            :         GEN rhs = mpqs_factorback(h, relations);
    2068                 :            :         GEN Qx_2 = remii(sqri(Y), h->N);
    2069                 :            :         if (!equalii(Qx_2, rhs))
    2070                 :            :         {
    2071                 :            :           PRINT_IF_VERBOSE("\b(!)\n");
    2072                 :            :           err_printf("MPQS: %Ps @ %Ps :%s\n", Y, Qx, relations);
    2073                 :            :           err_printf("\tQx_2 = %Ps\n", Qx_2);
    2074                 :            :           err_printf("\t rhs = %Ps\n", rhs);
    2075                 :            :           pari_err_BUG("MPQS: wrong full relation found");
    2076                 :            :         }
    2077                 :            :         else
    2078                 :            :           PRINT_IF_VERBOSE("\b(:)");
    2079                 :            :         avma = av1;
    2080                 :            :       }
    2081                 :            : #endif
    2082                 :            :     }
    2083         [ +  + ]:      31138 :     else if (cmpis(Qx, h->lp_bound) > 0)
    2084                 :            :     { /* TODO: check for double large prime */
    2085                 :            : #ifdef MPQS_USE_HISTOGRAMS
    2086                 :            :       /* bump useless-candidates counter at candidate's value */
    2087                 :            :       if (h->do_histograms) h->histo_drop[sa[x]-128]++;
    2088                 :            : #endif
    2089                 :            :       PRINT_IF_VERBOSE("\b.");
    2090                 :            :     }
    2091                 :            :     else
    2092                 :            :     { /* if (mpqs_isprime(itos(Qx))) */
    2093                 :      29310 :       mpqs_add_0(&relations_end);
    2094                 :     157793 :       fprintf(LPREL, "%s @ %s :%s\n", itostr(Qx), itostr(Y), relations);
    2095                 :            : #ifdef MPQS_USE_HISTOGRAMS
    2096                 :            :       /* bump LP relations counter at candidate's value */
    2097                 :            :       if (h->do_histograms) h->histo_lprl[sa[x]-128]++;
    2098                 :            : #endif
    2099                 :            : #ifdef MPQS_DEBUG
    2100                 :            :       {
    2101                 :            :         pari_sp av1 = avma;
    2102                 :            :         GEN rhs = mpqs_factorback(h, relations);
    2103                 :            :         GEN Qx_2 = remii(sqri(Y), h->N);
    2104                 :            : 
    2105                 :            :         rhs = modii(mulii(rhs, Qx), h->N);
    2106                 :            :         if (!equalii(Qx_2, rhs))
    2107                 :            :         {
    2108                 :            :           PRINT_IF_VERBOSE("\b(!)\n");
    2109                 :            :           err_printf("MPQS: %Ps @ %Ps :%s\n", Y, Qx, relations);
    2110                 :            :           err_printf("\tQx_2 = %Ps\n", Qx_2);
    2111                 :            :           err_printf("\t rhs = %Ps\n", rhs);
    2112                 :            :           pari_err_BUG("MPQS: wrong large prime relation found");
    2113                 :            :         }
    2114                 :            :         else
    2115                 :            :           PRINT_IF_VERBOSE("\b(;)");
    2116                 :            :         avma = av1;
    2117                 :            :       }
    2118                 :            : #endif
    2119                 :            :     }
    2120                 :            : 
    2121                 :            : #ifdef MPQS_DEBUG_AVMA
    2122                 :            :     err_printf("MPQS DEBUG: eval loop end, avma = 0x%lX\n", (ulong)avma);
    2123                 :            : #endif
    2124                 :            :   } /* for */
    2125                 :            :   PRINT_IF_VERBOSE("\n");
    2126                 :            : #ifdef MPQS_DEBUG_AVMA
    2127                 :            :   err_printf("MPQS DEBUG: leave eval cand, avma = 0x%lX\n", (ulong)avma);
    2128                 :            : #endif
    2129                 :      25185 :   return number_of_relations;
    2130                 :            : }
    2131                 :            : 
    2132                 :            : /*********************************************************************/
    2133                 :            : /**                                                                 **/
    2134                 :            : /**                      COMBINING RELATIONS                        **/
    2135                 :            : /**                                                                 **/
    2136                 :            : /*********************************************************************/
    2137                 :            : 
    2138                 :            : /* combines the large prime relations in COMB to full relations in FNEW.
    2139                 :            :  * FNEW is assumed to be open for writing / appending. */
    2140                 :            : 
    2141                 :            : typedef struct {
    2142                 :            :   long q;
    2143                 :            :   char Y[MPQS_STRING_LENGTH];
    2144                 :            :   char E[MPQS_STRING_LENGTH];
    2145                 :            : } mpqs_lp_entry;
    2146                 :            : 
    2147                 :            : static void
    2148                 :       5350 : mpqs_set_exponents(long *ei, char *r)
    2149                 :            : {
    2150                 :            :   char *s, b[MPQS_STRING_LENGTH], *tok;
    2151                 :            :   long e;
    2152                 :            : 
    2153                 :       5350 :   strcpy(b, r);
    2154                 :       5350 :   s = paristrtok_r(b, " \n", &tok);
    2155         [ +  - ]:      79270 :   while (s != NULL)
    2156                 :            :   {
    2157         [ +  + ]:      79270 :     e = atol(s); if (!e) break;
    2158                 :      73920 :     s = paristrtok_r(NULL, " \n", &tok);
    2159                 :      73920 :     ei[ atol(s) ] += e;
    2160                 :      73920 :     s = paristrtok_r(NULL, " \n", &tok);
    2161                 :            :   }
    2162                 :       5350 : }
    2163                 :            : 
    2164                 :            : static void
    2165                 :       5285 : set_lp_entry(mpqs_lp_entry *e, char *buf)
    2166                 :            : {
    2167                 :            :   char *s1, *s2;
    2168                 :       5285 :   s1 = buf; s2 = strchr(s1, ' '); *s2 = '\0';
    2169                 :       5285 :   e->q = atol(s1);
    2170                 :       5285 :   s1 = s2 + 3; s2 = strchr(s1, ' '); *s2 = '\0';
    2171                 :       5285 :   strcpy(e->Y, s1);
    2172                 :       5285 :   s1 = s2 + 3; s2 = strchr(s1, '\n'); *s2 = '\0';
    2173                 :       5285 :   strcpy(e->E, s1);
    2174                 :       5285 : }
    2175                 :            : 
    2176                 :            : static long
    2177                 :        272 : mpqs_combine_large_primes(mpqs_handle_t *h,
    2178                 :            :                           FILE *COMB, pariFILE *pFNEW, GEN *f)
    2179                 :            : {
    2180                 :        272 :   pari_sp av0 = avma, av, av2;
    2181                 :            :   char new_relation[MPQS_STRING_LENGTH], buf[MPQS_STRING_LENGTH];
    2182                 :            :   mpqs_lp_entry e[2]; /* we'll use the two alternatingly */
    2183                 :        272 :   long *ei, ei_size = h->size_of_FB + 2;
    2184                 :            :   long old_q;
    2185                 :            :   GEN inv_q, Y1, Y2, new_Y, new_Y1;
    2186                 :        272 :   long i, l, c = 0;
    2187                 :            : 
    2188                 :        272 :   *f = NULL;
    2189         [ -  + ]:        272 :   if (!fgets(buf, MPQS_STRING_LENGTH, COMB)) return 0; /* should not happen */
    2190                 :            : 
    2191                 :        272 :   ei = (long *) new_chunk(ei_size);
    2192                 :        272 :   av = avma;
    2193                 :            :   /* put first lp relation in row 0 of e */
    2194                 :        272 :   set_lp_entry(&e[0], buf);
    2195                 :            : 
    2196                 :        272 :   i = 1; /* second relation will go into row 1 */
    2197                 :        272 :   old_q = e[0].q;
    2198         [ -  + ]:        272 :   while (!invmod(utoipos(old_q), h->N, &inv_q)) /* can happen */
    2199                 :            :   {
    2200                 :          0 :     inv_q = gcdii(inv_q, h->N);
    2201                 :            :     /* inv_q can no longer be 1 here (it could while we were doing this mod
    2202                 :            :      * kN instead of mod N), but never mind - we're not in the fast path
    2203                 :            :      * at this point.  It could be N when N is quite small;  or we might
    2204                 :            :      * just have found a divisor by sheer luck. */
    2205 [ #  # ][ #  # ]:          0 :     if (is_pm1(inv_q) || equalii(inv_q, h->N)) /* pity */
    2206                 :            :     {
    2207                 :            : #ifdef MPQS_DEBUG
    2208                 :            :       err_printf("MPQS: skipping relation with non-invertible q\n");
    2209                 :            : #endif
    2210         [ #  # ]:          0 :       if (!fgets(buf, MPQS_STRING_LENGTH, COMB)) { avma = av0; return 0; }
    2211                 :          0 :       avma = av;
    2212                 :          0 :       set_lp_entry(&e[0], buf);
    2213                 :          0 :       old_q = e[0].q; continue;
    2214                 :            :     }
    2215                 :          0 :     *f = gerepileuptoint(av0, inv_q);
    2216                 :          0 :     return c;
    2217                 :            :   }
    2218                 :        272 :   Y1 = strtoi(e[0].Y);
    2219                 :        272 :   av2 = avma; /* preserve inv_q and Y1 */
    2220                 :            : 
    2221         [ +  + ]:       5285 :   while (fgets(buf, MPQS_STRING_LENGTH, COMB))
    2222                 :            :   {
    2223                 :       5013 :     set_lp_entry(&e[i], buf);
    2224         [ +  + ]:       5013 :     if (e[i].q != old_q)
    2225                 :            :     {
    2226                 :            :       /* switch to combining a new bunch, swapping the rows */
    2227                 :       2338 :       old_q = e[i].q;
    2228                 :       2338 :       avma = av; /* discard old inv_q and Y1 */
    2229         [ -  + ]:       2338 :       if (!invmod(utoipos(old_q), h->N, &inv_q)) /* can happen --GN */
    2230                 :            :       {
    2231                 :          0 :         inv_q = gcdii(inv_q, h->N);
    2232 [ #  # ][ #  # ]:          0 :         if (is_pm1(inv_q) || equalii(inv_q, h->N)) /* pity */
    2233                 :            :         {
    2234                 :            : #ifdef MPQS_DEBUG
    2235                 :            :           err_printf("MPQS: skipping relation with non-invertible q\n");
    2236                 :            : #endif
    2237                 :          0 :           old_q = -1; /* sentinel */
    2238                 :          0 :           av2 = avma = av;
    2239                 :          0 :           continue; /* discard this combination */
    2240                 :            :         }
    2241                 :          0 :         *f = gerepileuptoint(av0, inv_q);
    2242                 :          0 :         return c;
    2243                 :            :       }
    2244                 :       2338 :       Y1 = strtoi(e[i].Y);
    2245                 :       2338 :       i = 1 - i; /* subsequent relations go to other row */
    2246                 :       2338 :       av2 = avma; /* preserve inv_q and Y1 */
    2247                 :       2338 :       continue;
    2248                 :            :     }
    2249                 :            :     /* count and combine the two we've got, and continue in the same row */
    2250                 :       2675 :     c++;
    2251                 :       2675 :     memset((void *)ei, 0, ei_size * sizeof(long));
    2252                 :       2675 :     mpqs_set_exponents(ei, e[0].E);
    2253                 :       2675 :     mpqs_set_exponents(ei, e[1].E);
    2254                 :       2675 :     Y2 = strtoi(e[i].Y);
    2255                 :       2675 :     new_Y = modii(mulii(mulii(Y1, Y2), inv_q), h->N);
    2256                 :       2675 :     new_Y1 = subii(h->N, new_Y);
    2257         [ +  + ]:       2675 :     if (absi_cmp(new_Y1, new_Y) < 0) new_Y = new_Y1;
    2258                 :       2675 :     strcpy(new_relation, itostr(new_Y));
    2259                 :       2675 :     strcat(new_relation, " :");
    2260         [ +  + ]:       2675 :     if (ei[1] & 1) strcat(new_relation, " 1 1");
    2261         [ +  + ]:    3157755 :     for (l = 2; l < ei_size; l++)
    2262         [ +  + ]:    3155080 :       if (ei[l])
    2263                 :            :       {
    2264                 :      58575 :         sprintf(buf, " %ld %ld", ei[l], l);
    2265                 :      58575 :         strcat(new_relation, buf);
    2266                 :            :       }
    2267                 :       2675 :     strcat(new_relation, " 0");
    2268         [ -  + ]:       2675 :     if (DEBUGLEVEL >= 6)
    2269                 :            :     {
    2270                 :          0 :       err_printf("MPQS: combining\n");
    2271                 :          0 :       err_printf("    {%ld @ %s : %s}\n", old_q, e[1-i].Y, e[1-i].E);
    2272                 :          0 :       err_printf("  * {%ld @ %s : %s}\n", e[i].q, e[i].Y, e[i].E);
    2273                 :          0 :       err_printf(" == {%s}\n", new_relation);
    2274                 :            :     }
    2275                 :       2675 :     strcat(new_relation, "\n");
    2276                 :            : 
    2277                 :            : #ifdef MPQS_DEBUG
    2278                 :            :     {
    2279                 :            :       GEN Qx_2, prod;
    2280                 :            :       char *s = strchr(new_relation, ':') + 2;
    2281                 :            :       pari_sp av1 = avma;
    2282                 :            : 
    2283                 :            :       Qx_2 = modii(sqri(new_Y), h->N);
    2284                 :            :       prod = mpqs_factorback(h, s);
    2285                 :            :       if (!equalii(Qx_2, prod))
    2286                 :            :         pari_err_BUG("MPQS: combined large prime relation is false");
    2287                 :            :       avma = av1;
    2288                 :            :     }
    2289                 :            : #endif
    2290                 :            : 
    2291                 :       2675 :     pari_fputs(new_relation, pFNEW);
    2292                 :       2675 :     avma = av2;
    2293                 :            :   } /* while */
    2294                 :            : 
    2295         [ -  + ]:        272 :   if (DEBUGLEVEL >= 4)
    2296         [ #  # ]:          0 :     err_printf("MPQS: combined %ld full relation%s\n", c, (c!=1 ? "s" : ""));
    2297                 :        272 :   avma = av0; return c;
    2298                 :            : }
    2299                 :            : 
    2300                 :            : /*********************************************************************/
    2301                 :            : /**                                                                 **/
    2302                 :            : /**                    FROM RELATIONS TO DIVISORS                   **/
    2303                 :            : /**                                                                 **/
    2304                 :            : /*********************************************************************/
    2305                 :            : 
    2306                 :            : /* create and read an F2m from a relation file FREL.
    2307                 :            :  * Also record the position of each relation in the file for later use
    2308                 :            :  * rows = size_of_FB+1, cols = rel */
    2309                 :            : static GEN
    2310                 :        385 : stream_read_F2m(pariFILE *pFREL, long rows, long cols, long *fpos)
    2311                 :            : {
    2312                 :        385 :   FILE *FREL = pFREL->file;
    2313                 :            :   long i, e, p;
    2314                 :            :   char buf[MPQS_STRING_LENGTH], *s;
    2315                 :            :   GEN m;
    2316                 :        385 :   long space = 2*((nbits2nlong(rows)+3)*cols+1);
    2317         [ -  + ]:        385 :   if ((long)((GEN)avma - (GEN)pari_mainstack->bot) < space)
    2318                 :            :   {
    2319                 :          0 :     pari_sp av = avma;
    2320                 :          0 :     m = gclone(zero_F2m(rows, cols));
    2321         [ #  # ]:          0 :     if (DEBUGLEVEL>=4)
    2322                 :          0 :       err_printf("MPQS: allocating %ld words for Gauss\n",space);
    2323                 :          0 :     avma = av;
    2324                 :            :   }
    2325                 :            :   else
    2326                 :        385 :     m = zero_F2m_copy(rows, cols);
    2327                 :        385 :   for (i = 0;; i++)
    2328                 :            :   {
    2329                 :     125937 :     char *tok=NULL;
    2330 [ +  + ][ -  + ]:     125937 :     if (i < cols && (fpos[i] = ftell(FREL)) < 0)
    2331                 :          0 :       pari_err_FILE("full relations file [ftell]", pFREL->name);
    2332         [ +  + ]:     125937 :     if (!fgets(buf, MPQS_STRING_LENGTH, FREL)) break;
    2333                 :     125552 :     s = strchr(buf, ':');
    2334         [ -  + ]:     125552 :     if (!s) pari_err_FILE("full relations file [strchr]", pFREL->name);
    2335                 :     125552 :     s = paristrtok_r(s+2, " \n", &tok);
    2336         [ +  - ]:    1807212 :     while (s != NULL)
    2337                 :            :     {
    2338         [ +  + ]:    1807212 :       e = atol(s); if (!e) break;
    2339                 :    1681660 :       s = paristrtok_r(NULL, " \n", &tok);
    2340                 :    1681660 :       p = atol(s);
    2341         [ +  + ]:    1681660 :       if (e & 1) F2m_set(m, p, i+1);
    2342                 :    1681660 :       s = paristrtok_r(NULL, " \n", &tok);
    2343                 :            :     }
    2344                 :     125552 :   }
    2345         [ -  + ]:        385 :   if (i != cols)
    2346                 :            :   {
    2347         [ #  # ]:          0 :     err_printf("MPQS: full relations file %s than expected",
    2348                 :            :                i > cols ? "longer" : "shorter");
    2349                 :          0 :     pari_err(e_BUG, "MPQS [panicking]");
    2350                 :            :   }
    2351                 :        385 :   return m;
    2352                 :            : }
    2353                 :            : 
    2354                 :            : /* NB: overwrites rel */
    2355                 :            : static GEN
    2356                 :      78902 : mpqs_add_relation(GEN Y_prod, GEN N, long *ei, char *rel)
    2357                 :            : {
    2358                 :      78902 :   pari_sp av = avma;
    2359                 :            :   GEN res;
    2360                 :      78902 :   char *s, *tok=NULL;
    2361                 :            : 
    2362                 :      78902 :   s = strchr(rel, ':') - 1;
    2363                 :      78902 :   *s = '\0';
    2364                 :            : 
    2365                 :      78902 :   res = remii(mulii(Y_prod, strtoi(rel)), N);
    2366                 :            : 
    2367                 :      78902 :   s = paristrtok_r(s + 3, " \n", &tok);
    2368         [ +  - ]:    1141239 :   while (s != NULL)
    2369                 :            :   {
    2370                 :    1141239 :     long e = atol(s), i;
    2371         [ +  + ]:    1141239 :     if (!e) break;
    2372                 :    1062337 :     s = paristrtok_r(NULL, " \n", &tok);
    2373                 :    1062337 :     i = atol(s); /* bug in g++-3.4.1: miscompiles ei[ atol(s) ] */
    2374                 :    1062337 :     ei[i] += e;
    2375                 :    1062337 :     s = paristrtok_r(NULL, " \n", &tok);
    2376                 :            :   }
    2377                 :      78902 :   return gerepileuptoint(av, res);
    2378                 :            : }
    2379                 :            : 
    2380                 :            : static char*
    2381                 :      78902 : mpqs_get_relation(char *buf, long pos, pariFILE *pFREL)
    2382                 :            : {
    2383         [ -  + ]:      78902 :   if (fseek(pFREL->file, pos, SEEK_SET))
    2384                 :          0 :     pari_err_FILE("FREL file [fseek]", pFREL->name);
    2385         [ -  + ]:      78902 :   if (!fgets(buf, MPQS_STRING_LENGTH, pFREL->file))
    2386                 :          0 :     pari_err_FILE("FREL file [fgets]", pFREL->name);
    2387                 :      78902 :   return buf;
    2388                 :            : }
    2389                 :            : 
    2390                 :            : #define isprobableprime(n) (MR_Jaeschke((n),17))
    2391                 :            : static int
    2392                 :        880 : split(GEN N, GEN *e, GEN *res)
    2393                 :            : {
    2394                 :            :   ulong mask;
    2395                 :            :   long flag;
    2396                 :            :   GEN base;
    2397         [ +  + ]:        880 :   if (isprobableprime(N)) { *e = gen_1; return 1; }
    2398         [ -  + ]:         55 :   if (Z_issquareall(N, &base))
    2399                 :            :   { /* squares could cost us a lot of time */
    2400                 :            :     /* GN20050707: as used now, this is always called with res!=NULL */
    2401                 :          0 :     *res = base;
    2402                 :          0 :     *e = gen_2;
    2403         [ #  # ]:          0 :     if (DEBUGLEVEL >= 5) err_printf("MPQS: decomposed a square\n");
    2404                 :          0 :     return 1;
    2405                 :            :   }
    2406                 :         55 :   mask = 7;
    2407                 :            :   /* 5th/7th powers aren't worth the trouble. OTOH once we have the hooks for
    2408                 :            :    * dealing with cubes, higher powers can be handled essentially for free) */
    2409         [ -  + ]:         55 :   if ( (flag = is_357_power(N, &base, &mask)) )
    2410                 :            :   {
    2411                 :          0 :     *res = base;
    2412                 :          0 :     *e = utoipos(flag);
    2413         [ #  # ]:          0 :     if (DEBUGLEVEL >= 5)
    2414         [ #  # ]:          0 :       err_printf("MPQS: decomposed a %s\n",
    2415                 :            :                  (flag == 3 ? "cube" :
    2416         [ #  # ]:          0 :                   (flag == 5 ? "5th power" : "7th power")));
    2417                 :          0 :     return 1;
    2418                 :            :   }
    2419                 :        880 :   *e = gen_0; return 0; /* known composite */
    2420                 :            : }
    2421                 :            : 
    2422                 :            : static GEN
    2423                 :        385 : mpqs_solve_linear_system(mpqs_handle_t *h, pariFILE *pFREL, long rel)
    2424                 :            : {
    2425                 :        385 :   GEN N = h->N, X, Y_prod, X_plus_Y, D1, res, new_res;
    2426                 :        385 :   mpqs_FB_entry_t *FB = h->FB;
    2427                 :        385 :   pari_sp av=avma, av2, av3, lim, lim3;
    2428                 :            :   long *fpos, *ei;
    2429                 :            :   long i, j, H_cols, H_rows;
    2430                 :            :   long res_last, res_next, res_size, res_max;
    2431                 :            :   GEN  m, ker_m;
    2432                 :            :   long done, rank;
    2433                 :            :   char buf[MPQS_STRING_LENGTH];
    2434                 :            : 
    2435                 :        385 :   fpos = (long *) pari_malloc(rel * sizeof(long));
    2436                 :            : 
    2437                 :        385 :   m = stream_read_F2m(pFREL, h->size_of_FB+1, rel, fpos);
    2438         [ -  + ]:        385 :   if (DEBUGLEVEL >= 7)
    2439                 :          0 :     err_printf("\\\\ MATRIX READ BY MPQS\nFREL=%Ps\n",m);
    2440                 :            : 
    2441                 :        385 :   ker_m = F2m_ker_sp(m,0); rank = lg(ker_m)-1;
    2442         [ -  + ]:        385 :   if (isclone(m)) gunclone(m);
    2443                 :            : 
    2444         [ -  + ]:        385 :   if (DEBUGLEVEL >= 4)
    2445                 :            :   {
    2446         [ #  # ]:          0 :     if (DEBUGLEVEL >= 7)
    2447                 :            :     {
    2448                 :          0 :       err_printf("\\\\ KERNEL COMPUTED BY MPQS\n");
    2449                 :          0 :       err_printf("KERNEL=%Ps\n",ker_m);
    2450                 :            :     }
    2451                 :          0 :     err_printf("MPQS: Gauss done: kernel has rank %ld, taking gcds...\n", rank);
    2452                 :            :   }
    2453                 :            : 
    2454                 :        385 :   H_rows = rel;
    2455                 :        385 :   H_cols = rank;
    2456                 :            : 
    2457         [ -  + ]:        385 :   if (!H_cols)
    2458                 :            :   { /* trivial kernel. Fail gracefully: main loop may look for more relations */
    2459         [ #  # ]:          0 :     if (DEBUGLEVEL >= 3)
    2460                 :          0 :       pari_warn(warner, "MPQS: no solutions found from linear system solver");
    2461                 :          0 :     pari_free(fpos); /* ei not yet allocated */
    2462                 :          0 :     avma = av; return NULL; /* no factors found */
    2463                 :            :   }
    2464                 :            : 
    2465                 :            :   /* If the rank is r, we can expect up to 2^r pairwise coprime factors,
    2466                 :            :    * but it may happen that a kernel basis vector contributes nothing new to
    2467                 :            :    * the decomposition.  We allocate room for up to eight factors initially
    2468                 :            :    * (certainly adequate when one or two basis vectors work), adjusting this
    2469                 :            :    * down at the end to what we actually found, or up if we are very lucky and
    2470                 :            :    * find more factors.  In the upper half of our vector, we store information
    2471                 :            :    * about which factors we know to be composite (zero) or believe to be
    2472                 :            :    * composite (NULL) or suspect to be prime (one), or an exponent (two
    2473                 :            :    * or some t_INT) if it is a proper power */
    2474                 :        385 :   av2 = avma; lim = stack_lim(av2,1);
    2475         [ +  + ]:        385 :   if (rank > (long)BITS_IN_LONG - 2)
    2476                 :        269 :     res_max = LONG_MAX; /* the common case, unfortunately */
    2477                 :            :   else
    2478                 :        116 :     res_max = 1L<<rank; /* max number of factors we can hope for */
    2479                 :        385 :   res_size = 8; /* no. of factors we can store in this res */
    2480                 :        385 :   res = cgetg(2*res_size+1, t_VEC);
    2481         [ +  + ]:       6545 :   for (i=2*res_size; i; i--) res[i] = 0;
    2482                 :        385 :   res_next = res_last = 1;
    2483                 :            : 
    2484                 :        385 :   ei = (long *) pari_malloc((h->size_of_FB + 2) * sizeof(long));
    2485                 :            : 
    2486         [ +  - ]:       1017 :   for (i = 1; i <= H_cols; i++)
    2487                 :            :   { /* loop over kernel basis */
    2488                 :       1017 :     X = Y_prod = gen_1;
    2489                 :       1017 :     memset((void *)ei, 0, (h->size_of_FB + 2) * sizeof(long));
    2490                 :            : 
    2491                 :       1017 :     av3 = avma; lim3 = stack_lim(av3,1);
    2492         [ +  + ]:     460559 :     for (j = 1; j <= H_rows; j++)
    2493                 :            :     {
    2494         [ +  + ]:     459542 :       if (F2m_coeff(ker_m, j, i))
    2495                 :      78902 :         Y_prod = mpqs_add_relation(Y_prod, N, ei,
    2496                 :      78902 :                                    mpqs_get_relation(buf, fpos[j-1], pFREL));
    2497         [ -  + ]:     459542 :       if (low_stack(lim3, stack_lim(av3,1)))
    2498                 :            :       {
    2499         [ #  # ]:          0 :         if(DEBUGMEM>1) pari_warn(warnmem,"[1]: mpqs_solve_linear_system");
    2500                 :          0 :         Y_prod = gerepileuptoint(av3, Y_prod);
    2501                 :            :       }
    2502                 :            :     }
    2503                 :       1017 :     Y_prod = gerepileuptoint(av3, Y_prod);
    2504                 :            : 
    2505                 :       1017 :     av3 = avma; lim3 = stack_lim(av3,1);
    2506         [ +  + ]:     396999 :     for (j = 2; j <= h->size_of_FB + 1; j++)
    2507         [ +  + ]:     395982 :       if (ei[j])
    2508                 :            :       {
    2509         [ -  + ]:     128835 :         if (ei[j] & 1) pari_err_BUG("MPQS (relation is a nonsquare)");
    2510                 :     128835 :         X = remii(mulii(X,
    2511                 :     257670 :                         Fp_powu(utoipos(FB[j].fbe_p), (ulong)ei[j]>>1, N)),
    2512                 :            :                   N);
    2513         [ -  + ]:     128835 :         if (low_stack(lim3, stack_lim(av3,1)))
    2514                 :            :         {
    2515         [ #  # ]:          0 :           if(DEBUGMEM>1) pari_warn(warnmem,"[2]: mpqs_solve_linear_system");
    2516                 :          0 :           X = gerepileupto(av3, X);
    2517                 :            :         }
    2518                 :            :       }
    2519                 :       1017 :     X = gerepileuptoint(av3, X);
    2520         [ -  + ]:       1017 :     if (MPQS_DEBUGLEVEL >= 1)
    2521                 :            :     {
    2522         [ #  # ]:          0 :       if (signe(remii(subii(sqri(X), sqri(Y_prod)), N)))
    2523                 :            :       { /* shouldn't happen */
    2524                 :          0 :         err_printf("MPQS: X^2 - Y^2 != 0 mod N\n");
    2525                 :          0 :         err_printf("\tindex i = %ld\n", i);
    2526                 :          0 :         pari_warn(warner, "MPQS: wrong relation found after Gauss");
    2527                 :            :       }
    2528                 :            :     }
    2529                 :            :     /* At this point, X^2 == Y^2 mod N.  Indeed, something stronger is true:
    2530                 :            :      * We have gcd(X-Y, N) * gcd(X+Y, N) == N.  Why?
    2531                 :            :      * First, N divides X^2 - Y^2, so it divides the lefthand side.
    2532                 :            :      * Second, let P be any prime factor of N.  If P were to divide both
    2533                 :            :      * X-Y and X+Y, then it would divide their sum 2X.  But X (in the present
    2534                 :            :      * backwards notation!) is a product of powers of FB primes, and no FB
    2535                 :            :      * prime is a divisor of N, or we would have found out about it already
    2536                 :            :      * whilst constructing the FB.
    2537                 :            :      * Therefore in the following it is sufficient to work with gcd(X+Y, N)
    2538                 :            :      * alone, and never look at gcd(X-Y, N).
    2539                 :            :      */
    2540                 :       1017 :     done = 0; /* (re-)counts probably-prime factors (or powers whose bases we
    2541                 :            :                * don't want to handle any further) */
    2542                 :       1017 :     X_plus_Y = addii(X, Y_prod);
    2543         [ +  + ]:       1017 :     if (res_next < 3)
    2544                 :            :     { /* we still haven't decomposed the original N, and want both a gcd and
    2545                 :            :        * its cofactor. */
    2546                 :        906 :       D1 = gcdii(X_plus_Y, N);
    2547 [ +  + ][ +  + ]:        906 :       if (is_pm1(D1) || equalii(D1,N)) { avma = av3; continue; }
    2548                 :            :       /* got something that works */
    2549         [ -  + ]:        385 :       if (DEBUGLEVEL >= 5)
    2550         [ #  # ]:          0 :         err_printf("MPQS: splitting N after %ld kernel vector%s\n",
    2551                 :            :                    i+1, (i? "s" : ""));
    2552                 :            :       /* GN20050707 Fixed:
    2553                 :            :        * Don't divide N in place. We still need it for future X and Y_prod
    2554                 :            :        * computations! */
    2555                 :        385 :       gel(res,1) = diviiexact(N, D1);
    2556                 :        385 :       gel(res,2) = D1;
    2557                 :        385 :       res_last = res_next = 3;
    2558                 :            : 
    2559         [ +  + ]:        385 :       if ( split(gel(res,1),  &gel(res,res_size+1), &gel(res,1)) ) done++;
    2560         [ +  + ]:        385 :       if ( split(D1, &gel(res,res_size+2), &gel(res,2)) ) done++;
    2561         [ +  + ]:        385 :       if (done == 2) break;     /* both factors look prime or were powers */
    2562                 :            :       /* GN20050707: moved following line down to here, was before the
    2563                 :            :        * two split() invocations.  Very rare case anyway. */
    2564         [ -  + ]:         55 :       if (res_max == 2) break; /* two out of two possible factors seen */
    2565         [ -  + ]:         55 :       if (DEBUGLEVEL >= 5)
    2566                 :          0 :         err_printf("MPQS: got two factors, looking for more...\n");
    2567                 :            :     }
    2568                 :            :     else
    2569                 :            :     { /* we already have factors */
    2570         [ +  + ]:        388 :       for (j=1; j < res_next; j++)
    2571                 :            :       { /* loop over known-composite factors */
    2572 [ +  - ][ +  + ]:        277 :         if (gel(res,res_size+j) && gel(res,res_size+j) != gen_0)
    2573                 :            :         {
    2574                 :        166 :           done++; continue; /* skip probable primes etc */
    2575                 :            :         }
    2576                 :            :         /* actually, also skip square roots of squares etc.  They are a lot
    2577                 :            :          * smaller than the original N, and should be easy to deal with later */
    2578                 :        111 :         av3 = avma;
    2579                 :        111 :         D1 = gcdii(X_plus_Y, gel(res,j));
    2580 [ +  + ][ +  + ]:        111 :         if (is_pm1(D1) || equalii(D1, gel(res,j))) { avma = av3; continue; }
    2581                 :            :         /* got one which splits this factor */
    2582         [ -  + ]:         55 :         if (DEBUGLEVEL >= 5)
    2583                 :          0 :           err_printf("MPQS: resplitting a factor after %ld kernel vectors\n",
    2584                 :            :                      i+1);      /* always plural */
    2585                 :            :         /* first make sure there's room for another factor */
    2586         [ -  + ]:         55 :         if (res_next > res_size)
    2587                 :            :         { /* need to reallocate (_very_ rare case) */
    2588                 :          0 :           long i1, size = 2*res_size;
    2589                 :            :           GEN RES;
    2590         [ #  # ]:          0 :           if (size > res_max) size = res_max;
    2591                 :          0 :           RES = cgetg(2*size+1, t_VEC);
    2592         [ #  # ]:          0 :           for (i1=2*size; i1>=res_next; i1--) gel(RES,i1) = NULL;
    2593         [ #  # ]:          0 :           for (i1=1; i1<res_next; i1++)
    2594                 :            :           {
    2595                 :            :             /* GN20050707:
    2596                 :            :              * on-stack contents of RES must be rejuvenated */
    2597 [ #  # ][ #  # ]:          0 :             icopyifstack(gel(res,i1), gel(RES,i1)); /* factors */
    2598         [ #  # ]:          0 :             if ( gel(res,res_size+i1) )
    2599 [ #  # ][ #  # ]:          0 :               icopyifstack(gel(res,res_size+i1), gel(RES,size+i1));
    2600                 :            :             /* primality tags */
    2601                 :            :           }
    2602                 :          0 :           res = RES; res_size = size;   /* res_next unchanged */
    2603                 :            :         }
    2604                 :            :         /* now there is room; divide into existing factor and store the
    2605                 :            :            new gcd */
    2606                 :         55 :         diviiz(gel(res,j), D1, gel(res,j));
    2607                 :         55 :         gel(res,res_next) = D1;
    2608                 :            : 
    2609                 :            :         /* following overwrites the old known-composite indication at j */
    2610         [ +  - ]:         55 :         if (split( gel(res,j), &gel(res,res_size+j), &gel(res,j)) ) done++;
    2611                 :            :         /* GN20050707 Fixed:
    2612                 :            :          * Don't increment done when the newly stored factor seems to be
    2613                 :            :          * prime or otherwise devoid of interest - this happens later
    2614                 :            :          * when we routinely revisit it during the present inner loop. */
    2615                 :         55 :         (void)split(D1, &gel(res,res_size+res_next), &gel(res,res_next));
    2616                 :            : 
    2617                 :            :         /* GN20050707: Following line moved down to here, was before the
    2618                 :            :          * two split() invocations. */
    2619         [ -  + ]:         55 :         if (++res_next > res_max)
    2620                 :            :         {
    2621                 :            :           /* all possible factors seen, outer loop postprocessing will
    2622                 :            :            * proceed to break out of the outer loop below. */
    2623                 :          0 :           break;
    2624                 :            :         }
    2625                 :            :       }       /* loop over known composite factors */
    2626                 :            : 
    2627         [ +  + ]:        111 :       if (res_next > res_last)
    2628                 :            :       {
    2629                 :         55 :         res_last = res_next - 1; /* we might have resplit more than one */
    2630         [ -  + ]:         55 :         if (DEBUGLEVEL >= 5)
    2631         [ #  # ]:          0 :           err_printf("MPQS: got %ld factors%s\n", res_last,
    2632                 :            :                      (done < res_last ? ", looking for more..." : ""));
    2633                 :         55 :         res_last = res_next;
    2634                 :            :       }
    2635                 :            :       /* break out of the outer loop when we have seen res_max factors, and
    2636                 :            :        * also when all current factors are probable primes */
    2637 [ +  - ][ +  + ]:        111 :       if (res_next > res_max || done == res_next - 1) break;
    2638                 :            :     } /* end case of further splitting of existing factors */
    2639         [ -  + ]:        111 :     if (low_stack(lim, stack_lim(av2,1)))
    2640                 :            :     {
    2641                 :            :       long i1;
    2642         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"[3]: mpqs_solve_linear_system");
    2643                 :            :       /* gcopy would have a problem with our NULL pointers... */
    2644                 :          0 :       new_res = cgetg(lg(res), t_VEC);
    2645         [ #  # ]:          0 :       for (i1=2*res_size; i1>=res_next; i1--) new_res[i1] = 0;
    2646         [ #  # ]:          0 :       for (i1=1; i1<res_next; i1++)
    2647                 :            :       {
    2648 [ #  # ][ #  # ]:          0 :         icopyifstack(gel(res,i1), gel(new_res,i1));
    2649                 :            :         /* GN20050707: the marker GENs might need rejuvenating, too */
    2650         [ #  # ]:          0 :         if (gel(res,res_size+i1))
    2651 [ #  # ][ #  # ]:          0 :           icopyifstack(gel(res,res_size+i1), gel(new_res,res_size+i1));
    2652                 :            :       }
    2653                 :          0 :       res = gerepileupto(av2, new_res);
    2654                 :            :     }
    2655                 :            :   } /* for (loop over kernel basis) */
    2656                 :            : 
    2657                 :        385 :   pari_free(ei); pari_free(fpos);
    2658         [ -  + ]:        385 :   if (res_next < 3) { avma = av; return NULL; } /* no factors found */
    2659                 :            : 
    2660                 :            :   /* normal case:  convert internal format to ifac format as described in
    2661                 :            :    * src/basemath/ifactor1.c  (triples of components: value, exponent, class
    2662                 :            :    * [unknown, known composite, known prime]),  clean up and return result */
    2663                 :        385 :   res_last = res_next - 1; /* number of distinct factors found */
    2664                 :        385 :   new_res = cgetg(3*res_last + 1, t_VEC);
    2665         [ -  + ]:        385 :   if (DEBUGLEVEL >= 6)
    2666                 :          0 :     err_printf("MPQS: wrapping up vector of %ld factors\n", res_last);
    2667         [ +  + ]:       1210 :   for (i=1,j=1; i <= res_last; i++)
    2668                 :            :   {
    2669                 :        825 :     GEN F = gel(res, res_size+i);
    2670 [ +  - ][ +  - ]:        825 :     icopyifstack(gel(res,i), gel(new_res,j++)); /* factor */
    2671                 :       1650 :     gel(new_res,j++) = /* exponent */
    2672                 :        825 :       F ? (F == gen_0 ? gen_1
    2673 [ +  - ][ -  + ]:        825 :                       : (isonstack(F) ? icopy(F) : F))
    2674         [ +  - ]:       1650 :         : gen_1; /* F was NULL */
    2675                 :       1650 :     gel(new_res,j++) = /* class */
    2676         [ -  + ]:        825 :       F == gen_0 ? gen_0 :      /* known composite */
    2677                 :            :         NULL;           /* base of power or suspected prime --
    2678                 :            :                                    mark as `unknown' */
    2679         [ -  + ]:        825 :     if (DEBUGLEVEL >= 6)
    2680         [ #  # ]:          0 :       err_printf("\tpackaging %ld: %Ps ^%ld (%s)\n", i, res[i],
    2681                 :          0 :                  itos(gel(new_res,j-2)), (F == gen_0 ? "comp." : "unknown"));
    2682                 :            :   }
    2683                 :        385 :   return gerepileupto(av, new_res);
    2684                 :            : }
    2685                 :            : 
    2686                 :            : /*********************************************************************/
    2687                 :            : /**                                                                 **/
    2688                 :            : /**               MAIN ENTRY POINT AND DRIVER ROUTINE               **/
    2689                 :            : /**                                                                 **/
    2690                 :            : /*********************************************************************/
    2691                 :            : 
    2692                 :            : 
    2693                 :            : /* All percentages below are actually fixed-point quantities scaled by 10
    2694                 :            :  * (value of 1 means 0.1%, 1000 means 100%) */
    2695                 :            : 
    2696                 :            : /* Factors N using the self-initializing multipolynomial quadratic sieve
    2697                 :            :  * (SIMPQS).  Returns one of the two factors, or (usually) a vector of factors
    2698                 :            :  * and exponents and information about which ones are still composite, or NULL
    2699                 :            :  * when something goes wrong or when we can't seem to make any headway. */
    2700                 :            : 
    2701                 :            : /* TODO: this function to be renamed mpqs_main() with several extra parameters,
    2702                 :            :  * with mpqs() as a wrapper for the standard case, so we can do partial runs
    2703                 :            :  * across several machines etc.  (from gp or a dedicated C program). --GN */
    2704                 :            : static GEN
    2705                 :        385 : mpqs_i(mpqs_handle_t *handle)
    2706                 :            : {
    2707                 :        385 :   GEN N = handle->N, fact; /* will in the end hold our factor(s) */
    2708                 :            :   mpqs_int32_t size_of_FB; /* size of the factor base */
    2709                 :            :   mpqs_FB_entry_t *FB; /* factor base */
    2710                 :            :   mpqs_int32_t M;               /* sieve interval size [-M, M] */
    2711                 :            : 
    2712                 :            :   /* local loop / auxiliary vars */
    2713                 :            :   ulong p;
    2714                 :            : 
    2715                 :            :   /* already exists in the handle, keep for convenience */
    2716                 :            :   long lp_bound;                /* size limit for large primes */
    2717                 :            :   long lp_scale;                /* ...relative to largest FB prime */
    2718                 :            : 
    2719                 :            :   /* bookkeeping */
    2720                 :            :   long tc;                      /* # of candidates found in one iteration */
    2721                 :            :   long tp;                      /* # of recently sorted LP rels */
    2722                 :        385 :   long tff = 0;                 /* # recently found full rels from sieving */
    2723                 :            :   long tfc;                     /* # full rels recently combined from LPs */
    2724                 :        385 :   double tfc_ratio = 0;         /* recent (tfc + tff) / tff */
    2725                 :            :   ulong sort_interval;          /* determine when to sort and merge */
    2726                 :            :   ulong followup_sort_interval; /* temporary files (scaled percentages) */
    2727                 :        385 :   long percentage = 0;          /* scaled by 10, see comment above */
    2728                 :            :   double net_yield;
    2729                 :        385 :   long total_full_relations = 0, total_partial_relations = 0, total_no_cand = 0;
    2730                 :        385 :   long vain_iterations = 0, good_iterations = 0, iterations = 0;
    2731                 :            : #ifdef MPQS_USE_HISTOGRAMS
    2732                 :            :   long histo_checkpoint = MPQS_MIN_CANDS_FOR_HISTO;
    2733                 :            : #endif
    2734                 :            : 
    2735                 :            :   pariFILE *pFNEW, *pLPNEW, *pCOMB, *pFREL, *pLPREL;
    2736                 :            :   char *dir, *COMB_str, *FREL_str, *FNEW_str, *LPREL_str, *LPNEW_str, *TMP_str;
    2737                 :            :   pari_timer T;
    2738                 :            : 
    2739                 :            : /* END: global variables to disappear as soon as possible */
    2740                 :            : 
    2741                 :            : /******************************/
    2742                 :            : 
    2743                 :        385 :   pari_sp av = avma;
    2744                 :            : 
    2745         [ -  + ]:        385 :   if (DEBUGLEVEL >= 4)
    2746                 :            :   {
    2747                 :          0 :     timer_start(&T);
    2748                 :          0 :     err_printf("MPQS: number to factor N = %Ps\n", N);
    2749                 :            :   }
    2750                 :            : 
    2751                 :        385 :   handle->digit_size_N = decimal_len(N);
    2752         [ -  + ]:        385 :   if (handle->digit_size_N > MPQS_MAX_DIGIT_SIZE_KN)
    2753                 :            :   {
    2754                 :          0 :     pari_warn(warner, "MPQS: number too big to be factored with MPQS,\n\tgiving up");
    2755                 :          0 :     return NULL;
    2756                 :            :   }
    2757                 :            : 
    2758         [ -  + ]:        385 :   if (DEBUGLEVEL >= 4)
    2759                 :          0 :     err_printf("MPQS: factoring number of %ld decimal digits\n",
    2760                 :            :                handle->digit_size_N);
    2761                 :            : 
    2762                 :        385 :   p = mpqs_find_k(handle);
    2763         [ -  + ]:        385 :   if (p) { avma = av; return utoipos(p); }
    2764         [ -  + ]:        385 :   if (DEBUGLEVEL >= 5) err_printf("MPQS: found multiplier %ld for N\n",
    2765                 :          0 :                                   handle->_k->k);
    2766                 :        385 :   handle->kN = muliu(N, handle->_k->k);
    2767                 :            : 
    2768         [ -  + ]:        385 :   if (!mpqs_set_parameters(handle))
    2769                 :            :   {
    2770                 :          0 :     pari_warn(warner,
    2771                 :            :         "MPQS: number too big to be factored with MPQS,\n\tgiving up");
    2772                 :          0 :     return NULL;
    2773                 :            :   }
    2774                 :            : 
    2775                 :        385 :   size_of_FB = handle->size_of_FB;
    2776                 :        385 :   M = handle->M;
    2777                 :        385 :   sort_interval = handle->first_sort_point;
    2778                 :        385 :   followup_sort_interval = handle->sort_pt_interval;
    2779                 :            : 
    2780         [ -  + ]:        385 :   if (DEBUGLEVEL >= 5)
    2781                 :          0 :     err_printf("MPQS: creating factor base and allocating arrays...\n");
    2782                 :        385 :   FB = mpqs_create_FB(handle, &p);
    2783         [ -  + ]:        385 :   if (p) { avma = av; return utoipos(p); }
    2784                 :        385 :   mpqs_sieve_array_ctor(handle);
    2785                 :        385 :   mpqs_poly_ctor(handle);
    2786                 :            : 
    2787                 :        385 :   lp_bound = handle->largest_FB_p;
    2788         [ -  + ]:        385 :   if (lp_bound > MPQS_LP_BOUND) lp_bound = MPQS_LP_BOUND;
    2789                 :            :   /* don't allow large primes to have room for two factors both bigger than
    2790                 :            :    * what the FB contains (...yet!) */
    2791                 :        385 :   lp_scale = handle->lp_scale;
    2792         [ -  + ]:        385 :   if (lp_scale >= handle->largest_FB_p)
    2793                 :          0 :     lp_scale = handle->largest_FB_p - 1;
    2794                 :        385 :   lp_bound *= lp_scale;
    2795                 :        385 :   handle->lp_bound = lp_bound;
    2796                 :            : 
    2797                 :        385 :   handle->dkN = gtodouble(handle->kN);
    2798                 :            :   /* compute the threshold and fill in the byte-sized scaled logarithms */
    2799                 :        385 :   mpqs_set_sieve_threshold(handle);
    2800                 :            : 
    2801         [ -  + ]:        385 :   if (!mpqs_locate_A_range(handle)) return NULL;
    2802                 :            : 
    2803         [ -  + ]:        385 :   if (DEBUGLEVEL >= 4)
    2804                 :            :   {
    2805                 :          0 :     err_printf("MPQS: sieving interval = [%ld, %ld]\n", -(long)M, (long)M);
    2806                 :            :     /* that was a little white lie, we stop one position short at the top */
    2807                 :          0 :     err_printf("MPQS: size of factor base = %ld\n",
    2808                 :            :                (long)size_of_FB);
    2809                 :          0 :     err_printf("MPQS: striving for %ld relations\n",
    2810                 :          0 :                (long)handle->target_no_rels);
    2811                 :          0 :     err_printf("MPQS: coefficients A will be built from %ld primes each\n",
    2812                 :          0 :                (long)handle->omega_A);
    2813                 :          0 :     err_printf("MPQS: primes for A to be chosen near FB[%ld] = %ld\n",
    2814                 :          0 :                (long)handle->index2_FB,
    2815                 :          0 :                (long)FB[handle->index2_FB].fbe_p);
    2816                 :          0 :     err_printf("MPQS: smallest prime used for sieving FB[%ld] = %ld\n",
    2817                 :          0 :                (long)handle->index1_FB,
    2818                 :          0 :                (long)FB[handle->index1_FB].fbe_p);
    2819                 :          0 :     err_printf("MPQS: largest prime in FB = %ld\n",
    2820                 :          0 :                (long)handle->largest_FB_p);
    2821                 :          0 :     err_printf("MPQS: bound for `large primes' = %ld\n", (long)lp_bound);
    2822                 :            :   }
    2823                 :            : 
    2824         [ -  + ]:        385 :   if (DEBUGLEVEL >= 5)
    2825                 :            :   {
    2826                 :          0 :     err_printf("MPQS: sieve threshold = %u\n",
    2827                 :          0 :                (unsigned int)handle->sieve_threshold);
    2828                 :            :   }
    2829                 :            : 
    2830         [ -  + ]:        385 :   if (DEBUGLEVEL >= 4)
    2831                 :            :   {
    2832                 :          0 :     err_printf("MPQS: first sorting at %ld%%, then every %3.1f%% / %3.1f%%\n",
    2833                 :            :                sort_interval/10, followup_sort_interval/10.,
    2834                 :            :                followup_sort_interval/20.);
    2835                 :            :   }
    2836                 :            : 
    2837                 :            : 
    2838                 :            :   /* main loop which
    2839                 :            :    * - computes polynomials and their zeros (SI)
    2840                 :            :    * - does the sieving
    2841                 :            :    * - tests candidates of the sieve array */
    2842                 :            : 
    2843                 :            :   /* Let (A, B_i) the current pair of coeffs. If i == 0 a new A is generated */
    2844                 :        385 :   handle->index_j = (mpqs_uint32_t)-1;  /* increment below will have it start at 0 */
    2845                 :            : 
    2846         [ -  + ]:        385 :   if (DEBUGLEVEL >= 5) err_printf("MPQS: starting main loop\n");
    2847                 :            :   /* compute names for the temp files we'll need */
    2848                 :        385 :   dir = pari_unique_dir("MPQS");
    2849                 :        385 :   TMP_str   = mpqs_get_filename(dir, "LPTMP");
    2850                 :        385 :   FREL_str  = mpqs_get_filename(dir, "FREL");
    2851                 :        385 :   FNEW_str  = mpqs_get_filename(dir, "FNEW");
    2852                 :        385 :   LPREL_str = mpqs_get_filename(dir, "LPREL");
    2853                 :        385 :   LPNEW_str = mpqs_get_filename(dir, "LPNEW");
    2854                 :        385 :   COMB_str  = mpqs_get_filename(dir, "COMB");
    2855                 :            : #define unlink_all()\
    2856                 :            :       pari_unlink(FREL_str);\
    2857                 :            :       pari_unlink(FNEW_str);\
    2858                 :            :       pari_unlink(LPREL_str);\
    2859                 :            :       pari_unlink(LPNEW_str);\
    2860                 :            :       if (pCOMB) pari_unlink(COMB_str);\
    2861                 :            :       rmdir(dir); pari_free(dir);
    2862                 :            : 
    2863                 :        385 :   pFREL = pari_fopen_or_fail(FREL_str,  WRITE); pari_fclose(pFREL);
    2864                 :        385 :   pLPREL = pari_fopen_or_fail(LPREL_str,  WRITE); pari_fclose(pLPREL);
    2865                 :        385 :   pFNEW = pari_fopen_or_fail(FNEW_str,  WRITE);
    2866                 :        385 :   pLPNEW= pari_fopen_or_fail(LPNEW_str, WRITE);
    2867                 :        385 :   pCOMB = NULL;
    2868                 :            : 
    2869                 :            :   for(;;)
    2870                 :            :   { /* FNEW and LPNEW are open for writing */
    2871                 :      25375 :     iterations++;
    2872                 :            :     /* self initialization: compute polynomial and its zeros */
    2873                 :      25375 :     mpqs_self_init(handle);
    2874         [ -  + ]:      25375 :     if (handle->bin_index == 0)
    2875                 :            :     { /* have run out of primes for A */
    2876                 :            :       /* We might change some parameters.  For the moment, simply give up */
    2877         [ #  # ]:          0 :       if (DEBUGLEVEL >= 2)
    2878                 :          0 :         err_printf("MPQS: Ran out of primes for A, giving up.\n");
    2879                 :          0 :       pari_fclose(pFNEW);
    2880                 :          0 :       pari_fclose(pLPNEW);
    2881                 :            :       /* FREL, LPREL are closed at this point */
    2882         [ #  # ]:          0 :       unlink_all(); avma = av; return NULL;
    2883                 :            :     }
    2884                 :            : 
    2885                 :      25375 :     memset((void*)(handle->sieve_array), 0, (M << 1) * sizeof(unsigned char));
    2886                 :      25375 :     mpqs_sieve(handle);
    2887                 :            : 
    2888                 :      25375 :     tc = mpqs_eval_sieve(handle);
    2889                 :      25375 :     total_no_cand += tc;
    2890         [ -  + ]:      25375 :     if (DEBUGLEVEL >= 6)
    2891         [ #  # ]:          0 :       err_printf("MPQS: found %lu candidate%s\n", tc, (tc==1? "" : "s"));
    2892                 :            : 
    2893         [ +  + ]:      25375 :     if (tc)
    2894                 :            :     {
    2895                 :      25185 :       long t = mpqs_eval_cand(handle, tc, pFNEW->file, pLPNEW->file);
    2896                 :      25185 :       total_full_relations += t;
    2897                 :      25185 :       tff += t;
    2898                 :      25185 :       good_iterations++;
    2899                 :            :     }
    2900                 :            : 
    2901                 :            : #ifdef MPQS_USE_HISTOGRAMS
    2902                 :            :     if (handle->do_histograms && !handle->done_histograms &&
    2903                 :            :         total_no_cand >= histo_checkpoint)
    2904                 :            :     {
    2905                 :            :       int res = mpqs_eval_histograms(handle);
    2906                 :            :       if (res >= 0)
    2907                 :            :       {
    2908                 :            :         /* retry later */
    2909                 :            :         if (res > 0)
    2910                 :            :           /* histo_checkpoint *= 2.6; */
    2911                 :            :           handle->do_histograms = 0; /* no, don't retry later */
    2912                 :            :         else
    2913                 :            :           histo_checkpoint += (MPQS_MIN_CANDS_FOR_HISTO /* >> 1 */);
    2914                 :            :       }
    2915                 :            :       else
    2916                 :            :         handle->done_histograms = 1;
    2917                 :            :     }
    2918                 :            : #endif
    2919                 :            : 
    2920                 :      25375 :     percentage =
    2921                 :      25375 :       (long)((1000.0 * total_full_relations) / handle->target_no_rels);
    2922                 :            : 
    2923         [ +  + ]:      25375 :     if ((ulong)percentage < sort_interval) continue;
    2924                 :            :     /* most main loops continue here! */
    2925                 :            : 
    2926                 :            :     /* Extra processing when we have completed a sort interval: */
    2927         [ -  + ]:       3233 :     if (DEBUGLEVEL >= 3)
    2928                 :            :     {
    2929         [ #  # ]:          0 :       if (DEBUGLEVEL >= 4)
    2930                 :          0 :         err_printf("\nMPQS: passing the %3.1f%% sort point, time = %ld ms\n",
    2931                 :            :                    sort_interval/10., timer_delay(&T));
    2932                 :            :       else
    2933                 :          0 :         err_printf("\nMPQS: passing the %3.1f%% sort point\n",
    2934                 :            :                    sort_interval/10.);
    2935                 :          0 :       err_flush();
    2936                 :            :     }
    2937                 :            : 
    2938                 :            :     /* sort LPNEW and merge it into LPREL, diverting combinables into COMB */
    2939                 :       3233 :     pari_fclose(pLPNEW);
    2940                 :       3233 :     (void)mpqs_sort_lp_file(LPNEW_str);
    2941                 :       3233 :     pCOMB = pari_fopen_or_fail(COMB_str, WRITE);
    2942                 :       3233 :     tp = mpqs_mergesort_lp_file(LPREL_str, LPNEW_str, TMP_str, pCOMB);
    2943                 :       3233 :     pari_fclose(pCOMB);
    2944                 :       3233 :     pLPNEW = pari_fopen_or_fail(LPNEW_str, WRITE);
    2945                 :            : 
    2946                 :            :     /* combine whatever there is to be combined */
    2947                 :       3233 :     tfc = 0;
    2948         [ +  + ]:       3233 :     if (tp > 0)
    2949                 :            :     {
    2950                 :            :       /* build full relations out of large prime relations */
    2951                 :        272 :       pCOMB = pari_fopen_or_fail(COMB_str, READ);
    2952                 :        272 :       tfc = mpqs_combine_large_primes(handle, pCOMB->file, pFNEW, &fact);
    2953                 :        272 :       pari_fclose(pCOMB);
    2954                 :            :       /* now FREL, LPREL are closed and FNEW, LPNEW are still open */
    2955         [ -  + ]:        272 :       if (fact)
    2956                 :            :       { /* factor found during combining */
    2957         [ #  # ]:          0 :         if (DEBUGLEVEL >= 4)
    2958                 :            :         {
    2959                 :          0 :           err_printf("\nMPQS: split N whilst combining, time = %ld ms\n",
    2960                 :            :                      timer_delay(&T));
    2961                 :          0 :           err_printf("MPQS: found factor = %Ps\n", fact);
    2962                 :            :         }
    2963                 :          0 :         pari_fclose(pLPNEW);
    2964                 :          0 :         pari_fclose(pFNEW);
    2965         [ #  # ]:          0 :         unlink_all();
    2966                 :          0 :         return gerepileupto(av, fact);
    2967                 :            :       }
    2968                 :        272 :       total_partial_relations += tp;
    2969                 :            :     }
    2970                 :            : 
    2971                 :            :     /* sort FNEW and merge it into FREL */
    2972                 :       3233 :     pari_fclose(pFNEW);
    2973                 :       3233 :     (void)mpqs_sort_lp_file(FNEW_str);
    2974                 :            :     /* definitive count (combinables combined, and duplicates removed) */
    2975                 :       3233 :     total_full_relations = mpqs_mergesort_lp_file(FREL_str, FNEW_str, TMP_str, NULL);
    2976                 :            :     /* FNEW stays closed until we need to reopen it for another iteration */
    2977                 :            : 
    2978                 :            :     /* Due to the removal of duplicates, percentage may actually decrease at
    2979                 :            :      * this point.  Looks funny in the diagnostics but is nothing to worry
    2980                 :            :      * about: we _are_ making progress. */
    2981                 :       3233 :     percentage =
    2982                 :       3233 :       (long)((1000.0 * total_full_relations) / handle->target_no_rels);
    2983                 :       3233 :     net_yield =
    2984         [ +  - ]:       3233 :       (total_full_relations * 100.) / (total_no_cand ? total_no_cand : 1);
    2985                 :       3233 :     vain_iterations =
    2986                 :       3233 :       (long)((1000.0 * (iterations - good_iterations)) / iterations);
    2987                 :            : 
    2988                 :            :     /* Now estimate the current full relations yield rate:  we directly see
    2989                 :            :      * each time through the main loop how many full relations we're getting
    2990                 :            :      * as such from the sieve  (tff since the previous checkpoint),  but
    2991                 :            :      * only at checkpoints do we see how many we're typically combining
    2992                 :            :      * (tfc).  So we're really producing (tfc+tff)/tff as many full rels,
    2993                 :            :      * and when we get close to 100%, we should bias the next interval by
    2994                 :            :      * the inverse ratio.
    2995                 :            :      * Avoid drawing conclusions from too-small samples during very short
    2996                 :            :      * follow-on intervals  (in this case we'll just re-use an earlier
    2997                 :            :      * estimated ratio). */
    2998 [ +  + ][ +  - ]:       3233 :     if ((tfc >= 16) && (tff >= 20))
    2999                 :         65 :       tfc_ratio = (tfc + tff + 0.) / tff; /* floating-point division */
    3000                 :       3233 :     tff = 0;                    /* reset this count (tfc is always fresh) */
    3001                 :            : 
    3002         [ +  + ]:       3233 :     if (percentage >= 1000)     /* when Gauss had failed */
    3003                 :        385 :       sort_interval = percentage + 2;
    3004         [ +  + ]:       2848 :     else if (percentage >= 820)
    3005                 :            :     {
    3006         [ +  + ]:        979 :       if (tfc_ratio > 1.)
    3007                 :            :       {
    3008         [ +  + ]:         67 :         if (percentage + (followup_sort_interval >> 1) * tfc_ratio > 994)
    3009                 :            :         {
    3010                 :            :           /* aim for a _slight_ overshoot */
    3011                 :         27 :           sort_interval = (ulong)(percentage + 2 +
    3012                 :         27 :             (1000 - percentage) / tfc_ratio);
    3013                 :            :         }
    3014         [ -  + ]:         40 :         else if (percentage >= 980)
    3015                 :          0 :           sort_interval = percentage + 8;
    3016                 :            :         else
    3017                 :         40 :           sort_interval = percentage + (followup_sort_interval >> 1);
    3018                 :            :       }
    3019                 :            :       else
    3020                 :            :       {
    3021         [ +  + ]:        912 :         if (percentage >= 980)
    3022                 :        107 :           sort_interval = percentage + 10;
    3023                 :            :         else
    3024                 :        805 :           sort_interval = percentage + (followup_sort_interval >> 1);
    3025 [ +  + ][ +  - ]:        912 :         if (sort_interval >= 1000 && percentage < 1000)
    3026                 :        186 :           sort_interval = 1000;
    3027                 :            :       }
    3028                 :            :     }
    3029                 :            :     else
    3030                 :       1869 :       sort_interval = percentage + followup_sort_interval;
    3031                 :            : 
    3032         [ -  + ]:       3233 :     if (DEBUGLEVEL >= 4)
    3033                 :            :     {
    3034         [ #  # ]:          0 :       err_printf("MPQS: done sorting%s, time = %ld ms\n",
    3035                 :            :                  tp > 0 ? " and combining" : "", timer_delay(&T));
    3036                 :          0 :       err_printf("MPQS: found %3.1f%% of the required relations\n",
    3037                 :            :                  percentage/10.);
    3038         [ #  # ]:          0 :       if (DEBUGLEVEL >= 5)
    3039                 :            :       { /* total_full_relations are always plural */
    3040                 :            :         /* GN20050708: present code doesn't succeed in discarding all
    3041                 :            :          * dups, so don't lie about it... */
    3042                 :          0 :         err_printf("MPQS: found %ld full relations\n",
    3043                 :            :                    total_full_relations);
    3044         [ #  # ]:          0 :         if (lp_scale > 1)
    3045                 :          0 :         err_printf("MPQS:   (%ld of these from partial relations)\n",
    3046                 :            :                    total_partial_relations);
    3047                 :          0 :         err_printf("MPQS: Net yield: %4.3g full relations per 100 candidates\n",
    3048                 :            :                    net_yield);
    3049                 :          0 :         err_printf("MPQS:            %4.3g full relations per 100 polynomials\n",
    3050                 :          0 :                    (total_full_relations * 100.) / iterations);
    3051                 :          0 :         err_printf("MPQS: %4.1f%% of the polynomials yielded no candidates\n",
    3052                 :            :                    vain_iterations/10.);
    3053                 :          0 :         err_printf("MPQS: next sort point at %3.1f%%\n", sort_interval/10.);
    3054                 :            :       }
    3055                 :            :     }
    3056                 :            : 
    3057         [ +  + ]:       3233 :     if (percentage < 1000)
    3058                 :            :     {
    3059                 :       2848 :       pFNEW = pari_fopen_or_fail(FNEW_str, WRITE);
    3060                 :            :       /* LPNEW and FNEW are again open for writing */
    3061                 :       2848 :       continue; /* main loop */
    3062                 :            :     }
    3063                 :            : 
    3064                 :            :     /* percentage >= 1000, which implies total_full_relations > size_of_FB:
    3065                 :            :        try finishing it off */
    3066                 :            : 
    3067                 :            :     /* solve the system over F_2 */
    3068                 :            :     /* present code does NOT in fact guarantee absence of dup FRELs,
    3069                 :            :      * therefore removing the adjective "distinct" for the time being */
    3070         [ -  + ]:        385 :     if (DEBUGLEVEL >= 4)
    3071                 :          0 :       err_printf("\nMPQS: starting Gauss over F_2 on %ld relations\n",
    3072                 :            :                  total_full_relations);
    3073                 :        385 :     pFREL = pari_fopen_or_fail(FREL_str, READ);
    3074                 :        385 :     fact = mpqs_solve_linear_system(handle, pFREL, total_full_relations);
    3075                 :        385 :     pari_fclose(pFREL);
    3076                 :            : 
    3077         [ +  - ]:        385 :     if (fact)
    3078                 :            :     { /* solution found */
    3079         [ -  + ]:        385 :       if (DEBUGLEVEL >= 4)
    3080                 :            :       {
    3081                 :          0 :         err_printf("\nMPQS: time in Gauss and gcds = %ld ms\n", timer_delay(&T));
    3082         [ #  # ]:          0 :         if (typ(fact) == t_INT) err_printf("MPQS: found factor = %Ps\n", fact);
    3083                 :            :         else
    3084                 :            :         {
    3085                 :          0 :           long j, nf = (lg(fact)-1)/3;
    3086         [ #  # ]:          0 :           if (nf == 2)
    3087                 :            :             /* GN20050707: Changed the arrangement of the two factors,
    3088                 :            :              * to match the debug diagnostics in mpqs_solve_linear_system()
    3089                 :            :              * above */
    3090                 :          0 :             err_printf("MPQS: found factors = %Ps\n\tand %Ps\n",
    3091                 :          0 :                         fact[1], fact[4]);
    3092                 :            :           else
    3093                 :            :           {
    3094                 :            :             /* GN20050707: Changed loop to scan upwards instead of downwards,
    3095                 :            :              * to match the debug diagnostics in mpqs_solve_linear_system()
    3096                 :            :              * above */
    3097                 :          0 :             err_printf("MPQS: found %ld factors =\n", nf);
    3098         [ #  # ]:          0 :             for (j=1; j<=nf; j++)
    3099         [ #  # ]:          0 :               err_printf("\t%Ps%s\n", fact[3*j-2], (j<nf ? "," : ""));
    3100                 :            :           }
    3101                 :            :         }
    3102                 :            :       }
    3103                 :        385 :       pari_fclose(pLPNEW);
    3104         [ +  - ]:        385 :       unlink_all();
    3105                 :            :       /* fact not safe for a gerepilecopy(): segfaults on one of the NULL
    3106                 :            :        * markers. However, it is a nice connected object, and it resides
    3107                 :            :        * already the top of the stack, so... --GN */
    3108                 :        385 :       return gerepileupto(av, fact);
    3109                 :            :     }
    3110                 :            :     else
    3111                 :            :     {
    3112         [ #  # ]:          0 :       if (DEBUGLEVEL >= 4)
    3113                 :            :       {
    3114                 :          0 :         err_printf("\nMPQS: time in Gauss and gcds = %ld ms\n",timer_delay(&T));
    3115                 :          0 :         err_printf("MPQS: no factors found.\n");
    3116         [ #  # ]:          0 :         if (percentage <= MPQS_ADMIT_DEFEAT)
    3117                 :          0 :           err_printf("\nMPQS: restarting sieving ...\n");
    3118                 :            :         else
    3119                 :          0 :           err_printf("\nMPQS: giving up.\n");
    3120                 :            :       }
    3121         [ #  # ]:          0 :       if (percentage > MPQS_ADMIT_DEFEAT)
    3122                 :            :       {
    3123                 :          0 :         pari_fclose(pLPNEW);
    3124         [ #  # ]:          0 :         unlink_all(); avma = av; return NULL;
    3125                 :            :       }
    3126                 :          0 :       pFNEW = pari_fopen_or_fail(FNEW_str, WRITE);
    3127                 :            :     }
    3128                 :      25375 :   } /* main loop */
    3129                 :            : }
    3130                 :            : GEN
    3131                 :        385 : mpqs(GEN N)
    3132                 :            : {
    3133                 :        385 :   mpqs_handle_t *handle = mpqs_handle_ctor(N);
    3134                 :        385 :   GEN fact = mpqs_i(handle);
    3135                 :        385 :   mpqs_handle_dtor(handle); return fact;
    3136                 :            : }

Generated by: LCOV version 1.9