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

Generated by: LCOV version 1.13