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 12720 : static GEN rel_q(GEN c) { return gel(c,3); }
97 25440 : static GEN rel_Y(GEN c) { return gel(c,1); }
98 25440 : static GEN rel_p(GEN c) { return gel(c,2); }
99 :
100 : static void
101 382752 : frel_add(hashtable *frel, GEN R)
102 : {
103 382752 : ulong h = hash_GEN(R);
104 382752 : if (!hash_search2(frel, (void*)R, h))
105 382745 : hash_insert2(frel, (void*)R, (void*)1, h);
106 382752 : }
107 :
108 : /*********************************************************************/
109 : /** INITIAL SIZING **/
110 : /*********************************************************************/
111 : /* # of decimal digits of argument */
112 : static long
113 9877 : decimal_len(GEN N)
114 9877 : { 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 3290 : mpqs_set_parameters(mpqs_handle_t *h)
121 : {
122 : long s, D;
123 : const mpqs_parameterset_t *P;
124 :
125 3290 : h->digit_size_kN = D = decimal_len(h->kN);
126 3290 : if (D > MPQS_MAX_DIGIT_SIZE_KN) return 0;
127 3290 : P = &(mpqs_parameters[maxss(0, D - 9)]);
128 3290 : h->tolerance = P->tolerance;
129 3290 : h->lp_scale = P->lp_scale;
130 : /* make room for prime factors of k if any: */
131 3290 : 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 3290 : h->target_rels = (s >= 200 ? s + 10 : (mpqs_int32_t)(s * 1.05));
135 3290 : h->M = P->M;
136 3290 : h->omega_A = P->omega_A;
137 3290 : h->no_B = 1UL << (P->omega_A - 1);
138 3290 : h->pmin_index1 = P->pmin_index1;
139 : /* certain subscripts into h->FB should also be offset by omega_k: */
140 3290 : h->index0_FB = 3 + h->_k->omega_k;
141 3290 : 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 3290 : 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 3290 : 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 3290 : 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 3290 : long size_IAH_chunk = (h->size_of_FB + 2) * sizeof(mpqs_inv_A_H_t);
168 3290 : char *fbp = (char*)stack_malloc(size_FB_chunk + 64);
169 3290 : char *iahp = (char*)stack_malloc(size_IAH_chunk + 64);
170 : long fbl, iahl;
171 :
172 3290 : h->FB_chunk = (void *)fbp;
173 3290 : h->invAH_chunk = (void *)iahp;
174 : /* round up to next higher 64-bytes-aligned address */
175 3290 : fbl = (((long)fbp) + 64) & ~0x3FL;
176 : /* and put the actual array there */
177 3290 : h->FB = (mpqs_FB_entry_t *)fbl;
178 :
179 3290 : iahl = (((long)iahp) + 64) & ~0x3FL;
180 3290 : h->inv_A_H = (mpqs_inv_A_H_t *)iahl;
181 3290 : 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 3290 : mpqs_sieve_array_ctor(mpqs_handle_t *h)
188 : {
189 3290 : long size = (h->M << 1) + 1;
190 3290 : mpqs_int32_t size_of_FB = h->size_of_FB;
191 :
192 3290 : h->sieve_array = (unsigned char *) stack_calloc_align(size, sizeof(mpqs_mask));
193 3290 : h->sieve_array_end = h->sieve_array + size - 2;
194 3290 : h->sieve_array_end[1] = 255; /* sentinel */
195 3290 : 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 3290 : 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 3290 : h->relaprimes = (long *) stack_malloc((size_of_FB << 1) * sizeof(long));
201 3290 : }
202 :
203 : /* allocate GENs for current polynomial and self-initialization scratch data */
204 : static void
205 3290 : mpqs_poly_ctor(mpqs_handle_t *h)
206 : {
207 3290 : mpqs_int32_t i, w = h->omega_A;
208 3290 : h->per_A_pr = (mpqs_per_A_prime_t *)
209 3290 : 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 3290 : h->A = cgeti(w + 2);
214 3290 : h->B = cgeti(w + 3);
215 14016 : for (i = 0; i < w; i++) h->per_A_pr[i]._H = cgeti(w + 2);
216 3290 : }
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 3290 : mpqs_find_k(mpqs_handle_t *h)
231 : {
232 3290 : const pari_sp av = avma;
233 3290 : const long N_mod_8 = mod8(h->N), N_mod_4 = N_mod_8 & 3;
234 3290 : long dl = decimal_len(h->N);
235 3290 : long D = maxss(0, minss(dl,MPQS_MAX_DIGIT_SIZE_KN)-9);
236 3290 : 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 3290 : ulong MPQS_NB_MULTIPLIERS = dl < 40 ? 5 : MPQS_POSSIBLE_MULTIPLIERS;
244 : ulong p, i, nbk;
245 :
246 32773 : for (i = nbk = 0; i < numberof(cand_multipliers); i++)
247 : {
248 32773 : const mpqs_multiplier_t *cand_k = &cand_multipliers[i];
249 32773 : long k = cand_k->k;
250 : double v;
251 32773 : if ((k & 3) != N_mod_4) continue; /* want kN = 1 (mod 4) */
252 17120 : v = -log((double)k)/2;
253 17120 : if ((k & 7) == N_mod_8) v += M_LN2; /* kN = 1 (mod 8) */
254 17120 : cache[nbk].np = 0;
255 17120 : cache[nbk]._k = cand_k;
256 17120 : cache[nbk].value = v;
257 17120 : 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 3290 : if (nbk > MPQS_POSSIBLE_MULTIPLIERS) nbk = MPQS_POSSIBLE_MULTIPLIERS;
262 3290 : u_forprime_init(&S, 2, ULONG_MAX);
263 717982 : while ( (p = u_forprime_next(&S)) )
264 : {
265 717982 : long kroNp = kroiu(h->N, p), seen = 0;
266 717982 : if (!kroNp) return p;
267 5940312 : for (i = 0; i < nbk; i++)
268 : {
269 : long krokp;
270 5222330 : if (cache[i].np > MPQS_MULTIPLIER_SEARCH_DEPTH) continue;
271 4886983 : seen++;
272 4886983 : krokp = krouu(cache[i]._k->k % p, p);
273 4886983 : if (krokp == kroNp) /* kronecker(k*N, p)=1 */
274 : {
275 2435210 : cache[i].value += 2*log((double) p)/p;
276 2435210 : cache[i].np++;
277 2451773 : } else if (krokp == 0)
278 : {
279 18880 : cache[i].value += log((double) p)/p;
280 18880 : cache[i].np++;
281 : }
282 : }
283 717982 : if (!seen) break; /* we're gone through SEARCH_DEPTH primes for all k */
284 : }
285 3290 : if (!p) pari_err_OVERFLOW("mpqs_find_k [ran out of primes]");
286 : {
287 3290 : long best_i = 0;
288 3290 : double v = cache[0].value;
289 17120 : for (i = 1; i < nbk; i++)
290 13830 : if (cache[i].value > v) { best_i = i; v = cache[i].value; }
291 3290 : 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 3290 : mpqs_create_FB(mpqs_handle_t *h, ulong *f)
321 : {
322 3290 : mpqs_FB_entry_t *FB = mpqs_FB_ctor(h);
323 3290 : const pari_sp av = avma;
324 3290 : mpqs_int32_t size = h->size_of_FB;
325 : long i;
326 3290 : mpqs_uint32_t k = h->_k->k;
327 : forprime_t S;
328 :
329 3290 : h->largest_FB_p = 0; /* -Wall */
330 3290 : FB[2].fbe_p = 2;
331 : /* the fbe_logval and the fbe_sqrt_kN for 2 are never used */
332 3290 : FB[2].fbe_flags = MPQS_FBE_CLEAR;
333 5641 : for (i = 3; i < h->index0_FB; i++)
334 : { /* this loop executes h->_k->omega_k = 0, 1, or 2 times */
335 2351 : mpqs_uint32_t kp = (ulong)h->_k->kp[i-3];
336 2351 : if (MPQS_DEBUGLEVEL >= 7) err_printf(",<%lu>", (ulong)kp);
337 2351 : FB[i].fbe_p = kp;
338 : /* we could flag divisors of k here, but no need so far */
339 2351 : FB[i].fbe_flags = MPQS_FBE_CLEAR;
340 2351 : FB[i].fbe_flogp = (float)log2((double) kp);
341 2351 : FB[i].fbe_sqrt_kN = 0;
342 : }
343 3290 : (void)u_forprime_init(&S, 3, ULONG_MAX);
344 693203 : while (i < size + 2)
345 : {
346 689913 : ulong p = u_forprime_next(&S);
347 689913 : if (p > k || k % p)
348 : {
349 687562 : ulong kNp = umodiu(h->kN, p);
350 687562 : long kr = krouu(kNp, p);
351 687562 : if (kr != -1)
352 : {
353 351834 : if (kr == 0) { *f = p; return FB; }
354 351834 : FB[i].fbe_p = (mpqs_uint32_t) p;
355 351834 : FB[i].fbe_flags = MPQS_FBE_CLEAR;
356 : /* dyadic logarithm of p; single precision suffices */
357 351834 : 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 351834 : FB[i++].fbe_sqrt_kN = (mpqs_uint32_t)Fl_sqrt(kNp, p);
363 : }
364 : }
365 : }
366 3290 : set_avma(av);
367 3290 : 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 3290 : FB[i].fbe_p = 0; /* sentinel */
376 3290 : 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 6099 : for (i = h->index0_FB; FB[i].fbe_p != 0; i++)
380 6099 : if (FB[i].fbe_p >= h->pmin_index1) break;
381 3290 : h->index1_FB = i;
382 : /* with our parameters this will never fall off the end of the FB */
383 3290 : *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 3290 : mpqs_set_sieve_threshold(mpqs_handle_t *h)
418 : {
419 3290 : mpqs_FB_entry_t *FB = h->FB;
420 : double log_maxval, log_multiplier;
421 : long i;
422 :
423 3290 : h->l2sqrtkN = 0.5 * log2(h->dkN);
424 3290 : h->l2M = log2((double)h->M);
425 3290 : log_maxval = h->l2sqrtkN + h->l2M - MPQS_A_FUDGE;
426 3290 : log_multiplier = 232.0 / log_maxval;
427 3290 : h->sieve_threshold = (unsigned char) (log_multiplier *
428 3290 : (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 3290 : 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 3290 : 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 3290 : if (DEBUGLEVEL >= 5)
445 0 : err_printf("MPQS: computing logarithm approximations for p_i in FB\n");
446 355124 : for (i = h->index0_FB; i < h->size_of_FB + 2; i++)
447 351834 : FB[i].fbe_logval = (unsigned char) (log_multiplier * FB[i].fbe_flogp);
448 3290 : }
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 3290 : 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 3290 : long i = h->index0_FB + 2*(h->omega_A) - 4;
470 : double l2_target_pA;
471 3290 : mpqs_FB_entry_t *FB = h->FB;
472 :
473 3290 : h->l2_target_A = (h->l2sqrtkN - h->l2M - MPQS_A_FUDGE);
474 3290 : l2_target_pA = h->l2_target_A / h->omega_A;
475 :
476 : /* find the sweet spot, normally shouldn't take long */
477 40731 : 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 3290 : 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 3290 : 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 8763 : mpqs_increment(mpqs_uint32_t *x)
531 : {
532 8763 : mpqs_uint32_t r1_mask, r01_mask, slider=1UL;
533 :
534 8763 : switch (*x & 0x1F)
535 : { /* 32-way computed jump handles 22 out of 32 cases */
536 106 : case 29:
537 106 : (*x)++; break; /* shifts a single bit, but we postprocess this case */
538 0 : case 26:
539 0 : (*x) += 2; break; /* again */
540 5967 : case 1: case 3: case 6: case 9: case 11:
541 : case 17: case 19: case 22: case 25: case 27:
542 5967 : (*x) += 3; return;
543 70 : case 20:
544 70 : (*x) += 4; break; /* again */
545 238 : case 5: case 12: case 14: case 21:
546 238 : (*x) += 5; return;
547 1482 : case 2: case 7: case 13: case 18: case 23:
548 1482 : (*x) += 6; return;
549 0 : case 10:
550 0 : (*x) += 7; return;
551 0 : case 8:
552 0 : (*x) += 8; break; /* and again */
553 268 : case 4: case 15:
554 268 : (*x) += 12; return;
555 632 : default: /* 0, 16, 24, 28, 30, 31 */
556 : /* isolate rightmost 1 */
557 632 : r1_mask = ((*x ^ (*x - 1)) + 1) >> 1;
558 : /* isolate rightmost 1 which has a 0 to its left */
559 632 : 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 632 : if (r1_mask == r01_mask) { *x += r1_mask; break; }
563 632 : 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 501 : if (r1_mask == 2) { *x += (r01_mask>>1) + 1; return; }
567 741 : while (r01_mask > r1_mask && slider < r1_mask)
568 : {
569 494 : r01_mask >>= 1; slider <<= 1;
570 : }
571 247 : *x += r01_mask + slider - 1;
572 247 : return;
573 : }
574 : /* post-process cases which couldn't be finalized above */
575 307 : r1_mask = ((*x ^ (*x - 1)) + 1) >> 1;
576 307 : r01_mask = ((*x ^ (*x + r1_mask)) + r1_mask) >> 2;
577 307 : if (r1_mask == r01_mask) { *x += r1_mask; return; }
578 307 : if (r1_mask == 1) { *x += r01_mask; return; }
579 176 : if (r1_mask == 2) { *x += (r01_mask>>1) + 1; return; }
580 140 : while (r01_mask > r1_mask && slider < r1_mask)
581 : {
582 70 : r01_mask >>= 1; slider <<= 1;
583 : }
584 70 : *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 12123 : mpqs_si_choose_primes(mpqs_handle_t *h)
602 : {
603 12123 : mpqs_FB_entry_t *FB = h->FB;
604 12123 : mpqs_per_A_prime_t *per_A_pr = h->per_A_pr;
605 12123 : double l2_last_p = h->l2_target_A;
606 12123 : mpqs_int32_t omega_A = h->omega_A;
607 : int i, j, v2, prev_last_p_idx;
608 12123 : 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 12123 : 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 3360 : h->bin_index = (1UL << (omega_A - 1)) - 1;
629 3360 : prev_last_p_idx = 0;
630 : }
631 : else
632 : { /* clear out old flags */
633 44359 : for (i = 0; i < omega_A; i++) MPQS_FLG(i) = MPQS_FBE_CLEAR;
634 8763 : prev_last_p_idx = MPQS_I(omega_A-1);
635 :
636 8763 : if (room > 30) room = 30;
637 8763 : 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 8763 : mpqs_increment(&h->bin_index);
646 8763 : if (h->index2_moved)
647 : {
648 38 : 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 8763 : if ((h->bin_index & room_mask) != 0)
653 : { /* Yes. Turn on the index2_moved flag in the handle */
654 70 : h->index2_FB += 2; /* caller to check this isn't too large!!! */
655 70 : h->index2_moved = 1;
656 70 : h->bin_index = 0;
657 70 : 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 70 : 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 12053 : bits = h->bin_index;
667 12053 : 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 12053 : j = h->index2_FB;
673 12053 : v2 = vals((long)bits);
674 12053 : if (v2) { j -= v2; bits >>= v2; }
675 34269 : for (i = omega_A - 2; i >= 0; i--)
676 : {
677 34269 : MPQS_I(i) = j;
678 34269 : l2_last_p -= MPQS_LP(i);
679 34269 : MPQS_FLG(i) |= MPQS_FBE_DIVIDES_A;
680 34269 : bits &= ~1UL;
681 34269 : if (!bits) break; /* i = 0 */
682 22216 : v2 = vals((long)bits); /* > 0 */
683 22216 : bits >>= v2; j -= v2;
684 : }
685 : /* Choose the larger prime. Note we keep index2_FB <= size_of_FB - 3 */
686 150360 : for (j = h->index2_FB + 1; (p = FB[j].fbe_p); j++)
687 150255 : 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 12053 : if (p && j == prev_last_p_idx) { j++; p = FB[j].fbe_p; }
695 12053 : MPQS_I(omega_A - 1) = p? j: h->size_of_FB + 1;
696 12053 : MPQS_FLG(omega_A - 1) |= MPQS_FBE_DIVIDES_A;
697 :
698 12053 : 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 12053 : 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 143546 : mpqs_self_init(mpqs_handle_t *h)
735 : {
736 143546 : const ulong size_of_FB = h->size_of_FB + 1;
737 143546 : mpqs_FB_entry_t *FB = h->FB;
738 143546 : mpqs_inv_A_H_t *inv_A_H = h->inv_A_H;
739 143546 : const pari_sp av = avma;
740 143546 : GEN p1, A = h->A, B = h->B, mC;
741 143546 : 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 143546 : 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 8763 : h->index_j = 0;
751 8763 : h->index_i++; /* count finished A's */
752 : }
753 :
754 143546 : if (h->index_j == 0)
755 : { /* compute first polynomial with new A */
756 : GEN a, b, A2;
757 12053 : if (!mpqs_si_choose_primes(h))
758 : { /* Ran out of room towards small primes, and index2_FB was raised. */
759 70 : if (size_of_FB - h->index2_FB < 4) return 0; /* Fail */
760 70 : (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 12053 : a = b = NULL;
766 58375 : for (i = 0; i < h->omega_A; i++)
767 : {
768 46322 : ulong p = MPQS_AP(i);
769 46322 : a = a? muliu(a, p): utoipos(p);
770 : }
771 12053 : 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 58375 : for (i = 0; i < h->omega_A; i++)
777 : {
778 46322 : ulong p = MPQS_AP(i);
779 46322 : GEN t = divis(A, (long)p);
780 46322 : t = remii(mulii(t, muluu(Fl_inv(umodiu(t, p), p), MPQS_SQRT(i))), A);
781 46322 : affii(t, MPQS_H(i));
782 46322 : b = b? addii(b, t): t;
783 : }
784 : /* ensure b = 1 mod 4 */
785 12053 : if (mod2(b) == 0)
786 5821 : b = addii(b, mului(mod4(A), A)); /* b += (A % 4) * A; */
787 :
788 12053 : affii(b, B); set_avma(av);
789 :
790 12053 : 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 4140297 : for (j = 3; (ulong)j <= size_of_FB; j++)
798 : {
799 : ulong s, mb, t, m, p, iA2, iA;
800 4128244 : if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
801 4081922 : p = (ulong)FB[j].fbe_p;
802 4081922 : m = h->M % p;
803 4081922 : iA2 = Fl_inv(umodiu(A2, p), p); /* = 1/(2*A) mod p_j */
804 4081922 : iA = iA2 << 1; if (iA > p) iA -= p;
805 4081922 : mb = umodiu(B, p); if (mb) mb = p - mb; /* mb = -B mod p */
806 4081922 : s = FB[j].fbe_sqrt_kN;
807 4081922 : t = Fl_add(m, Fl_mul(Fl_sub(mb, s, p), iA2, p), p);
808 4081922 : FB[j].fbe_start1 = (mpqs_int32_t)t;
809 4081922 : FB[j].fbe_start2 = (mpqs_int32_t)Fl_add(t, Fl_mul(s, iA, p), p);
810 24685309 : for (i = 0; i < h->omega_A - 1; i++)
811 : {
812 20603387 : ulong h = umodiu(MPQS_H(i), p);
813 20603387 : 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 131493 : ulong p, v2 = vals(h->index_j); /* new starting positions for sieving */
832 131493 : j = h->index_j >> v2;
833 131493 : p1 = shifti(MPQS_H(v2), 1);
834 131493 : if (j & 2)
835 : { /* j = 3 mod 4 */
836 86607851 : for (j = 3; (ulong)j <= size_of_FB; j++)
837 : {
838 86548337 : if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
839 86190049 : p = (ulong)FB[j].fbe_p;
840 86190049 : FB[j].fbe_start1 = Fl_sub(FB[j].fbe_start1, MPQS_INV_A_H(v2,j), p);
841 86190049 : FB[j].fbe_start2 = Fl_sub(FB[j].fbe_start2, MPQS_INV_A_H(v2,j), p);
842 : }
843 59514 : p1 = addii(B, p1);
844 : }
845 : else
846 : { /* j = 1 mod 4 */
847 90950479 : for (j = 3; (ulong)j <= size_of_FB; j++)
848 : {
849 90878500 : if (FB[j].fbe_flags & MPQS_FBE_DIVIDES_A) continue;
850 90471784 : p = (ulong)FB[j].fbe_p;
851 90471784 : FB[j].fbe_start1 = Fl_add(FB[j].fbe_start1, MPQS_INV_A_H(v2,j), p);
852 90471784 : FB[j].fbe_start2 = Fl_add(FB[j].fbe_start2, MPQS_INV_A_H(v2,j), p);
853 : }
854 71979 : p1 = subii(B, p1);
855 : }
856 131493 : 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 143546 : mC = diviiexact(subii(h->kN, sqri(B)), shifti(A, 2)); /* coefficient -C */
864 954872 : for (i = 0; i < h->omega_A; i++)
865 : {
866 811326 : ulong p = MPQS_AP(i), s = h->M + Fl_div(umodiu(mC, p), umodiu(B, p), p);
867 811326 : 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, mC, FB[j].fbe_p, FB[j].fbe_start1);
873 : check_root(h, mC, FB[j].fbe_p, FB[j].fbe_start2);
874 : }
875 : #endif
876 143546 : 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 143546 : set_avma(av);
881 : #ifdef MPQS_DEBUG
882 : err_printf("MPQS DEBUG: leave self init, avma = 0x%lX\n", (ulong)avma);
883 : #endif
884 143546 : 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 810897 : mpqs_sieve_p(unsigned char *B, unsigned char *E, long p4, long p,
893 : unsigned char logp)
894 : {
895 810897 : 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 19439685 : while (e - B >= 0) /* signed comparison */
900 : {
901 18628788 : (*B) += logp, B += p;
902 18628788 : (*B) += logp, B += p;
903 18628788 : (*B) += logp, B += p;
904 18628788 : (*B) += logp, B += p;
905 : }
906 2766326 : while (E - B >= 0) (*B) += logp, B += p;
907 810897 : }
908 :
909 : INLINE void
910 126914333 : mpqs_sieve_p1(unsigned char *B, unsigned char *E, long s1, long s2,
911 : unsigned char logp)
912 : {
913 569033504 : while (E - B >= 0)
914 : {
915 482373242 : (*B) += logp, B += s1;
916 482373242 : if (E - B < 0) break;
917 442119171 : (*B) += logp, B += s2;
918 : }
919 126914333 : }
920 :
921 : INLINE void
922 53085237 : mpqs_sieve_p2(unsigned char *B, unsigned char *E, long p4, long s1, long s2,
923 : unsigned char logp)
924 : {
925 53085237 : 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 794093164 : while (e - B >= 0) /* signed comparison */
930 : {
931 741007927 : (*B) += logp, B += s1;
932 741007927 : (*B) += logp, B += s2;
933 741007927 : (*B) += logp, B += s1;
934 741007927 : (*B) += logp, B += s2;
935 741007927 : (*B) += logp, B += s1;
936 741007927 : (*B) += logp, B += s2;
937 741007927 : (*B) += logp, B += s1;
938 741007927 : (*B) += logp, B += s2;
939 : }
940 162969087 : while (E - B >= 0) {(*B) += logp, B += s1; if (E - B < 0) break; (*B) += logp, B += s2;}
941 53085237 : }
942 : static void
943 143546 : mpqs_sieve(mpqs_handle_t *h)
944 : {
945 143546 : long p, l = h->index1_FB;
946 143546 : mpqs_FB_entry_t *FB = &(h->FB[l]);
947 143546 : unsigned char *S = h->sieve_array, *Send = h->sieve_array_end;
948 143546 : long size = h->M << 1, size4 = size >> 3;
949 143546 : memset((void*)S, 0, size * sizeof(unsigned char));
950 54039680 : for ( ; (p = FB->fbe_p) && p <= size4; FB++) /* l++ */
951 : {
952 53896134 : unsigned char logp = FB->fbe_logval;
953 53896134 : 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 53896134 : if (s1 == s2) mpqs_sieve_p(S + s1, Send, p << 2, p, logp);
956 : else
957 : {
958 53085237 : if (s1>s2) lswap(s1,s2)
959 53085237 : mpqs_sieve_p2(S + s1, Send, p << 2, s2-s1,p+s1-s2, logp);
960 : }
961 : }
962 127057879 : for ( ; (p = FB->fbe_p) && p <= size; FB++) /* l++ */
963 : {
964 126914333 : unsigned char logp = FB->fbe_logval;
965 126914333 : 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 126914333 : if (s1 == s2) mpqs_sieve_p(S + s1, Send, p << 2, p, logp);
968 : else
969 : {
970 126914333 : if (s1>s2) lswap(s1,s2)
971 126914333 : mpqs_sieve_p1(S + s1, Send, s2-s1, p+s1-s2, logp);
972 : }
973 : }
974 143546 : 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 143546 : }
982 :
983 : /* Could use the fact that 4 | M, but let the compiler worry about unrolling. */
984 : static long
985 143546 : mpqs_eval_sieve(mpqs_handle_t *h)
986 : {
987 143546 : long x = 0, count = 0, M2 = h->M << 1;
988 143546 : unsigned char t = h->sieve_threshold;
989 143546 : unsigned char *S = h->sieve_array;
990 143546 : mpqs_bit_array * U = (mpqs_bit_array *) S;
991 143546 : long *cand = h->candidates;
992 143546 : 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 10664729 : while (count < MPQS_CANDIDATE_ARRAY_SIZE - 1)
997 : {
998 : long j, y;
999 506078324 : while (!TEST(U[x]&mpqs_mask)) x++;
1000 10664729 : y = x*sizemask;
1001 166178765 : for (j=0; j<sizemask; j++, y++)
1002 : {
1003 155657582 : if (y >= M2)
1004 143546 : { cand[count] = 0; return count; }
1005 155514036 : if (S[y]>=t)
1006 528706 : cand[count++] = y;
1007 : }
1008 10521183 : 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 366585 : rel_to_ei(GEN ei, GEN relp)
1071 : {
1072 366585 : long j, l = lg(relp);
1073 5317752 : for (j=1; j<l; j++)
1074 : {
1075 4951167 : long e = relp[j] >> REL_OFFSET, i = relp[j] & REL_MASK;
1076 4951167 : ei[i] += e;
1077 : }
1078 366585 : }
1079 : static void
1080 7558518 : mpqs_add_factor(GEN relp, long *i, ulong ei, ulong pi)
1081 7558518 : { relp[++*i] = pi | (ei << REL_OFFSET); }
1082 :
1083 : static int
1084 12720 : zv_is_even(GEN V)
1085 : {
1086 12720 : long i, l = lg(V);
1087 1098145 : for (i=1; i<l; i++)
1088 1097559 : if (odd(uel(V,i))) return 0;
1089 586 : return 1;
1090 : }
1091 :
1092 : static GEN
1093 12720 : combine_large_primes(mpqs_handle_t *h, GEN rel1, GEN rel2)
1094 : {
1095 12720 : GEN new_Y, new_Y1, Y1 = rel_Y(rel1), Y2 = rel_Y(rel2);
1096 12720 : long l, lei = h->size_of_FB + 1, nb = 0;
1097 12720 : GEN ei, relp, iq, q = rel_q(rel1);
1098 :
1099 12720 : if (!invmod(q, h->N, &iq)) return equalii(iq, h->N)? NULL: iq; /* rare */
1100 12720 : ei = zero_zv(lei);
1101 12720 : rel_to_ei(ei, rel_p(rel1));
1102 12720 : rel_to_ei(ei, rel_p(rel2));
1103 12720 : if (zv_is_even(ei)) return NULL;
1104 12134 : new_Y = modii(mulii(mulii(Y1, Y2), iq), h->N);
1105 12134 : new_Y1 = subii(h->N, new_Y);
1106 12134 : if (abscmpii(new_Y1, new_Y) < 0) new_Y = new_Y1;
1107 12134 : relp = cgetg(MAX_PE_PAIR+1,t_VECSMALL);
1108 12134 : if (odd(ei[1])) mpqs_add_factor(relp, &nb, 1, 1);
1109 20934649 : for (l = 2; l <= lei; l++)
1110 20922515 : if (ei[l]) mpqs_add_factor(relp, &nb, ei[l],l);
1111 12134 : setlg(relp, nb+1);
1112 12134 : 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 12134 : return mkvec2(new_Y, relp);
1132 : }
1133 :
1134 : /* nc candidates */
1135 : static GEN
1136 133863 : mpqs_eval_cand(mpqs_handle_t *h, long nc, hashtable *frel, hashtable *lprel)
1137 : {
1138 133863 : mpqs_FB_entry_t *FB = h->FB;
1139 133863 : GEN A = h->A, B = h->B;
1140 133863 : long *relaprimes = h->relaprimes, *candidates = h->candidates;
1141 : long pi, i;
1142 : int pii;
1143 133863 : mpqs_per_A_prime_t *per_A_pr = h->per_A_pr;
1144 :
1145 662569 : for (i = 0; i < nc; i++)
1146 : {
1147 528706 : pari_sp btop = avma;
1148 528706 : GEN Qx, Qx_part, Y, relp = cgetg(MAX_PE_PAIR+1,t_VECSMALL);
1149 528706 : long powers_of_2, p, x = candidates[i], nb = 0;
1150 528706 : int relaprpos = 0;
1151 : long k;
1152 528706 : unsigned char thr = h->sieve_array[x];
1153 : /* Y = 2*A*x + B, Qx = Y^2/(4*A) = Q(x) */
1154 528706 : Y = addii(mulis(A, 2 * (x - h->M)), B);
1155 528706 : Qx = subii(sqri(Y), h->kN); /* != 0 since N not a square and (N,k) = 1 */
1156 528706 : if (signe(Qx) < 0)
1157 : {
1158 283810 : setabssign(Qx);
1159 283810 : 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 528706 : powers_of_2 = vali(Qx);
1164 528706 : Qx = shifti(Qx, -powers_of_2);
1165 528706 : 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 2673467 : for (Qx_part = A, pi = 3; pi< h->index1_FB; pi++)
1191 : {
1192 2144761 : ulong p = FB[pi].fbe_p;
1193 2144761 : long xp = x % p;
1194 : /* Here we used that MPQS_FBE_DIVIDES_A = 1. */
1195 :
1196 2144761 : 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 538306 : ulong ei = FB[pi].fbe_flags & MPQS_FBE_DIVIDES_A;
1199 538306 : relaprimes[relaprpos++] = pi;
1200 538306 : relaprimes[relaprpos++] = 1 + ei;
1201 538306 : Qx_part = muliu(Qx_part, p);
1202 : }
1203 : }
1204 314246903 : for ( ; thr && (p = FB[pi].fbe_p); pi++)
1205 : {
1206 313718197 : long xp = x % p;
1207 : /* Here we used that MPQS_FBE_DIVIDES_A = 1. */
1208 :
1209 313718197 : 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 3387360 : ulong ei = FB[pi].fbe_flags & MPQS_FBE_DIVIDES_A;
1212 3387360 : relaprimes[relaprpos++] = pi;
1213 3387360 : relaprimes[relaprpos++] = 1 + ei;
1214 3387360 : Qx_part = muliu(Qx_part, p);
1215 3387360 : thr -= FB[pi].fbe_logval;
1216 : }
1217 : }
1218 3119958 : for (k = 0; k< h->omega_A; k++)
1219 : {
1220 2591252 : long pi = MPQS_I(k);
1221 2591252 : ulong p = FB[pi].fbe_p;
1222 2591252 : long xp = x % p;
1223 2591252 : 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 2514027 : relaprimes[relaprpos++] = pi;
1226 2514027 : 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 528706 : 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 6968399 : for (pii = 0; pii < relaprpos; pii += 2)
1243 : {
1244 6439693 : ulong r, ei = relaprimes[pii+1];
1245 : GEN q;
1246 :
1247 6439693 : pi = relaprimes[pii];
1248 : /* p | k (identified by its index before index0_FB)* or p | A (ei = 0) */
1249 6439693 : if ((mpqs_int32_t)pi < h->index0_FB || ei == 0)
1250 : {
1251 2572161 : mpqs_add_factor(relp, &nb, 1, pi);
1252 2572161 : continue;
1253 : }
1254 3867532 : p = FB[pi].fbe_p;
1255 : /* p might still divide the current adjusted Qx. Try it. */
1256 3867532 : switch(cmpiu(Qx, p))
1257 : {
1258 87803 : case 0: ei++; Qx = gen_1; break;
1259 1251532 : case 1:
1260 1251532 : q = absdiviu_rem(Qx, p, &r);
1261 1338309 : while (r == 0) { ei++; Qx = q; q = absdiviu_rem(Qx, p, &r); }
1262 1251532 : break;
1263 : }
1264 3867532 : 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 528706 : setlg(relp, nb+1);
1272 528706 : if (is_pm1(Qx))
1273 : {
1274 370618 : GEN rel = gerepilecopy(btop, mkvec2(absi_shallow(Y), relp));
1275 : #ifdef MPQS_DEBUG
1276 : mpqs_check_rel(h, rel);
1277 : #endif
1278 370618 : frel_add(frel, rel);
1279 : }
1280 158088 : else if (cmpiu(Qx, h->lp_bound) <= 0)
1281 : {
1282 145974 : ulong q = itou(Qx);
1283 145974 : GEN rel = mkvec3(absi_shallow(Y),relp,Qx);
1284 145974 : GEN col = hash_haskey_GEN(lprel, (void*)q);
1285 : #ifdef MPQS_DEBUG
1286 : mpqs_check_rel(h, rel);
1287 : #endif
1288 145974 : if (!col) /* relation up to large prime */
1289 133254 : hash_insert(lprel, (void*)q, (void*)gerepilecopy(btop,rel));
1290 12720 : else if ((rel = combine_large_primes(h, rel, col)))
1291 : {
1292 12134 : if (typ(rel) == t_INT) return rel; /* very unlikely */
1293 : #ifdef MPQS_DEBUG
1294 : mpqs_check_rel(h, rel);
1295 : #endif
1296 12134 : 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 12114 : set_avma(btop);
1305 : }
1306 : }
1307 : PRINT_IF_VERBOSE("\n");
1308 133863 : return NULL;
1309 : }
1310 :
1311 : /*********************************************************************/
1312 : /** FROM RELATIONS TO DIVISORS **/
1313 : /*********************************************************************/
1314 :
1315 : /* create an F2m from a relations list */
1316 : static GEN
1317 3380 : rels_to_F2Ms(GEN rel)
1318 : {
1319 3380 : long i, cols = lg(rel)-1;
1320 3380 : GEN m = cgetg(cols+1, t_VEC);
1321 390355 : for (i = 1; i <= cols; i++)
1322 : {
1323 386975 : GEN relp = gmael(rel,i,2), rel2;
1324 386975 : long j, l = lg(relp), o = 0, k;
1325 5536801 : for (j = 1; j < l; j++)
1326 5149826 : if (odd(relp[j] >> REL_OFFSET)) o++;
1327 386975 : rel2 = cgetg(o+1, t_VECSMALL);
1328 5536801 : for (j = 1, k = 1; j < l; j++)
1329 5149826 : if (odd(relp[j] >> REL_OFFSET))
1330 4738390 : rel2[k++] = relp[j] & REL_MASK;
1331 386975 : gel(m, i) = rel2;
1332 : }
1333 3380 : return m;
1334 : }
1335 :
1336 : static int
1337 7282 : split(GEN *D, long *e)
1338 : {
1339 : ulong mask;
1340 : long flag;
1341 7282 : if (MR_Jaeschke(*D)) { *e = 1; return 1; } /* probable prime */
1342 547 : if (Z_issquareall(*D, D))
1343 : { /* squares could cost us a lot of time */
1344 168 : if (DEBUGLEVEL >= 5) err_printf("MPQS: decomposed a square\n");
1345 168 : *e = 2; return 1;
1346 : }
1347 379 : 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 379 : 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 365 : *e = 0; return 0; /* known composite */
1357 : }
1358 :
1359 : /* return a GEN structure containing NULL but safe for gerepileupto */
1360 : static GEN
1361 3380 : mpqs_solve_linear_system(mpqs_handle_t *h, hashtable *frel)
1362 : {
1363 3380 : mpqs_FB_entry_t *FB = h->FB;
1364 3380 : pari_sp av = avma;
1365 3380 : GEN rels = hash_keys_GEN(frel), N = h->N, r, c, res, ei, M, Ker;
1366 : long i, j, nrows, rlast, rnext, rmax, rank;
1367 :
1368 3380 : M = rels_to_F2Ms(rels);
1369 3380 : Ker = F2Ms_ker(M, h->size_of_FB+1); rank = lg(Ker)-1;
1370 3380 : 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 3380 : 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 3380 : ei = cgetg(h->size_of_FB + 2, t_VECSMALL);
1387 3380 : rmax = logint(N, utoi(3));
1388 3380 : if (rank <= BITS_IN_LONG-2)
1389 3357 : rmax = minss(rmax, 1L<<rank); /* max # of factors we can hope for */
1390 3380 : r = cgetg(rmax+1, t_VEC);
1391 3380 : c = cgetg(rmax+1, t_VECSMALL);
1392 3380 : rnext = rlast = 1;
1393 3380 : nrows = lg(M)-1;
1394 9020 : for (i = 1; i <= rank; i++)
1395 : { /* loop over kernel basis */
1396 8916 : GEN X = gen_1, Y_prod = gen_1, X_plus_Y, D;
1397 8916 : pari_sp av2 = avma, av3;
1398 8916 : long done = 0; /* # probably-prime factors or powers whose bases we won't
1399 : * handle any further */
1400 8916 : memset((void *)(ei+1), 0, (h->size_of_FB + 1) * sizeof(long));
1401 1003088 : for (j = 1; j <= nrows; j++)
1402 994172 : if (F2m_coeff(Ker, j, i))
1403 : {
1404 341145 : GEN R = gel(rels,j);
1405 341145 : Y_prod = gerepileuptoint(av2, remii(mulii(Y_prod, gel(R,1)), N));
1406 341145 : rel_to_ei(ei, gel(R,2));
1407 : }
1408 8916 : av3 = avma;
1409 935236 : for (j = 2; j <= h->size_of_FB + 1; j++)
1410 926320 : if (ei[j])
1411 : {
1412 584603 : GEN q = utoipos(FB[j].fbe_p);
1413 584603 : if (ei[j] & 1) pari_err_BUG("MPQS (relation is a nonsquare)");
1414 584603 : X = remii(mulii(X, Fp_powu(q, (ulong)ei[j]>>1, N)), N);
1415 584603 : X = gerepileuptoint(av3, X);
1416 : }
1417 8916 : 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 8916 : X_plus_Y = addii(X, Y_prod);
1429 8916 : if (rnext == 1)
1430 : { /* we still haven't decomposed, and want both a gcd and its cofactor. */
1431 8284 : D = gcdii(X_plus_Y, N);
1432 8284 : if (is_pm1(D) || equalii(D,N)) { set_avma(av2); continue; }
1433 : /* got something that works */
1434 3290 : if (DEBUGLEVEL >= 5)
1435 0 : err_printf("MPQS: splitting N after %ld kernel vector%s\n",
1436 : i+1, (i? "s" : ""));
1437 3290 : gel(r,1) = diviiexact(N, D);
1438 3290 : gel(r,2) = D;
1439 3290 : rlast = rnext = 3;
1440 3290 : if (split(&gel(r,1), &c[1])) done++;
1441 3290 : if (split(&gel(r,2), &c[2])) done++;
1442 3290 : if (done == 2 || rmax == 2) break;
1443 365 : if (DEBUGLEVEL >= 5)
1444 0 : err_printf("MPQS: got two factors, looking for more...\n");
1445 : }
1446 : else
1447 : { /* we already have factors */
1448 2247 : 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 1615 : if (c[j]) { done++; continue; }
1453 632 : av3 = avma; D = gcdii(X_plus_Y, gel(r,j));
1454 632 : if (is_pm1(D) || equalii(D, gel(r,j))) { set_avma(av3); continue; }
1455 : /* got one which splits this factor */
1456 351 : if (DEBUGLEVEL >= 5)
1457 0 : err_printf("MPQS: resplitting a factor after %ld kernel vectors\n",
1458 : i+1);
1459 351 : gel(r,j) = diviiexact(gel(r,j), D);
1460 351 : gel(r,rnext) = D;
1461 351 : 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 351 : (void)split(&gel(r,rnext), &c[rnext]);
1465 351 : if (++rnext > rmax) break; /* all possible factors seen */
1466 : } /* loop over known composite factors */
1467 :
1468 632 : if (rnext > rlast)
1469 : {
1470 351 : if (DEBUGLEVEL >= 5)
1471 0 : err_printf("MPQS: got %ld factors%s\n", rlast - 1,
1472 : (done < rlast ? ", looking for more..." : ""));
1473 351 : 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 632 : if (rnext > rmax || done == rnext - 1) break;
1478 : }
1479 : }
1480 3380 : 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 3290 : rlast = rnext - 1; /* # of distinct factors found */
1485 3290 : res = cgetg(3*rlast + 1, t_VEC);
1486 3290 : if (DEBUGLEVEL >= 6) err_printf("MPQS: wrapping up %ld factors\n", rlast);
1487 10221 : for (i = j = 1; i <= rlast; i++, j += 3)
1488 : {
1489 6931 : long C = c[i];
1490 6931 : icopyifstack(gel(r,i), gel(res,j)); /* factor */
1491 6931 : gel(res,j+1) = C <= 1? gen_1: utoipos(C); /* exponent */
1492 6931 : gel(res,j+2) = C ? NULL: gen_0; /* unknown or known composite */
1493 6931 : 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 3290 : 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 3297 : mpqs(GEN N)
1513 : {
1514 3297 : 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 3297 : pari_sp av = avma;
1523 :
1524 3297 : if (DEBUGLEVEL >= 4) err_printf("MPQS: number to factor N = %Ps\n", N);
1525 3297 : if (size_N > MPQS_MAX_DIGIT_SIZE_KN) { toolarge(); return NULL; }
1526 3290 : 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 3290 : H.N = N;
1532 3290 : H.bin_index = 0;
1533 3290 : H.index_i = 0;
1534 3290 : H.index_j = 0;
1535 3290 : H.index2_moved = 0;
1536 3290 : p = mpqs_find_k(&H);
1537 3290 : if (p) return gc_utoipos(av,p);
1538 3290 : if (DEBUGLEVEL >= 5)
1539 0 : err_printf("MPQS: found multiplier %ld for N\n", H._k->k);
1540 3290 : H.kN = muliu(N, H._k->k);
1541 3290 : if (!mpqs_set_parameters(&H)) { toolarge(); return NULL; }
1542 :
1543 3290 : if (DEBUGLEVEL >= 5)
1544 0 : err_printf("MPQS: creating factor base and allocating arrays...\n");
1545 3290 : FB = mpqs_create_FB(&H, &p);
1546 3290 : if (p) return gc_utoipos(av, p);
1547 3290 : mpqs_sieve_array_ctor(&H);
1548 3290 : mpqs_poly_ctor(&H);
1549 :
1550 3290 : 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 3290 : H.lp_bound *= minss(H.lp_scale, H.largest_FB_p - 1);
1554 3290 : H.dkN = gtodouble(H.kN);
1555 : /* compute the threshold and fill in the byte-sized scaled logarithms */
1556 3290 : mpqs_set_sieve_threshold(&H);
1557 3290 : if (!mpqs_locate_A_range(&H)) return NULL;
1558 3290 : 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 3290 : H.index_j = (mpqs_uint32_t)-1; /* increment below will have it start at 0 */
1584 :
1585 3290 : dbg_target = H.target_rels / 100.;
1586 3290 : DEFEAT = H.target_rels * 1.5;
1587 3290 : hash_init_GEN(&frel, H.target_rels, gequal, 1);
1588 3290 : hash_init_ulong(&lprel,H.target_rels, 1);
1589 : for(;;)
1590 140256 : {
1591 : long tc;
1592 : /* self initialization: compute polynomial and its zeros */
1593 143546 : 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 143546 : mpqs_sieve(&H);
1600 143546 : tc = mpqs_eval_sieve(&H);
1601 143546 : if (DEBUGLEVEL >= 6)
1602 0 : err_printf("MPQS: found %lu candidate%s\n", tc, (tc==1? "" : "s"));
1603 143546 : if (tc)
1604 : {
1605 133863 : fact = mpqs_eval_cand(&H, tc, &frel, &lprel);
1606 133863 : 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 143546 : 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 143546 : if (frel.nb < (ulong)H.target_rels) continue; /* main loop */
1624 :
1625 3380 : 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 3380 : fact = mpqs_solve_linear_system(&H, &frel);
1632 3380 : if (fact)
1633 : { /* solution found */
1634 3290 : 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 3290 : return gerepileupto(av, fact);
1647 : }
1648 90 : 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 90 : if (frel.nb >= DEFEAT) return gc_NULL(av);
1658 90 : H.target_rels += 10;
1659 : }
1660 : }
|
