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-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 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.12.1 lcov report (development 24038-ebe36f6c4) Lines: 664 852 77.9 %
Date: 2019-07-23 05:53:17 Functions: 31 32 96.9 %
Legend: Lines: hit not hit

          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             : #include "pari.h"
      74             : #include "paripriv.h"
      75             : 
      76             : /** DEBUG **/
      77             : /* #define MPQS_DEBUG_VERBOSE 1 */
      78             : #include "mpqs.h"
      79             : 
      80             : #define REL_OFFSET 20
      81             : #define REL_MASK ((1UL<<REL_OFFSET)-1)
      82             : #define MAX_PE_PAIR 60
      83             : #define DEFAULT_VEC_LEN 17
      84             : 
      85      109255 : static GEN rel_q(GEN c) { return gel(c,1); }
      86       20640 : static GEN rel_Y(GEN c) { return gel(c,2); }
      87       20640 : static GEN rel_p(GEN c) { return gel(c,3); }
      88             : 
      89             : static void
      90      100471 : frel_add(hashtable *frel, GEN R)
      91             : {
      92      100471 :   ulong h = hash_GEN(R);
      93      100471 :   if (!hash_search2(frel, (void*)R, h))
      94       99799 :     hash_insert2(frel, (void*)R, (void*)1, h);
      95      100471 : }
      96             : static void
      97       24999 : vec_frel_add(hashtable *frel, GEN V)
      98             : {
      99       24999 :   long i, l = lg(V);
     100      115150 :   for (i = 1; i<l ; i++)
     101       90151 :     frel_add(frel, gel(V,i));
     102       24999 : }
     103             : 
     104             : static GEN
     105      224763 : vec_extend(GEN frel, GEN rel, long nfrel)
     106             : {
     107      224763 :   long lfrel = lg(frel)-1;
     108      224763 :   if (nfrel > lfrel)
     109             :   {
     110         103 :     lfrel *= 2;
     111         103 :     frel = vec_lengthen(frel, lfrel);
     112         103 :     if (DEBUGLEVEL >= 4) err_printf("MPQS: extending store to %ld\n",lfrel);
     113             :   }
     114      224763 :   gel(frel, nfrel) = rel;
     115      224763 :   return frel;
     116             : }
     117             : 
     118             : /*********************************************************************/
     119             : /**                                                                 **/
     120             : /**                         INITIAL SIZING                          **/
     121             : /**                                                                 **/
     122             : /*********************************************************************/
     123             : /* number of decimal digits of argument - for parameter choosing and for
     124             :  * diagnostics */
     125             : static long
     126         558 : decimal_len(GEN N)
     127         558 : { pari_sp av = avma; return gc_long(av, 1+logint(N,stoi(10))); }
     128             : 
     129             : /* To be called after choosing k and putting kN into the handle:
     130             :  * Pick up the requested parameter set for the given size of kN in decimal
     131             :  * digits and fill in various fields in the handle.  Return 0 when kN is
     132             :  * too large, 1 when we're ok. */
     133             : static int
     134         279 : mpqs_set_parameters(mpqs_handle_t *h)
     135             : {
     136             :   long i;
     137             :   const mpqs_parameterset_t *P;
     138             : 
     139         279 :   h->digit_size_kN = decimal_len(h->kN);
     140         279 :   if (h->digit_size_kN <= 9)
     141           0 :     i = 0;
     142         279 :   else if (h->digit_size_kN > MPQS_MAX_DIGIT_SIZE_KN)
     143           0 :     return 0;
     144             :   else
     145         279 :     i = h->digit_size_kN - 9;
     146             : 
     147             :   /* (cf PARI bug#235) the following has always been, and will remain,
     148             :    * a moving target... increased thresholds from 64, 80 to 79, 86
     149             :    * respectively --GN20050601.  Note that the new values correspond to
     150             :    * kN having >= 86 or >= 95 decimal digits, respectively.  Note also
     151             :    * that the current sizing parameters for 90 or more digits are based
     152             :    * on 100% theory and 0% practice. */
     153         279 :   if (i >= 79)
     154           0 :     pari_warn(warner, "MPQS: factoring this number will take %s hours:\nN = %Ps",
     155             :         i >= 86 ? "many": "several", h->N);
     156             : 
     157         279 :   if (DEBUGLEVEL >= 5)
     158             :   {
     159           0 :     err_printf("MPQS: kN = %Ps\n", h->kN);
     160           0 :     err_printf("MPQS: kN has %ld decimal digits\n", h->digit_size_kN);
     161             :   }
     162             : 
     163         279 :   P = &(mpqs_parameters[i]);
     164         279 :   h->tolerance        = P->tolerance;
     165         279 :   h->lp_scale         = P->lp_scale;
     166             :   /* make room for prime factors of k if any: */
     167         279 :   h->size_of_FB       = P->size_of_FB + h->_k->omega_k;
     168             :   /* for the purpose of Gauss elimination etc., prime factors of k behave
     169             :    * like real FB primes, so take them into account when setting the goal: */
     170         558 :   h->target_no_rels   = (h->size_of_FB >= 200 ?
     171         468 :                          h->size_of_FB + 70 :
     172         189 :                          (mpqs_int32_t)(h->size_of_FB * 1.35));
     173         279 :   h->M                = P->M;
     174         279 :   h->omega_A          = P->omega_A;
     175         279 :   h->no_B             = 1UL << (P->omega_A - 1);
     176         279 :   h->pmin_index1      = P->pmin_index1;
     177             :   /* certain subscripts into h->FB should also be offset by omega_k: */
     178         279 :   h->index0_FB        = 3 + h->_k->omega_k;
     179             :   /* following are converted from % to parts per thousand: */
     180         279 :   h->first_sort_point = 10 * P->first_sort_point;
     181         279 :   h->sort_pt_interval = 10 * P->sort_pt_interval;
     182             : 
     183         279 :   if (DEBUGLEVEL >= 5)
     184             :   {
     185           0 :     double mb = (h->size_of_FB + 1)/(8.*1048576.) * h->target_no_rels;
     186           0 :     err_printf("\t(estimated memory needed: %4.1fMBy)\n", mb);
     187             :   }
     188         279 :   return 1;
     189             : }
     190             : 
     191             : /*********************************************************************/
     192             : /**                                                                 **/
     193             : /**                       OBJECT HOUSEKEEPING                       **/
     194             : /**                                                                 **/
     195             : /*********************************************************************/
     196             : 
     197             : /* factor base constructor. Really a home-grown memalign(3c) underneath.
     198             :  * We don't want FB entries to straddle L1 cache line boundaries, and
     199             :  * malloc(3c) only guarantees alignment adequate for all primitive data
     200             :  * types of the platform ABI - typically to 8 or 16 byte boundaries.
     201             :  * Also allocate the inv_A_H array.
     202             :  * The FB array pointer is returned for convenience */
     203             : static mpqs_FB_entry_t *
     204         279 : mpqs_FB_ctor(mpqs_handle_t *h)
     205             : {
     206             :   /* leave room for slots 0, 1, and sentinel slot at the end of the array */
     207         279 :   long size_FB_chunk = (h->size_of_FB + 3) * sizeof(mpqs_FB_entry_t);
     208             :   /* like FB, except this one does not have a sentinel slot at the end */
     209         279 :   long size_IAH_chunk = (h->size_of_FB + 2) * sizeof(mpqs_inv_A_H_t);
     210         279 :   char *fbp = (char*)stack_malloc(size_FB_chunk + 64);
     211         279 :   char *iahp = (char*)stack_malloc(size_IAH_chunk + 64);
     212             :   long fbl, iahl;
     213             : 
     214         279 :   h->FB_chunk = (void *)fbp;
     215         279 :   h->invAH_chunk = (void *)iahp;
     216             :   /* round up to next higher 64-bytes-aligned address */
     217         279 :   fbl = (((long)fbp) + 64) & ~0x3FL;
     218             :   /* and put the actual array there */
     219         279 :   h->FB = (mpqs_FB_entry_t *)fbl;
     220             : 
     221         279 :   iahl = (((long)iahp) + 64) & ~0x3FL;
     222         279 :   h->inv_A_H = (mpqs_inv_A_H_t *)iahl;
     223         279 :   return (mpqs_FB_entry_t *)fbl;
     224             : }
     225             : 
     226             : /* sieve array constructor;  also allocates the candidates array
     227             :  * and temporary storage for relations under construction */
     228             : static void
     229         279 : mpqs_sieve_array_ctor(mpqs_handle_t *h)
     230             : {
     231         279 :   long size = (h->M << 1) + 1;
     232         279 :   mpqs_int32_t size_of_FB = h->size_of_FB;
     233             : 
     234         279 :   h->sieve_array = (unsigned char *) stack_malloc(size * sizeof(unsigned char));
     235         279 :   h->sieve_array_end = h->sieve_array + size - 2;
     236         279 :   h->sieve_array_end[1] = 255; /* sentinel */
     237         279 :   h->candidates = (long *)stack_malloc(MPQS_CANDIDATE_ARRAY_SIZE * sizeof(long));
     238             :   /* whereas mpqs_self_init() uses size_of_FB+1, we just use the size as
     239             :    * it is, not counting FB[1], to start off the following estimate */
     240         279 :   if (size_of_FB > MAX_PE_PAIR) size_of_FB = MAX_PE_PAIR;
     241             :   /* and for tracking which primes occur in the current relation: */
     242         279 :   h->relaprimes = (long *) stack_malloc((size_of_FB << 1) * sizeof(long));
     243         279 : }
     244             : 
     245             : /* mpqs() calls the following (after recording avma) to allocate GENs for
     246             :  * the current polynomial and self-initialization scratchpad data on the
     247             :  * PARI stack.  This space is released by mpqs() itself at the end. */
     248             : static void
     249         279 : mpqs_poly_ctor(mpqs_handle_t *h)
     250             : {
     251             :   mpqs_int32_t i;
     252         279 :   long size_per = h->omega_A * sizeof(mpqs_per_A_prime_t);
     253             : 
     254         279 :   h->per_A_pr = (mpqs_per_A_prime_t *) stack_calloc(size_per);
     255             :   /* Sizing:  A is the product of omega_A primes, each well below word
     256             :    * size.
     257             :    * |B| is bounded by (omega_A + 4) * A, so can have at most one word
     258             :    * more, and that's generous.
     259             :    * |C| is less than A*M^2, so can take at most two words more than A.
     260             :    * The array H holds residues modulo A, so the same size as used for A
     261             :    * is sufficient. */
     262         279 :   h->A = cgeti(h->omega_A + 2);
     263         279 :   h->B = cgeti(h->omega_A + 3);
     264             : #ifdef MPQS_DEBUG
     265             :   h->C = cgeti(h->omega_A + 4);
     266             : #endif
     267        1308 :   for (i = 0; i < h->omega_A; i++)
     268        1029 :     h->per_A_pr[i]._H = cgeti(h->omega_A + 2);
     269             :   /* the handle starts out all zero, so in particular bin_index and index_i
     270             :    * are initially 0.
     271             :    * TODO: index_j currently initialized in mqps() but this is going to
     272             :    * change. */
     273         279 : }
     274             : 
     275             : /* TODO: relationsdb handle */
     276             : 
     277             : /*********************************************************************/
     278             : /**                                                                 **/
     279             : /**                        FACTOR BASE SETUP                        **/
     280             : /**                                                                 **/
     281             : /*********************************************************************/
     282             : /* fill in the best-guess multiplier k for N. We force kN = 1 mod 4.
     283             :  * Caller should proceed to fill in kN */
     284             : static ulong
     285         279 : mpqs_find_k(mpqs_handle_t *h)
     286             : {
     287         279 :   const pari_sp av = avma;
     288         279 :   const long N_mod_8 = mod8(h->N), N_mod_4 = N_mod_8 & 3;
     289             :   forprime_t S;
     290             :   struct {
     291             :     const mpqs_multiplier_t *_k;
     292             :     long np; /* number of primes in factorbase so far for this k */
     293             :     double value; /* the larger, the better */
     294             :   } cache[MPQS_POSSIBLE_MULTIPLIERS];
     295             :   ulong p, i, nbk;
     296             : 
     297        2706 :   for (i = nbk = 0; i < numberof(cand_multipliers); i++)
     298             :   {
     299        2706 :     const mpqs_multiplier_t *cand_k = &cand_multipliers[i];
     300        2706 :     long k = cand_k->k;
     301             :     double v;
     302        2706 :     if ((k & 3) != N_mod_4) continue; /* want kN = 1 (mod 4) */
     303        1395 :     v = -0.35 * log2((double)k);
     304        1395 :     if ((k & 7) == N_mod_8) v += M_LN2; /* kN = 1 (mod 8) */
     305        1395 :     cache[nbk].np = 0;
     306        1395 :     cache[nbk]._k = cand_k;
     307        1395 :     cache[nbk].value = v;
     308        1395 :     if (++nbk == MPQS_POSSIBLE_MULTIPLIERS) break; /* enough */
     309             :   }
     310             :   /* next test is an impossible situation: kills spurious gcc-5.1 warnings
     311             :    * "array subscript is above array bounds" */
     312         279 :   if (nbk > MPQS_POSSIBLE_MULTIPLIERS) nbk = MPQS_POSSIBLE_MULTIPLIERS;
     313         279 :   u_forprime_init(&S, 2, ULONG_MAX);
     314         279 :   while ( (p = u_forprime_next(&S)) )
     315             :   {
     316        5147 :     ulong Np = umodiu(h->N, p);
     317        5147 :     long kroNp, seen = 0;
     318        5147 :     if (!Np) return p;
     319        5147 :     kroNp = krouu(Np, p);
     320       30882 :     for (i = 0; i < nbk; i++)
     321             :     {
     322       25735 :       if (cache[i].np > MPQS_MULTIPLIER_SEARCH_DEPTH) continue;
     323       17415 :       seen++;
     324       17415 :       if (krouu(cache[i]._k->k % p, p) == kroNp) /* kronecker(k*N, p)=1 */
     325             :       {
     326        8370 :         cache[i].value += log2((double) p)/p;
     327        8370 :         cache[i].np++;
     328             :       }
     329             :     }
     330        5147 :     if (!seen) break; /* we're gone through SEARCH_DEPTH primes for all k */
     331             :   }
     332         279 :   if (!p) pari_err_OVERFLOW("mpqs_find_k [ran out of primes]");
     333             :   {
     334         279 :     long best_i = 0;
     335         279 :     double v = cache[0].value;
     336        1395 :     for (i = 1; i < nbk; i++)
     337        1116 :       if (cache[i].value > v) { best_i = i; v = cache[i].value; }
     338         279 :     h->_k = cache[best_i]._k; return gc_ulong(av,0);
     339             :   }
     340             : }
     341             : 
     342             : /******************************/
     343             : /* Create a factor base of 'size' primes p_i such that legendre(k*N, p_i) != -1
     344             :  * We could have shifted subscripts down from their historical arrangement,
     345             :  * but this seems too risky for the tiny potential gain in memory economy.
     346             :  * The real constraint is that the subscripts of anything which later shows
     347             :  * up at the Gauss stage must be nonnegative, because the exponent vectors
     348             :  * there use the same subscripts to refer to the same FB entries.  Thus in
     349             :  * particular, the entry representing -1 could be put into FB[0], but could
     350             :  * not be moved to FB[-1] (although mpqs_FB_ctor() could be easily adapted
     351             :  * to support negative subscripts).-- The historically grown layout is:
     352             :  * FB[0] is unused.
     353             :  * FB[1] is not explicitly used but stands for -1.
     354             :  * FB[2] contains 2 (always).
     355             :  * Before we are called, the size_of_FB field in the handle will already have
     356             :  * been adjusted by _k->omega_k, so there's room for the primes dividing k,
     357             :  * which when present will occupy FB[3] and following.
     358             :  * The "real" odd FB primes begin at FB[h->index0_FB].
     359             :  * FB[size_of_FB+1] is the last prime p_i.
     360             :  * FB[size_of_FB+2] is a sentinel to simplify some of our loops.
     361             :  * Thus we allocate size_of_FB+3 slots for FB.
     362             :  *
     363             :  * If a prime factor of N is found during the construction, it is returned
     364             :  * in f, otherwise f = 0. */
     365             : 
     366             : /* returns the FB array pointer for convenience */
     367             : static mpqs_FB_entry_t *
     368         279 : mpqs_create_FB(mpqs_handle_t *h, ulong *f)
     369             : {
     370         279 :   mpqs_FB_entry_t *FB = mpqs_FB_ctor(h);
     371         279 :   const pari_sp av = avma;
     372         279 :   mpqs_int32_t size = h->size_of_FB;
     373             :   long i;
     374         279 :   mpqs_uint32_t k = h->_k->k;
     375             :   forprime_t S;
     376             : 
     377         279 :   FB[2].fbe_p = 2;
     378             :   /* the fbe_logval and the fbe_sqrt_kN for 2 are never used */
     379         279 :   FB[2].fbe_flags = MPQS_FBE_CLEAR;
     380         279 :   (void)u_forprime_init(&S, 3, ULONG_MAX);
     381             : 
     382             :   /* the first loop executes h->_k->omega_k = 0, 1, or 2 times */
     383         405 :   for (i = 3; i < h->index0_FB; i++)
     384             :   {
     385         126 :     mpqs_uint32_t kp = (ulong)h->_k->kp[i-3];
     386         126 :     if (MPQS_DEBUGLEVEL >= 7) err_printf(",<%lu>", (ulong)kp);
     387         126 :     FB[i].fbe_p = kp;
     388             :     /* we *could* flag divisors of k here, but so far I see no need,
     389             :      * and no flags bit has been assigned for the purpose */
     390         126 :     FB[i].fbe_flags = MPQS_FBE_CLEAR;
     391         126 :     FB[i].fbe_flogp = (float) log2((double) kp);
     392         126 :     FB[i].fbe_sqrt_kN = 0;
     393             :   }
     394             :   /* now i == h->index0_FB */
     395      170256 :   while (i < size + 2)
     396             :   {
     397      169698 :     ulong p = u_forprime_next(&S);
     398      169698 :     if (p > k || k % p)
     399             :     {
     400      169572 :       ulong kN_mod_p = umodiu(h->kN, p);
     401      169572 :       long kr = krouu(kN_mod_p, p);
     402      169572 :       if (kr != -1)
     403             :       {
     404       84746 :         if (kr == 0) { *f = p; return FB; }
     405       84746 :         FB[i].fbe_p = (mpqs_uint32_t) p;
     406       84746 :         FB[i].fbe_flags = MPQS_FBE_CLEAR;
     407             :         /* dyadic logarithm of p; single precision suffices */
     408       84746 :         FB[i].fbe_flogp = (float) log2((double)p);
     409             :         /* cannot yet fill in fbe_logval because the scaling multiplier
     410             :          * depends on the largest prime in FB, as yet unknown */
     411             : 
     412             :         /* x such that x^2 = kN (mod p_i) */
     413       84746 :         FB[i++].fbe_sqrt_kN = (mpqs_uint32_t)Fl_sqrt(kN_mod_p, p);
     414             :       }
     415             :     }
     416             :   }
     417         279 :   set_avma(av);
     418         279 :   if (MPQS_DEBUGLEVEL >= 7)
     419             :   {
     420           0 :     err_printf("MPQS: FB [-1,2");
     421           0 :     for (i = 3; i < h->index0_FB; i++) err_printf(",<%lu>", FB[i].fbe_p);
     422           0 :     for (; i < size + 2; i++) err_printf(",%lu", FB[i].fbe_p);
     423           0 :     err_printf("]\n");
     424             :   }
     425             : 
     426         279 :   FB[i].fbe_p = 0;              /* sentinel */
     427         279 :   h->largest_FB_p = FB[i-1].fbe_p; /* at subscript size_of_FB + 1 */
     428             : 
     429             :   /* locate the smallest prime that will be used for sieving */
     430         710 :   for (i = h->index0_FB; FB[i].fbe_p != 0; i++)
     431         710 :     if (FB[i].fbe_p >= h->pmin_index1) break;
     432         279 :   h->index1_FB = i;
     433             :   /* with our parameters this will never fall of the end of the FB */
     434         279 :   *f = 0; return FB;
     435             : }
     436             : 
     437             : /*********************************************************************/
     438             : /**                                                                 **/
     439             : /**                      MISC HELPER FUNCTIONS                      **/
     440             : /**                                                                 **/
     441             : /*********************************************************************/
     442             : 
     443             : /* Effect of the following:  multiplying the base-2 logarithm of some
     444             :  * quantity by log_multiplier will rescale something of size
     445             :  *    log2 ( sqrt(kN) * M / (largest_FB_prime)^tolerance )
     446             :  * to 232.  Note that sqrt(kN) * M is just A*M^2, the value our polynomials
     447             :  * take at the outer edges of the sieve interval.  The scale here leaves
     448             :  * a little wiggle room for accumulated rounding errors from the approximate
     449             :  * byte-sized scaled logarithms for the factor base primes which we add up
     450             :  * in the sieving phase.-- The threshold is then chosen so that a point in
     451             :  * the sieve has to reach a result which, under the same scaling, represents
     452             :  *    log2 ( sqrt(kN) * M / (largest_FB_prime)^tolerance )
     453             :  * in order to be accepted as a candidate. */
     454             : /* The old formula was...
     455             :  *   log_multiplier =
     456             :  *      127.0 / (0.5 * log2 (handle->dkN)
     457             :  *               + log2((double)M)
     458             :  *               - tolerance * log2((double)handle->largest_FB_p)
     459             :  *               );
     460             :  * and we used to use this with a constant threshold of 128. */
     461             : 
     462             : /* NOTE: We used to divide log_multiplier by an extra factor 2, and in
     463             :  * compensation we were multiplying by 2 when the fbe_logp fields were being
     464             :  * filled in, making all those bytes even.  Tradeoff: the extra bit of
     465             :  * precision is helpful, but interferes with a possible sieving optimization
     466             :  * (artifically shift right the logp's of primes in A, and just run over both
     467             :  * arithmetical progressions  (which coincide in this case)  instead of
     468             :  * skipping the second one, to avoid the conditional branch in the
     469             :  * mpqs_sieve() loops).  We could still do this, but might lose a little bit
     470             :  * accuracy for those primes.  Probably no big deal. */
     471             : static void
     472         279 : mpqs_set_sieve_threshold(mpqs_handle_t *h)
     473             : {
     474         279 :   mpqs_FB_entry_t *FB = h->FB;
     475             :   long i;
     476             :   double log_maxval;
     477             :   double log_multiplier;
     478             : 
     479         279 :   h->l2sqrtkN = 0.5 * log2(h->dkN);
     480         279 :   h->l2M = log2((double)h->M);
     481         279 :   log_maxval = h->l2sqrtkN + h->l2M - MPQS_A_FUDGE;
     482         279 :   log_multiplier = 232.0 / log_maxval;
     483         279 :   h->sieve_threshold =
     484         279 :     (unsigned char) (log_multiplier *
     485             :                      (log_maxval
     486         279 :                       - h->tolerance * log2((double)h->largest_FB_p)
     487             :                       )
     488         279 :                      ) + 1;
     489             :   /* That "+ 1" really helps - we may want to tune towards somewhat smaller
     490             :    * tolerances  (or introduce self-tuning one day)... */
     491             : 
     492             :   /* If this turns out to be <128, scream loudly.
     493             :    * That means that the FB or the tolerance or both are way too
     494             :    * large for the size of kN.  (Normally, the threshold should end
     495             :    * up in the 150...170 range.) */
     496         279 :   if (h->sieve_threshold < 128) {
     497           0 :     h->sieve_threshold = 128;
     498           0 :     pari_warn(warner,
     499             :         "MPQS: sizing out of tune, FB size or tolerance\n\ttoo large");
     500             :   }
     501             : 
     502             :   /* Now fill in the byte-sized approximate scaled logarithms of p_i */
     503         279 :   if (DEBUGLEVEL >= 5)
     504             :   {
     505           0 :     err_printf("MPQS: computing logarithm approximations for p_i in FB\n");
     506             :   }
     507       85025 :   for (i = h->index0_FB; i < h->size_of_FB + 2; i++)
     508             :   {
     509      169492 :     FB[i].fbe_logval =
     510       84746 :       (unsigned char) (log_multiplier * FB[i].fbe_flogp);
     511             :   }
     512         279 : }
     513             : 
     514             : /* Given the partially populated handle, find the optimum place in the FB
     515             :  * to pick prime factors for A from.  The lowest admissible subscript is
     516             :  * index0_FB, but unless kN is very small, we stay away a bit from that.
     517             :  * The highest admissible is size_of_FB + 1, where the largest FB prime
     518             :  * resides.  The ideal corner is about (sqrt(kN)/M) ^ (1/omega_A),
     519             :  * so that A will end up of size comparable to sqrt(kN)/M;  experimentally
     520             :  * it seems desirable to stay slightly below this.  Moreover, the selection
     521             :  * of the individual primes happens to err on the large side, for which we
     522             :  * compensate a bit, using the (small positive) quantity MPQS_A_FUDGE.
     523             :  * We rely on a few auxiliary fields in the handle to be already set by
     524             :  * mqps_set_sieve_threshold() before we are called.
     525             :  * Return 1 on success, and 0 otherwise. */
     526             : static int
     527         279 : mpqs_locate_A_range(mpqs_handle_t *h)
     528             : {
     529             :   /* i will be counted up to the desirable index2_FB + 1, and omega_A is never
     530             :    * less than 3, and we want
     531             :    *   index2_FB - (omega_A - 1) + 1 >= index0_FB + omega_A - 3,
     532             :    * so: */
     533         279 :   long i = h->index0_FB + 2*(h->omega_A) - 4;
     534             :   double l2_target_pA;
     535         279 :   mpqs_FB_entry_t *FB = h->FB;
     536             : 
     537         279 :   h->l2_target_A = (h->l2sqrtkN - h->l2M - MPQS_A_FUDGE);
     538         279 :   l2_target_pA = h->l2_target_A / h->omega_A;
     539             : 
     540             :   /* find the sweet spot, normally shouldn't take long */
     541         279 :   while ((FB[i].fbe_p != 0) && (FB[i].fbe_flogp <= l2_target_pA)) i++;
     542             : 
     543             : #ifdef MPQS_DEBUG_LOCATE_A_RANGE
     544             :   err_printf("MPQS DEBUG: omega_A=%ld, index0=%ld, i=%ld\n",
     545             :              (long) h->omega_A, (long) h->index0_FB, i);
     546             : #endif
     547             : 
     548             :   /* check whether this hasn't walked off the top end... */
     549             :   /* The following should actually NEVER happen. */
     550         279 :   if (i > h->size_of_FB - 3)
     551             :   { /* this isn't going to work at all. */
     552           0 :     pari_warn(warner,
     553             :         "MPQS: sizing out of tune, FB too small or\n\tway too few primes in A");
     554           0 :     return 0;
     555             :   }
     556         279 :   h->index2_FB = i - 1;
     557             : #ifdef MPQS_DEBUG_LOCATE_A_RANGE
     558             :   err_printf("MPQS DEBUG: index2_FB = %ld\n", i - 1);
     559             : #endif
     560             :   /* GN20050723
     561             :    * assert: index0_FB + (omega_A - 3) [the lowest FB subscript eligible to
     562             :    * be used in picking primes for A]  plus  (omega_A - 2)  does not exceed
     563             :    * index2_FB  [the subscript from which the choice of primes for A starts,
     564             :    * putting omega_A - 1 of them at or below index2_FB, and the last and
     565             :    * largest one above, cf. mpqs_si_choose_primes() below].
     566             :    * Moreover, index2_FB indicates the last prime below the ideal size, unless
     567             :    * (when kN is very small) the ideal size was too small to use. */
     568             : 
     569         279 :   return 1;
     570             : }
     571             : 
     572             : /*********************************************************************/
     573             : /**                                                                 **/
     574             : /**                       SELF-INITIALIZATION                       **/
     575             : /**                                                                 **/
     576             : /*********************************************************************/
     577             : 
     578             : #ifdef MPQS_DEBUG
     579             : /* Debug-only helper routine: check correctness of the root z mod p_i
     580             :  * by evaluting A * z^2 + B * z + C mod p_i  (which should be 0).
     581             :  * C is written as (B*B-kN)/(4*A) */
     582             : static void
     583             : check_root(mpqs_handle_t *h, long p, long start)
     584             : {
     585             :   long z = start - ((long)(h->M) % p);
     586             :   if (smodis(addii(h->C, mulsi(z, addii(h->B, mulsi(z, h->A)))), p))
     587             :   {
     588             :     err_printf("MPQS: p = %ld\n", p);
     589             :     err_printf("MPQS: A = %Ps\n", h->A);
     590             :     err_printf("MPQS: B = %Ps\n", h->B);
     591             :     err_printf("MPQS: C = %Ps\n", h->C);
     592             :     err_printf("MPQS: z = %ld\n", z);
     593             :     pari_err_BUG("MPQS: self_init: found wrong polynomial");
     594             :   }
     595             : }
     596             : #endif
     597             : 
     598             : /* There are four parts to self-initialization, which are exercised at
     599             :  * somewhat different times:
     600             :  * - choosing a new coefficient A  (selecting the prime factors to go into it,
     601             :  *   and doing the required bookkeeping in the FB entries, including clearing
     602             :  *   out the flags from the previous cohort), together with:
     603             :  * - doing the actual computations attached to a new A
     604             :  * - choosing a new B keeping the same A (very much simpler and quicker)
     605             :  * - and a small common bit that needs to happen in both cases.
     606             :  * As to the first item, the new scheme works as follows:
     607             :  * We pick omega_A - 1 prime factors for A below the index2_FB point which
     608             :  * marks their ideal size, and one prime above this point, choosing the
     609             :  * latter so as to get log2(A) as close as possible to l2_target_A.
     610             :  * The lower prime factors are chosen using bit patterns of constant weight,
     611             :  * gradually moving away from index2_FB towards smaller FB subscripts.
     612             :  * If this bumps into index0_FB  (might happen for very small input),  we
     613             :  * back up by increasing index2_FB by two, and from then on choosing only
     614             :  * bit patterns with either or both of their bottom bits set, so at least
     615             :  * one of the omega_A - 1 smaller prime factor will be beyond the original
     616             :  * index2_FB point.  In this way we avoid re-using A's which had already
     617             :  * been done.
     618             :  * (The choice of the upper "flyer" prime is of course constrained by the
     619             :  * size of the FB, which normally should never be anywhere close to becoming
     620             :  * a problem.  In unfortunate cases, e.g. for very small kN, we might have
     621             :  * to live with a rather non-optimal choice.
     622             :  * Then again, MPQS as such is surprisingly robust.  One day, I had got the
     623             :  * order of entries in mpqs_parameterset_t mixed up, and was running on a
     624             :  * smallish N with a huuuuge factor base and extremely tiny sieve interval,
     625             :  * and it still managed to factor it fairly quickly...)
     626             :  *
     627             :  * Mathematically, there isn't much more to this than the usual formula for
     628             :  * solving a quadratic  (over the field of p elements for each prime p in
     629             :  * the FB which doesn't divide A),  solving a linear equation for each of
     630             :  * the primes which do divide A, and precomputing differences between roots
     631             :  * mod p so we can adjust the roots quickly when we change B.
     632             :  * See Thomas Sosnowski's diploma thesis.
     633             :  */
     634             : 
     635             : /* Helper function:
     636             :  * Increment *x (!=0) to a larger value which has the same number of 1s in its
     637             :  * binary representation.  Wraparound can be detected by the caller as long as
     638             :  * we keep total_no_of_primes_for_A strictly less than BITS_IN_LONG.
     639             :  *
     640             :  * Changed switch to increment *x in all cases to the next larger number
     641             :  * which (a) has the same count of 1 bits and (b) does not arise from the
     642             :  * old value by moving a single 1 bit one position to the left  (which was
     643             :  * undesirable for the sieve). --GN based on discussion with TP */
     644             : INLINE void
     645        1854 : mpqs_increment(mpqs_uint32_t *x)
     646             : {
     647        1854 :   mpqs_uint32_t r1_mask, r01_mask, slider=1UL;
     648             : 
     649             :   /* 32-way computed jump handles 22 out of 32 cases */
     650        1854 :   switch (*x & 0x1F)
     651             :   {
     652             :   case 29:
     653          39 :     (*x)++; break; /* shifts a single bit, but we postprocess this case */
     654             :   case 26:
     655           0 :     (*x) += 2; break; /* again */
     656             :   case 1: case 3: case 6: case 9: case 11:
     657             :   case 17: case 19: case 22: case 25: case 27:
     658        1080 :     (*x) += 3; return;
     659             :   case 20:
     660          34 :     (*x) += 4; break; /* again */
     661             :   case 5: case 12: case 14: case 21:
     662          55 :     (*x) += 5; return;
     663             :   case 2: case 7: case 13: case 18: case 23:
     664         347 :     (*x) += 6; return;
     665             :   case 10:
     666           0 :     (*x) += 7; return;
     667             :   case 8:
     668           0 :     (*x) += 8; break; /* and again */
     669             :   case 4: case 15:
     670          85 :     (*x) += 12; return;
     671             :   default: /* 0, 16, 24, 28, 30, 31 */
     672             :     /* isolate rightmost 1 */
     673         214 :     r1_mask = ((*x ^ (*x - 1)) + 1) >> 1;
     674             :     /* isolate rightmost 1 which has a 0 to its left */
     675         214 :     r01_mask = ((*x ^ (*x + r1_mask)) + r1_mask) >> 2;
     676             :     /* simple cases.  Both of these shift a single bit one position to the
     677             :        left, and will need postprocessing */
     678         214 :     if (r1_mask == r01_mask) { *x += r1_mask; break; }
     679         214 :     if (r1_mask == 1) { *x += r01_mask; break; }
     680             :     /* general case -- idea: add r01_mask, kill off as many 1 bits as possible
     681             :      * to its right while at the same time filling in 1 bits from the LSB. */
     682         169 :     if (r1_mask == 2) { *x += (r01_mask>>1) + 1; return; }
     683         364 :     while (r01_mask > r1_mask && slider < r1_mask)
     684             :     {
     685         182 :       r01_mask >>= 1; slider <<= 1;
     686             :     }
     687          91 :     *x += r01_mask + slider - 1;
     688          91 :     return;
     689             :   }
     690             :   /* post-process all cases which couldn't be finalized above.  If we get
     691             :      here, slider still has its original value. */
     692         118 :   r1_mask = ((*x ^ (*x - 1)) + 1) >> 1;
     693         118 :   r01_mask = ((*x ^ (*x + r1_mask)) + r1_mask) >> 2;
     694         118 :   if (r1_mask == r01_mask) { *x += r1_mask; return; }
     695         118 :   if (r1_mask == 1) { *x += r01_mask; return; }
     696          73 :   if (r1_mask == 2) { *x += (r01_mask>>1) + 1; return; }
     697         121 :   while (r01_mask > r1_mask && slider < r1_mask)
     698             :   {
     699          53 :     r01_mask >>= 1; slider <<= 1;
     700             :   }
     701          34 :   *x += r01_mask + slider - 1;
     702          34 :   return;
     703             : }
     704             : 
     705             : /* self-init (1): advancing the bit pattern, and choice of primes for A.
     706             :  * The first time this is called, it finds h->bin_index == 0, which tells us
     707             :  * to initialize things from scratch.  On later occasions, we need to begin
     708             :  * by clearing the MPQS_FBE_DIVIDES_A bit in the fbe_flags of the former
     709             :  * prime factors of A.  We use, of course, the per_A_pr array for finding
     710             :  * them.  Upon successful return, that array will have been filled in, and
     711             :  * the flag bits will have been turned on again in the right places.
     712             :  * We return 1 (true) when we could set things up successfully, and 0 when
     713             :  * we found we'd be using more bits to the left in bin_index than we have
     714             :  * matching primes for in the FB.  In the latter case, bin_index will be
     715             :  * zeroed out, index2_FB will be incremented by 2, index2_moved will be
     716             :  * turned on, and the caller, after checking that index2_FB has not become
     717             :  * too large, should just call us again, which then is guaranteed to succeed:
     718             :  * we'll start again with a right-justified sequence of 1 bits in bin_index,
     719             :  * now interpreted as selecting primes relative to the new index2_FB. */
     720             : #ifndef MPQS_DEBUG_SI_CHOOSE_PRIMES
     721             : #  define MPQS_DEBUG_SI_CHOOSE_PRIMES 0
     722             : #endif
     723             : INLINE int
     724        2133 : mpqs_si_choose_primes(mpqs_handle_t *h)
     725             : {
     726        2133 :   mpqs_FB_entry_t *FB = h->FB;
     727        2133 :   mpqs_per_A_prime_t *per_A_pr = h->per_A_pr;
     728        2133 :   double l2_last_p = h->l2_target_A;
     729        2133 :   mpqs_int32_t omega_A = h->omega_A;
     730             :   int i, j, v2, prev_last_p_idx;
     731        2133 :   int room = h->index2_FB - h->index0_FB - omega_A + 4;
     732             :   /* GN 20050723:  I.e., index2_FB minus (index0_FB + omega_A - 3) plus 1
     733             :    * The notion of room here (cf mpqs_locate_A_range() above) is the number
     734             :    * of primes at or below index2_FB which are eligible for A.
     735             :    * At the very least, we need omega_A - 1 of them, and it is guaranteed
     736             :    * by mpqs_locate_A_range() that at least this many are available when we
     737             :    * first get here.  The logic here ensures that the lowest FB slot used
     738             :    * for A is never less than index0_FB + omega_A - 3.  In other words, when
     739             :    * omega_A == 3 (very small kN), we allow ourselves to reach all the way
     740             :    * down to index0_FB;  otherwise, we keep away from it by at least one
     741             :    * position.  For omega_A >= 4 this avoids situations where the selection
     742             :    * of the smaller primes here has advanced to a lot of very small ones, and
     743             :    * the single last larger one has soared away to bump into the top end of
     744             :    * the FB. */
     745             :   mpqs_uint32_t room_mask;
     746             :   mpqs_int32_t p;
     747             :   ulong bits;
     748             : 
     749             :   /* XXX also clear the index_j field here? */
     750        2133 :   if (h->bin_index == 0)
     751             :   { /* first time here, or after increasing index2_FB, initialize to a pattern
     752             :      * of omega_A - 1 consecutive right-justified 1 bits.
     753             :      * Caller will have ensured that there are enough primes for this in the
     754             :      * FB below index2_FB. */
     755         279 :     h->bin_index = (1UL << (omega_A - 1)) - 1;
     756         279 :     prev_last_p_idx = 0;
     757             :   }
     758             :   else
     759             :   { /* clear out the old flags */
     760       10730 :     for (i = 0; i < omega_A; i++)
     761        8876 :       MPQS_FLG(i) &= ~MPQS_FBE_DIVIDES_A;
     762        1854 :     prev_last_p_idx = MPQS_I(omega_A-1);
     763             : 
     764             :     /* find out how much maneuvering room we have before we're up against
     765             :      * the index0_FB wall */
     766        1854 :     if (room > 30) room = 30;
     767        1854 :     room_mask = ~((1UL << room) - 1);
     768             : 
     769             :     /* bump bin_index to the next acceptable value.  If index2_moved is off,
     770             :      * call mpqs_increment() just once;  otherwise, repeat until there's
     771             :      * something in the least significant 2 bits - this to ensure that we
     772             :      * never re-use an A which we'd used before increasing index2_FB - but
     773             :      * also stop if something shows up in the forbidden bits on the left
     774             :      * where we'd run out of bits or out of subscripts  (i.e. walk beyond
     775             :      * index0_FB + omega_A - 3). */
     776        1854 :     mpqs_increment(&h->bin_index);
     777        1854 :     if (h->index2_moved)
     778             :     {
     779           0 :       while ((h->bin_index & (room_mask | 0x3)) == 0)
     780           0 :         mpqs_increment(&h->bin_index);
     781             :     }
     782             :     /* ok so did we fall off the edge on the left? */
     783        1854 :     if ((h->bin_index & room_mask) != 0)
     784             :     {
     785             :       /* Yes.  Turn on the index2_moved flag in the handle */
     786           0 :       h->index2_FB += 2;        /* caller to check this isn't too large!!! */
     787           0 :       h->index2_moved = 1;
     788           0 :       h->bin_index = 0;
     789           0 :       if (MPQS_DEBUG_SI_CHOOSE_PRIMES || (MPQS_DEBUGLEVEL >= 5))
     790           0 :         err_printf("MPQS: wrapping, more primes for A now chosen near FB[%ld] = %ld\n",
     791           0 :                    (long)h->index2_FB,
     792           0 :                    (long)FB[h->index2_FB].fbe_p);
     793           0 :       return 0;                 /* back off - caller should retry */
     794             :     }
     795             :   }
     796             :   /* assert: we aren't occupying any of the room_mask bits now, and if
     797             :    * index2_moved had already been on, at least one of the two LSBs is on */
     798        2133 :   bits = h->bin_index;
     799        2133 :   if (MPQS_DEBUG_SI_CHOOSE_PRIMES || (MPQS_DEBUGLEVEL >= 6))
     800           0 :     err_printf("MPQS: new bit pattern for primes for A: 0x%lX\n", bits);
     801             : 
     802             :   /* map bits to FB subscripts, counting downward with bit 0 corresponding
     803             :    * to index2_FB, and accumulate logarithms against l2_last_p */
     804        2133 :   j = h->index2_FB;
     805        2133 :   v2 = vals((long)bits);
     806        2133 :   if (v2) { j -= v2; bits >>= v2; }
     807        7772 :   for (i = omega_A - 2; i >= 0; i--)
     808             :   {
     809        7772 :     MPQS_I(i) = j;
     810        7772 :     l2_last_p -= MPQS_LP(i);
     811        7772 :     MPQS_FLG(i) |= MPQS_FBE_DIVIDES_A;
     812        7772 :     bits &= ~1UL;
     813        7772 :     if (!bits) break;           /* that was the i=0 iteration */
     814        5639 :     v2 = vals((long)bits);
     815        5639 :     j -= v2;
     816        5639 :     bits >>= v2;
     817             :   }
     818             :   /* Choose the larger prime.  Note we keep index2_FB <= size_of_FB - 3 */
     819       37052 :   for (j = h->index2_FB + 1; (p = FB[j].fbe_p) != 0; j++)
     820       37052 :     if (FB[j].fbe_flogp > l2_last_p) break;
     821             :   /* GN 20050724: The following trick avoids generating a relatively large
     822             :    * proportion of duplicate relations when the last prime happens to fall
     823             :    * into an area where there are large gaps from one FB prime to the next,
     824             :    * and would otherwise often be repeated  (so that successive A's would
     825             :    * wind up too similar to each other).  While this trick isn't perfect,
     826             :    * it seems to get rid of a major part of the potential duplication. */
     827        2133 :   if (p != 0 && j == prev_last_p_idx) { j++; p = FB[j].fbe_p; }
     828        4266 :   MPQS_I(omega_A - 1) = (p == 0 ? /* did we fall off the end of the FB? */
     829        2133 :                          h->size_of_FB + 1 : /* then improvise */
     830             :                          j);
     831        2133 :   MPQS_FLG(omega_A - 1) |= MPQS_FBE_DIVIDES_A;
     832             : 
     833        2133 :   if (MPQS_DEBUG_SI_CHOOSE_PRIMES || (MPQS_DEBUGLEVEL >= 6))
     834             :   {
     835           0 :     err_printf("MPQS: chose primes for A");
     836           0 :     for (i = 0; i < omega_A; i++)
     837           0 :       err_printf(" FB[%ld]=%ld%s", (long) MPQS_I(i), (long) MPQS_AP(i),
     838           0 :                  i < omega_A - 1 ? "," : "\n");
     839             :   }
     840        2133 :   return 1;
     841             : }
     842             : 
     843             : /* compute coefficients of the sieving polynomial for self initializing
     844             :  * variant. Coefficients A and B are returned and several tables are
     845             :  * updated. A and B are kept on the PARI stack in preallocated GENs; so is C
     846             :  * when MPQS_DEBUG */
     847             : static void
     848       46985 : mpqs_self_init(mpqs_handle_t *h)
     849             : {
     850       46985 :   const ulong size_of_FB = h->size_of_FB + 1;
     851       46985 :   mpqs_FB_entry_t *FB = h->FB;
     852       46985 :   mpqs_inv_A_H_t *inv_A_H = h->inv_A_H;
     853       46985 :   const pari_sp av = avma;
     854             :   GEN p1, p2;
     855       46985 :   GEN A = h->A;
     856       46985 :   GEN B = h->B;
     857       46985 :   mpqs_per_A_prime_t *per_A_pr = h->per_A_pr;
     858             :   long i, j;
     859             :   long inv_A2;
     860             :   ulong p;
     861             : 
     862             : #ifdef MPQS_DEBUG_AVMA
     863             :   err_printf("MPQS DEBUG: enter self init, avma = 0x%lX\n", (ulong)avma);
     864             : #endif
     865             : 
     866             :   /* when all of the B's have already been used, choose new A ;
     867             :      this is indicated by setting index_j to 0 */
     868       46985 :   if (++h->index_j == (mpqs_uint32_t)h->no_B)
     869             :   {
     870        1854 :     h->index_j = 0;
     871        1854 :     h->index_i++; /* count finished A's */
     872             :   }
     873             : 
     874       46985 :   if (h->index_j == 0)
     875             :   /* "mpqs_self_init_A()" */
     876             :   { /* compute first polynomial with new A */
     877        2133 :     if (!mpqs_si_choose_primes(h))
     878             :     {
     879             :       /* We ran out of room towards small primes, and index2_FB was raised.
     880             :        * Check that we're still ok in that direction before re-trying the
     881             :        * operation, which then is guaranteed to succeed. The invariant we
     882             :        * maintain towards the top end is that h->size_of_FB - h->index2_FB >= 3,
     883             :        * but note that our size_of_FB is one larger. */
     884             : 
     885             :       /* "throw exception up to caller." ( bin_index set to impossible value
     886             :        * 0 by mpqs_si_choose_primes() */
     887           0 :       if (size_of_FB - h->index2_FB < 4) return;
     888           0 :       (void) mpqs_si_choose_primes(h);
     889             :     }
     890             :     /* assert: bin_index and per_A_pr are now populated with consistent
     891             :      * values */
     892             : 
     893             :     /* compute A = product of omega_A primes given by bin_index */
     894        2133 :     p1 = NULL;
     895       12038 :     for (i = 0; i < h->omega_A; i++)
     896             :     {
     897        9905 :       p = (ulong) MPQS_AP(i);
     898        9905 :       p1 = p1 ? muliu(p1, p): utoipos(p);
     899             :     }
     900        2133 :     affii(p1, A); set_avma(av);
     901             : 
     902             :     /* Compute H[i], 0 <= i < omega_A.  Also compute the initial
     903             :      * B = sum(v_i*H[i]), by taking all v_i = +1 */
     904             :     /* TODO: following needs to be changed later for segmented FB and sieve
     905             :      * interval, where we'll want to precompute several B's. */
     906        2133 :     p2 = NULL;
     907       12038 :     for (i = 0; i < h->omega_A; i++)
     908             :     {
     909        9905 :       p = (ulong) MPQS_AP(i);
     910        9905 :       p1 = divis(A, (long)p);
     911        9905 :       p1 = muliu(p1, Fl_inv(umodiu(p1, p), p));
     912        9905 :       p1 = muliu(p1, MPQS_SQRT(i));
     913        9905 :       affii(remii(p1, A), MPQS_H(i));
     914        9905 :       p2 = p2 ? addii(p2, MPQS_H(i)) : MPQS_H(i);
     915             :     }
     916        2133 :     affii(p2, B);
     917        2133 :     set_avma(av);                  /* must happen outside the loop */
     918             : 
     919             :     /* ensure B = 1 mod 4 */
     920        2133 :     if (mod2(B) == 0)
     921        1056 :       affii(addii(B, mului(mod4(A), A)), B); /* B += (A % 4) * A; */
     922             : 
     923        2133 :     p1 = shifti(A, 1);
     924             :     /* compute the roots z1, z2, of the polynomial Q(x) mod p_j and
     925             :      * initialize start1[i] with the first value p_i | Q(z1 + i p_j)
     926             :      * initialize start2[i] with the first value p_i | Q(z2 + i p_j)
     927             :      * The following loop "miraculously" does The Right Thing for the
     928             :      * primes dividing k (where sqrt_kN is 0 mod p).  Primes dividing A
     929             :      * are skipped here, and are handled further down in the common part
     930             :      * of SI. */
     931     2044489 :     for (j = 3; (ulong)j <= size_of_FB; j++)
     932             :     {
     933             :       ulong mb, tmp1, tmp2, m;
     934     2042356 :       if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
     935     2032451 :       p = (ulong)FB[j].fbe_p; m = h->M % p;
     936     2032451 :       inv_A2 = Fl_inv(umodiu(p1, p), p); /* = 1/(2*A) mod p_j */
     937     2032451 :       mb = umodiu(B, p); if (mb) mb = p - mb;
     938             :       /* mb = -B mod p */
     939     2032451 :       tmp1 = Fl_sub(mb, FB[j].fbe_sqrt_kN, p);
     940     2032451 :       tmp1 = Fl_mul(tmp1, inv_A2, p);
     941     2032451 :       FB[j].fbe_start1 = (mpqs_int32_t)Fl_add(tmp1, m, p);
     942             : 
     943     2032451 :       tmp2 = Fl_add(mb, FB[j].fbe_sqrt_kN, p);
     944     2032451 :       tmp2 = Fl_mul(tmp2, inv_A2, p);
     945     2032451 :       FB[j].fbe_start2 = (mpqs_int32_t)Fl_add(tmp2, m, p);
     946    13054755 :       for (i = 0; i < h->omega_A - 1; i++)
     947             :       {
     948    11022304 :         ulong h = umodiu(MPQS_H(i), p) << 1; if (h > p) h -= p;
     949    11022304 :         MPQS_INV_A_H(i,j) = Fl_mul(h, inv_A2, p); /* 1/A * H[i] mod p_j */
     950             :       }
     951             :     }
     952             :   }
     953             :   else
     954             :   /* "mpqs_self_init_B()" */
     955             :   { /* no "real" computation -- use recursive formula */
     956             :     /* The following exploits that B is the sum of omega_A terms +-H[i].
     957             :      * Each time we switch to a new B, we choose a new pattern of signs;
     958             :      * the precomputation of the inv_A_H array allows us to change the
     959             :      * two arithmetic progressions equally fast.  The choice of sign
     960             :      * patterns does *not* follow the bit pattern of the ordinal number
     961             :      * of B in the current cohort;  rather, we use a Gray code, changing
     962             :      * only one sign each time.  When the i-th rightmost bit of the new
     963             :      * ordinal number index_j of B is 1, the sign of H[i] is changed;
     964             :      * the next bit to the left tells us whether we should be adding or
     965             :      * subtracting the difference term.  We never need to change the sign
     966             :      * of H[omega_A-1] (the topmost one), because that would just give us
     967             :      * the same sieve items Q(x) again with the opposite sign of x.  This
     968             :      * is why we only precomputed inv_A_H up to i = omega_A - 2. */
     969             : 
     970       44852 :     ulong v2 = 0;               /* 2-valuation of h->index_j */
     971             : 
     972       44852 :     j = h->index_j;
     973             :     /* could use vals() here, but we need to right shift the bit pattern
     974             :      * anyway in order to find out which inv_A_H entries must be added to or
     975             :      * subtracted from the modular roots */
     976       44852 :     while ((j & 1) == 0) { v2++; j >>= 1; }
     977             : 
     978             :     /* v2 = v_2(index_j), determine new starting positions for sieving */
     979       44852 :     p1 = shifti(MPQS_H(v2), 1);
     980       44852 :     if (j & 2)
     981             :     { /* j = 3 mod 4 */
     982    49647810 :       for (j = 3; (ulong)j <= size_of_FB; j++)
     983             :       {
     984    49626526 :         if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
     985    49489795 :         p = (ulong)FB[j].fbe_p;
     986    49489795 :         FB[j].fbe_start1 = Fl_sub(FB[j].fbe_start1, MPQS_INV_A_H(v2,j), p);
     987    49489795 :         FB[j].fbe_start2 = Fl_sub(FB[j].fbe_start2, MPQS_INV_A_H(v2,j), p);
     988             :       }
     989       21284 :       p1 = addii(B, p1);
     990             :     }
     991             :     else
     992             :     { /* j = 1 mod 4 */
     993    51827324 :       for (j = 3; (ulong)j <= size_of_FB; j++)
     994             :       {
     995    51803756 :         if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
     996    51656357 :         p = (ulong)FB[j].fbe_p;
     997    51656357 :         FB[j].fbe_start1 = Fl_add(FB[j].fbe_start1, MPQS_INV_A_H(v2,j), p);
     998    51656357 :         FB[j].fbe_start2 = Fl_add(FB[j].fbe_start2, MPQS_INV_A_H(v2,j), p);
     999             :       }
    1000       23568 :       p1 = subii(B, p1);
    1001             :     }
    1002       44852 :     affii(p1, B);
    1003             :   }
    1004       46985 :   set_avma(av);
    1005             : 
    1006             :   /* p=2 is a special case.  start1[2], start2[2] are never looked at,
    1007             :    * so don't bother setting them. */
    1008             : 
    1009             :   /* now compute zeros of polynomials that have only one zero mod p
    1010             :      because p divides the coefficient A */
    1011             : 
    1012       46985 :   p1 = diviiexact(subii(h->kN, sqri(B)), shifti(A, 2)); /* coefficient -C */
    1013      341020 :   for (i = 0; i < h->omega_A; i++)
    1014             :   {
    1015             :     ulong tmp, s;
    1016      294035 :     p = (ulong) MPQS_AP(i);
    1017      294035 :     tmp = Fl_div(umodiu(p1, p), umodiu(B, p), p); s = (tmp + h->M) % p;
    1018      294035 :     FB[MPQS_I(i)].fbe_start1 = (mpqs_int32_t)s;
    1019      294035 :     FB[MPQS_I(i)].fbe_start2 = (mpqs_int32_t)s;
    1020             :   }
    1021             : 
    1022       46985 :   if (MPQS_DEBUGLEVEL >= 6)
    1023           0 :     err_printf("MPQS: chose Q_%ld(x) = %Ps x^2 %c %Ps x + C\n",
    1024           0 :                (long) h->index_j, h->A,
    1025           0 :                signe(h->B) < 0? '-': '+', absi_shallow(h->B));
    1026             : 
    1027       46985 :   set_avma(av);
    1028             : 
    1029             : #ifdef MPQS_DEBUG
    1030             :   /* stash C into the handle.  Since check_root() is the only thing which
    1031             :    * uses it, and only for debugging, C doesn't even exist as a field in
    1032             :    * the handle unless we're built with MPQS_DEBUG. */
    1033             :   affii(negi(p1), h->C);
    1034             :   for (j = 3; j <= size_of_FB; j++)
    1035             :   {
    1036             :     check_root(h, FB[j].fbe_p, FB[j].fbe_start1);
    1037             :     check_root(h, FB[j].fbe_p, FB[j].fbe_start2); set_avma(av);
    1038             :   }
    1039             :   if (DEBUGLEVEL >= 6)
    1040             :     PRINT_IF_VERBOSE("MPQS: checking of roots of Q(x) was successful\n");
    1041             : #endif
    1042             : 
    1043             : #ifdef MPQS_DEBUG_AVMA
    1044             :   err_printf("MPQS DEBUG: leave self init, avma = 0x%lX\n", (ulong)avma);
    1045             : #endif
    1046             : }
    1047             : 
    1048             : /*********************************************************************/
    1049             : /**                                                                 **/
    1050             : /**                           THE SIEVE                             **/
    1051             : /**                                                                 **/
    1052             : /*********************************************************************/
    1053             : 
    1054             : /* Main sieving routine:
    1055             :  * p4 = 4*p (used for loop unrolling)
    1056             :  * log_p: approximation for log(p)
    1057             :  * begin: points to a sieve array
    1058             :  * end: points to the end of the sieve array
    1059             :  * starting_sieving_index: marks the first FB element used for sieving */
    1060             : INLINE void
    1061   206330201 : mpqs_sieve_p(unsigned char *begin, unsigned char *end,
    1062             :              long p4, long p, unsigned char log_p)
    1063             : {
    1064   206330201 :   register unsigned char *e = end - p4;
    1065             :   /* Loop unrolled some time ago. It might be better to let the compiler worry
    1066             :    * about *this* kind of optimization, based on its knowledge of whatever
    1067             :    * useful tricks the machine instruction set architecture is offering
    1068             :    * ("speculative loads" being the buzzword). --GN */
    1069   917149128 :   while (e - begin >= 0) /* signed comparison */
    1070             :   {
    1071   504488726 :     (*begin) += log_p, begin += p;
    1072   504488726 :     (*begin) += log_p, begin += p;
    1073   504488726 :     (*begin) += log_p, begin += p;
    1074   504488726 :     (*begin) += log_p, begin += p;
    1075             :   }
    1076   886308685 :   while (end - begin >= 0) /* again */
    1077   473648283 :     (*begin) += log_p, begin += p;
    1078   206330201 : }
    1079             : 
    1080             : static void
    1081       46985 : mpqs_sieve(mpqs_handle_t *h)
    1082             : {
    1083       46985 :   long p, l = h->index1_FB;
    1084             :   mpqs_FB_entry_t *ptr_FB;
    1085       46985 :   unsigned char *sieve_array = h->sieve_array;
    1086       46985 :   unsigned char *sieve_array_end = h->sieve_array_end;
    1087             : 
    1088   103359103 :   for (ptr_FB = &(h->FB[l]); (p = ptr_FB->fbe_p) != 0; ptr_FB++, l++)
    1089             :   {
    1090   103312118 :     unsigned char log_p = ptr_FB->fbe_logval;
    1091   103312118 :     long start1 = ptr_FB->fbe_start1;
    1092   103312118 :     long start2 = ptr_FB->fbe_start2;
    1093             : 
    1094             :     /* sieve with FB[l] from start_1[l] */
    1095             :     /* if start1 != start2 sieve with FB[l] from start_2[l] */
    1096             :     /* Maybe it is more efficient not to have a conditional branch in
    1097             :      * the present loop body, and instead to right-shift log_p one bit
    1098             :      * based on a flag bit telling us that we're on a one-root prime?
    1099             :      * And instead roll the two invocations of mpqs_sieve_p into one. */
    1100   103312118 :     mpqs_sieve_p(sieve_array + start1, sieve_array_end, p << 2, p, log_p);
    1101   103312118 :     if (start1 != start2)
    1102   103018083 :       mpqs_sieve_p(sieve_array + start2, sieve_array_end, p << 2, p, log_p);
    1103             :   }
    1104       46985 : }
    1105             : 
    1106             : /******************************/
    1107             : 
    1108             : /* Could make shameless use of the fact that M is divisible by 4, but
    1109             :  * let the compiler worry about loop unrolling.  Indeed I wonder whether
    1110             :  * modern compilers woudln't have an easier time optimizing this if it
    1111             :  * were written as array accesses.  Doing so. */
    1112             : static long
    1113       46985 : mpqs_eval_sieve(mpqs_handle_t *h)
    1114             : {
    1115       46985 :   long x = 0, count = 0, M_2 = h->M << 1;
    1116       46985 :   unsigned char th = h->sieve_threshold;
    1117       46985 :   unsigned char *sieve_array = h->sieve_array;
    1118       46985 :   long *candidates = h->candidates;
    1119             : 
    1120             :   /* The following variation on the original is marginally faster with a
    1121             :    * good optimizing compiler.  Exploiting the sentinel, we don't need to
    1122             :    * check for x < M_2 in the inner while loop - this more than makes up
    1123             :    * for the "lack" of explicit unrolling.  Optimizations like loop
    1124             :    * unrolling are best left to the compiler anyway... */
    1125      301541 :   while (count < MPQS_CANDIDATE_ARRAY_SIZE - 1)
    1126             :   {
    1127      254556 :     while (sieve_array[x] < th) x++;
    1128      254556 :     if (x >= M_2) break;
    1129      207571 :     candidates[count++] = x++;
    1130             :   }
    1131             : 
    1132       46985 :   candidates[count] = 0; return count;
    1133             : }
    1134             : 
    1135             : /*********************************************************************/
    1136             : /**                                                                 **/
    1137             : /**                     CONSTRUCTING RELATIONS                      **/
    1138             : /**                                                                 **/
    1139             : /*********************************************************************/
    1140             : 
    1141             : /* Main relation routine */
    1142             : static void
    1143     3230554 : mpqs_add_factor(GEN relp, long *i, ulong ei, ulong pi)
    1144             : {
    1145     3230554 :   ++*i;
    1146     3230554 :   relp[*i]=pi | (ei<<REL_OFFSET);
    1147     3230554 : }
    1148             : 
    1149             : /* only used for debugging */
    1150             : static void
    1151           0 : split_relp(GEN rel, GEN *pt_relp, GEN *pt_relc)
    1152             : {
    1153           0 :   long j, l = lg(rel);
    1154           0 :   GEN relp = cgetg(l, t_VECSMALL);
    1155           0 :   GEN relc = cgetg(l, t_VECSMALL);
    1156           0 :   for (j=1; j<l; j++)
    1157             :   {
    1158           0 :     relc[j] = rel[j]>>REL_OFFSET;
    1159           0 :     relp[j] = rel[j]&REL_MASK;
    1160             :   }
    1161           0 :   *pt_relp = relp;
    1162           0 :   *pt_relc = relc;
    1163           0 : }
    1164             : 
    1165             : #ifdef MPQS_DEBUG
    1166             : static GEN
    1167             : mpqs_factorback(mpqs_handle_t *h, GEN relp)
    1168             : {
    1169             :   GEN N = h->N, prod = gen_1;
    1170             :   long i, j, l = lg(relp);
    1171             :   mpqs_FB_entry_t *FB = h->FB;
    1172             : 
    1173             :   for (j=1; j<l; j++)
    1174             :   {
    1175             :     long e = relp[j]>>REL_OFFSET;
    1176             :     i = relp[j]&REL_MASK;
    1177             :     /* special case -1 */
    1178             :     if (i == 1) prod = Fp_neg(prod,N);
    1179             :     else prod = Fp_mul(prod, Fp_powu(utoipos(FB[i].fbe_p), e, N), N);
    1180             :   }
    1181             :   return prod;
    1182             : }
    1183             : #endif
    1184             : 
    1185             : /* NB FREL, LPREL are actually FNEW, LPNEW when we get called */
    1186             : static long
    1187       45725 : mpqs_eval_cand(mpqs_handle_t *h, long number_of_cand,
    1188             :                GEN *FREL, GEN *LPREL)
    1189             : {
    1190       45725 :   pari_sp av = avma;
    1191       45725 :   long number_of_relations = 0;
    1192       45725 :   long *relaprimes = h->relaprimes;
    1193             :   ulong i, pi;
    1194       45725 :   mpqs_FB_entry_t *FB = h->FB;
    1195       45725 :   GEN A = h->A;
    1196       45725 :   GEN B = h->B;                 /* we don't need coefficient C here */
    1197             :   int pii;
    1198       45725 :   long *candidates = h->candidates;
    1199             :   long nb;
    1200             :   GEN frel, lprel;
    1201       45725 :   long nfrel = 1, nlprel = 1;
    1202       45725 :   frel = cgetg(DEFAULT_VEC_LEN, t_VEC);
    1203       45725 :   lprel = cgetg(DEFAULT_VEC_LEN, t_VEC);
    1204      253296 :   for (i = 0; i < (ulong)number_of_cand; i++)
    1205             :   {
    1206      207571 :     pari_sp btop = avma;
    1207             :     GEN Qx, Qx_part, A_2x_plus_B, Y;
    1208             :     long powers_of_2, p;
    1209      207571 :     long x = candidates[i];
    1210      207571 :     long x_minus_M = x - h->M;
    1211      207571 :     int relaprpos = 0;
    1212      207571 :     GEN relp = cgetg(MAX_PE_PAIR+1,t_VECSMALL);
    1213             : 
    1214      207571 :     nb=0;
    1215             : 
    1216             : #ifdef MPQS_DEBUG_VERYVERBOSE
    1217             :     err_printf("%c", (char)('0' + i%10));
    1218             : #endif
    1219             : 
    1220             :     /* A_2x_plus_B = (A*(2x)+B), Qx = (A*(2x)+B)^2/(4*A) = Q(x) */
    1221      207571 :     A_2x_plus_B = addii(mulis(A, 2 * x_minus_M), B);
    1222      207571 :     Y = absi_shallow(A_2x_plus_B);
    1223             : 
    1224      207571 :     Qx = subii(sqri(A_2x_plus_B), h->kN);
    1225             : 
    1226             :     /* When N is relatively small, it may happen that Qx is outright
    1227             :      * divisible by N at this point.  In any case, when no extensive prior
    1228             :      * trial division / Rho / ECM had been attempted, gcd(Qx,N) may turn
    1229             :      * out to be a nontrivial factor of N  (larger than what the FB contains
    1230             :      * or we'd have found it already, but possibly smaller than the large-
    1231             :      * prime bound).  This is too rare to check for here in the inner loop,
    1232             :      * but it will be caught if such an LP relation is ever combined with
    1233             :      * another. */
    1234             : 
    1235             :     /* Qx cannot possibly vanish here */
    1236      207571 :     if (!signe(Qx)) { PRINT_IF_VERBOSE("<+>"); continue; }
    1237      207571 :     else if (signe(Qx) < 0) {
    1238      113524 :       setabssign(Qx);
    1239      113524 :       mpqs_add_factor(relp, &nb, 1, 1); /* i = 1, ei = 1, pi */
    1240             :     }
    1241             : 
    1242             :     /* divide by powers of 2;  we're really dealing with 4*A*Q(x), so we
    1243             :      * always have at least 2^2 here, and at least 2^3 when kN is 1 mod 4 */
    1244      207571 :     powers_of_2 = vali(Qx);
    1245      207571 :     Qx = shifti(Qx, -powers_of_2);
    1246      207571 :     mpqs_add_factor(relp, &nb, powers_of_2, 2); /* i = 1, ei = 1, pi */
    1247             : 
    1248             :     /* That has dealt with a possible -1 and the power of 2.  First pass
    1249             :      * over odd primes in FB: pick up all possible divisors of Qx including
    1250             :      * those sitting in k or in A, and remember them in relaprimes. Do not
    1251             :      * yet worry about possible repeated factors, these will be found in the
    1252             :      * second pass. */
    1253      207571 :     Qx_part = A;
    1254             : 
    1255             :     /* The first pass recognizes divisors of A by their corresponding flags
    1256             :      * bit in the FB entry.  (Divisors of k require no special treatment at
    1257             :      * this stage.)  We construct a preliminary table of FB subscripts and
    1258             :      * "exponents" of the FB primes which divide Qx.  (We store subscripts
    1259             :      * rather than the primes themselves because the string representation
    1260             :      * of a relation is in terms of the subscripts.)
    1261             :      * We must distinguish three cases so we can do the right thing in the
    1262             :      * 2nd pass: prime not in A which divides Qx, prime in A which does not
    1263             :      * divide Qx/A, prime in A which does divide Qx/A. The first and third
    1264             :      * kinds need checking for repeated factors, the second kind doesn't. The
    1265             :      * first and second kinds contribute 1 to the exponent in the relation,
    1266             :      * the 3rd kind contributes 2. We store 1,0,2 respectively in these three
    1267             :      * cases.
    1268             :      * Factors in common with k are much simpler - if they occur, they occur
    1269             :      * exactly to the first power, and this makes no difference in the first
    1270             :      * pass - here they behave just like every normal odd factor base prime.
    1271             :      */
    1272             : 
    1273   381335500 :     for (pi = 3; (p = FB[pi].fbe_p); pi++)
    1274             :     {
    1275   381127929 :       long tmp_p = x % p;
    1276   381127929 :       ulong ei = 0;
    1277             : 
    1278             :       /* Here we use that MPQS_FBE_DIVIDES_A equals 1. */
    1279   381127929 :       ei = FB[pi].fbe_flags & MPQS_FBE_DIVIDES_A;
    1280             : 
    1281   381127929 :       if (tmp_p == FB[pi].fbe_start1 || tmp_p == FB[pi].fbe_start2)
    1282             :       { /* p divides Q(x)/A (and possibly A), 1st or 3rd case */
    1283     1452174 :         relaprimes[relaprpos++] = pi;
    1284     1452174 :         relaprimes[relaprpos++] = 1 + ei;
    1285     1452174 :         Qx_part = muliu(Qx_part, p);
    1286             :       }
    1287   379675755 :       else if (ei)
    1288             :       { /* p divides A but does not divide Q(x)/A, 2nd case */
    1289     1207591 :         relaprimes[relaprpos++] = pi;
    1290     1207591 :         relaprimes[relaprpos++] = 0;
    1291             :       }
    1292             :     }
    1293             : 
    1294             :     /* We have now accumulated the known factors of Qx except for possible
    1295             :      * repeated factors and for possible large primes.  Divide off what we
    1296             :      * have.  (This is faster than dividing off A and each prime separately.)
    1297             :      */
    1298      207571 :     Qx = diviiexact(Qx, Qx_part);
    1299             :     /* (ToDo: MPQS_DEBUG sanity check...) */
    1300             : 
    1301             : #ifdef MPQS_DEBUG_AVMA
    1302             :     err_printf("MPQS DEBUG: eval loop 3, avma = 0x%lX\n", (ulong)avma);
    1303             : #endif
    1304             : 
    1305             :     /* second pass - deal with any repeated factors, and write out the string
    1306             :      * representation of the tentative relation. At this point, the only
    1307             :      * primes which can occur again in the adjusted Qx are those in relaprimes
    1308             :      * which are followed by 1 or 2.  We must pick up those followed by a 0,
    1309             :      * too, though. */
    1310             :     PRINT_IF_VERBOSE("a");
    1311     2867336 :     for (pii = 0; pii < relaprpos; pii+=2)
    1312             :     {
    1313             :       long remd_p;
    1314     2659765 :       ulong ei = relaprimes[pii+1];
    1315             :       GEN Qx_div_p;
    1316     2659765 :       pi = relaprimes[pii];
    1317             : 
    1318             :       /* p | k (identified by its position before index0_FB)
    1319             :        * or p | A (ei = 0) */
    1320     2659765 :       if ((mpqs_int32_t)pi < h->index0_FB || ei == 0) {
    1321     1224034 :         mpqs_add_factor(relp, &nb, 1, pi);
    1322     1224034 :         continue;
    1323             :       }
    1324     1435731 :       p = FB[pi].fbe_p;
    1325             : #ifdef MPQS_DEBUG_CANDIDATE_EVALUATION
    1326             :       err_printf("MPQS DEBUG: Qx=%Ps p=%ld\n", Qx, p);
    1327             : #endif
    1328             :       /* otherwise p might still divide the current adjusted Qx. Try it... */
    1329             :       /* NOTE: break out of loop when remaining Qx is 1 ?  Or rather, suppress
    1330             :        * the trial divisions, since we still need to write our string.
    1331             :        * Actually instead of testing for 1, test whether Qx is smaller than p;
    1332             :        * cf Karim's mail from 20050124.  If it is, without being 1, then it has
    1333             :        * a common factor with k.  But those factors are soon going to have
    1334             :        * disappeared before we get here.  However, inserting
    1335             :        * an explicit if (!is_pm1(Qx)) here did not help any. */
    1336     1435731 :       Qx_div_p = divis_rem(Qx, p, &remd_p);
    1337     2915762 :       while (remd_p == 0) {
    1338       44300 :         ei++; Qx = Qx_div_p;
    1339       44300 :         Qx_div_p = divis_rem(Qx, p, &remd_p);
    1340             :       }
    1341     1435731 :       mpqs_add_factor(relp, &nb, ei, pi);
    1342             :     }
    1343             : 
    1344             : #ifdef MPQS_DEBUG_AVMA
    1345             :     err_printf("MPQS DEBUG: eval loop 4, avma = 0x%lX\n", (ulong)avma);
    1346             : #endif
    1347             :     PRINT_IF_VERBOSE("\bb");
    1348      207571 :     if (is_pm1(Qx))
    1349             :     {
    1350             :       GEN rel;
    1351       90151 :       setlg(relp, nb+1);
    1352       90151 :       rel = gerepilecopy(btop, mkvec2(Y, relp));
    1353       90151 :       frel = vec_extend(frel, rel, nfrel++);
    1354       90151 :       number_of_relations++;
    1355             : 
    1356             : #ifdef MPQS_DEBUG
    1357             :       {
    1358             :         pari_sp av1 = avma;
    1359             :         GEN rhs = mpqs_factorback(h, relp);
    1360             :         GEN Y = gel(rel,1);
    1361             :         GEN Qx_2 = remii(sqri(Y), h->N);
    1362             :         GEN relpp, relpc;
    1363             :         split_relp(relp,&relpp,&relpc);
    1364             :         if (!equalii(Qx_2, rhs))
    1365             :         {
    1366             :           PRINT_IF_VERBOSE("\b(!)\n");
    1367             :           err_printf("MPQS: %Ps @ %Ps : %Ps %Ps\n", Y, Qx,relpp,relpc);
    1368             :           err_printf("\tQx_2 = %Ps\n", Qx_2);
    1369             :           err_printf("\t rhs = %Ps\n", rhs);
    1370             :           pari_err_BUG("MPQS: wrong full relation found");
    1371             :         }
    1372             :         else
    1373             :           PRINT_IF_VERBOSE("\b(:)");
    1374             :         set_avma(av1);
    1375             :       }
    1376             : #endif
    1377             :     }
    1378      117420 :     else if (cmpis(Qx, h->lp_bound) > 0)
    1379             :     { /* TODO: check for double large prime */
    1380             :       PRINT_IF_VERBOSE("\b.");
    1381       18485 :       set_avma(btop);
    1382             :     }
    1383             :     else
    1384             :     { /* if (mpqs_isprime(itos(Qx))) */
    1385             :       GEN rel;
    1386       98935 :       setlg(relp, nb+1);
    1387       98935 :       rel = gerepilecopy(btop, mkvec3(Qx,Y,relp));
    1388       98935 :       lprel = vec_extend(lprel, rel, nlprel++);
    1389             : #ifdef MPQS_DEBUG
    1390             :       {
    1391             :         pari_sp av1 = avma;
    1392             :         GEN rhs = mpqs_factorback(h, relp);
    1393             :         GEN Qx = gel(rel,1), Y = gel(rel,2);
    1394             :         GEN Qx_2 = remii(sqri(Y), h->N);
    1395             :         split_relp(relp,&relpp,&relpc);
    1396             : 
    1397             :         rhs = modii(mulii(rhs, Qx), h->N);
    1398             :         if (!equalii(Qx_2, rhs))
    1399             :         {
    1400             :           PRINT_IF_VERBOSE("\b(!)\n");
    1401             :           err_printf("MPQS: %Ps @ %Ps :%s %Ps %Ps\n", Y, Qx, relpp, relpc);
    1402             :           err_printf("\tQx_2 = %Ps\n", Qx_2);
    1403             :           err_printf("\t rhs = %Ps\n", rhs);
    1404             :           pari_err_BUG("MPQS: wrong large prime relation found");
    1405             :         }
    1406             :         else
    1407             :           PRINT_IF_VERBOSE("\b(;)");
    1408             :         set_avma(av1);
    1409             :       }
    1410             : #endif
    1411             :     }
    1412             : 
    1413             :   } /* for */
    1414             :   PRINT_IF_VERBOSE("\n");
    1415       45725 :   if (nfrel==1) frel=NULL;   else setlg(frel, nfrel);
    1416       45725 :   if (nlprel==1) lprel=NULL; else setlg(lprel, nlprel);
    1417       45725 :   *FREL = frel; *LPREL = lprel;
    1418       45725 :   if (!frel && !lprel) { set_avma(av); return 0; }
    1419       44798 :   if (!frel) *LPREL = gerepilecopy(av, lprel);
    1420       24999 :   else gerepileall(av, lprel ? 2: 1, FREL, LPREL);
    1421       44798 :   return number_of_relations;
    1422             : }
    1423             : 
    1424             : /*********************************************************************/
    1425             : /**                                                                 **/
    1426             : /**                      COMBINING RELATIONS                        **/
    1427             : /**                                                                 **/
    1428             : /*********************************************************************/
    1429             : 
    1430             : /* combines the large prime relations in COMB to full relations in FNEW.*/
    1431             : 
    1432             : static void
    1433      113410 : rel_to_ei(GEN ei, GEN relp)
    1434             : {
    1435      113410 :   long j, l = lg(relp);
    1436     1929744 :   for (j=1; j<l; j++)
    1437             :   {
    1438     1816334 :     long e = relp[j]>>REL_OFFSET;
    1439     1816334 :     long i = relp[j]&REL_MASK;
    1440     1816334 :     ei[i] += e;
    1441             :   }
    1442      113410 : }
    1443             : 
    1444             : static GEN
    1445       10320 : combine_large_primes(mpqs_handle_t *h, GEN rel1, GEN rel2)
    1446             : {
    1447       10320 :   pari_sp av = avma;
    1448       10320 :   GEN new_Y, new_Y1, Y1 = rel_Y(rel1), Y2 = rel_Y(rel2);
    1449       10320 :   long l, lei = h->size_of_FB + 1, nb = 0;
    1450       10320 :   GEN ei, relp, inv_q, q = rel_q(rel1);
    1451             : 
    1452       10320 :   if (!invmod(q, h->N, &inv_q)) /* can happen --GN */
    1453             :   {
    1454           0 :     inv_q = gcdii(inv_q, h->N);
    1455           0 :     if (is_pm1(inv_q) || equalii(inv_q, h->N)) /* pity */
    1456             :     {
    1457             : #ifdef MPQS_DEBUG
    1458             :       err_printf("MPQS: skipping relation with non-invertible q\n");
    1459             : #endif
    1460           0 :       set_avma(av); return NULL;
    1461             :     }
    1462           0 :     return inv_q;
    1463             :   }
    1464       10320 :   ei = zero_zv(lei);
    1465       10320 :   relp = cgetg(MAX_PE_PAIR+1,t_VECSMALL);
    1466             : 
    1467       10320 :   rel_to_ei(ei, rel_p(rel1));
    1468       10320 :   rel_to_ei(ei, rel_p(rel2));
    1469       10320 :   new_Y = modii(mulii(mulii(Y1, Y2), inv_q), h->N);
    1470       10320 :   new_Y1 = subii(h->N, new_Y);
    1471       10320 :   if (abscmpii(new_Y1, new_Y) < 0) new_Y = new_Y1;
    1472       10320 :   if (odd(ei[1]))
    1473        4919 :     mpqs_add_factor(relp, &nb, 1, 1);
    1474    27515799 :   for (l = 2; l <= lei; l++)
    1475    27505479 :     if (ei[l])
    1476      244775 :       mpqs_add_factor(relp, &nb, ei[l],l);
    1477       10320 :   if (DEBUGLEVEL >= 6)
    1478             :   {
    1479             :     GEN relpp, relpc;
    1480             :     GEN rel1p, rel1c;
    1481             :     GEN rel2p, rel2c;
    1482           0 :     split_relp(relp,&relpp,&relpc);
    1483           0 :     split_relp(rel1,&rel1p,&rel1c);
    1484           0 :     split_relp(rel2,&rel2p,&rel2c);
    1485           0 :     err_printf("MPQS: combining\n");
    1486           0 :     err_printf("    {%Ps @ %Ps : %Ps}\n", rel_q(rel1), Y1, rel1p, rel1c);
    1487           0 :     err_printf("  * {%Ps @ %Ps : %Ps}\n", rel_q(rel2), Y2, rel2p, rel2c);
    1488           0 :     err_printf(" == {%Ps, %Ps}\n", relpp, relpc);
    1489             :   }
    1490       10320 :   setlg(relp, nb+1);
    1491             : 
    1492             : #ifdef MPQS_DEBUG
    1493             :   {
    1494             :     GEN Qx_2, prod;
    1495             :     pari_sp av1 = avma;
    1496             : 
    1497             :     Qx_2 = modii(sqri(new_Y), h->N);
    1498             :     prod = mpqs_factorback(h, relp);
    1499             :     if (!equalii(Qx_2, prod))
    1500             :       pari_err_BUG("MPQS: combined large prime relation is false");
    1501             :     set_avma(av1);
    1502             :   }
    1503             : #endif
    1504       10320 :   return gerepilecopy(av, mkvec2(new_Y, relp));
    1505             : }
    1506             : 
    1507             : static GEN
    1508       35677 : mpqs_combine_large_primes(mpqs_handle_t *h, hashtable *lprel, GEN LPNEW, hashtable *frel)
    1509             : {
    1510       35677 :   long j, lpnew = lg(LPNEW);
    1511      134612 :   for (j = 1; j < lpnew; j++)
    1512             :   {
    1513       98935 :     GEN rel = gel(LPNEW,j);
    1514       98935 :     ulong q = itou(rel_q(rel));
    1515       98935 :     GEN col = hash_haskey_GEN(lprel, (void*)q);
    1516       98935 :     if (!col)
    1517       88615 :       hash_insert(lprel, (void*)q, (void*)rel);
    1518             :     else
    1519             :     {
    1520       10320 :       GEN f = combine_large_primes(h, rel, col);
    1521       10320 :       if (f)
    1522             :       {
    1523       10320 :         if (typ(f) == t_INT) return f;
    1524       10320 :         frel_add(frel, f);
    1525             :       }
    1526             :     }
    1527             :   }
    1528       35677 :   return NULL;
    1529             : }
    1530             : 
    1531             : /*********************************************************************/
    1532             : /**                                                                 **/
    1533             : /**                    FROM RELATIONS TO DIVISORS                   **/
    1534             : /**                                                                 **/
    1535             : /*********************************************************************/
    1536             : 
    1537             : /* create and read an F2m from a relation file FREL.
    1538             :  * Also record the position of each relation in the file for later use
    1539             :  * rows = size_of_FB+1, cols = rel */
    1540             : static GEN
    1541         279 : stream_read_F2m(GEN rel, long rows, long cols)
    1542             : {
    1543         279 :   long i, lr = lg(rel);
    1544             :   GEN m;
    1545         279 :   long space = 2*((nbits2nlong(rows)+3)*cols+1);
    1546         279 :   if ((long)((GEN)avma - (GEN)pari_mainstack->bot) < space)
    1547             :   {
    1548           0 :     pari_sp av = avma;
    1549           0 :     m = gclone(zero_F2m(rows, cols));
    1550           0 :     if (DEBUGLEVEL>=4)
    1551           0 :       err_printf("MPQS: allocating %ld words for Gauss\n",space);
    1552           0 :     set_avma(av);
    1553             :   }
    1554             :   else
    1555         279 :     m = zero_F2m_copy(rows, cols);
    1556      100078 :   for (i = 1; i < lr; i++)
    1557             :   {
    1558       99799 :     GEN relp = gmael(rel,i,2);
    1559       99799 :     long j, l = lg(relp);
    1560     1587335 :     for (j = 1; j < l; j++)
    1561     1487536 :       if (odd(relp[j]>>REL_OFFSET))
    1562     1358597 :          F2m_set(m, relp[j]&REL_MASK, i);
    1563             :   }
    1564         279 :   if (i-1 != cols)
    1565           0 :     pari_err_BUG(stack_sprintf("MPQS: full relations %s than expected",
    1566             :                i > cols ? "longer" : "shorter"));
    1567         279 :   return m;
    1568             : }
    1569             : 
    1570             : static GEN
    1571       92770 : mpqs_add_relation(GEN Y_prod, GEN N, GEN ei, GEN r)
    1572             : {
    1573       92770 :   pari_sp av = avma;
    1574       92770 :   GEN res = remii(mulii(Y_prod, gel(r,1)), N);
    1575       92770 :   rel_to_ei(ei, gel(r,2));
    1576       92770 :   return gerepileuptoint(av, res);
    1577             : }
    1578             : 
    1579             : static int
    1580         586 : split(GEN N, GEN *e, GEN *res)
    1581             : {
    1582             :   ulong mask;
    1583             :   long flag;
    1584             :   GEN base;
    1585         586 :   if (MR_Jaeschke(N)) { *e = gen_1; return 1; } /* probable prime */
    1586          14 :   if (Z_issquareall(N, &base))
    1587             :   { /* squares could cost us a lot of time */
    1588             :     /* GN20050707: as used now, this is always called with res!=NULL */
    1589           0 :     *res = base;
    1590           0 :     *e = gen_2;
    1591           0 :     if (DEBUGLEVEL >= 5) err_printf("MPQS: decomposed a square\n");
    1592           0 :     return 1;
    1593             :   }
    1594          14 :   mask = 7;
    1595             :   /* 5th/7th powers aren't worth the trouble. OTOH once we have the hooks for
    1596             :    * dealing with cubes, higher powers can be handled essentially for free) */
    1597          14 :   if ( (flag = is_357_power(N, &base, &mask)) )
    1598             :   {
    1599           0 :     *res = base;
    1600           0 :     *e = utoipos(flag);
    1601           0 :     if (DEBUGLEVEL >= 5)
    1602           0 :       err_printf("MPQS: decomposed a %s\n",
    1603             :                  (flag == 3 ? "cube" :
    1604           0 :                   (flag == 5 ? "5th power" : "7th power")));
    1605           0 :     return 1;
    1606             :   }
    1607          14 :   *e = gen_0; return 0; /* known composite */
    1608             : }
    1609             : 
    1610             : static GEN
    1611         279 : mpqs_solve_linear_system(mpqs_handle_t *h, GEN frel, long rel)
    1612             : {
    1613         279 :   GEN N = h->N, X, Y_prod, X_plus_Y, D1, res, new_res;
    1614         279 :   mpqs_FB_entry_t *FB = h->FB;
    1615         279 :   pari_sp av=avma, av2, av3;
    1616             :   GEN ei;
    1617             :   long i, j, H_cols, H_rows;
    1618             :   long res_last, res_next, res_size, res_max;
    1619             :   GEN  m, ker_m;
    1620             :   long done, rank;
    1621             : 
    1622         279 :   m = stream_read_F2m(frel, h->size_of_FB+1, rel);
    1623         279 :   if (DEBUGLEVEL >= 7)
    1624           0 :     err_printf("\\\\ MATRIX READ BY MPQS\nFREL=%Ps\n",m);
    1625             : 
    1626         279 :   ker_m = F2m_ker_sp(m,0); rank = lg(ker_m)-1;
    1627         279 :   clone_unlock(m);
    1628             : 
    1629         279 :   if (DEBUGLEVEL >= 4)
    1630             :   {
    1631           0 :     if (DEBUGLEVEL >= 7)
    1632             :     {
    1633           0 :       err_printf("\\\\ KERNEL COMPUTED BY MPQS\n");
    1634           0 :       err_printf("KERNEL=%Ps\n",ker_m);
    1635             :     }
    1636           0 :     err_printf("MPQS: Gauss done: kernel has rank %ld, taking gcds...\n", rank);
    1637             :   }
    1638             : 
    1639         279 :   H_rows = rel;
    1640         279 :   H_cols = rank;
    1641             : 
    1642         279 :   if (!H_cols)
    1643             :   { /* trivial kernel. Fail gracefully: main loop may look for more relations */
    1644           0 :     if (DEBUGLEVEL >= 3)
    1645           0 :       pari_warn(warner, "MPQS: no solutions found from linear system solver");
    1646           0 :     return gc_NULL(av); /* no factors found */
    1647             :   }
    1648             : 
    1649             :   /* If the rank is r, we can expect up to 2^r pairwise coprime factors,
    1650             :    * but it may happen that a kernel basis vector contributes nothing new to
    1651             :    * the decomposition.  We allocate room for up to eight factors initially
    1652             :    * (certainly adequate when one or two basis vectors work), adjusting this
    1653             :    * down at the end to what we actually found, or up if we are very lucky and
    1654             :    * find more factors.  In the upper half of our vector, we store information
    1655             :    * about which factors we know to be composite (zero) or believe to be
    1656             :    * composite (NULL) or suspect to be prime (one), or an exponent (two
    1657             :    * or some t_INT) if it is a proper power */
    1658         279 :   ei = cgetg(h->size_of_FB + 2, t_VECSMALL);
    1659         279 :   av2 = avma;
    1660         279 :   if (rank > (long)BITS_IN_LONG - 2)
    1661         133 :     res_max = LONG_MAX; /* the common case, unfortunately */
    1662             :   else
    1663         146 :     res_max = 1L<<rank; /* max number of factors we can hope for */
    1664         279 :   res_size = 8; /* no. of factors we can store in this res */
    1665         279 :   res = cgetg(2*res_size+1, t_VEC);
    1666         279 :   for (i=2*res_size; i; i--) res[i] = 0;
    1667         279 :   res_next = res_last = 1;
    1668             : 
    1669         970 :   for (i = 1; i <= H_cols; i++)
    1670             :   { /* loop over kernel basis */
    1671         970 :     X = Y_prod = gen_1;
    1672         970 :     memset((void *)(ei+1), 0, (h->size_of_FB + 1) * sizeof(long));
    1673             : 
    1674         970 :     av3 = avma;
    1675     1224335 :     for (j = 1; j <= H_rows; j++)
    1676             :     {
    1677     1223365 :       if (F2m_coeff(ker_m, j, i))
    1678       92770 :         Y_prod = mpqs_add_relation(Y_prod, N, ei, gel(frel,j));
    1679     1223365 :       if (gc_needed(av3,1))
    1680             :       {
    1681           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"[1]: mpqs_solve_linear_system");
    1682           0 :         Y_prod = gerepileuptoint(av3, Y_prod);
    1683             :       }
    1684             :     }
    1685         970 :     Y_prod = gerepileuptoint(av3, Y_prod);
    1686             : 
    1687         970 :     av3 = avma;
    1688     1157797 :     for (j = 2; j <= h->size_of_FB + 1; j++)
    1689     1156827 :       if (ei[j])
    1690             :       {
    1691      158794 :         if (ei[j] & 1) pari_err_BUG("MPQS (relation is a nonsquare)");
    1692      317588 :         X = remii(mulii(X,
    1693      317588 :                         Fp_powu(utoipos(FB[j].fbe_p), (ulong)ei[j]>>1, N)),
    1694             :                   N);
    1695      158794 :         if (gc_needed(av3,1))
    1696             :         {
    1697           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"[2]: mpqs_solve_linear_system");
    1698           0 :           X = gerepileupto(av3, X);
    1699             :         }
    1700             :       }
    1701         970 :     X = gerepileuptoint(av3, X);
    1702         970 :     if (MPQS_DEBUGLEVEL >= 1)
    1703             :     {
    1704           0 :       if (!dvdii(subii(sqri(X), sqri(Y_prod)), N))
    1705             :       { /* shouldn't happen */
    1706           0 :         err_printf("MPQS: X^2 - Y^2 != 0 mod N\n");
    1707           0 :         err_printf("\tindex i = %ld\n", i);
    1708           0 :         pari_warn(warner, "MPQS: wrong relation found after Gauss");
    1709             :       }
    1710             :     }
    1711             :     /* At this point, X^2 == Y^2 mod N.  Indeed, something stronger is true:
    1712             :      * We have gcd(X-Y, N) * gcd(X+Y, N) == N.  Why?
    1713             :      * First, N divides X^2 - Y^2, so it divides the lefthand side.
    1714             :      * Second, let P be any prime factor of N.  If P were to divide both
    1715             :      * X-Y and X+Y, then it would divide their sum 2X.  But X (in the present
    1716             :      * backwards notation!) is a product of powers of FB primes, and no FB
    1717             :      * prime is a divisor of N, or we would have found out about it already
    1718             :      * whilst constructing the FB.
    1719             :      * Therefore in the following it is sufficient to work with gcd(X+Y, N)
    1720             :      * alone, and never look at gcd(X-Y, N).
    1721             :      */
    1722         970 :     done = 0; /* (re-)counts probably-prime factors (or powers whose bases we
    1723             :                * don't want to handle any further) */
    1724         970 :     X_plus_Y = addii(X, Y_prod);
    1725         970 :     if (res_next < 3)
    1726             :     { /* we still haven't decomposed the original N, and want both a gcd and
    1727             :        * its cofactor. */
    1728         949 :       D1 = gcdii(X_plus_Y, N);
    1729         949 :       if (is_pm1(D1) || equalii(D1,N)) { set_avma(av3); continue; }
    1730             :       /* got something that works */
    1731         279 :       if (DEBUGLEVEL >= 5)
    1732           0 :         err_printf("MPQS: splitting N after %ld kernel vector%s\n",
    1733             :                    i+1, (i? "s" : ""));
    1734             :       /* GN20050707 Fixed:
    1735             :        * Don't divide N in place. We still need it for future X and Y_prod
    1736             :        * computations! */
    1737         279 :       gel(res,1) = diviiexact(N, D1);
    1738         279 :       gel(res,2) = D1;
    1739         279 :       res_last = res_next = 3;
    1740             : 
    1741         279 :       if ( split(gel(res,1),  &gel(res,res_size+1), &gel(res,1)) ) done++;
    1742         279 :       if ( split(D1, &gel(res,res_size+2), &gel(res,2)) ) done++;
    1743         279 :       if (done == 2) break;     /* both factors look prime or were powers */
    1744             :       /* GN20050707: moved following line down to here, was before the
    1745             :        * two split() invocations.  Very rare case anyway. */
    1746          14 :       if (res_max == 2) break; /* two out of two possible factors seen */
    1747          14 :       if (DEBUGLEVEL >= 5)
    1748           0 :         err_printf("MPQS: got two factors, looking for more...\n");
    1749             :     }
    1750             :     else
    1751             :     { /* we already have factors */
    1752          77 :       for (j=1; j < res_next; j++)
    1753             :       { /* loop over known-composite factors */
    1754          56 :         if (gel(res,res_size+j) && gel(res,res_size+j) != gen_0)
    1755             :         {
    1756          35 :           done++; continue; /* skip probable primes etc */
    1757             :         }
    1758             :         /* actually, also skip square roots of squares etc.  They are a lot
    1759             :          * smaller than the original N, and should be easy to deal with later */
    1760          21 :         av3 = avma;
    1761          21 :         D1 = gcdii(X_plus_Y, gel(res,j));
    1762          21 :         if (is_pm1(D1) || equalii(D1, gel(res,j))) { set_avma(av3); continue; }
    1763             :         /* got one which splits this factor */
    1764          14 :         if (DEBUGLEVEL >= 5)
    1765           0 :           err_printf("MPQS: resplitting a factor after %ld kernel vectors\n",
    1766             :                      i+1);      /* always plural */
    1767             :         /* first make sure there's room for another factor */
    1768          14 :         if (res_next > res_size)
    1769             :         { /* need to reallocate (_very_ rare case) */
    1770           0 :           long i1, size = 2*res_size;
    1771             :           GEN RES;
    1772           0 :           if (size > res_max) size = res_max;
    1773           0 :           RES = cgetg(2*size+1, t_VEC);
    1774           0 :           for (i1=2*size; i1>=res_next; i1--) gel(RES,i1) = NULL;
    1775           0 :           for (i1=1; i1<res_next; i1++)
    1776             :           {
    1777             :             /* GN20050707:
    1778             :              * on-stack contents of RES must be rejuvenated */
    1779           0 :             icopyifstack(gel(res,i1), gel(RES,i1)); /* factors */
    1780           0 :             if ( gel(res,res_size+i1) )
    1781           0 :               icopyifstack(gel(res,res_size+i1), gel(RES,size+i1));
    1782             :             /* primality tags */
    1783             :           }
    1784           0 :           res = RES; res_size = size;   /* res_next unchanged */
    1785             :         }
    1786             :         /* now there is room; divide into existing factor and store the
    1787             :            new gcd */
    1788          14 :         diviiz(gel(res,j), D1, gel(res,j));
    1789          14 :         gel(res,res_next) = D1;
    1790             : 
    1791             :         /* following overwrites the old known-composite indication at j */
    1792          14 :         if (split( gel(res,j), &gel(res,res_size+j), &gel(res,j)) ) done++;
    1793             :         /* GN20050707 Fixed:
    1794             :          * Don't increment done when the newly stored factor seems to be
    1795             :          * prime or otherwise devoid of interest - this happens later
    1796             :          * when we routinely revisit it during the present inner loop. */
    1797          14 :         (void)split(D1, &gel(res,res_size+res_next), &gel(res,res_next));
    1798             : 
    1799             :         /* GN20050707: Following line moved down to here, was before the
    1800             :          * two split() invocations. */
    1801          14 :         if (++res_next > res_max)
    1802             :         {
    1803             :           /* all possible factors seen, outer loop postprocessing will
    1804             :            * proceed to break out of the outer loop below. */
    1805           0 :           break;
    1806             :         }
    1807             :       }       /* loop over known composite factors */
    1808             : 
    1809          21 :       if (res_next > res_last)
    1810             :       {
    1811          14 :         res_last = res_next - 1; /* we might have resplit more than one */
    1812          14 :         if (DEBUGLEVEL >= 5)
    1813           0 :           err_printf("MPQS: got %ld factors%s\n", res_last,
    1814             :                      (done < res_last ? ", looking for more..." : ""));
    1815          14 :         res_last = res_next;
    1816             :       }
    1817             :       /* break out of the outer loop when we have seen res_max factors, and
    1818             :        * also when all current factors are probable primes */
    1819          21 :       if (res_next > res_max || done == res_next - 1) break;
    1820             :     } /* end case of further splitting of existing factors */
    1821          21 :     if (gc_needed(av2,1))
    1822             :     {
    1823             :       long i1;
    1824           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"[3]: mpqs_solve_linear_system");
    1825             :       /* gcopy would have a problem with our NULL pointers... */
    1826           0 :       new_res = cgetg(lg(res), t_VEC);
    1827           0 :       for (i1=2*res_size; i1>=res_next; i1--) new_res[i1] = 0;
    1828           0 :       for (i1=1; i1<res_next; i1++)
    1829             :       {
    1830           0 :         icopyifstack(gel(res,i1), gel(new_res,i1));
    1831             :         /* GN20050707: the marker GENs might need rejuvenating, too */
    1832           0 :         if (gel(res,res_size+i1))
    1833           0 :           icopyifstack(gel(res,res_size+i1), gel(new_res,res_size+i1));
    1834             :       }
    1835           0 :       res = gerepileupto(av2, new_res);
    1836             :     }
    1837             :   } /* for (loop over kernel basis) */
    1838             : 
    1839         279 :   if (res_next < 3) return gc_NULL(av); /* no factors found */
    1840             : 
    1841             :   /* normal case:  convert internal format to ifac format as described in
    1842             :    * src/basemath/ifactor1.c  (triples of components: value, exponent, class
    1843             :    * [unknown, known composite, known prime]),  clean up and return result */
    1844         279 :   res_last = res_next - 1; /* number of distinct factors found */
    1845         279 :   new_res = cgetg(3*res_last + 1, t_VEC);
    1846         279 :   if (DEBUGLEVEL >= 6)
    1847           0 :     err_printf("MPQS: wrapping up vector of %ld factors\n", res_last);
    1848         851 :   for (i=1,j=1; i <= res_last; i++)
    1849             :   {
    1850         572 :     GEN F = gel(res, res_size+i);
    1851         572 :     icopyifstack(gel(res,i), gel(new_res,j++)); /* factor */
    1852        1144 :     gel(new_res,j++) = /* exponent */
    1853         572 :       F ? (F == gen_0 ? gen_1
    1854         572 :                       : (isonstack(F) ? icopy(F) : F))
    1855        1144 :         : gen_1; /* F was NULL */
    1856        1144 :     gel(new_res,j++) = /* class */
    1857         572 :       F == gen_0 ? gen_0 :      /* known composite */
    1858             :         NULL;           /* base of power or suspected prime --
    1859             :                                    mark as `unknown' */
    1860         572 :     if (DEBUGLEVEL >= 6)
    1861           0 :       err_printf("\tpackaging %ld: %Ps ^%ld (%s)\n", i, res[i],
    1862           0 :                  itos(gel(new_res,j-2)), (F == gen_0 ? "comp." : "unknown"));
    1863             :   }
    1864         279 :   return gerepileupto(av, new_res);
    1865             : }
    1866             : 
    1867             : /*********************************************************************/
    1868             : /**                                                                 **/
    1869             : /**               MAIN ENTRY POINT AND DRIVER ROUTINE               **/
    1870             : /**                                                                 **/
    1871             : /*********************************************************************/
    1872             : /* All percentages below are actually fixed-point quantities scaled by 10
    1873             :  * (value of 1 means 0.1%, 1000 means 100%) */
    1874             : 
    1875             : /* Factors N using the self-initializing multipolynomial quadratic sieve
    1876             :  * (SIMPQS).  Returns one of the two factors, or (usually) a vector of factors
    1877             :  * and exponents and information about which ones are still composite, or NULL
    1878             :  * when something goes wrong or when we can't seem to make any headway. */
    1879             : GEN
    1880         279 : mpqs(GEN N)
    1881             : {
    1882             :   mpqs_handle_t H;
    1883             :   GEN fact; /* will in the end hold our factor(s) */
    1884             :   mpqs_int32_t size_of_FB; /* size of the factor base */
    1885             :   mpqs_FB_entry_t *FB; /* factor base */
    1886             : 
    1887             :   /* local loop / auxiliary vars */
    1888             :   ulong p;
    1889             :   /* bookkeeping */
    1890             :   long size_N; /* ~ log_10(N) */
    1891             :   long tc;                      /* # of candidates found in one iteration */
    1892         279 :   long tff = 0;                 /* # recently found full rels from sieving */
    1893             :   long tfc;                     /* # full rels recently combined from LPs */
    1894         279 :   double tfc_ratio = 0;         /* recent (tfc + tff) / tff */
    1895             :   ulong sort_interval;          /* determine when to sort and merge */
    1896         279 :   long percentage = 0;          /* scaled by 10, see comment above */
    1897             :   double net_yield;
    1898         279 :   long total_full_relations = 0, total_partial_relations = 0, total_no_cand = 0;
    1899         279 :   long vain_iterations = 0, good_iterations = 0, iterations = 0;
    1900             : 
    1901             :   pari_timer T;
    1902             :   GEN fnew, vnew;
    1903             :   long nvnew;
    1904             :   hashtable lprel, frel;
    1905         279 :   pari_sp av = avma;
    1906             : 
    1907         279 :   if (DEBUGLEVEL >= 4)
    1908             :   {
    1909           0 :     timer_start(&T);
    1910           0 :     err_printf("MPQS: number to factor N = %Ps\n", N);
    1911             :   }
    1912             : 
    1913         279 :   H.N = N;
    1914         279 :   H.bin_index = 0;
    1915         279 :   H.index_i = 0;
    1916         279 :   H.index_j = 0;
    1917         279 :   H.index2_moved = 0;
    1918             : 
    1919         279 :   size_N = decimal_len(N);
    1920         279 :   if (size_N > MPQS_MAX_DIGIT_SIZE_KN)
    1921             :   {
    1922           0 :     pari_warn(warner, "MPQS: number too big to be factored with MPQS,\n\tgiving up");
    1923           0 :     return NULL;
    1924             :   }
    1925             : 
    1926         279 :   if (DEBUGLEVEL >= 4)
    1927           0 :     err_printf("MPQS: factoring number of %ld decimal digits\n", size_N);
    1928             : 
    1929         279 :   p = mpqs_find_k(&H);
    1930         279 :   if (p) { set_avma(av); return utoipos(p); }
    1931         279 :   if (DEBUGLEVEL >= 5) err_printf("MPQS: found multiplier %ld for N\n",
    1932           0 :                                   H._k->k);
    1933         279 :   H.kN = muliu(N, H._k->k);
    1934             : 
    1935         279 :   if (!mpqs_set_parameters(&H))
    1936             :   {
    1937           0 :     pari_warn(warner,
    1938             :         "MPQS: number too big to be factored with MPQS,\n\tgiving up");
    1939           0 :     return NULL;
    1940             :   }
    1941             : 
    1942         279 :   size_of_FB = H.size_of_FB;
    1943         279 :   sort_interval = H.first_sort_point;
    1944             : 
    1945         279 :   if (DEBUGLEVEL >= 5)
    1946           0 :     err_printf("MPQS: creating factor base and allocating arrays...\n");
    1947         279 :   FB = mpqs_create_FB(&H, &p);
    1948         279 :   if (p) { set_avma(av); return utoipos(p); }
    1949         279 :   mpqs_sieve_array_ctor(&H);
    1950         279 :   mpqs_poly_ctor(&H);
    1951             : 
    1952         279 :   H.lp_bound = minss(H.largest_FB_p, MPQS_LP_BOUND);
    1953             :   /* don't allow large primes to have room for two factors both bigger than
    1954             :    * what the FB contains (...yet!) */
    1955         279 :   H.lp_bound *= minss(H.lp_scale, H.largest_FB_p - 1);
    1956             : 
    1957         279 :   H.dkN = gtodouble(H.kN);
    1958             :   /* compute the threshold and fill in the byte-sized scaled logarithms */
    1959         279 :   mpqs_set_sieve_threshold(&H);
    1960             : 
    1961         279 :   if (!mpqs_locate_A_range(&H)) return NULL;
    1962             : 
    1963         279 :   if (DEBUGLEVEL >= 4)
    1964             :   {
    1965           0 :     err_printf("MPQS: sieving interval = [%ld, %ld]\n", -(long)H.M, (long)H.M);
    1966             :     /* that was a little white lie, we stop one position short at the top */
    1967           0 :     err_printf("MPQS: size of factor base = %ld\n",
    1968             :                (long)size_of_FB);
    1969           0 :     err_printf("MPQS: striving for %ld relations\n",
    1970           0 :                (long)H.target_no_rels);
    1971           0 :     err_printf("MPQS: coefficients A will be built from %ld primes each\n",
    1972           0 :                (long)H.omega_A);
    1973           0 :     err_printf("MPQS: primes for A to be chosen near FB[%ld] = %ld\n",
    1974           0 :                (long)H.index2_FB,
    1975           0 :                (long)FB[H.index2_FB].fbe_p);
    1976           0 :     err_printf("MPQS: smallest prime used for sieving FB[%ld] = %ld\n",
    1977           0 :                (long)H.index1_FB,
    1978           0 :                (long)FB[H.index1_FB].fbe_p);
    1979           0 :     err_printf("MPQS: largest prime in FB = %ld\n",
    1980           0 :                (long)H.largest_FB_p);
    1981           0 :     err_printf("MPQS: bound for `large primes' = %ld\n", (long)H.lp_bound);
    1982             :   }
    1983             : 
    1984         279 :   if (DEBUGLEVEL >= 5)
    1985           0 :     err_printf("MPQS: sieve threshold = %u\n",
    1986           0 :                (unsigned int)H.sieve_threshold);
    1987             : 
    1988         279 :   if (DEBUGLEVEL >= 4)
    1989           0 :     err_printf("MPQS: first sorting at %ld%%, then every %3.1f%% / %3.1f%%\n",
    1990           0 :                sort_interval/10, H.sort_pt_interval/10.,
    1991           0 :                H.sort_pt_interval/20.);
    1992             : 
    1993             :   /* main loop which
    1994             :    * - computes polynomials and their zeros (SI)
    1995             :    * - does the sieving
    1996             :    * - tests candidates of the sieve array */
    1997             : 
    1998             :   /* Let (A, B_i) the current pair of coeffs. If i == 0 a new A is generated */
    1999         279 :   H.index_j = (mpqs_uint32_t)-1;  /* increment below will have it start at 0 */
    2000             : 
    2001         279 :   if (DEBUGLEVEL >= 5) err_printf("MPQS: starting main loop\n");
    2002             : 
    2003         279 :   hash_init_GEN(&frel, H.target_no_rels, gequal, 1);
    2004         279 :   hash_init_ulong(&lprel,H.target_no_rels, 1);
    2005         279 :   vnew = cgetg((long)(sort_interval * (H.target_no_rels/1000.))+2, t_VEC);
    2006         279 :   nvnew = 1;
    2007             :   for(;;)
    2008       46706 :   {
    2009             :     long i, fnb;
    2010       46985 :     iterations++;
    2011             :     /* self initialization: compute polynomial and its zeros */
    2012       46985 :     mpqs_self_init(&H);
    2013       46985 :     if (H.bin_index == 0)
    2014             :     { /* have run out of primes for A */
    2015             :       /* We might change some parameters.  For the moment, simply give up */
    2016           0 :       if (DEBUGLEVEL >= 2)
    2017           0 :         err_printf("MPQS: Ran out of primes for A, giving up.\n");
    2018           0 :       return gc_NULL(av);
    2019             :     }
    2020             : 
    2021       46985 :     memset((void*)(H.sieve_array), 0, (H.M << 1) * sizeof(unsigned char));
    2022       46985 :     mpqs_sieve(&H);
    2023             : 
    2024       46985 :     tc = mpqs_eval_sieve(&H);
    2025       46985 :     total_no_cand += tc;
    2026       46985 :     if (DEBUGLEVEL >= 6)
    2027           0 :       err_printf("MPQS: found %lu candidate%s\n", tc, (tc==1? "" : "s"));
    2028             : 
    2029       46985 :     if (tc)
    2030             :     {
    2031             :       GEN lpnew;
    2032       45725 :       long t = mpqs_eval_cand(&H, tc, &fnew, &lpnew);
    2033       45725 :       total_full_relations += t;
    2034       45725 :       tff += t;
    2035       45725 :       good_iterations++;
    2036       45725 :       if (fnew) vec_frel_add(&frel, fnew);
    2037       45725 :       if (lpnew) vnew = vec_extend(vnew, lpnew, nvnew++);
    2038             :     }
    2039             : 
    2040       46985 :     percentage =
    2041       46985 :       (long)((1000.0 * total_full_relations) / H.target_no_rels);
    2042             : 
    2043       46985 :     if ((ulong)percentage < sort_interval) continue;
    2044             :     /* most main loops continue here! */
    2045             : 
    2046             :     /* Extra processing when we have completed a sort interval: */
    2047        2058 :     if (DEBUGLEVEL >= 3)
    2048             :     {
    2049           0 :       if (DEBUGLEVEL >= 4)
    2050           0 :         err_printf("\nMPQS: passing the %3.1f%% sort point, time = %ld ms\n",
    2051             :                    sort_interval/10., timer_delay(&T));
    2052             :       else
    2053           0 :         err_printf("\nMPQS: passing the %3.1f%% sort point\n",
    2054             :                    sort_interval/10.);
    2055             :     }
    2056             : 
    2057             :     /* combine whatever there is to be combined */
    2058        2058 :     tfc = 0;
    2059             :     /* build full relations out of large prime relations */
    2060        2058 :     fnb = frel.nb;
    2061       37735 :     for (i=1; i<nvnew; i++)
    2062             :     {
    2063       35677 :       GEN fact = mpqs_combine_large_primes(&H, &lprel, gel(vnew,i) , &frel);
    2064       35677 :       if (fact)
    2065             :       { /* factor found during combining */
    2066           0 :         if (DEBUGLEVEL >= 4)
    2067             :         {
    2068           0 :           err_printf("\nMPQS: split N whilst combining, time = %ld ms\n",
    2069             :               timer_delay(&T));
    2070           0 :           err_printf("MPQS: found factor = %Ps\n", fact);
    2071             :         }
    2072           0 :         return gerepileupto(av, fact);
    2073             :       }
    2074             :     }
    2075        2058 :     nvnew = 1;
    2076        2058 :     tfc = frel.nb - fnb;
    2077        2058 :     if (DEBUGLEVEL >= 4)
    2078           0 :       err_printf("MPQS: combined %ld full relation%s\n", tfc, (tfc!=1 ? "s" : ""));
    2079        2058 :     total_partial_relations += tfc;
    2080             : 
    2081             :     /* sort FNEW and merge it into frel */
    2082        2058 :     total_full_relations = frel.nb;
    2083             : 
    2084             :     /* Due to the removal of duplicates, percentage may decrease at
    2085             :      * Nothing to worry about: we _are_ making progress. */
    2086        2058 :     percentage = (long)((1000.0 * total_full_relations) / H.target_no_rels);
    2087        2058 :     net_yield = (total_full_relations * 100.) / (total_no_cand? total_no_cand: 1);
    2088        2058 :     vain_iterations = (long)((1000.0 * (iterations - good_iterations)) / iterations);
    2089             : 
    2090             :     /* Now estimate the current full relations yield rate:  we directly see
    2091             :      * each time through the main loop how many full relations we're getting
    2092             :      * as such from the sieve  (tff since the previous checkpoint),  but
    2093             :      * only at checkpoints do we see how many we're typically combining
    2094             :      * (tfc).  So we're really producing (tfc+tff)/tff as many full rels,
    2095             :      * and when we get close to 100%, we should bias the next interval by
    2096             :      * the inverse ratio.
    2097             :      * Avoid drawing conclusions from too-small samples during very short
    2098             :      * follow-on intervals  (in this case we'll just re-use an earlier
    2099             :      * estimated ratio). */
    2100        2058 :     if (tfc >= 16 && tff >= 20)
    2101          91 :       tfc_ratio = (tfc + tff + 0.) / tff; /* floating-point division */
    2102        2058 :     tff = 0; /* reset this count (tfc is always fresh) */
    2103             : 
    2104        2058 :     if (percentage >= 1000) /* when Gauss had failed */
    2105         279 :       sort_interval = percentage + 2;
    2106        1779 :     else if (percentage >= 820)
    2107             :     {
    2108         593 :       if (tfc_ratio > 1.)
    2109             :       {
    2110          54 :         if (percentage + (H.sort_pt_interval>> 1) * tfc_ratio > 994)
    2111          40 :           sort_interval = (ulong)(percentage + 2 +
    2112          20 :             (1000 - percentage) / tfc_ratio); /* aim for a slight overshoot */
    2113          34 :         else if (percentage >= 980)
    2114           0 :           sort_interval = percentage + 8;
    2115             :         else
    2116          34 :           sort_interval = percentage + (H.sort_pt_interval >> 1);
    2117             :       }
    2118             :       else
    2119             :       {
    2120         539 :         if (percentage >= 980)
    2121          49 :           sort_interval = percentage + 10;
    2122             :         else
    2123         490 :           sort_interval = percentage + (H.sort_pt_interval >> 1);
    2124         539 :         if (sort_interval >= 1000 && percentage < 1000)
    2125         161 :           sort_interval = 1000;
    2126             :       }
    2127             :     }
    2128             :     else
    2129        1186 :       sort_interval = percentage + H.sort_pt_interval;
    2130             : 
    2131        2058 :     if (DEBUGLEVEL >= 4)
    2132             :     {
    2133           0 :       err_printf("MPQS: done sorting%s, time = %ld ms\n",
    2134             :                  " and combining", timer_delay(&T));
    2135           0 :       err_printf("MPQS: found %3.1f%% of the required relations\n",
    2136             :                  percentage/10.);
    2137           0 :       if (DEBUGLEVEL >= 5)
    2138             :       { /* total_full_relations are always plural */
    2139             :         /* GN20050708: present code doesn't succeed in discarding all
    2140             :          * dups, so don't lie about it... */
    2141           0 :         err_printf("MPQS: found %ld full relations\n",
    2142             :                    total_full_relations);
    2143           0 :         if (H.lp_scale > 1)
    2144           0 :           err_printf("MPQS:   (%ld of these from partial relations)\n",
    2145             :                      total_partial_relations);
    2146           0 :         err_printf("MPQS: Net yield: %4.3g full relations per 100 candidates\n",
    2147             :                    net_yield);
    2148           0 :         err_printf("MPQS:            %4.3g full relations per 100 polynomials\n",
    2149           0 :                    (total_full_relations * 100.) / iterations);
    2150           0 :         err_printf("MPQS: %4.1f%% of the polynomials yielded no candidates\n",
    2151             :                    vain_iterations/10.);
    2152           0 :         err_printf("MPQS: next sort point at %3.1f%%\n", sort_interval/10.);
    2153             :       }
    2154             :     }
    2155        2058 :     if (percentage < 1000) continue; /* main loop */
    2156             : 
    2157             :     /* percentage >= 1000, which implies total_full_relations > size_of_FB:
    2158             :        try finishing it off */
    2159             : 
    2160         279 :     if (DEBUGLEVEL >= 4)
    2161           0 :       err_printf("\nMPQS: starting Gauss over F_2 on %ld distinct relations\n",
    2162             :                  total_full_relations);
    2163         279 :     fact = mpqs_solve_linear_system(&H, hash_keys(&frel), total_full_relations);
    2164         279 :     if (fact)
    2165             :     { /* solution found */
    2166         279 :       if (DEBUGLEVEL >= 4)
    2167             :       {
    2168           0 :         err_printf("\nMPQS: time in Gauss and gcds = %ld ms\n", timer_delay(&T));
    2169           0 :         if (typ(fact) == t_INT) err_printf("MPQS: found factor = %Ps\n", fact);
    2170             :         else
    2171             :         {
    2172           0 :           long j, nf = (lg(fact)-1)/3;
    2173           0 :           if (nf == 2)
    2174             :             /* GN20050707: Changed the arrangement of the two factors,
    2175             :              * to match the debug diagnostics in mpqs_solve_linear_system()
    2176             :              * above */
    2177           0 :             err_printf("MPQS: found factors = %Ps\n\tand %Ps\n",
    2178           0 :                         fact[1], fact[4]);
    2179             :           else
    2180             :           {
    2181             :             /* GN20050707: Changed loop to scan upwards instead of downwards,
    2182             :              * to match the debug diagnostics in mpqs_solve_linear_system()
    2183             :              * above */
    2184           0 :             err_printf("MPQS: found %ld factors =\n", nf);
    2185           0 :             for (j=1; j<=nf; j++)
    2186           0 :               err_printf("\t%Ps%s\n", fact[3*j-2], (j<nf ? "," : ""));
    2187             :           }
    2188             :         }
    2189             :       }
    2190             :       /* fact not safe for a gerepilecopy(): segfaults on one of the NULL
    2191             :        * markers. However, it is a nice connected object, and it resides
    2192             :        * already the top of the stack, so... --GN */
    2193         279 :       return gerepileupto(av, fact);
    2194             :     }
    2195             :     else
    2196             :     {
    2197           0 :       if (DEBUGLEVEL >= 4)
    2198             :       {
    2199           0 :         err_printf("\nMPQS: time in Gauss and gcds = %ld ms\n",timer_delay(&T));
    2200           0 :         err_printf("MPQS: no factors found.\n");
    2201           0 :         if (percentage <= MPQS_ADMIT_DEFEAT)
    2202           0 :           err_printf("\nMPQS: restarting sieving ...\n");
    2203             :         else
    2204           0 :           err_printf("\nMPQS: giving up.\n");
    2205             :       }
    2206           0 :       if (percentage > MPQS_ADMIT_DEFEAT)
    2207           0 :         return gc_NULL(av);
    2208             :     }
    2209             :   } /* main loop */
    2210             : }

Generated by: LCOV version 1.13