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 24988-2584e74448) Lines: 519 629 82.5 %
Date: 2020-01-26 05:57:03 Functions: 28 29 96.6 %
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             : /* Self-Initializing Multi-Polynomial Quadratic Sieve, based on code developed
      15             :  * as part of the LiDIA project.
      16             :  *
      17             :  * Original version: Thomas Papanikolaou and Xavier Roblot
      18             :  * Extensively modified by The PARI group. */
      19             : /* Notation commonly used in this file, and sketch of algorithm:
      20             :  *
      21             :  * Given an odd integer N > 1 to be factored, we throw in a small odd squarefree
      22             :  * multiplier k so as to make kN = 1 mod 4 and to have many small primes over
      23             :  * which X^2 - kN splits.  We compute a factor base FB of such primes then
      24             :  * look for values x0 such that Q0(x0) = x0^2 - kN can be decomposed over FB,
      25             :  * up to a possible factor dividing k and a possible "large prime". Relations
      26             :  * involving the latter can be combined into full relations which don't; full
      27             :  * relations, by Gaussian elimination over F2 for the exponent vectors lead us
      28             :  * to an expression X^2 - Y^2 divisible by N and hopefully to a nontrivial
      29             :  * splitting when we compute gcd(X + Y, N).  Note that this can never
      30             :  * split prime powers.
      31             :  *
      32             :  * Candidates x0 are found by sieving along arithmetic progressions modulo the
      33             :  * small primes in FB and evaluation of candidates picks out those x0 where
      34             :  * many of these progressions coincide, resulting in a highly divisible Q0(x0).
      35             :  *
      36             :  * The Multi-Polynomial version improves this by choosing a modest subset of
      37             :  * FB primes (let A be their product) and forcing these to divide Q0(x).
      38             :  * Write Q(x) = Q0(2Ax + B) = (2Ax + B)^2 - kN = 4A(Ax^2 + Bx + C), where B is
      39             :  * suitably chosen.  For each A, there are 2^omega_A possible values for B
      40             :  * but we'll use only half of these, since the other half is easily covered by
      41             :  * exploiting the symmetry x -> -x of the original Q0. The "Self-Initializating"
      42             :  * bit refers to the fact that switching from one B to the next is fast, whereas
      43             :  * switching to the next A involves some recomputation (C is never needed).
      44             :  * Thus we quickly run through many polynomials sharing the same A.
      45             :  *
      46             :  * The sieve ranges over values x0 such that |x0| < M  (we use x = x0 + M
      47             :  * as array subscript).  The coefficients A are chosen so that A*M ~ sqrt(kN).
      48             :  * Then |B| is bounded by ~ (j+4)*A, and |C| = -C ~ (M/4)*sqrt(kN), so
      49             :  * Q(x0)/(4A) takes values roughly between -|C| and 3|C|.
      50             :  *
      51             :  * Refinements. We do not use the smallest FB primes for sieving, incorporating
      52             :  * them only after selecting candidates).  The substition of 2Ax+B into
      53             :  * X^2 - kN, with odd B, forces 2 to occur; when kN is 1 mod 8, it occurs at
      54             :  * least to the 3rd power; when kN = 5 mod 8, it occurs exactly to the 2nd
      55             :  * power.  We never sieve on 2 and always pull out the power of 2 directly. The
      56             :  * prime factors of k show up whenever 2Ax + B has a factor in common with k;
      57             :  * we don't sieve on these either but easily recognize them in a candidate. */
      58             : #include "pari.h"
      59             : #include "paripriv.h"
      60             : 
      61             : /** DEBUG **/
      62             : /* #define MPQS_DEBUG_VERBOSE 1 */
      63             : #include "mpqs.h"
      64             : 
      65             : #define REL_OFFSET 20
      66             : #define REL_MASK ((1UL<<REL_OFFSET)-1)
      67             : #define MAX_PE_PAIR 60
      68             : 
      69       17779 : static GEN rel_q(GEN c) { return gel(c,3); }
      70       35558 : static GEN rel_Y(GEN c) { return gel(c,1); }
      71       35558 : static GEN rel_p(GEN c) { return gel(c,2); }
      72             : 
      73             : static void
      74      144053 : frel_add(hashtable *frel, GEN R)
      75             : {
      76      144053 :   ulong h = hash_GEN(R);
      77      144053 :   if (!hash_search2(frel, (void*)R, h))
      78      143409 :     hash_insert2(frel, (void*)R, (void*)1, h);
      79      144053 : }
      80             : 
      81             : /*********************************************************************/
      82             : /**                         INITIAL SIZING                          **/
      83             : /*********************************************************************/
      84             : /* # of decimal digits of argument */
      85             : static long
      86         677 : decimal_len(GEN N)
      87         677 : { pari_sp av = avma; return gc_long(av, 1+logint(N, utoipos(10))); }
      88             : 
      89             : /* To be called after choosing k and putting kN into the handle:
      90             :  * Pick up the parameters for given size of kN in decimal digits and fill in
      91             :  * the handle. Return 0 when kN is too large, 1 when we're ok. */
      92             : static int
      93         335 : mpqs_set_parameters(mpqs_handle_t *h)
      94             : {
      95             :   long s, D;
      96             :   const mpqs_parameterset_t *P;
      97             : 
      98         335 :   h->digit_size_kN = D = decimal_len(h->kN);
      99         335 :   if (D > MPQS_MAX_DIGIT_SIZE_KN) return 0;
     100         335 :   P = &(mpqs_parameters[maxss(0, D - 9)]);
     101         335 :   h->tolerance   = P->tolerance;
     102         335 :   h->lp_scale    = P->lp_scale;
     103             :   /* make room for prime factors of k if any: */
     104         335 :   h->size_of_FB  = s = P->size_of_FB + h->_k->omega_k;
     105             :   /* for the purpose of Gauss elimination etc., prime factors of k behave
     106             :    * like real FB primes, so take them into account when setting the goal: */
     107         335 :   h->target_rels = (s >= 200 ? s + 10 : (mpqs_int32_t)(s * 1.05));
     108         335 :   h->M           = P->M;
     109         335 :   h->omega_A     = P->omega_A;
     110         335 :   h->no_B        = 1UL << (P->omega_A - 1);
     111         335 :   h->pmin_index1 = P->pmin_index1;
     112             :   /* certain subscripts into h->FB should also be offset by omega_k: */
     113         335 :   h->index0_FB   = 3 + h->_k->omega_k;
     114         335 :   if (DEBUGLEVEL >= 5)
     115             :   {
     116           0 :     err_printf("MPQS: kN = %Ps\n", h->kN);
     117           0 :     err_printf("MPQS: kN has %ld decimal digits\n", D);
     118           0 :     err_printf("\t(estimated memory needed: %4.1fMBy)\n",
     119           0 :                (s + 1)/8388608. * h->target_rels);
     120             :   }
     121         335 :   return 1;
     122             : }
     123             : 
     124             : /*********************************************************************/
     125             : /**                       OBJECT HOUSEKEEPING                       **/
     126             : /*********************************************************************/
     127             : 
     128             : /* factor base constructor. Really a home-grown memalign(3c) underneath.
     129             :  * We don't want FB entries to straddle L1 cache line boundaries, and
     130             :  * malloc(3c) only guarantees alignment adequate for all primitive data
     131             :  * types of the platform ABI - typically to 8 or 16 byte boundaries.
     132             :  * Also allocate the inv_A_H array.
     133             :  * The FB array pointer is returned for convenience */
     134             : static mpqs_FB_entry_t *
     135         335 : mpqs_FB_ctor(mpqs_handle_t *h)
     136             : {
     137             :   /* leave room for slots 0, 1, and sentinel slot at the end of the array */
     138         335 :   long size_FB_chunk = (h->size_of_FB + 3) * sizeof(mpqs_FB_entry_t);
     139             :   /* like FB, except this one does not have a sentinel slot at the end */
     140         335 :   long size_IAH_chunk = (h->size_of_FB + 2) * sizeof(mpqs_inv_A_H_t);
     141         335 :   char *fbp = (char*)stack_malloc(size_FB_chunk + 64);
     142         335 :   char *iahp = (char*)stack_malloc(size_IAH_chunk + 64);
     143             :   long fbl, iahl;
     144             : 
     145         335 :   h->FB_chunk = (void *)fbp;
     146         335 :   h->invAH_chunk = (void *)iahp;
     147             :   /* round up to next higher 64-bytes-aligned address */
     148         335 :   fbl = (((long)fbp) + 64) & ~0x3FL;
     149             :   /* and put the actual array there */
     150         335 :   h->FB = (mpqs_FB_entry_t *)fbl;
     151             : 
     152         335 :   iahl = (((long)iahp) + 64) & ~0x3FL;
     153         335 :   h->inv_A_H = (mpqs_inv_A_H_t *)iahl;
     154         335 :   return (mpqs_FB_entry_t *)fbl;
     155             : }
     156             : 
     157             : /* sieve array constructor;  also allocates the candidates array
     158             :  * and temporary storage for relations under construction */
     159             : static void
     160         335 : mpqs_sieve_array_ctor(mpqs_handle_t *h)
     161             : {
     162         335 :   long size = (h->M << 1) + 1;
     163         335 :   mpqs_int32_t size_of_FB = h->size_of_FB;
     164             : 
     165         335 :   h->sieve_array = (unsigned char *) stack_malloc(size * sizeof(unsigned char));
     166         335 :   h->sieve_array_end = h->sieve_array + size - 2;
     167         335 :   h->sieve_array_end[1] = 255; /* sentinel */
     168         335 :   h->candidates = (long *)stack_malloc(MPQS_CANDIDATE_ARRAY_SIZE * sizeof(long));
     169             :   /* whereas mpqs_self_init() uses size_of_FB+1, we just use the size as
     170             :    * it is, not counting FB[1], to start off the following estimate */
     171         335 :   if (size_of_FB > MAX_PE_PAIR) size_of_FB = MAX_PE_PAIR;
     172             :   /* and for tracking which primes occur in the current relation: */
     173         335 :   h->relaprimes = (long *) stack_malloc((size_of_FB << 1) * sizeof(long));
     174         335 : }
     175             : 
     176             : /* allocate GENs for current polynomial and self-initialization scratch data */
     177             : static void
     178         335 : mpqs_poly_ctor(mpqs_handle_t *h)
     179             : {
     180         335 :   mpqs_int32_t i, w = h->omega_A;
     181         335 :   h->per_A_pr = (mpqs_per_A_prime_t *)
     182         335 :                 stack_calloc(w * sizeof(mpqs_per_A_prime_t));
     183             :   /* A is the product of w primes, each below word size.
     184             :    * |B| <= (w + 4) * A, so can have at most one word more
     185             :    * H holds residues modulo A: the same size as used for A is sufficient. */
     186         335 :   h->A = cgeti(w + 2);
     187         335 :   h->B = cgeti(w + 3);
     188         335 :   for (i = 0; i < w; i++) h->per_A_pr[i]._H = cgeti(w + 2);
     189         335 : }
     190             : 
     191             : /*********************************************************************/
     192             : /**                        FACTOR BASE SETUP                        **/
     193             : /*********************************************************************/
     194             : /* fill in the best-guess multiplier k for N. We force kN = 1 mod 4.
     195             :  * Caller should proceed to fill in kN */
     196             : static ulong
     197         335 : mpqs_find_k(mpqs_handle_t *h)
     198             : {
     199         335 :   const pari_sp av = avma;
     200         335 :   const long N_mod_8 = mod8(h->N), N_mod_4 = N_mod_8 & 3;
     201             :   forprime_t S;
     202             :   struct {
     203             :     const mpqs_multiplier_t *_k;
     204             :     long np; /* number of primes in factorbase so far for this k */
     205             :     double value; /* the larger, the better */
     206             :   } cache[MPQS_POSSIBLE_MULTIPLIERS];
     207             :   ulong p, i, nbk;
     208             : 
     209        3245 :   for (i = nbk = 0; i < numberof(cand_multipliers); i++)
     210             :   {
     211        3245 :     const mpqs_multiplier_t *cand_k = &cand_multipliers[i];
     212        3245 :     long k = cand_k->k;
     213             :     double v;
     214        3245 :     if ((k & 3) != N_mod_4) continue; /* want kN = 1 (mod 4) */
     215        1675 :     v = -0.35 * log2((double)k);
     216        1675 :     if ((k & 7) == N_mod_8) v += M_LN2; /* kN = 1 (mod 8) */
     217        1675 :     cache[nbk].np = 0;
     218        1675 :     cache[nbk]._k = cand_k;
     219        1675 :     cache[nbk].value = v;
     220        1675 :     if (++nbk == MPQS_POSSIBLE_MULTIPLIERS) break; /* enough */
     221             :   }
     222             :   /* next test is an impossible situation: kills spurious gcc-5.1 warnings
     223             :    * "array subscript is above array bounds" */
     224         335 :   if (nbk > MPQS_POSSIBLE_MULTIPLIERS) nbk = MPQS_POSSIBLE_MULTIPLIERS;
     225         335 :   u_forprime_init(&S, 2, ULONG_MAX);
     226         335 :   while ( (p = u_forprime_next(&S)) )
     227             :   {
     228        6309 :     long kroNp = kroiu(h->N, p), seen = 0;
     229        6309 :     if (!kroNp) return p;
     230       37854 :     for (i = 0; i < nbk; i++)
     231             :     {
     232       31545 :       if (cache[i].np > MPQS_MULTIPLIER_SEARCH_DEPTH) continue;
     233       21979 :       seen++;
     234       21979 :       if (krouu(cache[i]._k->k % p, p) == kroNp) /* kronecker(k*N, p)=1 */
     235             :       {
     236       10050 :         cache[i].value += log2((double) p)/p;
     237       10050 :         cache[i].np++;
     238             :       }
     239             :     }
     240        6309 :     if (!seen) break; /* we're gone through SEARCH_DEPTH primes for all k */
     241             :   }
     242         335 :   if (!p) pari_err_OVERFLOW("mpqs_find_k [ran out of primes]");
     243             :   {
     244         335 :     long best_i = 0;
     245         335 :     double v = cache[0].value;
     246        1675 :     for (i = 1; i < nbk; i++)
     247        1340 :       if (cache[i].value > v) { best_i = i; v = cache[i].value; }
     248         335 :     h->_k = cache[best_i]._k; return gc_ulong(av,0);
     249             :   }
     250             : }
     251             : 
     252             : /* Create a factor base of 'size' primes p_i such that legendre(k*N, p_i) != -1
     253             :  * We could have shifted subscripts down from their historical arrangement,
     254             :  * but this seems too risky for the tiny potential gain in memory economy.
     255             :  * The real constraint is that the subscripts of anything which later shows
     256             :  * up at the Gauss stage must be nonnegative, because the exponent vectors
     257             :  * there use the same subscripts to refer to the same FB entries.  Thus in
     258             :  * particular, the entry representing -1 could be put into FB[0], but could
     259             :  * not be moved to FB[-1] (although mpqs_FB_ctor() could be easily adapted
     260             :  * to support negative subscripts).-- The historically grown layout is:
     261             :  * FB[0] is unused.
     262             :  * FB[1] is not explicitly used but stands for -1.
     263             :  * FB[2] contains 2 (always).
     264             :  * Before we are called, the size_of_FB field in the handle will already have
     265             :  * been adjusted by _k->omega_k, so there's room for the primes dividing k,
     266             :  * which when present will occupy FB[3] and following.
     267             :  * The "real" odd FB primes begin at FB[h->index0_FB].
     268             :  * FB[size_of_FB+1] is the last prime p_i.
     269             :  * FB[size_of_FB+2] is a sentinel to simplify some of our loops.
     270             :  * Thus we allocate size_of_FB+3 slots for FB.
     271             :  *
     272             :  * If a prime factor of N is found during the construction, it is returned
     273             :  * in f, otherwise f = 0. */
     274             : 
     275             : /* returns the FB array pointer for convenience */
     276             : static mpqs_FB_entry_t *
     277         335 : mpqs_create_FB(mpqs_handle_t *h, ulong *f)
     278             : {
     279         335 :   mpqs_FB_entry_t *FB = mpqs_FB_ctor(h);
     280         335 :   const pari_sp av = avma;
     281         335 :   mpqs_int32_t size = h->size_of_FB;
     282             :   long i;
     283         335 :   mpqs_uint32_t k = h->_k->k;
     284             :   forprime_t S;
     285             : 
     286         335 :   FB[2].fbe_p = 2;
     287             :   /* the fbe_logval and the fbe_sqrt_kN for 2 are never used */
     288         335 :   FB[2].fbe_flags = MPQS_FBE_CLEAR;
     289         489 :   for (i = 3; i < h->index0_FB; i++)
     290             :   { /* this loop executes h->_k->omega_k = 0, 1, or 2 times */
     291         154 :     mpqs_uint32_t kp = (ulong)h->_k->kp[i-3];
     292         154 :     if (MPQS_DEBUGLEVEL >= 7) err_printf(",<%lu>", (ulong)kp);
     293         154 :     FB[i].fbe_p = kp;
     294             :     /* we could flag divisors of k here, but no need so far */
     295         154 :     FB[i].fbe_flags = MPQS_FBE_CLEAR;
     296         154 :     FB[i].fbe_flogp = (float)log2((double) kp);
     297         154 :     FB[i].fbe_sqrt_kN = 0;
     298             :   }
     299         335 :   (void)u_forprime_init(&S, 3, ULONG_MAX);
     300      280961 :   while (i < size + 2)
     301             :   {
     302      280291 :     ulong p = u_forprime_next(&S);
     303      280291 :     if (p > k || k % p)
     304             :     {
     305      280137 :       ulong kNp = umodiu(h->kN, p);
     306      280137 :       long kr = krouu(kNp, p);
     307      280137 :       if (kr != -1)
     308             :       {
     309      139430 :         if (kr == 0) { *f = p; return FB; }
     310      139430 :         FB[i].fbe_p = (mpqs_uint32_t) p;
     311      139430 :         FB[i].fbe_flags = MPQS_FBE_CLEAR;
     312             :         /* dyadic logarithm of p; single precision suffices */
     313      139430 :         FB[i].fbe_flogp = (float)log2((double)p);
     314             :         /* cannot yet fill in fbe_logval because the scaling multiplier
     315             :          * depends on the largest prime in FB, as yet unknown */
     316             : 
     317             :         /* x such that x^2 = kN (mod p_i) */
     318      139430 :         FB[i++].fbe_sqrt_kN = (mpqs_uint32_t)Fl_sqrt(kNp, p);
     319             :       }
     320             :     }
     321             :   }
     322         335 :   set_avma(av);
     323         335 :   if (MPQS_DEBUGLEVEL >= 7)
     324             :   {
     325           0 :     err_printf("MPQS: FB [-1,2");
     326           0 :     for (i = 3; i < h->index0_FB; i++) err_printf(",<%lu>", FB[i].fbe_p);
     327           0 :     for (; i < size + 2; i++) err_printf(",%lu", FB[i].fbe_p);
     328           0 :     err_printf("]\n");
     329             :   }
     330             : 
     331         335 :   FB[i].fbe_p = 0;              /* sentinel */
     332         335 :   h->largest_FB_p = FB[i-1].fbe_p; /* at subscript size_of_FB + 1 */
     333             : 
     334             :   /* locate the smallest prime that will be used for sieving */
     335         913 :   for (i = h->index0_FB; FB[i].fbe_p != 0; i++)
     336         913 :     if (FB[i].fbe_p >= h->pmin_index1) break;
     337         335 :   h->index1_FB = i;
     338             :   /* with our parameters this will never fall off the end of the FB */
     339         335 :   *f = 0; return FB;
     340             : }
     341             : 
     342             : /*********************************************************************/
     343             : /**                      MISC HELPER FUNCTIONS                      **/
     344             : /*********************************************************************/
     345             : 
     346             : /* Effect of the following:  multiplying the base-2 logarithm of some
     347             :  * quantity by log_multiplier will rescale something of size
     348             :  *    log2 ( sqrt(kN) * M / (largest_FB_prime)^tolerance )
     349             :  * to 232.  Note that sqrt(kN) * M is just A*M^2, the value our polynomials
     350             :  * take at the outer edges of the sieve interval.  The scale here leaves
     351             :  * a little wiggle room for accumulated rounding errors from the approximate
     352             :  * byte-sized scaled logarithms for the factor base primes which we add up
     353             :  * in the sieving phase.-- The threshold is then chosen so that a point in
     354             :  * the sieve has to reach a result which, under the same scaling, represents
     355             :  *    log2 ( sqrt(kN) * M / (largest_FB_prime)^tolerance )
     356             :  * in order to be accepted as a candidate. */
     357             : /* The old formula was...
     358             :  *   log_multiplier =
     359             :  *      127.0 / (0.5 * log2 (handle->dkN) + log2((double)M)
     360             :  *               - tolerance * log2((double)handle->largest_FB_p));
     361             :  * and we used to use this with a constant threshold of 128. */
     362             : 
     363             : /* NOTE: We used to divide log_multiplier by an extra factor 2, and in
     364             :  * compensation we were multiplying by 2 when the fbe_logp fields were being
     365             :  * filled in, making all those bytes even.  Tradeoff: the extra bit of
     366             :  * precision is helpful, but interferes with a possible sieving optimization
     367             :  * (artifically shift right the logp's of primes in A, and just run over both
     368             :  * arithmetical progressions  (which coincide in this case)  instead of
     369             :  * skipping the second one, to avoid the conditional branch in the
     370             :  * mpqs_sieve() loops).  We could still do this, but might lose a little bit
     371             :  * accuracy for those primes.  Probably no big deal. */
     372             : static void
     373         335 : mpqs_set_sieve_threshold(mpqs_handle_t *h)
     374             : {
     375         335 :   mpqs_FB_entry_t *FB = h->FB;
     376             :   double log_maxval, log_multiplier;
     377             :   long i;
     378             : 
     379         335 :   h->l2sqrtkN = 0.5 * log2(h->dkN);
     380         335 :   h->l2M = log2((double)h->M);
     381         335 :   log_maxval = h->l2sqrtkN + h->l2M - MPQS_A_FUDGE;
     382         335 :   log_multiplier = 232.0 / log_maxval;
     383         670 :   h->sieve_threshold = (unsigned char) (log_multiplier *
     384         670 :     (log_maxval - h->tolerance * log2((double)h->largest_FB_p))) + 1;
     385             :   /* That "+ 1" really helps - we may want to tune towards somewhat smaller
     386             :    * tolerances  (or introduce self-tuning one day)... */
     387             : 
     388             :   /* If this turns out to be <128, scream loudly.
     389             :    * That means that the FB or the tolerance or both are way too
     390             :    * large for the size of kN.  (Normally, the threshold should end
     391             :    * up in the 150...170 range.) */
     392         335 :   if (h->sieve_threshold < 128) {
     393           0 :     h->sieve_threshold = 128;
     394           0 :     pari_warn(warner,
     395             :         "MPQS: sizing out of tune, FB size or tolerance\n\ttoo large");
     396             :   }
     397             : 
     398             :   /* Now fill in the byte-sized approximate scaled logarithms of p_i */
     399         335 :   if (DEBUGLEVEL >= 5)
     400           0 :     err_printf("MPQS: computing logarithm approximations for p_i in FB\n");
     401      139765 :   for (i = h->index0_FB; i < h->size_of_FB + 2; i++)
     402      139430 :     FB[i].fbe_logval = (unsigned char) (log_multiplier * FB[i].fbe_flogp);
     403         335 : }
     404             : 
     405             : /* Given the partially populated handle, find the optimum place in the FB
     406             :  * to pick prime factors for A from.  The lowest admissible subscript is
     407             :  * index0_FB, but unless kN is very small, we stay away a bit from that.
     408             :  * The highest admissible is size_of_FB + 1, where the largest FB prime
     409             :  * resides.  The ideal corner is about (sqrt(kN)/M) ^ (1/omega_A),
     410             :  * so that A will end up of size comparable to sqrt(kN)/M;  experimentally
     411             :  * it seems desirable to stay slightly below this.  Moreover, the selection
     412             :  * of the individual primes happens to err on the large side, for which we
     413             :  * compensate a bit, using the (small positive) quantity MPQS_A_FUDGE.
     414             :  * We rely on a few auxiliary fields in the handle to be already set by
     415             :  * mqps_set_sieve_threshold() before we are called.
     416             :  * Return 1 on success, and 0 otherwise. */
     417             : static int
     418         335 : mpqs_locate_A_range(mpqs_handle_t *h)
     419             : {
     420             :   /* i will be counted up to the desirable index2_FB + 1, and omega_A is never
     421             :    * less than 3, and we want
     422             :    *   index2_FB - (omega_A - 1) + 1 >= index0_FB + omega_A - 3,
     423             :    * so: */
     424         335 :   long i = h->index0_FB + 2*(h->omega_A) - 4;
     425             :   double l2_target_pA;
     426         335 :   mpqs_FB_entry_t *FB = h->FB;
     427             : 
     428         335 :   h->l2_target_A = (h->l2sqrtkN - h->l2M - MPQS_A_FUDGE);
     429         335 :   l2_target_pA = h->l2_target_A / h->omega_A;
     430             : 
     431             :   /* find the sweet spot, normally shouldn't take long */
     432         335 :   while (FB[i].fbe_p && FB[i].fbe_flogp <= l2_target_pA) i++;
     433             : 
     434             :   /* check whether this hasn't walked off the top end... */
     435             :   /* The following should actually NEVER happen. */
     436         335 :   if (i > h->size_of_FB - 3)
     437             :   { /* this isn't going to work at all. */
     438           0 :     pari_warn(warner,
     439             :         "MPQS: sizing out of tune, FB too small or\n\tway too few primes in A");
     440           0 :     return 0;
     441             :   }
     442         335 :   h->index2_FB = i - 1; return 1;
     443             :   /* assert: index0_FB + (omega_A - 3) [the lowest FB subscript used in primes
     444             :    * for A]  + (omega_A - 2) <= index2_FB  [the subscript from which the choice
     445             :    * of primes for A starts, putting omega_A - 1 of them at or below index2_FB,
     446             :    * and the last and largest one above, cf. mpqs_si_choose_primes]. Moreover,
     447             :    * index2_FB indicates the last prime below the ideal size, unless (when kN
     448             :    * is tiny) the ideal size was too small to use. */
     449             : }
     450             : 
     451             : /*********************************************************************/
     452             : /**                       SELF-INITIALIZATION                       **/
     453             : /*********************************************************************/
     454             : 
     455             : #ifdef MPQS_DEBUG
     456             : /* Debug-only helper routine: check correctness of the root z mod p_i
     457             :  * by evaluting A * z^2 + B * z + C mod p_i  (which should be 0). */
     458             : static void
     459             : check_root(mpqs_handle_t *h, GEN mC, long p, long start)
     460             : {
     461             :   pari_sp av = avma;
     462             :   long z = start - ((long)(h->M) % p);
     463             :   if (umodiu(subii(mulsi(z, addii(h->B, mulsi(z, h->A))), mC), p))
     464             :   {
     465             :     err_printf("MPQS: p = %ld\n", p);
     466             :     err_printf("MPQS: A = %Ps\n", h->A);
     467             :     err_printf("MPQS: B = %Ps\n", h->B);
     468             :     err_printf("MPQS: C = %Ps\n", negi(mC));
     469             :     err_printf("MPQS: z = %ld\n", z);
     470             :     pari_err_BUG("MPQS: self_init: found wrong polynomial");
     471             :   }
     472             :   avma = av;
     473             : }
     474             : #endif
     475             : 
     476             : /* Increment *x > 0 to a larger value which has the same number of 1s in its
     477             :  * binary representation.  Wraparound can be detected by the caller as long as
     478             :  * we keep total_no_of_primes_for_A strictly less than BITS_IN_LONG.
     479             :  *
     480             :  * Changed switch to increment *x in all cases to the next larger number
     481             :  * which (a) has the same count of 1 bits and (b) does not arise from the
     482             :  * old value by moving a single 1 bit one position to the left  (which was
     483             :  * undesirable for the sieve). --GN based on discussion with TP */
     484             : INLINE void
     485        2170 : mpqs_increment(mpqs_uint32_t *x)
     486             : {
     487        2170 :   mpqs_uint32_t r1_mask, r01_mask, slider=1UL;
     488             : 
     489        2170 :   switch (*x & 0x1F)
     490             :   { /* 32-way computed jump handles 22 out of 32 cases */
     491             :   case 29:
     492          81 :     (*x)++; break; /* shifts a single bit, but we postprocess this case */
     493             :   case 26:
     494           0 :     (*x) += 2; break; /* again */
     495             :   case 1: case 3: case 6: case 9: case 11:
     496             :   case 17: case 19: case 22: case 25: case 27:
     497        1067 :     (*x) += 3; return;
     498             :   case 20:
     499           6 :     (*x) += 4; break; /* again */
     500             :   case 5: case 12: case 14: case 21:
     501          27 :     (*x) += 5; return;
     502             :   case 2: case 7: case 13: case 18: case 23:
     503         501 :     (*x) += 6; return;
     504             :   case 10:
     505           0 :     (*x) += 7; return;
     506             :   case 8:
     507           0 :     (*x) += 8; break; /* and again */
     508             :   case 4: case 15:
     509         148 :     (*x) += 12; return;
     510             :   default: /* 0, 16, 24, 28, 30, 31 */
     511             :     /* isolate rightmost 1 */
     512         340 :     r1_mask = ((*x ^ (*x - 1)) + 1) >> 1;
     513             :     /* isolate rightmost 1 which has a 0 to its left */
     514         340 :     r01_mask = ((*x ^ (*x + r1_mask)) + r1_mask) >> 2;
     515             :     /* simple cases.  Both of these shift a single bit one position to the
     516             :        left, and will need postprocessing */
     517         340 :     if (r1_mask == r01_mask) { *x += r1_mask; break; }
     518         340 :     if (r1_mask == 1) { *x += r01_mask; break; }
     519             :     /* General case: add r01_mask, kill off as many 1 bits as possible to its
     520             :      * right while at the same time filling in 1 bits from the LSB. */
     521         246 :     if (r1_mask == 2) { *x += (r01_mask>>1) + 1; return; }
     522         420 :     while (r01_mask > r1_mask && slider < r1_mask)
     523             :     {
     524         210 :       r01_mask >>= 1; slider <<= 1;
     525             :     }
     526         105 :     *x += r01_mask + slider - 1;
     527         105 :     return;
     528             :   }
     529             :   /* post-process cases which couldn't be finalized above */
     530         181 :   r1_mask = ((*x ^ (*x - 1)) + 1) >> 1;
     531         181 :   r01_mask = ((*x ^ (*x + r1_mask)) + r1_mask) >> 2;
     532         181 :   if (r1_mask == r01_mask) { *x += r1_mask; return; }
     533         181 :   if (r1_mask == 1) { *x += r01_mask; return; }
     534          87 :   if (r1_mask == 2) { *x += (r01_mask>>1) + 1; return; }
     535          30 :   while (r01_mask > r1_mask && slider < r1_mask)
     536             :   {
     537          18 :     r01_mask >>= 1; slider <<= 1;
     538             :   }
     539           6 :   *x += r01_mask + slider - 1;
     540             : }
     541             : 
     542             : /* self-init (1): advancing the bit pattern, and choice of primes for A.
     543             :  * On first call, h->bin_index = 0. On later occasions, we need to begin
     544             :  * by clearing the MPQS_FBE_DIVIDES_A bit in the fbe_flags of the former
     545             :  * prime factors of A (use per_A_pr to find them). Upon successful return, that
     546             :  * array will have been filled in, and the flag bits will have been turned on
     547             :  * again in the right places.
     548             :  * Return 1 when all is fine and 0 when we found we'd be using more bits to
     549             :  * the left in bin_index than we have matching primes in the FB. In the latter
     550             :  * case, bin_index will be zeroed out, index2_FB will be incremented by 2,
     551             :  * index2_moved will be turned on; the caller, after checking that index2_FB
     552             :  * has not become too large, should just call us again, which then succeeds:
     553             :  * we'll start again with a right-justified sequence of 1 bits in bin_index,
     554             :  * now interpreted as selecting primes relative to the new index2_FB. */
     555             : INLINE int
     556        2505 : mpqs_si_choose_primes(mpqs_handle_t *h)
     557             : {
     558        2505 :   mpqs_FB_entry_t *FB = h->FB;
     559        2505 :   mpqs_per_A_prime_t *per_A_pr = h->per_A_pr;
     560        2505 :   double l2_last_p = h->l2_target_A;
     561        2505 :   mpqs_int32_t omega_A = h->omega_A;
     562             :   int i, j, v2, prev_last_p_idx;
     563        2505 :   int room = h->index2_FB - h->index0_FB - omega_A + 4;
     564             :   /* The notion of room here (cf mpqs_locate_A_range() above) is the number
     565             :    * of primes at or below index2_FB which are eligible for A. We need
     566             :    * >= omega_A - 1 of them, and it is guaranteed by mpqs_locate_A_range() that
     567             :    * this many are available: the lowest FB slot used for A is never less than
     568             :    * index0_FB + omega_A - 3. When omega_A = 3 (very small kN), we allow
     569             :    * ourselves to reach all the way down to index0_FB; otherwise, we keep away
     570             :    * from it by at least one position.  For omega_A >= 4 this avoids situations
     571             :    * where the selection of the smaller primes here has advanced to a lot of
     572             :    * very small ones, and the single last larger one has soared away to bump
     573             :    * into the top end of the FB. */
     574             :   mpqs_uint32_t room_mask;
     575             :   mpqs_int32_t p;
     576             :   ulong bits;
     577             : 
     578             :   /* XXX also clear the index_j field here? */
     579        2505 :   if (h->bin_index == 0)
     580             :   { /* first time here, or after increasing index2_FB, initialize to a pattern
     581             :      * of omega_A - 1 consecutive 1 bits. Caller has ensured that there are
     582             :      * enough primes for this in the FB below index2_FB. */
     583         335 :     h->bin_index = (1UL << (omega_A - 1)) - 1;
     584         335 :     prev_last_p_idx = 0;
     585             :   }
     586             :   else
     587             :   { /* clear out old flags */
     588        2170 :     for (i = 0; i < omega_A; i++) MPQS_FLG(i) = MPQS_FBE_CLEAR;
     589        2170 :     prev_last_p_idx = MPQS_I(omega_A-1);
     590             : 
     591        2170 :     if (room > 30) room = 30;
     592        2170 :     room_mask = ~((1UL << room) - 1);
     593             : 
     594             :     /* bump bin_index to next acceptable value. If index2_moved is off, call
     595             :      * mpqs_increment() once; otherwise, repeat until there's something in the
     596             :      * least significant 2 bits - to ensure that we never re-use an A which
     597             :      * we'd used before increasing index2_FB - but also stop if something shows
     598             :      * up in the forbidden bits on the left where we'd run out of bits or walk
     599             :      * beyond index0_FB + omega_A - 3. */
     600        2170 :     mpqs_increment(&h->bin_index);
     601        2170 :     if (h->index2_moved)
     602             :     {
     603           0 :       while ((h->bin_index & (room_mask | 0x3)) == 0)
     604           0 :         mpqs_increment(&h->bin_index);
     605             :     }
     606             :     /* did we fall off the edge on the left? */
     607        2170 :     if ((h->bin_index & room_mask) != 0)
     608             :     { /* Yes. Turn on the index2_moved flag in the handle */
     609           0 :       h->index2_FB += 2; /* caller to check this isn't too large!!! */
     610           0 :       h->index2_moved = 1;
     611           0 :       h->bin_index = 0;
     612           0 :       if (MPQS_DEBUGLEVEL >= 5)
     613           0 :         err_printf("MPQS: wrapping, more primes for A now chosen near FB[%ld] = %ld\n",
     614           0 :                    (long)h->index2_FB,
     615           0 :                    (long)FB[h->index2_FB].fbe_p);
     616           0 :       return 0; /* back off - caller should retry */
     617             :     }
     618             :   }
     619             :   /* assert: we aren't occupying any of the room_mask bits now, and if
     620             :    * index2_moved had already been on, at least one of the two LSBs is on */
     621        2505 :   bits = h->bin_index;
     622        2505 :   if (MPQS_DEBUGLEVEL >= 6)
     623           0 :     err_printf("MPQS: new bit pattern for primes for A: 0x%lX\n", bits);
     624             : 
     625             :   /* map bits to FB subscripts, counting downward with bit 0 corresponding
     626             :    * to index2_FB, and accumulate logarithms against l2_last_p */
     627        2505 :   j = h->index2_FB;
     628        2505 :   v2 = vals((long)bits);
     629        2505 :   if (v2) { j -= v2; bits >>= v2; }
     630       10277 :   for (i = omega_A - 2; i >= 0; i--)
     631             :   {
     632       10277 :     MPQS_I(i) = j;
     633       10277 :     l2_last_p -= MPQS_LP(i);
     634       10277 :     MPQS_FLG(i) |= MPQS_FBE_DIVIDES_A;
     635       10277 :     bits &= ~1UL;
     636       10277 :     if (!bits) break; /* i = 0 */
     637        7772 :     v2 = vals((long)bits); /* > 0 */
     638        7772 :     bits >>= v2; j -= v2;
     639             :   }
     640             :   /* Choose the larger prime.  Note we keep index2_FB <= size_of_FB - 3 */
     641       49217 :   for (j = h->index2_FB + 1; (p = FB[j].fbe_p); j++)
     642       49217 :     if (FB[j].fbe_flogp > l2_last_p) break;
     643             :   /* The following trick avoids generating a large proportion of duplicate
     644             :    * relations when the last prime falls into an area where there are large
     645             :    * gaps from one FB prime to the next, and would otherwise often be repeated
     646             :    * (so that successive A's would wind up too similar to each other). While
     647             :    * this trick isn't perfect, it gets rid of a major part of the potential
     648             :    * duplication. */
     649        2505 :   if (p && j == prev_last_p_idx) { j++; p = FB[j].fbe_p; }
     650        2505 :   MPQS_I(omega_A - 1) = p? j: h->size_of_FB + 1;
     651        2505 :   MPQS_FLG(omega_A - 1) |= MPQS_FBE_DIVIDES_A;
     652             : 
     653        2505 :   if (MPQS_DEBUGLEVEL >= 6)
     654             :   {
     655           0 :     err_printf("MPQS: chose primes for A");
     656           0 :     for (i = 0; i < omega_A; i++)
     657           0 :       err_printf(" FB[%ld]=%ld%s", (long)MPQS_I(i), (long)MPQS_AP(i),
     658           0 :                  i < omega_A - 1 ? "," : "\n");
     659             :   }
     660        2505 :   return 1;
     661             : }
     662             : 
     663             : /* There are 4 parts to self-initialization, exercised at different times:
     664             :  * - choosing a new sqfree coef. A (selecting its prime factors, FB bookkeeping)
     665             :  * - doing the actual computations attached to a new A
     666             :  * - choosing a new B keeping the same A (much simpler)
     667             :  * - a small common bit that needs to happen in both cases.
     668             :  * As to the first item, the scheme works as follows: pick omega_A - 1 prime
     669             :  * factors for A below the index2_FB point which marks their ideal size, and
     670             :  * one prime above this point, choosing the latter so log2(A) ~ l2_target_A.
     671             :  * Lower prime factors are chosen using bit patterns of constant weight,
     672             :  * gradually moving away from index2_FB towards smaller FB subscripts.
     673             :  * If this bumps into index0_FB (for very small input), back up by increasing
     674             :  * index2_FB by two, and from then on choosing only bit patterns with either or
     675             :  * both of their bottom bits set, so at least one of the omega_A - 1 smaller
     676             :  * prime factor will be beyond the original index2_FB point. In this way we
     677             :  * avoid re-using the same A. (The choice of the upper "flyer" prime is
     678             :  * constrained by the size of the FB, which normally should never a problem.
     679             :  * For tiny kN, we might have to live with a non-optimal choice.)
     680             :  *
     681             :  * Mathematically, we solve a quadratic (over F_p for each prime p in the FB
     682             :  * which doesn't divide A), a linear equation for each prime p | A, and
     683             :  * precompute differences between roots mod p so we can adjust the roots
     684             :  * quickly when we change B. See Thomas Sosnowski's Diplomarbeit. */
     685             : /* compute coefficients of sieving polynomial for self initializing variant.
     686             :  * Coefficients A and B are set (preallocated GENs) and several tables are
     687             :  * updated. */
     688             : static int
     689       66426 : mpqs_self_init(mpqs_handle_t *h)
     690             : {
     691       66426 :   const ulong size_of_FB = h->size_of_FB + 1;
     692       66426 :   mpqs_FB_entry_t *FB = h->FB;
     693       66426 :   mpqs_inv_A_H_t *inv_A_H = h->inv_A_H;
     694       66426 :   const pari_sp av = avma;
     695       66426 :   GEN p1, A = h->A, B = h->B;
     696       66426 :   mpqs_per_A_prime_t *per_A_pr = h->per_A_pr;
     697             :   long i, j;
     698             : 
     699             : #ifdef MPQS_DEBUG
     700             :   err_printf("MPQS DEBUG: enter self init, avma = 0x%lX\n", (ulong)avma);
     701             : #endif
     702       66426 :   if (++h->index_j == (mpqs_uint32_t)h->no_B)
     703             :   { /* all the B's have been used, choose new A; this is indicated by setting
     704             :      * index_j to 0 */
     705        2170 :     h->index_j = 0;
     706        2170 :     h->index_i++; /* count finished A's */
     707             :   }
     708             : 
     709       66426 :   if (h->index_j == 0)
     710             :   { /* compute first polynomial with new A */
     711             :     GEN a, b, A2;
     712        2505 :     if (!mpqs_si_choose_primes(h))
     713             :     { /* Ran out of room towards small primes, and index2_FB was raised. */
     714           0 :       if (size_of_FB - h->index2_FB < 4) return 0; /* Fail */
     715           0 :       (void)mpqs_si_choose_primes(h); /* now guaranteed to succeed */
     716             :     }
     717             :     /* bin_index and per_A_pr now populated with consistent values */
     718             : 
     719             :     /* compute A = product of omega_A primes given by bin_index */
     720        2505 :     a = b = NULL;
     721       15287 :     for (i = 0; i < h->omega_A; i++)
     722             :     {
     723       12782 :       ulong p = MPQS_AP(i);
     724       12782 :       a = a? muliu(a, p): utoipos(p);
     725             :     }
     726        2505 :     affii(a, A);
     727             :     /* Compute H[i], 0 <= i < omega_A.  Also compute the initial
     728             :      * B = sum(v_i*H[i]), by taking all v_i = +1
     729             :      * TODO: following needs to be changed later for segmented FB and sieve
     730             :      * interval, where we'll want to precompute several B's. */
     731       15287 :     for (i = 0; i < h->omega_A; i++)
     732             :     {
     733       12782 :       ulong p = MPQS_AP(i);
     734       12782 :       GEN t = divis(A, (long)p);
     735       12782 :       t = remii(mulii(t, muluu(Fl_inv(umodiu(t, p), p), MPQS_SQRT(i))), A);
     736       12782 :       affii(t, MPQS_H(i));
     737       12782 :       b = b? addii(b, t): t;
     738             :     }
     739        2505 :     affii(b, B); set_avma(av);
     740             : 
     741             :     /* ensure B = 1 mod 4 */
     742        2505 :     if (mod2(B) == 0)
     743        1217 :       affii(addii(B, mului(mod4(A), A)), B); /* B += (A % 4) * A; */
     744             : 
     745        2505 :     A2 = shifti(A, 1);
     746             :     /* compute the roots z1, z2, of the polynomial Q(x) mod p_j and
     747             :      * initialize start1[i] with the first value p_i | Q(z1 + i p_j)
     748             :      * initialize start2[i] with the first value p_i | Q(z2 + i p_j)
     749             :      * The following loop does The Right Thing for primes dividing k (where
     750             :      * sqrt_kN is 0 mod p). Primes dividing A are skipped here, and are handled
     751             :      * further down in the common part of SI. */
     752     2729338 :     for (j = 3; (ulong)j <= size_of_FB; j++)
     753             :     {
     754             :       ulong s, mb, t, m, p, iA2, iA;
     755     2726833 :       if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
     756     2714051 :       p = (ulong)FB[j].fbe_p;
     757     2714051 :       m = h->M % p;
     758     2714051 :       iA2 = Fl_inv(umodiu(A2, p), p); /* = 1/(2*A) mod p_j */
     759     2714051 :       iA = iA2 << 1; if (iA > p) iA -= p;
     760     2714051 :       mb = umodiu(B, p); if (mb) mb = p - mb; /* mb = -B mod p */
     761     2714051 :       s = FB[j].fbe_sqrt_kN;
     762     2714051 :       t = Fl_add(m, Fl_mul(Fl_sub(mb, s, p), iA2, p), p);
     763     2714051 :       FB[j].fbe_start1 = (mpqs_int32_t)t;
     764     2714051 :       FB[j].fbe_start2 = (mpqs_int32_t)Fl_add(t, Fl_mul(s, iA, p), p);
     765    17489885 :       for (i = 0; i < h->omega_A - 1; i++)
     766             :       {
     767    14775834 :         ulong h = umodiu(MPQS_H(i), p);
     768    14775834 :         MPQS_INV_A_H(i,j) = Fl_mul(h, iA, p); /* 1/A * H[i] mod p_j */
     769             :       }
     770             :     }
     771             :   }
     772             :   else
     773             :   { /* no "real" computation -- use recursive formula */
     774             :     /* The following exploits that B is the sum of omega_A terms +-H[i]. Each
     775             :      * time we switch to a new B, we choose a new pattern of signs; the
     776             :      * precomputation of the inv_A_H array allows us to change the two
     777             :      * arithmetic progressions equally fast. The choice of sign patterns does
     778             :      * not follow the bit pattern of the ordinal number of B in the current
     779             :      * cohort; rather, we use a Gray code, changing only one sign each time.
     780             :      * When the i-th rightmost bit of the new ordinal number index_j of B is 1,
     781             :      * the sign of H[i] is changed; the next bit to the left tells us whether
     782             :      * we should be adding or subtracting the difference term. We never need to
     783             :      * change the sign of H[omega_A-1] (the topmost one), because that would
     784             :      * just give us the same sieve items Q(x) again with the opposite sign
     785             :      * of x.  This is why we only precomputed inv_A_H up to i = omega_A - 2. */
     786       63921 :     ulong p, v2 = vals(h->index_j); /* new starting positions for sieving */
     787       63921 :     j = h->index_j >> v2;
     788       63921 :     p1 = shifti(MPQS_H(v2), 1);
     789       63921 :     if (j & 2)
     790             :     { /* j = 3 mod 4 */
     791    64087375 :       for (j = 3; (ulong)j <= size_of_FB; j++)
     792             :       {
     793    64056732 :         if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
     794    63859696 :         p = (ulong)FB[j].fbe_p;
     795    63859696 :         FB[j].fbe_start1 = Fl_sub(FB[j].fbe_start1, MPQS_INV_A_H(v2,j), p);
     796    63859696 :         FB[j].fbe_start2 = Fl_sub(FB[j].fbe_start2, MPQS_INV_A_H(v2,j), p);
     797             :       }
     798       30643 :       p1 = addii(B, p1);
     799             :     }
     800             :     else
     801             :     { /* j = 1 mod 4 */
     802    67020632 :       for (j = 3; (ulong)j <= size_of_FB; j++)
     803             :       {
     804    66987354 :         if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
     805    66776710 :         p = (ulong)FB[j].fbe_p;
     806    66776710 :         FB[j].fbe_start1 = Fl_add(FB[j].fbe_start1, MPQS_INV_A_H(v2,j), p);
     807    66776710 :         FB[j].fbe_start2 = Fl_add(FB[j].fbe_start2, MPQS_INV_A_H(v2,j), p);
     808             :       }
     809       33278 :       p1 = subii(B, p1);
     810             :     }
     811       63921 :     affii(p1, B);
     812             :   }
     813             : 
     814             :   /* p=2 is a special case.  start1[2], start2[2] are never looked at,
     815             :    * so don't bother setting them. */
     816             : 
     817             :   /* compute zeros of polynomials that have only one zero mod p since p | A */
     818       66426 :   p1 = diviiexact(subii(h->kN, sqri(B)), shifti(A, 2)); /* coefficient -C */
     819      486888 :   for (i = 0; i < h->omega_A; i++)
     820             :   {
     821      420462 :     ulong p = MPQS_AP(i), s = h->M + Fl_div(umodiu(p1, p), umodiu(B, p), p);
     822      420462 :     FB[MPQS_I(i)].fbe_start1 = FB[MPQS_I(i)].fbe_start2 = (mpqs_int32_t)(s % p);
     823             :   }
     824             : #ifdef MPQS_DEBUG
     825             :   for (j = 3; j <= size_of_FB; j++)
     826             :   {
     827             :     check_root(h, p1, FB[j].fbe_p, FB[j].fbe_start1);
     828             :     check_root(h, p1, FB[j].fbe_p, FB[j].fbe_start2);
     829             :   }
     830             : #endif
     831       66426 :   if (MPQS_DEBUGLEVEL >= 6)
     832           0 :     err_printf("MPQS: chose Q_%ld(x) = %Ps x^2 %c %Ps x + C\n",
     833           0 :                (long) h->index_j, h->A,
     834           0 :                signe(h->B) < 0? '-': '+', absi_shallow(h->B));
     835       66426 :   set_avma(av);
     836             : #ifdef MPQS_DEBUG
     837             :   err_printf("MPQS DEBUG: leave self init, avma = 0x%lX\n", (ulong)avma);
     838             : #endif
     839       66426 :   return 1;
     840             : }
     841             : 
     842             : /*********************************************************************/
     843             : /**                           THE SIEVE                             **/
     844             : /*********************************************************************/
     845             : /* p4 = 4*p, logp ~ log(p), B/E point to the beginning/end of a sieve array */
     846             : INLINE void
     847   266661458 : mpqs_sieve_p(unsigned char *B, unsigned char *E, long p4, long p,
     848             :              unsigned char logp)
     849             : {
     850   266661458 :   register unsigned char *e = E - p4;
     851             :   /* Unrolled loop. It might be better to let the compiler worry about this
     852             :    * kind of optimization, based on its knowledge of whatever useful tricks the
     853             :    * machine instruction set architecture is offering */
     854  1180070562 :   while (e - B >= 0) /* signed comparison */
     855             :   {
     856   646747646 :     (*B) += logp, B += p;
     857   646747646 :     (*B) += logp, B += p;
     858   646747646 :     (*B) += logp, B += p;
     859   646747646 :     (*B) += logp, B += p;
     860             :   }
     861   266661458 :   while (E - B >= 0) (*B) += logp, B += p;
     862   266661458 : }
     863             : 
     864             : static void
     865       66426 : mpqs_sieve(mpqs_handle_t *h)
     866             : {
     867       66426 :   long p, l = h->index1_FB;
     868             :   mpqs_FB_entry_t *FB;
     869       66426 :   unsigned char *S = h->sieve_array, *Send = h->sieve_array_end;
     870             : 
     871       66426 :   memset((void*)S, 0, (h->M << 1) * sizeof(unsigned char));
     872   133607386 :   for (FB = &(h->FB[l]); (p = FB->fbe_p); FB++) /* l++ */
     873             :   {
     874   133540960 :     unsigned char logp = FB->fbe_logval;
     875   133540960 :     long s1 = FB->fbe_start1, s2 = FB->fbe_start2;
     876             :     /* sieve with FB[l] from start1[l], and from start2[l] if s1 != s2 */
     877   133540960 :     mpqs_sieve_p(S + s1, Send, p << 2, p, logp);
     878   133540960 :     if (s1 != s2) mpqs_sieve_p(S + s2, Send, p << 2, p, logp);
     879             :   }
     880       66426 : }
     881             : 
     882             : /* Could use the fact that 4 | M, but let the compiler worry about unrolling. */
     883             : static long
     884       66426 : mpqs_eval_sieve(mpqs_handle_t *h)
     885             : {
     886       66426 :   long x = 0, count = 0, M2 = h->M << 1;
     887       66426 :   unsigned char t = h->sieve_threshold;
     888       66426 :   unsigned char *S = h->sieve_array;
     889       66426 :   long *cand = h->candidates;
     890             : 
     891             :   /* Exploiting the sentinel, we don't need to check for x < M2 in the inner
     892             :    * while loop; more than makes up for the lack of explicit unrolling. */
     893      455266 :   while (count < MPQS_CANDIDATE_ARRAY_SIZE - 1)
     894             :   {
     895      388840 :     while (S[x] < t) x++;
     896      388840 :     if (x >= M2) break;
     897      322414 :     cand[count++] = x++;
     898             :   }
     899       66426 :   cand[count] = 0; return count;
     900             : }
     901             : 
     902             : /*********************************************************************/
     903             : /**                     CONSTRUCTING RELATIONS                      **/
     904             : /*********************************************************************/
     905             : 
     906             : /* only used for debugging */
     907             : static void
     908           0 : split_relp(GEN rel, GEN *prelp, GEN *prelc)
     909             : {
     910           0 :   long j, l = lg(rel);
     911             :   GEN relp, relc;
     912           0 :   *prelp = relp = cgetg(l, t_VECSMALL);
     913           0 :   *prelc = relc = cgetg(l, t_VECSMALL);
     914           0 :   for (j=1; j<l; j++)
     915             :   {
     916           0 :     relc[j] = rel[j] >> REL_OFFSET;
     917           0 :     relp[j] = rel[j] & REL_MASK;
     918             :   }
     919           0 : }
     920             : 
     921             : #ifdef MPQS_DEBUG
     922             : static GEN
     923             : mpqs_factorback(mpqs_handle_t *h, GEN relp)
     924             : {
     925             :   GEN N = h->N, Q = gen_1;
     926             :   long j, l = lg(relp);
     927             :   for (j = 1; j < l; j++)
     928             :   {
     929             :     long e = relp[j] >> REL_OFFSET, i = relp[j] & REL_MASK;
     930             :     if (i == 1) Q = Fp_neg(Q,N); /* special case -1 */
     931             :     else Q = Fp_mul(Q, Fp_powu(utoipos(h->FB[i].fbe_p), e, N), N);
     932             :   }
     933             :   return Q;
     934             : }
     935             : static void
     936             : mpqs_check_rel(mpqs_handle_t *h, GEN c)
     937             : {
     938             :   pari_sp av = avma;
     939             :   int LP = (lg(c) == 4);
     940             :   GEN rhs = mpqs_factorback(h, rel_p(c));
     941             :   GEN Y = rel_Y(c), Qx_2 = remii(sqri(Y), h->N);
     942             :   if (LP) rhs = modii(mulii(rhs, rel_q(c)), h->N);
     943             :   if (!equalii(Qx_2, rhs))
     944             :   {
     945             :     GEN relpp, relpc;
     946             :     split_relp(rel_p(c), &relpp, &relpc);
     947             :     err_printf("MPQS: %Ps : %Ps %Ps\n", Y, relpp,relpc);
     948             :     err_printf("\tQx_2 = %Ps\n", Qx_2);
     949             :     err_printf("\t rhs = %Ps\n", rhs);
     950             :     pari_err_BUG(LP? "MPQS: wrong large prime relation found"
     951             :                    : "MPQS: wrong full relation found");
     952             :   }
     953             :   PRINT_IF_VERBOSE(LP? "\b(;)": "\b(:)");
     954             :   set_avma(av);
     955             : }
     956             : #endif
     957             : 
     958             : static void
     959      147696 : rel_to_ei(GEN ei, GEN relp)
     960             : {
     961      147696 :   long j, l = lg(relp);
     962     2383871 :   for (j=1; j<l; j++)
     963             :   {
     964     2236175 :     long e = relp[j] >> REL_OFFSET, i = relp[j] & REL_MASK;
     965     2236175 :     ei[i] += e;
     966             :   }
     967      147696 : }
     968             : static void
     969     5107169 : mpqs_add_factor(GEN relp, long *i, ulong ei, ulong pi)
     970     5107169 : { relp[++*i] = pi | (ei << REL_OFFSET); }
     971             : 
     972             : static GEN
     973       17779 : combine_large_primes(mpqs_handle_t *h, GEN rel1, GEN rel2)
     974             : {
     975       17779 :   GEN new_Y, new_Y1, Y1 = rel_Y(rel1), Y2 = rel_Y(rel2);
     976       17779 :   long l, lei = h->size_of_FB + 1, nb = 0;
     977       17779 :   GEN ei, relp, iq, q = rel_q(rel1);
     978             : 
     979       17779 :   if (!invmod(q, h->N, &iq)) return equalii(iq, h->N)? NULL: iq; /* rare */
     980       17779 :   ei = zero_zv(lei);
     981       17779 :   rel_to_ei(ei, rel_p(rel1));
     982       17779 :   rel_to_ei(ei, rel_p(rel2));
     983       17779 :   new_Y = modii(mulii(mulii(Y1, Y2), iq), h->N);
     984       17779 :   new_Y1 = subii(h->N, new_Y);
     985       17779 :   if (abscmpii(new_Y1, new_Y) < 0) new_Y = new_Y1;
     986       17779 :   relp = cgetg(MAX_PE_PAIR+1,t_VECSMALL);
     987       17779 :   if (odd(ei[1])) mpqs_add_factor(relp, &nb, 1, 1);
     988    37955416 :   for (l = 2; l <= lei; l++)
     989    37937637 :     if (ei[l]) mpqs_add_factor(relp, &nb, ei[l],l);
     990       17779 :   setlg(relp, nb+1);
     991       17779 :   if (DEBUGLEVEL >= 6)
     992             :   {
     993             :     GEN relpp, relpc, rel1p, rel1c, rel2p, rel2c;
     994           0 :     split_relp(relp,&relpp,&relpc);
     995           0 :     split_relp(rel1,&rel1p,&rel1c);
     996           0 :     split_relp(rel2,&rel2p,&rel2c);
     997           0 :     err_printf("MPQS: combining\n");
     998           0 :     err_printf("    {%Ps @ %Ps : %Ps}\n", q, Y1, rel1p, rel1c);
     999           0 :     err_printf("  * {%Ps @ %Ps : %Ps}\n", q, Y2, rel2p, rel2c);
    1000           0 :     err_printf(" == {%Ps, %Ps}\n", relpp, relpc);
    1001             :   }
    1002             : #ifdef MPQS_DEBUG
    1003             :   {
    1004             :     pari_sp av1 = avma;
    1005             :     if (!equalii(modii(sqri(new_Y), h->N), mpqs_factorback(h, relp)))
    1006             :       pari_err_BUG("MPQS: combined large prime relation is false");
    1007             :     set_avma(av1);
    1008             :   }
    1009             : #endif
    1010       17779 :   return mkvec2(new_Y, relp);
    1011             : }
    1012             : 
    1013             : /* nc candidates */
    1014             : static GEN
    1015       65092 : mpqs_eval_cand(mpqs_handle_t *h, long nc, hashtable *frel, hashtable *lprel)
    1016             : {
    1017       65092 :   mpqs_FB_entry_t *FB = h->FB;
    1018       65092 :   GEN A = h->A, B = h->B;
    1019       65092 :   long *relaprimes = h->relaprimes, *candidates = h->candidates;
    1020             :   long pi, i;
    1021             :   int pii;
    1022             : 
    1023      387506 :   for (i = 0; i < nc; i++)
    1024             :   {
    1025      322414 :     pari_sp btop = avma;
    1026      322414 :     GEN Qx, Qx_part, Y, relp = cgetg(MAX_PE_PAIR+1,t_VECSMALL);
    1027      322414 :     long powers_of_2, p, x = candidates[i], nb = 0;
    1028      322414 :     int relaprpos = 0;
    1029             : 
    1030             :     /* Y = 2*A*x + B, Qx = Y^2/(4*A) = Q(x) */
    1031      322414 :     Y = addii(mulis(A, 2 * (x - h->M)), B);
    1032      322414 :     Qx = subii(sqri(Y), h->kN); /* != 0 since N not a square and (N,k) = 1 */
    1033      322414 :     if (signe(Qx) < 0)
    1034             :     {
    1035      175752 :       setabssign(Qx);
    1036      175752 :       mpqs_add_factor(relp, &nb, 1, 1); /* i = 1, ei = 1, pi */
    1037             :     }
    1038             :     /* Qx > 0, divide by powers of 2; we're really dealing with 4*A*Q(x), so we
    1039             :      * always have at least 2^2 here, and at least 2^3 when kN = 1 mod 4 */
    1040      322414 :     powers_of_2 = vali(Qx);
    1041      322414 :     Qx = shifti(Qx, -powers_of_2);
    1042      322414 :     mpqs_add_factor(relp, &nb, powers_of_2, 2); /* i = 1, ei = 1, pi */
    1043             :     /* When N is small, it may happen that N | Qx outright. In any case, when
    1044             :      * no extensive prior trial division / Rho / ECM was attempted, gcd(Qx,N)
    1045             :      * may turn out to be a nontrivial factor of N (not in FB or we'd have
    1046             :      * found it already, but possibly smaller than the large prime bound). This
    1047             :      * is too rare to check for here in the inner loop, but it will be caught
    1048             :      * if such an LP relation is ever combined with another. */
    1049             : 
    1050             :     /* Pass 1 over odd primes in FB: pick up all possible divisors of Qx
    1051             :      * including those sitting in k or in A, and remember them in relaprimes.
    1052             :      * Do not yet worry about possible repeated factors, these will be found in
    1053             :      * the Pass 2. Pass 1 recognizes divisors of A by their corresponding flags
    1054             :      * bit in the FB entry. (Divisors of k are ignored at this stage.)
    1055             :      * We construct a preliminary table of FB subscripts and "exponents" of FB
    1056             :      * primes which divide Qx. (We store subscripts, not the primes themselves.)
    1057             :      * We distinguish three cases:
    1058             :      * 0) prime in A which does not divide Qx/A,
    1059             :      * 1) prime not in A which divides Qx/A,
    1060             :      * 2) prime in A which divides Qx/A.
    1061             :      * Cases 1 and 2 need checking for repeated factors, kind 0 doesn't.
    1062             :      * Cases 0 and 1 contribute 1 to the exponent in the relation, case 2
    1063             :      * contributes 2.
    1064             :      * Factors in common with k are simpler: if they occur, they occur
    1065             :      * exactly to the first power, and this makes no difference in Pass 1,
    1066             :      * so they behave just like every normal odd FB prime. */
    1067   550981701 :     for (Qx_part = A, pi = 3; (p = FB[pi].fbe_p); pi++)
    1068             :     {
    1069   550659287 :       long xp = x % p;
    1070   550659287 :       ulong ei = FB[pi].fbe_flags & MPQS_FBE_DIVIDES_A;
    1071             :       /* Here we used that MPQS_FBE_DIVIDES_A = 1. */
    1072             : 
    1073   550659287 :       if (xp == FB[pi].fbe_start1 || xp == FB[pi].fbe_start2)
    1074             :       { /* p divides Q(x)/A and possibly A, case 2 or 3 */
    1075     2251865 :         relaprimes[relaprpos++] = pi;
    1076     2251865 :         relaprimes[relaprpos++] = 1 + ei;
    1077     2251865 :         Qx_part = muliu(Qx_part, p);
    1078             :       }
    1079   548407422 :       else if (ei)
    1080             :       { /* p divides A but does not divide Q(x)/A, case 1 */
    1081     1939931 :         relaprimes[relaprpos++] = pi;
    1082     1939931 :         relaprimes[relaprpos++] = 0;
    1083             :       }
    1084             :     }
    1085             :     /* We have accumulated the known factors of Qx except for possible repeated
    1086             :      * factors and for possible large primes.  Divide off what we have.
    1087             :      * This is faster than dividing off A and each prime separately. */
    1088      322414 :     Qx = diviiexact(Qx, Qx_part);
    1089             : 
    1090             : #ifdef MPQS_DEBUG
    1091             :     err_printf("MPQS DEBUG: eval loop 3, avma = 0x%lX\n", (ulong)avma);
    1092             : #endif
    1093             :     /* Pass 2: deal with repeated factors and store tentative relation. At this
    1094             :      * point, the only primes which can occur again in the adjusted Qx are
    1095             :      * those in relaprimes which are followed by 1 or 2. We must pick up those
    1096             :      * followed by a 0, too. */
    1097             :     PRINT_IF_VERBOSE("a");
    1098     4514210 :     for (pii = 0; pii < relaprpos; pii += 2)
    1099             :     {
    1100     4191796 :       ulong r, ei = relaprimes[pii+1];
    1101             :       GEN q;
    1102             : 
    1103     4191796 :       pi = relaprimes[pii];
    1104             :       /* p | k (identified by its index before index0_FB)* or p | A (ei = 0) */
    1105     4191796 :       if ((mpqs_int32_t)pi < h->index0_FB || ei == 0)
    1106             :       {
    1107     1962555 :         mpqs_add_factor(relp, &nb, 1, pi);
    1108     1962555 :         continue;
    1109             :       }
    1110     2229241 :       p = FB[pi].fbe_p;
    1111             :       /* p might still divide the current adjusted Qx. Try it. */
    1112     2229241 :       switch(cmpiu(Qx, p))
    1113             :       {
    1114       32273 :         case 0: ei++; Qx = gen_1; break;
    1115             :         case 1:
    1116     1291253 :           q = absdiviu_rem(Qx, p, &r);
    1117     1291253 :           while (r == 0) { ei++; Qx = q; q = absdiviu_rem(Qx, p, &r); }
    1118     1291253 :           break;
    1119             :       }
    1120     2229241 :       mpqs_add_factor(relp, &nb, ei, pi);
    1121             :     }
    1122             : 
    1123             : #ifdef MPQS_DEBUG
    1124             :     err_printf("MPQS DEBUG: eval loop 4, avma = 0x%lX\n", (ulong)avma);
    1125             : #endif
    1126             :     PRINT_IF_VERBOSE("\bb");
    1127      322414 :     setlg(relp, nb+1);
    1128      322414 :     if (is_pm1(Qx))
    1129             :     {
    1130      126274 :       GEN rel = gerepilecopy(btop, mkvec2(absi_shallow(Y), relp));
    1131             : #ifdef MPQS_DEBUG
    1132             :       mpqs_check_rel(h, rel);
    1133             : #endif
    1134      126274 :       frel_add(frel, rel);
    1135             :     }
    1136      196140 :     else if (cmpiu(Qx, h->lp_bound) <= 0)
    1137             :     {
    1138      175990 :       ulong q = itou(Qx);
    1139      175990 :       GEN rel = mkvec3(absi_shallow(Y),relp,Qx);
    1140      175990 :       GEN col = hash_haskey_GEN(lprel, (void*)q);
    1141             : #ifdef MPQS_DEBUG
    1142             :       mpqs_check_rel(h, rel);
    1143             : #endif
    1144      175990 :       if (!col) /* relation up to large prime */
    1145      158211 :         hash_insert(lprel, (void*)q, (void*)gerepilecopy(btop,rel));
    1146       17779 :       else if ((rel = combine_large_primes(h, rel, col)))
    1147             :       {
    1148       17779 :         if (typ(rel) == t_INT) return rel; /* very unlikely */
    1149             : #ifdef MPQS_DEBUG
    1150             :         mpqs_check_rel(h, rel);
    1151             : #endif
    1152       17779 :         frel_add(frel, gerepilecopy(btop,rel));
    1153             :       }
    1154             :       else
    1155           0 :         set_avma(btop);
    1156             :     }
    1157             :     else
    1158             :     { /* TODO: check for double large prime */
    1159             :       PRINT_IF_VERBOSE("\b.");
    1160       20150 :       set_avma(btop);
    1161             :     }
    1162             :   }
    1163             :   PRINT_IF_VERBOSE("\n");
    1164       65092 :   return NULL;
    1165             : }
    1166             : 
    1167             : /*********************************************************************/
    1168             : /**                    FROM RELATIONS TO DIVISORS                   **/
    1169             : /*********************************************************************/
    1170             : 
    1171             : /* create an F2m from a relations list; rows = size_of_FB+1 */
    1172             : static GEN
    1173         335 : rels_to_F2m(GEN rel, long rows)
    1174             : {
    1175         335 :   long i, cols = lg(rel)-1;
    1176         335 :   GEN m = zero_F2m_copy(rows, cols);
    1177      143744 :   for (i = 1; i <= cols; i++)
    1178             :   {
    1179      143409 :     GEN relp = gmael(rel,i,2);
    1180      143409 :     long j, l = lg(relp);
    1181     2399428 :     for (j = 1; j < l; j++)
    1182     2256019 :       if (odd(relp[j] >> REL_OFFSET)) F2m_set(m, relp[j] & REL_MASK, i);
    1183             :   }
    1184         335 :   return m;
    1185             : }
    1186             : 
    1187             : static int
    1188         740 : split(GEN *D, long *e)
    1189             : {
    1190             :   ulong mask;
    1191             :   long flag;
    1192         740 :   if (MR_Jaeschke(*D)) { *e = 1; return 1; } /* probable prime */
    1193          84 :   if (Z_issquareall(*D, D))
    1194             :   { /* squares could cost us a lot of time */
    1195          35 :     if (DEBUGLEVEL >= 5) err_printf("MPQS: decomposed a square\n");
    1196          35 :     *e = 2; return 1;
    1197             :   }
    1198          49 :   mask = 7;
    1199             :   /* 5th/7th powers aren't worth the trouble. OTOH once we have the hooks for
    1200             :    * dealing with cubes, higher powers can be handled essentially for free) */
    1201          49 :   if ((flag = is_357_power(*D, D, &mask)))
    1202             :   {
    1203          14 :     if (DEBUGLEVEL >= 5)
    1204           0 :       err_printf("MPQS: decomposed a %s power\n", uordinal(flag));
    1205          14 :     *e = flag; return 1;
    1206             :   }
    1207          35 :   *e = 0; return 0; /* known composite */
    1208             : }
    1209             : 
    1210             : /* return a GEN structure containing NULL but safe for gerepileupto */
    1211             : static GEN
    1212         335 : mpqs_solve_linear_system(mpqs_handle_t *h, hashtable *frel)
    1213             : {
    1214         335 :   mpqs_FB_entry_t *FB = h->FB;
    1215         335 :   pari_sp av = avma;
    1216         335 :   GEN rels = hash_keys(frel), N = h->N, r, c, res, ei, M, Ker;
    1217             :   long i, j, nrows, rlast, rnext, rmax, rank;
    1218             : 
    1219         335 :   M = rels_to_F2m(rels, h->size_of_FB+1);
    1220         335 :   Ker = F2m_ker_sp(M,0); rank = lg(Ker)-1;
    1221         335 :   if (DEBUGLEVEL >= 4)
    1222             :   {
    1223           0 :     if (DEBUGLEVEL >= 7)
    1224           0 :       err_printf("\\\\ MPQS RELATION MATRIX\nFREL=%Ps\nKERNEL=%Ps\n",M, Ker);
    1225           0 :     err_printf("MPQS: Gauss done: kernel has rank %ld, taking gcds...\n", rank);
    1226             :   }
    1227         335 :   if (!rank)
    1228             :   { /* trivial kernel; main loop may look for more relations */
    1229           0 :     if (DEBUGLEVEL >= 3)
    1230           0 :       pari_warn(warner, "MPQS: no solutions found from linear system solver");
    1231           0 :     return gc_NULL(av); /* no factors found */
    1232             :   }
    1233             : 
    1234             :   /* Expect up to 2^rank pairwise coprime factors, but a kernel basis vector
    1235             :    * may not contribute to the decomposition; r stores the factors and c what
    1236             :    * we know about them (0: composite, 1: probably prime, >= 2: proper power) */
    1237         335 :   ei = cgetg(h->size_of_FB + 2, t_VECSMALL);
    1238         335 :   rmax = logint(N, utoi(3));
    1239         335 :   if (rank <= BITS_IN_LONG-2)
    1240         285 :     rmax = minss(rmax, 1L<<rank); /* max # of factors we can hope for */
    1241         335 :   r = cgetg(rmax+1, t_VEC);
    1242         335 :   c = cgetg(rmax+1, t_VECSMALL);
    1243         335 :   rnext = rlast = 1;
    1244         335 :   nrows = lg(M)-1;
    1245        1208 :   for (i = 1; i <= rank; i++)
    1246             :   { /* loop over kernel basis */
    1247        1208 :     GEN X = gen_1, Y_prod = gen_1, X_plus_Y, D;
    1248        1208 :     pari_sp av2 = avma, av3;
    1249        1208 :     long done = 0; /* # probably-prime factors or powers whose bases we won't
    1250             :                     * handle any further */
    1251        1208 :     memset((void *)(ei+1), 0, (h->size_of_FB + 1) * sizeof(long));
    1252     1403074 :     for (j = 1; j <= nrows; j++)
    1253     1401866 :       if (F2m_coeff(Ker, j, i))
    1254             :       {
    1255      112138 :         GEN R = gel(rels,j);
    1256      112138 :         Y_prod = gerepileuptoint(av2, remii(mulii(Y_prod, gel(R,1)), N));
    1257      112138 :         rel_to_ei(ei, gel(R,2));
    1258             :       }
    1259        1208 :     av3 = avma;
    1260     1390037 :     for (j = 2; j <= h->size_of_FB + 1; j++)
    1261     1388829 :       if (ei[j])
    1262             :       {
    1263      188086 :         GEN q = utoipos(FB[j].fbe_p);
    1264      188086 :         if (ei[j] & 1) pari_err_BUG("MPQS (relation is a nonsquare)");
    1265      188086 :         X = remii(mulii(X, Fp_powu(q, (ulong)ei[j]>>1, N)), N);
    1266      188086 :         X = gerepileuptoint(av3, X);
    1267             :       }
    1268        1208 :     if (MPQS_DEBUGLEVEL >= 1 && !dvdii(subii(sqri(X), sqri(Y_prod)), N))
    1269             :     {
    1270           0 :       err_printf("MPQS: X^2 - Y^2 != 0 mod N\n");
    1271           0 :       err_printf("\tindex i = %ld\n", i);
    1272           0 :       pari_warn(warner, "MPQS: wrong relation found after Gauss");
    1273             :     }
    1274             :     /* At this point, gcd(X-Y, N) * gcd(X+Y, N) = N:
    1275             :      * 1) N | X^2 - Y^2, so it divides the LHS;
    1276             :      * 2) let P be any prime factor of N. If P | X-Y and P | X+Y, then P | 2X
    1277             :      * But X is a product of powers of FB primes => coprime to N.
    1278             :      * Hence we work with gcd(X+Y, N) alone. */
    1279        1208 :     X_plus_Y = addii(X, Y_prod);
    1280        1208 :     if (rnext == 1)
    1281             :     { /* we still haven't decomposed, and want both a gcd and its cofactor. */
    1282        1113 :       D = gcdii(X_plus_Y, N);
    1283        1113 :       if (is_pm1(D) || equalii(D,N)) { set_avma(av2); continue; }
    1284             :       /* got something that works */
    1285         335 :       if (DEBUGLEVEL >= 5)
    1286           0 :         err_printf("MPQS: splitting N after %ld kernel vector%s\n",
    1287             :                    i+1, (i? "s" : ""));
    1288         335 :       gel(r,1) = diviiexact(N, D);
    1289         335 :       gel(r,2) = D;
    1290         335 :       rlast = rnext = 3;
    1291         335 :       if (split(&gel(r,1), &c[1])) done++;
    1292         335 :       if (split(&gel(r,2), &c[2])) done++;
    1293         335 :       if (done == 2 || rmax == 2) break;
    1294          35 :       if (DEBUGLEVEL >= 5)
    1295           0 :         err_printf("MPQS: got two factors, looking for more...\n");
    1296             :     }
    1297             :     else
    1298             :     { /* we already have factors */
    1299         320 :       for (j = 1; j < rnext; j++)
    1300             :       { /* loop over known-composite factors */
    1301             :         /* skip probable primes and also roots of pure powers: they are a lot
    1302             :          * smaller than N and should be easy to deal with later */
    1303         225 :         if (c[j]) { done++; continue; }
    1304          95 :         av3 = avma; D = gcdii(X_plus_Y, gel(r,j));
    1305          95 :         if (is_pm1(D) || equalii(D, gel(r,j))) { set_avma(av3); continue; }
    1306             :         /* got one which splits this factor */
    1307          35 :         if (DEBUGLEVEL >= 5)
    1308           0 :           err_printf("MPQS: resplitting a factor after %ld kernel vectors\n",
    1309             :                      i+1);
    1310          35 :         gel(r,j) = diviiexact(gel(r,j), D);
    1311          35 :         gel(r,rnext) = D;
    1312          35 :         if (split(&gel(r,j), &c[j])) done++;
    1313             :         /* Don't increment done: happens later when we revisit c[rnext] during
    1314             :          * the present inner loop. */
    1315          35 :         (void)split(&gel(r,rnext), &c[rnext]);
    1316          35 :         if (++rnext > rmax) break; /* all possible factors seen */
    1317             :       } /* loop over known composite factors */
    1318             : 
    1319          95 :       if (rnext > rlast)
    1320             :       {
    1321          35 :         if (DEBUGLEVEL >= 5)
    1322           0 :           err_printf("MPQS: got %ld factors%s\n", rlast - 1,
    1323             :                      (done < rlast ? ", looking for more..." : ""));
    1324          35 :         rlast = rnext;
    1325             :       }
    1326             :       /* break out if we have rmax factors or all current factors are probable
    1327             :        * primes or tiny roots from pure powers */
    1328          95 :       if (rnext > rmax || done == rnext - 1) break;
    1329             :     }
    1330             :   }
    1331         335 :   if (rnext == 1) return gc_NULL(av); /* no factors found */
    1332             : 
    1333             :   /* normal case: convert to ifac format as described in ifactor1.c (value,
    1334             :    * exponent, class [unknown, known composite, known prime]) */
    1335         335 :   rlast = rnext - 1; /* # of distinct factors found */
    1336         335 :   res = cgetg(3*rlast + 1, t_VEC);
    1337         335 :   if (DEBUGLEVEL >= 6) err_printf("MPQS: wrapping up %ld factors\n", rlast);
    1338        1040 :   for (i = j = 1; i <= rlast; i++, j += 3)
    1339             :   {
    1340         705 :     long C  = c[i];
    1341         705 :     icopyifstack(gel(r,i), gel(res,j)); /* factor */
    1342         705 :     gel(res,j+1) = C <= 1? gen_1: utoipos(C); /* exponent */
    1343         705 :     gel(res,j+2) = C ? NULL: gen_0; /* unknown or known composite */
    1344         705 :     if (DEBUGLEVEL >= 6)
    1345           0 :       err_printf("\tpackaging %ld: %Ps ^%ld (%s)\n", i, gel(r,i),
    1346           0 :                  itos(gel(res,j+1)), (C? "unknown": "composite"));
    1347             :   }
    1348         335 :   return res;
    1349             : }
    1350             : 
    1351             : /*********************************************************************/
    1352             : /**               MAIN ENTRY POINT AND DRIVER ROUTINE               **/
    1353             : /*********************************************************************/
    1354             : static void
    1355           7 : toolarge()
    1356           7 : { pari_warn(warner, "MPQS: number too big to be factored with MPQS,\n\tgiving up"); }
    1357             : 
    1358             : /* Factors N using the self-initializing multipolynomial quadratic sieve
    1359             :  * (SIMPQS).  Returns one of the two factors, or (usually) a vector of factors
    1360             :  * and exponents and information about which ones are still composite, or NULL
    1361             :  * when we can't seem to make any headway. */
    1362             : GEN
    1363         342 : mpqs(GEN N)
    1364             : {
    1365         342 :   const long size_N = decimal_len(N);
    1366             :   mpqs_handle_t H;
    1367             :   GEN fact; /* will in the end hold our factor(s) */
    1368             :   mpqs_FB_entry_t *FB; /* factor base */
    1369             :   double dbg_target, DEFEAT;
    1370             :   ulong p;
    1371             :   pari_timer T;
    1372             :   hashtable lprel, frel;
    1373         342 :   pari_sp av = avma;
    1374             : 
    1375         342 :   if (DEBUGLEVEL >= 4) err_printf("MPQS: number to factor N = %Ps\n", N);
    1376         342 :   if (size_N > MPQS_MAX_DIGIT_SIZE_KN) { toolarge(); return NULL; }
    1377         335 :   if (DEBUGLEVEL >= 4)
    1378             :   {
    1379           0 :     timer_start(&T);
    1380           0 :     err_printf("MPQS: factoring number of %ld decimal digits\n", size_N);
    1381             :   }
    1382         335 :   H.N = N;
    1383         335 :   H.bin_index = 0;
    1384         335 :   H.index_i = 0;
    1385         335 :   H.index_j = 0;
    1386         335 :   H.index2_moved = 0;
    1387         335 :   p = mpqs_find_k(&H);
    1388         335 :   if (p) { set_avma(av); return utoipos(p); }
    1389         335 :   if (DEBUGLEVEL >= 5)
    1390           0 :     err_printf("MPQS: found multiplier %ld for N\n", H._k->k);
    1391         335 :   H.kN = muliu(N, H._k->k);
    1392         335 :   if (!mpqs_set_parameters(&H)) { toolarge(); return NULL; }
    1393             : 
    1394         335 :   if (DEBUGLEVEL >= 5)
    1395           0 :     err_printf("MPQS: creating factor base and allocating arrays...\n");
    1396         335 :   FB = mpqs_create_FB(&H, &p);
    1397         335 :   if (p) { set_avma(av); return utoipos(p); }
    1398         335 :   mpqs_sieve_array_ctor(&H);
    1399         335 :   mpqs_poly_ctor(&H);
    1400             : 
    1401         335 :   H.lp_bound = minss(H.largest_FB_p, MPQS_LP_BOUND);
    1402             :   /* don't allow large primes to have room for two factors both bigger than
    1403             :    * what the FB contains (...yet!) */
    1404         335 :   H.lp_bound *= minss(H.lp_scale, H.largest_FB_p - 1);
    1405         335 :   H.dkN = gtodouble(H.kN);
    1406             :   /* compute the threshold and fill in the byte-sized scaled logarithms */
    1407         335 :   mpqs_set_sieve_threshold(&H);
    1408         335 :   if (!mpqs_locate_A_range(&H)) return NULL;
    1409         335 :   if (DEBUGLEVEL >= 4)
    1410             :   {
    1411           0 :     err_printf("MPQS: sieving interval = [%ld, %ld]\n", -(long)H.M, (long)H.M);
    1412             :     /* that was a little white lie, we stop one position short at the top */
    1413           0 :     err_printf("MPQS: size of factor base = %ld\n", (long)H.size_of_FB);
    1414           0 :     err_printf("MPQS: striving for %ld relations\n", (long)H.target_rels);
    1415           0 :     err_printf("MPQS: coefficients A will be built from %ld primes each\n",
    1416           0 :                (long)H.omega_A);
    1417           0 :     err_printf("MPQS: primes for A to be chosen near FB[%ld] = %ld\n",
    1418           0 :                (long)H.index2_FB, (long)FB[H.index2_FB].fbe_p);
    1419           0 :     err_printf("MPQS: smallest prime used for sieving FB[%ld] = %ld\n",
    1420           0 :                (long)H.index1_FB, (long)FB[H.index1_FB].fbe_p);
    1421           0 :     err_printf("MPQS: largest prime in FB = %ld\n", (long)H.largest_FB_p);
    1422           0 :     err_printf("MPQS: bound for `large primes' = %ld\n", (long)H.lp_bound);
    1423             :   }
    1424         335 :   if (DEBUGLEVEL >= 5)
    1425             :   {
    1426           0 :     err_printf("MPQS: sieve threshold = %u\n", (unsigned int)H.sieve_threshold);
    1427           0 :     err_printf("MPQS: starting main loop\n");
    1428             :   }
    1429             : 
    1430             :   /* main loop which
    1431             :    * - computes polynomials and their zeros (SI)
    1432             :    * - does the sieving
    1433             :    * - tests candidates of the sieve array */
    1434             : 
    1435             :   /* Let (A, B_i) the current pair of coeffs. If i == 0 a new A is generated */
    1436         335 :   H.index_j = (mpqs_uint32_t)-1;  /* increment below will have it start at 0 */
    1437             : 
    1438         335 :   dbg_target = H.target_rels / 100.;
    1439         335 :   DEFEAT = H.target_rels * 1.5;
    1440         335 :   hash_init_GEN(&frel, H.target_rels, gequal, 1);
    1441         335 :   hash_init_ulong(&lprel,H.target_rels, 1);
    1442             :   for(;;)
    1443       66091 :   {
    1444             :     long tc;
    1445             :     /* self initialization: compute polynomial and its zeros */
    1446       66426 :     if (!mpqs_self_init(&H))
    1447             :     { /* have run out of primes for A; give up */
    1448           0 :       if (DEBUGLEVEL >= 2)
    1449           0 :         err_printf("MPQS: Ran out of primes for A, giving up.\n");
    1450           0 :       return gc_NULL(av);
    1451             :     }
    1452       66426 :     mpqs_sieve(&H);
    1453       66426 :     tc = mpqs_eval_sieve(&H);
    1454       66426 :     if (DEBUGLEVEL >= 6)
    1455           0 :       err_printf("MPQS: found %lu candidate%s\n", tc, (tc==1? "" : "s"));
    1456       66426 :     if (tc)
    1457             :     {
    1458       65092 :       fact = mpqs_eval_cand(&H, tc, &frel, &lprel);
    1459       65092 :       if (fact)
    1460             :       { /* factor found during combining */
    1461           0 :         if (DEBUGLEVEL >= 4)
    1462             :         {
    1463           0 :           err_printf("\nMPQS: split N whilst combining, time = %ld ms\n",
    1464             :                      timer_delay(&T));
    1465           0 :           err_printf("MPQS: found factor = %Ps\n", fact);
    1466             :         }
    1467           0 :         return gerepileupto(av, fact);
    1468             :       }
    1469             :     }
    1470       66426 :     if (DEBUGLEVEL >= 4 && frel.nb > dbg_target)
    1471             :     {
    1472           0 :       err_printf("MPQS: found %ld / %ld required relations, time = %ld ms\n",
    1473             :                  frel.nb, H.target_rels, timer_delay(&T));
    1474           0 :       dbg_target += H.target_rels / 100.;
    1475             :     }
    1476       66426 :     if (frel.nb < (ulong)H.target_rels) continue; /* main loop */
    1477             : 
    1478         335 :     if (DEBUGLEVEL >= 4)
    1479           0 :       err_printf("\nMPQS: starting Gauss over F_2 on %ld distinct relations\n",
    1480             :                  frel.nb);
    1481         335 :     fact = mpqs_solve_linear_system(&H, &frel);
    1482         335 :     if (fact)
    1483             :     { /* solution found */
    1484         335 :       if (DEBUGLEVEL >= 4)
    1485             :       {
    1486           0 :         err_printf("\nMPQS: time in Gauss and gcds = %ld ms\n",timer_delay(&T));
    1487           0 :         if (typ(fact) == t_INT) err_printf("MPQS: found factor = %Ps\n", fact);
    1488             :         else
    1489             :         {
    1490           0 :           long j, nf = (lg(fact)-1)/3;
    1491           0 :           err_printf("MPQS: found %ld factors =\n", nf);
    1492           0 :           for (j = 1; j <= nf; j++)
    1493           0 :             err_printf("\t%Ps%s\n", gel(fact,3*j-2), (j < nf)? ",": "");
    1494             :         }
    1495             :       }
    1496         335 :       return gerepileupto(av, fact);
    1497             :     }
    1498           0 :     if (DEBUGLEVEL >= 4)
    1499             :     {
    1500           0 :       err_printf("\nMPQS: time in Gauss and gcds = %ld ms\n",timer_delay(&T));
    1501           0 :       err_printf("MPQS: no factors found.\n");
    1502           0 :       if (frel.nb < DEFEAT)
    1503           0 :         err_printf("\nMPQS: restarting sieving ...\n");
    1504             :       else
    1505           0 :         err_printf("\nMPQS: giving up.\n");
    1506             :     }
    1507           0 :     if (frel.nb >= DEFEAT) return gc_NULL(av);
    1508           0 :     H.target_rels += 10;
    1509             :   }
    1510             : }

Generated by: LCOV version 1.13