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 25819-e703fe1174) Lines: 595 718 82.9 %
Date: 2020-09-18 06:10:04 Functions: 31 32 96.9 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13