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 : #include "pari.h"
15 : #include "paripriv.h"
16 :
17 : #define DEBUGLEVEL DEBUGLEVEL_factorint
18 :
19 : /***********************************************************************/
20 : /** PRIMES IN SUCCESSION **/
21 : /***********************************************************************/
22 :
23 : /* map from prime residue classes mod 210 to their numbers in {0...47}.
24 : * Subscripts into this array take the form ((k-1)%210)/2, ranging from
25 : * 0 to 104. Unused entries are */
26 : #define NPRC 128 /* nonprime residue class */
27 :
28 : static unsigned char prc210_no[] = {
29 : 0, NPRC, NPRC, NPRC, NPRC, 1, 2, NPRC, 3, 4, NPRC, /* 21 */
30 : 5, NPRC, NPRC, 6, 7, NPRC, NPRC, 8, NPRC, 9, /* 41 */
31 : 10, NPRC, 11, NPRC, NPRC, 12, NPRC, NPRC, 13, 14, NPRC, /* 63 */
32 : NPRC, 15, NPRC, 16, 17, NPRC, NPRC, 18, NPRC, 19, /* 83 */
33 : NPRC, NPRC, 20, NPRC, NPRC, NPRC, 21, NPRC, 22, 23, NPRC, /* 105 */
34 : 24, 25, NPRC, 26, NPRC, NPRC, NPRC, 27, NPRC, NPRC, /* 125 */
35 : 28, NPRC, 29, NPRC, NPRC, 30, 31, NPRC, 32, NPRC, NPRC, /* 147 */
36 : 33, 34, NPRC, NPRC, 35, NPRC, NPRC, 36, NPRC, 37, /* 167 */
37 : 38, NPRC, 39, NPRC, NPRC, 40, 41, NPRC, NPRC, 42, NPRC, /* 189 */
38 : 43, 44, NPRC, 45, 46, NPRC, NPRC, NPRC, NPRC, 47, /* 209 */
39 : };
40 :
41 : /* first differences of the preceding */
42 : static unsigned char prc210_d1[] = {
43 : 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6,
44 : 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6,
45 : 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2,
46 : };
47 :
48 : static int
49 984000 : unextprime_overflow(ulong n)
50 : {
51 : #ifdef LONG_IS_64BIT
52 983414 : return (n > (ulong)-59);
53 : #else
54 586 : return (n > (ulong)-5);
55 : #endif
56 : }
57 :
58 : /* return 0 for overflow */
59 : ulong
60 1121103 : unextprime(ulong n)
61 : {
62 : long rc, rc0, rcn;
63 :
64 1121103 : switch(n) {
65 6858 : case 0: case 1: case 2: return 2;
66 2434 : case 3: return 3;
67 1668 : case 4: case 5: return 5;
68 1162 : case 6: case 7: return 7;
69 : }
70 1108981 : if (n <= maxprime())
71 : {
72 124964 : long i = PRIMES_search(n);
73 124964 : return i > 0? n: pari_PRIMES[-i];
74 : }
75 984005 : if (unextprime_overflow(n)) return 0;
76 : /* here n > 7 */
77 983971 : n |= 1; /* make it odd */
78 983971 : rc = rc0 = n % 210;
79 : /* find next prime residue class mod 210 */
80 : for(;;)
81 : {
82 2166691 : rcn = (long)(prc210_no[rc>>1]);
83 2166691 : if (rcn != NPRC) break;
84 1182720 : rc += 2; /* cannot wrap since 209 is coprime and rc odd */
85 : }
86 983971 : if (rc > rc0) n += rc - rc0;
87 : /* now find an actual (pseudo)prime */
88 : for(;;)
89 : {
90 11530693 : if (uisprime(n)) break;
91 10546722 : n += prc210_d1[rcn];
92 10546722 : if (++rcn > 47) rcn = 0;
93 : }
94 983998 : return n;
95 : }
96 :
97 : GEN
98 126847 : nextprime(GEN n)
99 : {
100 : long rc, rc0, rcn;
101 126847 : pari_sp av = avma;
102 :
103 126847 : if (typ(n) != t_INT)
104 : {
105 14 : n = gceil(n);
106 14 : if (typ(n) != t_INT) pari_err_TYPE("nextprime",n);
107 : }
108 126840 : if (signe(n) <= 0) { set_avma(av); return gen_2; }
109 126840 : if (lgefint(n) == 3)
110 : {
111 118710 : ulong k = unextprime(uel(n,2));
112 118710 : set_avma(av);
113 118710 : if (k) return utoipos(k);
114 : #ifdef LONG_IS_64BIT
115 6 : return uutoi(1,13);
116 : #else
117 1 : return uutoi(1,15);
118 : #endif
119 : }
120 : /* here n > 7 */
121 8130 : if (!mod2(n)) n = addui(1,n);
122 8130 : rc = rc0 = umodiu(n, 210);
123 : /* find next prime residue class mod 210 */
124 : for(;;)
125 : {
126 17701 : rcn = (long)(prc210_no[rc>>1]);
127 17701 : if (rcn != NPRC) break;
128 9571 : rc += 2; /* cannot wrap since 209 is coprime and rc odd */
129 : }
130 8130 : if (rc > rc0) n = addui(rc - rc0, n);
131 : /* now find an actual (pseudo)prime */
132 : for(;;)
133 : {
134 84466 : if (BPSW_psp(n)) break;
135 76336 : n = addui(prc210_d1[rcn], n);
136 76336 : if (++rcn > 47) rcn = 0;
137 : }
138 8130 : if (avma == av) return icopy(n);
139 8130 : return gerepileuptoint(av, n);
140 : }
141 :
142 : ulong
143 32 : uprecprime(ulong n)
144 : {
145 : long rc, rc0, rcn;
146 : { /* check if n <= 10 */
147 32 : if (n <= 1) return 0;
148 25 : if (n == 2) return 2;
149 18 : if (n <= 4) return 3;
150 18 : if (n <= 6) return 5;
151 18 : if (n <= 10) return 7;
152 : }
153 18 : if (n <= maxprimelim())
154 : {
155 0 : long i = PRIMES_search(n);
156 0 : return i > 0? n: pari_PRIMES[-i-1];
157 : }
158 : /* here n >= 11 */
159 18 : if (!(n % 2)) n--;
160 18 : rc = rc0 = n % 210;
161 : /* find previous prime residue class mod 210 */
162 : for(;;)
163 : {
164 36 : rcn = (long)(prc210_no[rc>>1]);
165 36 : if (rcn != NPRC) break;
166 18 : rc -= 2; /* cannot wrap since 1 is coprime and rc odd */
167 : }
168 18 : if (rc < rc0) n += rc - rc0;
169 : /* now find an actual (pseudo)prime */
170 : for(;;)
171 : {
172 36 : if (uisprime(n)) break;
173 18 : if (--rcn < 0) rcn = 47;
174 18 : n -= prc210_d1[rcn];
175 : }
176 18 : return n;
177 : }
178 :
179 : GEN
180 49 : precprime(GEN n)
181 : {
182 : long rc, rc0, rcn;
183 49 : pari_sp av = avma;
184 :
185 49 : if (typ(n) != t_INT)
186 : {
187 14 : n = gfloor(n);
188 14 : if (typ(n) != t_INT) pari_err_TYPE("nextprime",n);
189 : }
190 42 : if (signe(n) <= 0) { set_avma(av); return gen_0; }
191 42 : if (lgefint(n) <= 3)
192 : {
193 32 : ulong k = uel(n,2);
194 32 : return gc_utoi(av, uprecprime(k));
195 : }
196 10 : if (!mod2(n)) n = subiu(n,1);
197 10 : rc = rc0 = umodiu(n, 210);
198 : /* find previous prime residue class mod 210 */
199 : for(;;)
200 : {
201 20 : rcn = (long)(prc210_no[rc>>1]);
202 20 : if (rcn != NPRC) break;
203 10 : rc -= 2; /* cannot wrap since 1 is coprime and rc odd */
204 : }
205 10 : if (rc0 > rc) n = subiu(n, rc0 - rc);
206 : /* now find an actual (pseudo)prime */
207 : for(;;)
208 : {
209 48 : if (BPSW_psp(n)) break;
210 38 : if (--rcn < 0) rcn = 47;
211 38 : n = subiu(n, prc210_d1[rcn]);
212 : }
213 10 : if (avma == av) return icopy(n);
214 10 : return gerepileuptoint(av, n);
215 : }
216 :
217 : /* Find next single-word prime strictly larger than p.
218 : * If *n < pari_PRIMES[0], p is *n-th prime, otherwise imitate nextprime().
219 : * *rcn = NPRC or the correct residue class for the current p; we'll use this
220 : * to track the current prime residue class mod 210 once we're out of range of
221 : * the prime table, and we'll update it before that if it isn't NPRC.
222 : *
223 : * *q is incremented whenever q!=NULL and we wrap from 209 mod 210 to
224 : * 1 mod 210 */
225 : static ulong
226 4531317 : snextpr(ulong p, long *n, long *rcn, long *q, int (*ispsp)(ulong))
227 : {
228 4531317 : if (*n < pari_PRIMES[0])
229 : {
230 4531317 : ulong t, p1 = t = pari_PRIMES[++*n]; /* nextprime(p + 1) */
231 4531317 : if (*rcn != NPRC)
232 : {
233 15888894 : while (t > p)
234 : {
235 11373946 : t -= prc210_d1[*rcn];
236 11373946 : if (++*rcn > 47) { *rcn = 0; if (q) (*q)++; }
237 : }
238 : /* assert(d1 == p) */
239 : }
240 4531317 : return p1;
241 : }
242 0 : if (unextprime_overflow(p)) pari_err_OVERFLOW("snextpr");
243 : /* we are beyond the prime table, initialize if needed */
244 0 : if (*rcn == NPRC) *rcn = prc210_no[(p % 210) >> 1]; /* != NPRC */
245 : /* look for the next one */
246 : do {
247 0 : p += prc210_d1[*rcn];
248 0 : if (++*rcn > 47) { *rcn = 0; if (q) (*q)++; }
249 0 : } while (!ispsp(p));
250 0 : return p;
251 : }
252 :
253 : /********************************************************************/
254 : /** **/
255 : /** INTEGER FACTORIZATION **/
256 : /** **/
257 : /********************************************************************/
258 : int factor_add_primes = 0, factor_proven = 0;
259 :
260 : /***********************************************************************/
261 : /** **/
262 : /** FACTORIZATION (ECM) -- GN Jul-Aug 1998 **/
263 : /** Integer factorization using the elliptic curves method (ECM). **/
264 : /** ellfacteur() returns a non trivial factor of N, assuming N>0, **/
265 : /** is composite, and has no prime divisor below tridiv_bound(N) **/
266 : /** Thanks to Paul Zimmermann for much helpful advice and to **/
267 : /** Guillaume Hanrot and Igor Schein for intensive testing **/
268 : /** **/
269 : /***********************************************************************/
270 : #define nbcmax 64 /* max number of simultaneous curves */
271 :
272 : static const ulong TB1[] = {
273 : 142,172,208,252,305,370,450,545,661,801,972,1180,1430,
274 : 1735,2100,2550,3090,3745,4540,5505,6675,8090,9810,11900,
275 : 14420,17490,21200,25700,31160,37780UL,45810UL,55550UL,67350UL,
276 : 81660UL,99010UL,120050UL,145550UL,176475UL,213970UL,259430UL,
277 : 314550UL,381380UL,462415UL,560660UL,679780UL,824220UL,999340UL,
278 : 1211670UL,1469110UL,1781250UL,2159700UL,2618600UL,3175000UL,
279 : 3849600UL,4667500UL,5659200UL,6861600UL,8319500UL,10087100UL,
280 : 12230300UL,14828900UL,17979600UL,21799700UL,26431500UL,
281 : 32047300UL,38856400UL, /* 110 times that still fits into 32bits */
282 : #ifdef LONG_IS_64BIT
283 : 47112200UL,57122100UL,69258800UL,83974200UL,101816200UL,
284 : 123449000UL,149678200UL,181480300UL,220039400UL,266791100UL,
285 : 323476100UL,392204900UL,475536500UL,576573500UL,699077800UL,
286 : 847610500UL,1027701900UL,1246057200UL,1510806400UL,1831806700UL,
287 : 2221009800UL,2692906700UL,3265067200UL,3958794400UL,4799917500UL
288 : #endif
289 : };
290 : static const ulong TB1_for_stage[] = {
291 : /* Start below the optimal B1 for finding factors which would just have been
292 : * missed by pollardbrent(), and escalate, changing curves to give good
293 : * coverage of the small factor ranges. Entries grow faster than what would
294 : * be optimal but a table instead of a 2D array keeps the code simple */
295 : 500,520,560,620,700,800,900,1000,1150,1300,1450,1600,1800,2000,
296 : 2200,2450,2700,2950,3250,3600,4000,4400,4850,5300,5800,6400,
297 : 7100,7850,8700,9600,10600,11700,12900,14200,15700,17300,
298 : 19000,21000,23200,25500,28000,31000,34500UL,38500UL,43000UL,
299 : 48000UL,53800UL,60400UL,67750UL,76000UL,85300UL,95700UL,
300 : 107400UL,120500UL,135400UL,152000UL,170800UL,191800UL,215400UL,
301 : 241800UL,271400UL,304500UL,341500UL,383100UL,429700UL,481900UL,
302 : 540400UL,606000UL,679500UL,761800UL,854100UL,957500UL,1073500UL
303 : };
304 :
305 : /* addition/doubling/multiplication of a point on an 'elliptic curve mod N'
306 : * may result in one of three things:
307 : * - a new bona fide point
308 : * - a point at infinity (denominator divisible by N)
309 : * - a point at infinity mod some p | N but finite mod q | N betraying itself
310 : * by a denominator which has nontrivial gcd with N.
311 : *
312 : * In the second case, addition/doubling aborts, copying one of the summands
313 : * to the destination array of points unless they coincide.
314 : * Multiplication will stop at some unpredictable intermediate stage: The
315 : * destination will contain _some_ multiple of the input point, but not
316 : * necessarily the desired one, which doesn't matter. As long as we're
317 : * multiplying (B1 phase) we simply carry on with the next multiplier.
318 : * During the B2 phase, the only additions are the giant steps, and the
319 : * worst that can happen here is that we lose one residue class mod 210
320 : * of prime multipliers on 4 of the curves, so again, we ignore the problem
321 : * and just carry on.)
322 : *
323 : * Idea: select nbc curves mod N and one point P on each of them. For each
324 : * such P, compute [M]P = Q where M is the product of all powers <= B2 of
325 : * primes <= nextprime(B1). Then check whether [p]Q for p < nextprime(B2)
326 : * betrays a factor. This second stage looks separately at the primes in
327 : * each residue class mod 210, four curves at a time, and steps additively
328 : * to ever larger multipliers, by comparing X coordinates of points which we
329 : * would need to add in order to reach another prime multiplier in the same
330 : * residue class. 'Comparing' means that we accumulate a product of
331 : * differences of X coordinates, and from time to time take a gcd of this
332 : * product with N. Montgomery's multi-inverse trick is used heavily. */
333 :
334 : /* *** auxiliary functions for ellfacteur: *** */
335 : /* (Rx,Ry) <- (Px,Py)+(Qx,Qy) over Z/NZ, z=1/(Px-Qx). If Ry = NULL, don't set */
336 : static void
337 8291496 : FpE_add_i(GEN N, GEN z, GEN Px, GEN Py, GEN Qx, GEN Qy, GEN *Rx, GEN *Ry)
338 : {
339 8291496 : GEN slope = modii(mulii(subii(Py, Qy), z), N);
340 8291496 : GEN t = subii(sqri(slope), addii(Qx, Px));
341 8291496 : affii(modii(t, N), *Rx);
342 8291496 : if (Ry) {
343 8219188 : t = subii(mulii(slope, subii(Px, *Rx)), Py);
344 8219188 : affii(modii(t, N), *Ry);
345 : }
346 8291496 : }
347 : /* X -> Z; cannot add on one of the curves: make sure Z contains
348 : * something useful before letting caller proceed */
349 : static void
350 25650 : ZV_aff(long n, GEN *X, GEN *Z)
351 : {
352 25650 : if (X != Z) {
353 : long k;
354 1569330 : for (k = n; k--; ) affii(X[k],Z[k]);
355 : }
356 25650 : }
357 :
358 : /* Parallel addition on nbc curves, assigning the result to locations at and
359 : * following *X3, *Y3. (If Y-coords of result not desired, set Y=NULL.)
360 : * Safe even if (X3,Y3) = (X2,Y2), _not_ if (X1,Y1). It is also safe to
361 : * overwrite Y2 with X3. If nbc1 < nbc, the first summand is
362 : * assumed to hold only nbc1 distinct points, repeated as often as we need
363 : * them (to add one point on each of a few curves to several other points on
364 : * the same curves): only used with nbc1 = nbc or nbc1 = 4 | nbc.
365 : *
366 : * Return 0 [SUCCESS], 1 [N | den], 2 [gcd(den, N) is a factor of N, preserved
367 : * in gl.
368 : * Stack space is bounded by a constant multiple of lgefint(N)*nbc:
369 : * - Phase 2 creates 12 items on the stack per iteration, of which 4 are twice
370 : * as long and 1 is thrice as long as N, i.e. 18 units per iteration.
371 : * - Phase 1 creates 4 units.
372 : * Total can be as large as 4*nbcmax + 18*8 units; ecm_elladd2() is
373 : * just as bad, and elldouble() comes to 3*nbcmax + 29*8 units. */
374 : static int
375 240431 : ecm_elladd0(GEN N, GEN *gl, long nbc, long nbc1,
376 : GEN *X1, GEN *Y1, GEN *X2, GEN *Y2, GEN *X3, GEN *Y3)
377 : {
378 240431 : const ulong mask = (nbc1 == 4)? 3: ~0UL; /*nbc1 = 4 or nbc*/
379 240431 : GEN W[2*nbcmax], *A = W+nbc; /* W[0],A[0] unused */
380 : long i;
381 240431 : pari_sp av = avma;
382 :
383 240431 : W[1] = subii(X1[0], X2[0]);
384 7825896 : for (i=1; i<nbc; i++)
385 : { /*prepare for multi-inverse*/
386 7585465 : A[i] = subii(X1[i&mask], X2[i]); /* don't waste time reducing mod N */
387 7585465 : W[i+1] = modii(mulii(A[i], W[i]), N);
388 : }
389 240431 : if (!invmod(W[nbc], N, gl))
390 : {
391 18 : if (!equalii(N,*gl)) return 2;
392 0 : ZV_aff(nbc, X2,X3);
393 0 : if (Y3) ZV_aff(nbc, Y2,Y3);
394 0 : return gc_int(av,1);
395 : }
396 :
397 7825032 : while (i--) /* nbc times */
398 : {
399 7825032 : pari_sp av2 = avma;
400 7825032 : GEN Px = X1[i&mask], Py = Y1[i&mask], Qx = X2[i], Qy = Y2[i];
401 7825032 : GEN z = i? mulii(*gl,W[i]): *gl; /*1/(Px-Qx)*/
402 7825032 : FpE_add_i(N,z, Px,Py,Qx,Qy, X3+i, Y3? Y3+i: NULL);
403 7825032 : if (!i) break;
404 7584619 : set_avma(av2); *gl = modii(mulii(*gl, A[i]), N);
405 : }
406 240413 : return gc_int(av,0);
407 : }
408 :
409 : /* Shortcut, for use in cases where Y coordinates follow their corresponding
410 : * X coordinates, and first summand doesn't need to be repeated */
411 : static int
412 234487 : ecm_elladd(GEN N, GEN *gl, long nbc, GEN *X1, GEN *X2, GEN *X3) {
413 234487 : return ecm_elladd0(N, gl, nbc, nbc, X1, X1+nbc, X2, X2+nbc, X3, X3+nbc);
414 : }
415 :
416 : /* As ecm_elladd except it does twice as many additions (and hides even more
417 : * of the cost of the modular inverse); the net effect is the same as
418 : * ecm_elladd(nbc,X1,X2,X3) && ecm_elladd(nbc,X4,X5,X6). Safe to
419 : * have X2=X3, X5=X6, or X1,X2 coincide with X4,X5 in any order. */
420 : static int
421 7194 : ecm_elladd2(GEN N, GEN *gl, long nbc,
422 : GEN *X1, GEN *X2, GEN *X3, GEN *X4, GEN *X5, GEN *X6)
423 : {
424 7194 : GEN *Y1 = X1+nbc, *Y2 = X2+nbc, *Y3 = X3+nbc;
425 7194 : GEN *Y4 = X4+nbc, *Y5 = X5+nbc, *Y6 = X6+nbc;
426 7194 : GEN W[4*nbcmax], *A = W+2*nbc; /* W[0],A[0] unused */
427 : long i, j;
428 7194 : pari_sp av = avma;
429 :
430 7194 : W[1] = subii(X1[0], X2[0]);
431 233232 : for (i=1; i<nbc; i++)
432 : {
433 226038 : A[i] = subii(X1[i], X2[i]); /* don't waste time reducing mod N here */
434 226038 : W[i+1] = modii(mulii(A[i], W[i]), N);
435 : }
436 240426 : for (j=0; j<nbc; i++,j++)
437 : {
438 233232 : A[i] = subii(X4[j], X5[j]);
439 233232 : W[i+1] = modii(mulii(A[i], W[i]), N);
440 : }
441 7194 : if (!invmod(W[2*nbc], N, gl))
442 : {
443 0 : if (!equalii(N,*gl)) return 2;
444 0 : ZV_aff(2*nbc, X2,X3); /* hack: 2*nbc => copy Y2->Y3 */
445 0 : ZV_aff(2*nbc, X5,X6); /* also copy Y5->Y6 */
446 0 : return gc_int(av,1);
447 : }
448 :
449 240426 : while (j--) /* nbc times */
450 : {
451 233232 : pari_sp av2 = avma;
452 233232 : GEN Px = X4[j], Py = Y4[j], Qx = X5[j], Qy = Y5[j];
453 233232 : GEN z = mulii(*gl,W[--i]); /*1/(Px-Qx)*/
454 233232 : FpE_add_i(N,z, Px,Py, Qx,Qy, X6+j,Y6+j);
455 233232 : set_avma(av2); *gl = modii(mulii(*gl, A[i]), N);
456 : }
457 233232 : while (i--) /* nbc times */
458 : {
459 233232 : pari_sp av2 = avma;
460 233232 : GEN Px = X1[i], Py = Y1[i], Qx = X2[i], Qy = Y2[i];
461 233232 : GEN z = i? mulii(*gl, W[i]): *gl; /*1/(Px-Qx)*/
462 233232 : FpE_add_i(N,z, Px,Py, Qx,Qy, X3+i,Y3+i);
463 233232 : if (!i) break;
464 226038 : set_avma(av2); *gl = modii(mulii(*gl, A[i]), N);
465 : }
466 7194 : return gc_int(av,0);
467 : }
468 :
469 : /* Parallel doubling on nbc curves, assigning the result to locations at
470 : * and following *X2. Safe to be called with X2 equal to X1. Return
471 : * value as for ecm_elladd. If we find a point at infinity mod N,
472 : * and if X1 != X2, we copy the points at X1 to X2. */
473 : static int
474 40073 : elldouble(GEN N, GEN *gl, long nbc, GEN *X1, GEN *X2)
475 : {
476 40073 : GEN *Y1 = X1+nbc, *Y2 = X2+nbc;
477 : GEN W[nbcmax+1]; /* W[0] unused */
478 : long i;
479 40073 : pari_sp av = avma;
480 40073 : /*W[0] = gen_1;*/ W[1] = Y1[0];
481 1233544 : for (i=1; i<nbc; i++) W[i+1] = modii(mulii(Y1[i], W[i]), N);
482 40073 : if (!invmod(W[nbc], N, gl))
483 : {
484 0 : if (!equalii(N,*gl)) return 2;
485 0 : ZV_aff(2*nbc,X1,X2); /* also copies Y1->Y2 */
486 0 : return gc_int(av,1);
487 : }
488 1273617 : while (i--) /* nbc times */
489 : {
490 : pari_sp av2;
491 1233544 : GEN v, w, L, z = i? mulii(*gl,W[i]): *gl;
492 1233544 : if (i) *gl = modii(mulii(*gl, Y1[i]), N);
493 1233544 : av2 = avma;
494 1233544 : L = modii(mulii(addui(1, mului(3, Fp_sqr(X1[i],N))), z), N);
495 1233544 : if (signe(L)) /* half of zero is still zero */
496 1233544 : L = shifti(mod2(L)? addii(L, N): L, -1);
497 1233544 : v = modii(subii(sqri(L), shifti(X1[i],1)), N);
498 1233544 : w = modii(subii(mulii(L, subii(X1[i], v)), Y1[i]), N);
499 1233544 : affii(v, X2[i]);
500 1233544 : affii(w, Y2[i]);
501 1233544 : set_avma(av2);
502 : }
503 40073 : return gc_int(av,0);
504 : }
505 :
506 : /* Parallel multiplication by an odd prime k on nbc curves, storing the
507 : * result to locations at and following *X2. Safe to be called with X2 = X1.
508 : * Return values as ecm_elladd. Uses (a simplified variant of) Montgomery's
509 : * PRAC algorithm; see ftp://ftp.cwi.nl/pub/pmontgom/Lucas.ps.gz .
510 : * With thanks to Paul Zimmermann for the reference. --GN1998Aug13 */
511 : static int
512 208727 : get_rule(ulong d, ulong e)
513 : {
514 208727 : if (d <= e + (e>>2)) /* floor(1.25*e) */
515 : {
516 16630 : if ((d+e)%3 == 0) return 0; /* rule 1 */
517 9928 : if ((d-e)%6 == 0) return 1; /* rule 2 */
518 : }
519 : /* d <= 4*e but no ofl */
520 201971 : if ((d+3)>>2 <= e) return 2; /* rule 3, common case */
521 12148 : if ((d&1)==(e&1)) return 1; /* rule 4 = rule 2 */
522 6344 : if (!(d&1)) return 3; /* rule 5 */
523 1769 : if (d%3 == 0) return 4; /* rule 6 */
524 417 : if ((d+e)%3 == 0) return 5; /* rule 7 */
525 0 : if ((d-e)%3 == 0) return 6; /* rule 8 */
526 : /* when we get here, e is even, otherwise one of rules 4,5 would apply */
527 0 : return 7; /* rule 9 */
528 : }
529 :
530 : /* PRAC implementation notes - main changes against the paper version:
531 : * (1) The general function [m+n]P = f([m]P,[n]P,[m-n]P) collapses (for m!=n)
532 : * to an ecm_elladd() which does not depend on the third argument; thus
533 : * references to the third variable (C in the paper) can be eliminated.
534 : * (2) Since our multipliers are prime, the outer loop of the paper
535 : * version executes only once, and thus is invisible above.
536 : * (3) The first step in the inner loop of the paper version will always be
537 : * rule 3, but the addition requested by this rule amounts to a doubling, and
538 : * will always be followed by a swap, so we have unrolled this first iteration.
539 : * (4) Simplifications in rules 6 and 7 are possible given the above, and we
540 : * save one addition in each of the two cases. NB none of the other
541 : * ecm_elladd()s in the loop can ever degenerate into an elldouble.
542 : * (5) I tried to optimize for rule 3, which is used more frequently than all
543 : * others together, but it didn't improve things, so I removed the nested
544 : * tight loop again. --GN */
545 : /* The main loop body of ellfacteur() runs _slower_ under PRAC than under a
546 : * straightforward left-shift binary multiplication when N has <30 digits and
547 : * B1 is small; PRAC wins when N and B1 get larger. Weird. --GN */
548 : /* k>2 assumed prime, XAUX = scratchpad */
549 : static int
550 25650 : ellmult(GEN N, GEN *gl, long nbc, ulong k, GEN *X1, GEN *X2, GEN *XAUX)
551 : {
552 : ulong r, d, e, e1;
553 : int res;
554 25650 : GEN *A = X2, *B = XAUX, *T = XAUX + 2*nbc;
555 :
556 25650 : ZV_aff(2*nbc,X1,XAUX);
557 : /* first doubling picks up X1; after this we'll be working in XAUX and
558 : * X2 only, mostly via A and B and T */
559 25650 : if ((res = elldouble(N, gl, nbc, X1, X2)) != 0) return res;
560 :
561 : /* split the work at the golden ratio */
562 25650 : r = (ulong)(k*0.61803398875 + .5);
563 25650 : d = k - r;
564 25650 : e = r - d; /* d+e == r, so no danger of ofl below */
565 234377 : while (d != e)
566 : { /* apply one of the nine transformations from PM's Table 4. */
567 208727 : switch(get_rule(d,e))
568 : {
569 6702 : case 0: /* rule 1 */
570 6702 : if ( (res = ecm_elladd(N, gl, nbc, A, B, T)) ) return res;
571 6702 : if ( (res = ecm_elladd2(N, gl, nbc, T, A, A, T, B, B)) != 0) return res;
572 6702 : e1 = d - e; d = (d + e1)/3; e = (e - e1)/3; break;
573 5858 : case 1: /* rules 2 and 4 */
574 5858 : if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
575 5858 : if ( (res = elldouble(N, gl, nbc, A, A)) ) return res;
576 5858 : d = (d-e)>>1; break;
577 4575 : case 3: /* rule 5 */
578 4575 : if ( (res = elldouble(N, gl, nbc, A, A)) ) return res;
579 4575 : d >>= 1; break;
580 1352 : case 4: /* rule 6 */
581 1352 : if ( (res = elldouble(N, gl, nbc, A, T)) ) return res;
582 1352 : if ( (res = ecm_elladd(N, gl, nbc, T, A, A)) ) return res;
583 1352 : if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
584 1352 : d = d/3 - e; break;
585 189823 : case 2: /* rule 3 */
586 189823 : if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
587 189823 : d -= e; break;
588 417 : case 5: /* rule 7 */
589 417 : if ( (res = elldouble(N, gl, nbc, A, T)) ) return res;
590 417 : if ( (res = ecm_elladd2(N, gl, nbc, T, A, A, T, B, B)) != 0) return res;
591 417 : d = (d - 2*e)/3; break;
592 0 : case 6: /* rule 8 */
593 0 : if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
594 0 : if ( (res = elldouble(N, gl, nbc, A, T)) ) return res;
595 0 : if ( (res = ecm_elladd(N, gl, nbc, T, A, A)) ) return res;
596 0 : d = (d - e)/3; break;
597 0 : case 7: /* rule 9 */
598 0 : if ( (res = elldouble(N, gl, nbc, B, B)) ) return res;
599 0 : e >>= 1; break;
600 : }
601 : /* swap d <-> e and A <-> B if necessary */
602 208727 : if (d < e) { lswap(d,e); pswap(A,B); }
603 : }
604 25650 : return ecm_elladd(N, gl, nbc, XAUX, X2, X2);
605 : }
606 :
607 : struct ECM {
608 : pari_timer T;
609 : long nbc, nbc2, seed;
610 : GEN *X, *XAUX, *XT, *XD, *XB, *XB2, *XH, *Xh, *Yh;
611 : };
612 :
613 : /* memory layout in ellfacteur(): a large array of GEN pointers, and one
614 : * huge chunk of memory containing all the actual GEN (t_INT) objects.
615 : * nbc is constant throughout the invocation:
616 : * - The B1 stage of each iteration through the main loop needs little
617 : * space: enough for the X and Y coordinates of the current points,
618 : * and twice as much again as scratchpad for ellmult().
619 : * - The B2 stage, starting from some current set of points Q, needs, in
620 : * succession:
621 : * + space for [2]Q, [4]Q, ..., [10]Q, and [p]Q for building the helix;
622 : * + space for 48*nbc X and Y coordinates to hold the helix. This could
623 : * re-use [2]Q,...,[8]Q, but only with difficulty, since we don't
624 : * know in advance which residue class mod 210 our p is going to be in.
625 : * It can and should re-use [p]Q, though;
626 : * + space for (temporarily [30]Q and then) [210]Q, [420]Q, and several
627 : * further doublings until the giant step multiplier is reached. This
628 : * can re-use the remaining cells from above. The computation of [210]Q
629 : * will have been the last call to ellmult() within this iteration of the
630 : * main loop, so the scratchpad is now also free to be re-used. We also
631 : * compute [630]Q by a parallel addition; we'll need it later to get the
632 : * baby-step table bootstrapped a little faster.
633 : * + Finally, for no more than 4 curves at a time, room for up to 1024 X
634 : * coordinates only: the Y coordinates needed whilst setting up this baby
635 : * step table are temporarily stored in the upper half, and overwritten
636 : * during the last series of additions.
637 : *
638 : * Graphically: after end of B1 stage (X,Y are the coords of Q):
639 : * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
640 : * | X Y | scratch | [2]Q| [4]Q| [6]Q| [8]Q|[10]Q| ... | ...
641 : * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
642 : * *X *XAUX *XT *XD *XB
643 : *
644 : * [30]Q is computed from [10]Q. [210]Q can go into XY, etc:
645 : * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
646 : * |[210]|[420]|[630]|[840]|[1680,3360,6720,...,2048*210] |bstp table...
647 : * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
648 : * *X *XAUX *XT *XD [*XG, somewhere here] *XB .... *XH
649 : *
650 : * So we need (13 + 48) * 2 * nbc slots here + 4096 slots for the baby step
651 : * table (not all of which will be used when we start with a small B1, but
652 : * better to allocate and initialize ahead of time all the slots that might
653 : * be needed later).
654 : *
655 : * Note on memory locality: During the B2 phase, accesses to the helix
656 : * (once it is set up) will be clustered by curves (4 out of nbc at a time).
657 : * Accesses to the baby steps table will wander from one end of the array to
658 : * the other and back, one such cycle per giant step, and during a full cycle
659 : * we would expect on the order of 2E4 accesses when using the largest giant
660 : * step size. Thus we shouldn't be doing too bad with respect to thrashing
661 : * a 512KBy L2 cache. However, we don't want the baby step table to grow
662 : * larger than this, even if it would reduce the number of EC operations by a
663 : * few more per cent for very large B2, lest cache thrashing slow down
664 : * everything disproportionally. --GN */
665 : /* Auxiliary routines need < (3*nbc+240)*lN words on the PARI stack, in
666 : * addition to the spc*(lN+1) words occupied by our main table. */
667 : static void
668 57 : ECM_alloc(struct ECM *E, long lN)
669 : {
670 57 : const long bstpmax = 1024; /* max number of baby step table entries */
671 57 : long spc = (13 + 48) * E->nbc2 + bstpmax * 4;
672 57 : long len = spc + 385 + spc*lN;
673 57 : long i, tw = _evallg(lN) | evaltyp(t_INT);
674 57 : GEN w, *X = (GEN*)new_chunk(len);
675 : /* hack for X[i] = cgeti(lN). X = current point in B1 phase */
676 57 : w = (GEN)(X + spc + 385);
677 377001 : for (i = spc-1; i >= 0; i--) { X[i] = w; *w = tw; w += lN; }
678 57 : E->X = X;
679 57 : E->XAUX = E->X + E->nbc2; /* scratchpad for ellmult() */
680 57 : E->XT = E->XAUX + E->nbc2; /* ditto, will later hold [3*210]Q */
681 57 : E->XD = E->XT + E->nbc2; /* room for various multiples */
682 57 : E->XB = E->XD + 10*E->nbc2; /* start of baby steps table */
683 57 : E->XB2 = E->XB + 2 * bstpmax; /* middle of baby steps table */
684 57 : E->XH = E->XB2 + 2 * bstpmax; /* end of bstps table, start of helix */
685 57 : E->Xh = E->XH + 48*E->nbc2; /* little helix, X coords */
686 57 : E->Yh = E->XH + 192; /* ditto, Y coords */
687 : /* XG,YG set inside the main loop, since they depend on B2 */
688 : /* E.Xh range of 384 pointers not set; these will later duplicate the pointers
689 : * in the E.XH range, 4 curves at a time. Some of the cells reserved here for
690 : * the E.XB range will never be used, instead, we'll warp the pointers to
691 : * connect to (read-only) GENs in the X/E.XD range */
692 57 : }
693 : /* N.B. E->seed is not initialized here */
694 : static void
695 57 : ECM_init(struct ECM *E, GEN N, long nbc)
696 : {
697 57 : if (nbc < 0)
698 : { /* choose a sensible default */
699 57 : const long size = expi(N) + 1;
700 57 : nbc = ((size >> 3) << 2) - 80;
701 57 : if (nbc < 8) nbc = 8;
702 : }
703 57 : if (nbc > nbcmax) nbc = nbcmax;
704 57 : E->nbc = nbc;
705 57 : E->nbc2 = nbc << 1;
706 57 : ECM_alloc(E, lgefint(N));
707 57 : }
708 :
709 : static GEN
710 93 : ECM_loop(struct ECM *E, GEN N, ulong B1)
711 : {
712 93 : const ulong B2 = 110 * B1, B2_rt = usqrt(B2);
713 93 : const ulong nbc = E->nbc, nbc2 = E->nbc2;
714 : pari_sp av1, avtmp;
715 : long i, np, np0, gse, gss, bstp, bstp0, rcn0, rcn;
716 : ulong B2_p, m, p, p0;
717 : GEN g, *XG, *YG;
718 93 : GEN *X = E->X, *XAUX = E->XAUX, *XT = E->XT, *XD = E->XD;
719 93 : GEN *XB = E->XB, *XB2 = E->XB2, *XH = E->XH, *Xh = E->Xh, *Yh = E->Yh;
720 : /* pick curves */
721 4461 : for (i = nbc2; i--; ) affui(E->seed++, X[i]);
722 : /* pick giant step exponent and size */
723 93 : gse = B1 < 656
724 : ? (B1 < 200? 5: 6)
725 93 : : (B1 < 10500
726 : ? (B1 < 2625? 7: 8)
727 : : (B1 < 42000? 9: 10));
728 93 : gss = 1UL << gse;
729 : /* With 32 baby steps, a giant step corresponds to 32*420 = 13440,
730 : * appropriate for the smallest B2s. With 1024, a giant step will be 430080;
731 : * appropriate for B1 >~ 42000, where 512 baby steps would imply roughly
732 : * the same number of curve additions. */
733 93 : XG = XT + gse*nbc2; /* will later hold [2^(gse+1)*210]Q */
734 93 : YG = XG + nbc;
735 :
736 93 : if (DEBUGLEVEL >= 4) {
737 0 : err_printf("ECM: time = %6ld ms\nECM: B1 = %4lu,", timer_delay(&E->T), B1);
738 0 : err_printf("\tB2 = %6lu,\tgss = %4ld*420\n", B2, gss);
739 : }
740 93 : p = 2; np = 1; /* p is np-th prime */
741 :
742 : /* ---B1 PHASE--- */
743 : /* treat p=2 separately */
744 93 : B2_p = B2 >> 1;
745 1603 : for (m=1; m<=B2_p; m<<=1)
746 : {
747 1510 : int fl = elldouble(N, &g, nbc, X, X);
748 1510 : if (fl > 1) return g; else if (fl) break;
749 : }
750 93 : rcn = NPRC; /* multipliers begin at the beginning */
751 : /* p=3,...,nextprime(B1) */
752 6538 : while (p < B1 && p <= B2_rt)
753 : {
754 6445 : pari_sp av2 = avma;
755 6445 : p = snextpr(p, &np, &rcn, NULL, uisprime);
756 6445 : B2_p = B2/p; /* beware integer overflow on 32-bit CPUs */
757 22021 : for (m=1; m<=B2_p; m*=p)
758 : {
759 15576 : int fl = ellmult(N, &g, nbc, p, X, X, XAUX);
760 15576 : if (fl > 1) return g; else if (fl) break;
761 15576 : set_avma(av2);
762 : }
763 6445 : set_avma(av2);
764 : }
765 : /* primes p larger than sqrt(B2) appear only to the 1st power */
766 9924 : while (p < B1)
767 : {
768 9849 : pari_sp av2 = avma;
769 9849 : p = snextpr(p, &np, &rcn, NULL, uisprime);
770 9849 : if (ellmult(N, &g, nbc, p, X, X, XAUX) > 1) return g;
771 9831 : set_avma(av2);
772 : }
773 75 : if (DEBUGLEVEL >= 4) {
774 0 : err_printf("ECM: time = %6ld ms, B1 phase done, ", timer_delay(&E->T));
775 0 : err_printf("p = %lu, setting up for B2\n", p);
776 : }
777 :
778 : /* ---B2 PHASE--- */
779 : /* compute [2]Q,...,[10]Q, needed to build the helix */
780 75 : if (elldouble(N, &g, nbc, X, XD) > 1) return g; /*[2]Q*/
781 75 : if (elldouble(N, &g, nbc, XD, XD + nbc2) > 1) return g; /*[4]Q*/
782 75 : if (ecm_elladd(N, &g, nbc,
783 75 : XD, XD + nbc2, XD + (nbc<<2)) > 1) return g; /* [6]Q */
784 75 : if (ecm_elladd2(N, &g, nbc,
785 : XD, XD + (nbc<<2), XT + (nbc<<3),
786 75 : XD + nbc2, XD + (nbc<<2), XD + (nbc<<3)) > 1)
787 0 : return g; /* [8]Q and [10]Q */
788 75 : if (DEBUGLEVEL >= 7) err_printf("\t(got [2]Q...[10]Q)\n");
789 :
790 : /* get next prime (still using the foolproof test) */
791 75 : p = snextpr(p, &np, &rcn, NULL, uisprime);
792 : /* make sure we have the residue class number (mod 210) */
793 75 : if (rcn == NPRC)
794 : {
795 75 : rcn = prc210_no[(p % 210) >> 1];
796 75 : if (rcn == NPRC)
797 : {
798 0 : err_printf("ECM: %lu should have been prime but isn\'t\n", p);
799 0 : pari_err_BUG("ellfacteur");
800 : }
801 : }
802 :
803 : /* compute [p]Q and put it into its place in the helix */
804 75 : if (ellmult(N, &g, nbc, p, X, XH + rcn*nbc2, XAUX) > 1)
805 0 : return g;
806 75 : if (DEBUGLEVEL >= 7)
807 0 : err_printf("\t(got [p]Q, p = %lu = prc210_rp[%ld] mod 210)\n", p, rcn);
808 :
809 : /* save current p, np, and rcn; we'll need them more than once below */
810 75 : p0 = p; np0 = np; rcn0 = rcn;
811 75 : bstp0 = 0; /* p is at baby-step offset 0 from itself */
812 :
813 : /* fill up the helix, stepping forward through the prime residue classes
814 : * mod 210 until we're back at the r'class of p0. Keep updating p so
815 : * that we can print meaningful diagnostics if a factor shows up; don't
816 : * bother checking which of these p's are in fact prime */
817 3600 : for (i = 47; i; i--) /* 47 iterations */
818 : {
819 3525 : ulong dp = (ulong)prc210_d1[rcn];
820 3525 : p += dp;
821 3525 : if (rcn == 47)
822 : { /* wrap mod 210 */
823 75 : if (ecm_elladd(N, &g, nbc, XT+dp*nbc, XH+rcn*nbc2, XH) > 1) return g;
824 75 : rcn = 0; continue;
825 : }
826 3450 : if (ecm_elladd(N, &g, nbc, XT+dp*nbc, XH+rcn*nbc2, XH+rcn*nbc2+nbc2) > 1)
827 0 : return g;
828 3450 : rcn++;
829 : }
830 75 : if (DEBUGLEVEL >= 7) err_printf("\t(got initial helix)\n");
831 : /* compute [210]Q etc, needed for the baby step table */
832 75 : if (ellmult(N, &g, nbc, 3, XD + (nbc<<3), X, XAUX) > 1) return g;
833 75 : if (ellmult(N, &g, nbc, 7, X, X, XAUX) > 1) return g; /* [210]Q */
834 : /* this was the last call to ellmult() in the main loop body; may now
835 : * overwrite XAUX and slots XD and following */
836 75 : if (elldouble(N, &g, nbc, X, XAUX) > 1) return g; /* [420]Q */
837 75 : if (ecm_elladd(N, &g, nbc, X, XAUX, XT) > 1) return g;/*[630]Q*/
838 75 : if (ecm_elladd(N, &g, nbc, X, XT, XD) > 1) return g; /*[840]Q*/
839 561 : for (i=1; i <= gse; i++)
840 486 : if (elldouble(N, &g, nbc, XT + i*nbc2, XD + i*nbc2) > 1) return g;
841 : /* (the last iteration has initialized XG to [210*2^(gse+1)]Q) */
842 :
843 75 : if (DEBUGLEVEL >= 4)
844 0 : err_printf("ECM: time = %6ld ms, entering B2 phase, p = %lu\n",
845 : timer_delay(&E->T), p);
846 :
847 391 : for (i = nbc - 4; i >= 0; i -= 4)
848 : { /* loop over small sets of 4 curves at a time */
849 : GEN *Xb;
850 : long j, k;
851 323 : if (DEBUGLEVEL >= 6)
852 0 : err_printf("ECM: finishing curves %ld...%ld\n", i, i+3);
853 : /* Copy relevant pointers from XH to Xh. Memory layout in XH:
854 : * nbc X coordinates, nbc Y coordinates for residue class
855 : * 1 mod 210, then the same for r.c. 11 mod 210, etc. Memory layout for
856 : * Xh is: four X coords for 1 mod 210, four for 11 mod 210, ..., four
857 : * for 209 mod 210, then the corresponding Y coordinates in the same
858 : * order. This allows a giant step on Xh using just three calls to
859 : * ecm_elladd0() each acting on 64 points in parallel */
860 15827 : for (j = 48; j--; )
861 : {
862 15504 : k = nbc2*j + i;
863 15504 : m = j << 2; /* X coordinates */
864 15504 : Xh[m] = XH[k]; Xh[m+1] = XH[k+1];
865 15504 : Xh[m+2] = XH[k+2]; Xh[m+3] = XH[k+3];
866 15504 : k += nbc; /* Y coordinates */
867 15504 : Yh[m] = XH[k]; Yh[m+1] = XH[k+1];
868 15504 : Yh[m+2] = XH[k+2]; Yh[m+3] = XH[k+3];
869 : }
870 : /* Build baby step table of X coords of multiples of [210]Q. XB[4*j]
871 : * will point at X coords on four curves from [(j+1)*210]Q. Until
872 : * we're done, we need some Y coords as well, which we keep in the
873 : * second half of the table, overwriting them at the end when gse=10.
874 : * Multiples which we already have (by 1,2,3,4,8,16,...,2^gse) are
875 : * entered simply by copying the pointers, ignoring the few slots in w
876 : * that were initially reserved for them. Here are the initial entries */
877 969 : for (Xb=XB,k=2,j=i; k--; Xb=XB2,j+=nbc) /* first X, then Y coords */
878 : {
879 646 : Xb[0] = X[j]; Xb[1] = X[j+1]; /* [210]Q */
880 646 : Xb[2] = X[j+2]; Xb[3] = X[j+3];
881 646 : Xb[4] = XAUX[j]; Xb[5] = XAUX[j+1]; /* [420]Q */
882 646 : Xb[6] = XAUX[j+2]; Xb[7] = XAUX[j+3];
883 646 : Xb[8] = XT[j]; Xb[9] = XT[j+1]; /* [630]Q */
884 646 : Xb[10] = XT[j+2]; Xb[11] = XT[j+3];
885 646 : Xb += 4; /* points at [420]Q */
886 : /* ... entries at powers of 2 times 210 .... */
887 4057 : for (m = 2; m < (ulong)gse+k; m++) /* omit Y coords of [2^gse*210]Q */
888 : {
889 3411 : long m2 = m*nbc2 + j;
890 3411 : Xb += (2UL<<m); /* points at [2^m*210]Q */
891 3411 : Xb[0] = XAUX[m2]; Xb[1] = XAUX[m2+1];
892 3411 : Xb[2] = XAUX[m2+2]; Xb[3] = XAUX[m2+3];
893 : }
894 : }
895 323 : if (DEBUGLEVEL >= 7)
896 0 : err_printf("\t(extracted precomputed helix / baby step entries)\n");
897 : /* ... glue in between, up to 16*210 ... */
898 323 : if (ecm_elladd0(N, &g, 12, 4, /* 12 pts + (4 pts replicated thrice) */
899 : XB + 12, XB2 + 12,
900 : XB, XB2,
901 0 : XB + 16, XB2 + 16) > 1) return g; /*4+{1,2,3} = {5,6,7}*/
902 323 : if (ecm_elladd0(N, &g, 28, 4, /* 28 pts + (4 pts replicated 7fold) */
903 : XB + 28, XB2 + 28,
904 : XB, XB2,
905 0 : XB + 32, XB2 + 32) > 1) return g;/*8+{1...7} = {9...15}*/
906 : /* ... and the remainder of the lot */
907 1221 : for (m = 5; m <= (ulong)gse; m++)
908 : { /* fill in from 2^(m-1)+1 to 2^m-1 in chunks of 64 and 60 points */
909 898 : ulong m2 = 2UL << m; /* will point at 2^(m-1)+1 */
910 1977 : for (j = 0; (ulong)j < m2-64; j+=64) /* executed 0 times when m = 5 */
911 : {
912 1906 : if (ecm_elladd0(N, &g, 64, 4,
913 1079 : XB + m2-4, XB2 + m2-4,
914 1079 : XB + j, XB2 + j,
915 1906 : XB + m2+j, (m<(ulong)gse? XB2+m2+j: NULL)) > 1)
916 0 : return g;
917 : } /* j = m2-64 here, 60 points left */
918 1221 : if (ecm_elladd0(N, &g, 60, 4,
919 898 : XB + m2-4, XB2 + m2-4,
920 898 : XB + j, XB2 + j,
921 1221 : XB + m2+j, (m<(ulong)gse? XB2+m2+j: NULL)) > 1)
922 0 : return g;
923 : /* when m=gse, drop Y coords of result, and when both equal 1024,
924 : * overwrite Y coords of second argument with X coords of result */
925 : }
926 323 : if (DEBUGLEVEL >= 7) err_printf("\t(baby step table complete)\n");
927 : /* initialize a few other things */
928 323 : bstp = bstp0; p = p0; np = np0; rcn = rcn0;
929 323 : g = gen_1; av1 = avma;
930 : /* scratchspace for prod (x_i-x_j) */
931 323 : avtmp = (pari_sp)new_chunk(8 * lgefint(N));
932 : /* The correct entry in XB to use depends on bstp and on where we are
933 : * on the helix. As we skip from prime to prime, bstp is incremented
934 : * by snextpr each time we wrap around through residue class number 0
935 : * (1 mod 210), but the baby step should not be taken until rcn>=rcn0,
936 : * i.e. until we pass again the residue class of p0.
937 : *
938 : * The correct signed multiplier is thus k = bstp - (rcn < rcn0),
939 : * and the offset from XB is four times (|k| - 1). When k=0, we ignore
940 : * the current prime: if it had led to a factorization, this
941 : * would have been noted during the last giant step, or -- when we
942 : * first get here -- whilst initializing the helix. When k > gss,
943 : * we must do a giant step and bump bstp back by -2*gss.
944 : *
945 : * The gcd of the product of X coord differences against N is taken just
946 : * before we do a giant step. */
947 4515264 : while (p < B2)
948 : {/* loop over probable primes p0 < p <= nextprime(B2), inserting giant
949 : * steps as necessary */
950 4514948 : p = snextpr(p, &np, &rcn, &bstp, uis2psp); /* next probable prime */
951 : /* work out the corresponding baby-step multiplier */
952 4514948 : k = bstp - (rcn < rcn0 ? 1 : 0);
953 4514948 : if (k > gss)
954 : { /* giant-step time, take gcd */
955 1114 : g = gcdii(g, N);
956 1114 : if (!is_pm1(g) && !equalii(g, N)) return g;
957 1107 : g = gen_1; set_avma(av1);
958 2214 : while (k > gss)
959 : { /* giant step */
960 1107 : if (DEBUGLEVEL >= 7) err_printf("\t(giant step at p = %lu)\n", p);
961 1107 : if (ecm_elladd0(N, &g, 64, 4, XG + i, YG + i,
962 0 : Xh, Yh, Xh, Yh) > 1) return g;
963 1107 : if (ecm_elladd0(N, &g, 64, 4, XG + i, YG + i,
964 : Xh + 64, Yh + 64, Xh + 64, Yh + 64) > 1)
965 0 : return g;
966 1107 : if (ecm_elladd0(N, &g, 64, 4, XG + i, YG + i,
967 : Xh + 128, Yh + 128, Xh + 128, Yh + 128) > 1)
968 0 : return g;
969 1107 : bstp -= (gss << 1);
970 1107 : k = bstp - (rcn < rcn0? 1: 0); /* recompute multiplier */
971 : }
972 : }
973 4514941 : if (!k) continue; /* point of interest is already in Xh */
974 4490276 : if (k < 0) k = -k;
975 4490276 : m = ((ulong)k - 1) << 2;
976 : /* accumulate product of differences of X coordinates */
977 4490276 : j = rcn<<2;
978 4490276 : avma = avtmp; /* go to garbage zone; don't use set_avma */
979 4490276 : g = modii(mulii(g, subii(XB[m], Xh[j])), N);
980 4490276 : g = modii(mulii(g, subii(XB[m+1], Xh[j+1])), N);
981 4490276 : g = modii(mulii(g, subii(XB[m+2], Xh[j+2])), N);
982 4490276 : g = mulii(g, subii(XB[m+3], Xh[j+3]));
983 4490276 : set_avma(av1);
984 4490276 : g = modii(g, N);
985 : }
986 316 : set_avma(av1);
987 : }
988 68 : return NULL;
989 : }
990 :
991 : /* ellfacteur() tuned to be useful as a first stage before MPQS, especially for
992 : * large arguments, when 'insist' is false, and now also for the case when
993 : * 'insist' is true, vaguely following suggestions by Paul Zimmermann
994 : * (http://www.loria.fr/~zimmerma/records/ecmnet.html). --GN 1998Jul,Aug */
995 : static GEN
996 3267 : ellfacteur(GEN N, int insist)
997 : {
998 3267 : const long size = expi(N) + 1;
999 3267 : pari_sp av = avma;
1000 : struct ECM E;
1001 3267 : long nbc, dsn, dsnmax, rep = 0;
1002 3267 : if (insist)
1003 : {
1004 0 : const long DSNMAX = numberof(TB1)-1;
1005 0 : dsnmax = (size >> 2) - 10;
1006 0 : if (dsnmax < 0) dsnmax = 0;
1007 0 : else if (dsnmax > DSNMAX) dsnmax = DSNMAX;
1008 0 : E.seed = 1 + (nbcmax<<7)*(size&0xffff); /* seed for choice of curves */
1009 :
1010 0 : dsn = (size >> 3) - 5;
1011 0 : if (dsn < 0) dsn = 0; else if (dsn > 47) dsn = 47;
1012 : /* pick up the torch where noninsistent stage would have given up */
1013 0 : nbc = dsn + (dsn >> 2) + 9; /* 8 or more curves in parallel */
1014 0 : nbc &= ~3; /* 4 | nbc */
1015 : }
1016 : else
1017 : {
1018 3267 : dsn = (size - 140) >> 3;
1019 3267 : if (dsn < 0)
1020 : {
1021 : #ifndef __EMX__ /* unless DOS/EMX: MPQS's disk access is abysmally slow */
1022 3210 : if (DEBUGLEVEL >= 4)
1023 0 : err_printf("ECM: number too small to justify this stage\n");
1024 3210 : return NULL; /* too small, decline the task */
1025 : #endif
1026 : dsn = 0;
1027 57 : } else if (dsn > 12) dsn = 12;
1028 57 : rep = (size <= 248 ?
1029 57 : (size <= 176 ? (size - 124) >> 4 : (size - 148) >> 3) :
1030 18 : (size - 224) >> 1);
1031 : #ifdef __EMX__ /* DOS/EMX: extra rounds (shun MPQS) */
1032 : rep += 20;
1033 : #endif
1034 57 : dsnmax = 72;
1035 : /* Use disjoint sets of curves for non-insist and insist phases; moreover,
1036 : * repeated calls acting on factors of the same original number should try
1037 : * to use fresh curves. The following achieves this */
1038 57 : E.seed = 1 + (nbcmax<<3)*(size & 0xf);
1039 57 : nbc = -1;
1040 : }
1041 57 : ECM_init(&E, N, nbc);
1042 57 : if (DEBUGLEVEL >= 4)
1043 : {
1044 0 : timer_start(&E.T);
1045 0 : err_printf("ECM: working on %ld curves at a time; initializing", E.nbc);
1046 0 : if (!insist)
1047 : {
1048 0 : if (rep == 1) err_printf(" for one round");
1049 0 : else err_printf(" for up to %ld rounds", rep);
1050 : }
1051 0 : err_printf("...\n");
1052 : }
1053 57 : if (dsn > dsnmax) dsn = dsnmax;
1054 : for(;;)
1055 36 : {
1056 93 : ulong B1 = insist? TB1[dsn]: TB1_for_stage[dsn];
1057 93 : GEN g = ECM_loop(&E, N, B1);
1058 93 : if (g)
1059 : {
1060 25 : if (DEBUGLEVEL >= 4)
1061 0 : err_printf("ECM: time = %6ld ms\n\tfound factor = %Ps\n",
1062 : timer_delay(&E.T), g);
1063 25 : return gerepilecopy(av, g);
1064 : }
1065 68 : if (dsn < dsnmax)
1066 : {
1067 68 : if (insist) dsn++;
1068 68 : else { dsn += 2; if (dsn > dsnmax) dsn = dsnmax; }
1069 : }
1070 68 : if (!insist && !--rep)
1071 : {
1072 32 : if (DEBUGLEVEL >= 4)
1073 0 : err_printf("ECM: time = %6ld ms,\tellfacteur giving up.\n",
1074 : timer_delay(&E.T));
1075 32 : return gc_NULL(av);
1076 : }
1077 : }
1078 : }
1079 : /* assume rounds >= 1, seed >= 1, B1 <= ULONG_MAX / 110 */
1080 : GEN
1081 0 : Z_ECM(GEN N, long rounds, long seed, ulong B1)
1082 : {
1083 0 : pari_sp av = avma;
1084 : struct ECM E;
1085 : long i;
1086 0 : E.seed = seed;
1087 0 : ECM_init(&E, N, -1);
1088 0 : if (DEBUGLEVEL >= 4) timer_start(&E.T);
1089 0 : for (i = rounds; i--; )
1090 : {
1091 0 : GEN g = ECM_loop(&E, N, B1);
1092 0 : if (g) return gerepilecopy(av, g);
1093 : }
1094 0 : return gc_NULL(av);
1095 : }
1096 :
1097 : /***********************************************************************/
1098 : /** **/
1099 : /** FACTORIZATION (Pollard-Brent rho) --GN1998Jun18-26 **/
1100 : /** pollardbrent() returns a nontrivial factor of n, assuming n is **/
1101 : /** composite and has no small prime divisor, or NULL if going on **/
1102 : /** would take more time than we want to spend. Sometimes it finds **/
1103 : /** more than one factor, and returns a structure suitable for **/
1104 : /** interpretation by ifac_crack. (Cf Algo 8.5.2 in ACiCNT) **/
1105 : /** **/
1106 : /***********************************************************************/
1107 : #define VALUE(x) gel(x,0)
1108 : #define EXPON(x) gel(x,1)
1109 : #define CLASS(x) gel(x,2)
1110 :
1111 : INLINE void
1112 51817 : INIT(GEN x, GEN v, GEN e, GEN c) {
1113 51817 : VALUE(x) = v;
1114 51817 : EXPON(x) = e;
1115 51817 : CLASS(x) = c;
1116 51817 : }
1117 : static void
1118 45877 : ifac_delete(GEN x) { INIT(x,NULL,NULL,NULL); }
1119 :
1120 : static void
1121 0 : rho_dbg(pari_timer *T, long c, long msg_mask)
1122 : {
1123 0 : if (c & msg_mask) return;
1124 0 : err_printf("Rho: time = %6ld ms,\t%3ld round%s\n",
1125 : timer_delay(T), c, (c==1?"":"s"));
1126 : }
1127 :
1128 : static void
1129 28246075 : one_iter(GEN *x, GEN *P, GEN x1, GEN n, long delta)
1130 : {
1131 28246075 : *x = addis(remii(sqri(*x), n), delta);
1132 28197686 : *P = modii(mulii(*P, subii(x1, *x)), n);
1133 28250691 : }
1134 : /* Return NULL when we run out of time, or a single t_INT containing a
1135 : * nontrivial factor of n, or a vector of t_INTs, each triple of successive
1136 : * entries containing a factor, an exponent (equal to one), and a factor
1137 : * class (NULL for unknown or zero for known composite), matching the
1138 : * internal representation used by the ifac_*() routines below. Repeated
1139 : * factors may arise; the caller will sort the factors anyway. Result
1140 : * is not gerepile-able (contains NULL) */
1141 : static GEN
1142 802 : pollardbrent_i(GEN n, long size, long c0, long retries)
1143 : {
1144 802 : long tf = lgefint(n), delta, msg_mask, c, k, k1, l;
1145 : pari_sp av;
1146 : GEN x, x1, y, P, g, g1, res;
1147 : pari_timer T;
1148 :
1149 802 : if (DEBUGLEVEL >= 4) timer_start(&T);
1150 802 : c = c0 << 5; /* 2^5 iterations per round */
1151 1604 : msg_mask = (size >= 448? 0x1fff:
1152 802 : (size >= 192? (256L<<((size-128)>>6))-1: 0xff));
1153 802 : y = cgeti(tf);
1154 802 : x1= cgeti(tf);
1155 802 : av = avma;
1156 :
1157 802 : PB_RETRY:
1158 : /* trick to make a 'random' choice determined by n. Don't use x^2+0 or
1159 : * x^2-2, ever. Don't use x^2-3 or x^2-7 with a starting value of 2.
1160 : * x^2+4, x^2+9 are affine conjugate to x^2+1, so don't use them either.
1161 : *
1162 : * (the point being that when we get called again on a composite cofactor
1163 : * of something we've already seen, we had better avoid the same delta) */
1164 802 : switch ((size + retries) & 7)
1165 : {
1166 107 : case 0: delta= 1; break;
1167 177 : case 1: delta= -1; break;
1168 96 : case 2: delta= 3; break;
1169 73 : case 3: delta= 5; break;
1170 72 : case 4: delta= -5; break;
1171 56 : case 5: delta= 7; break;
1172 137 : case 6: delta= 11; break;
1173 : /* case 7: */
1174 84 : default: delta=-11; break;
1175 : }
1176 802 : if (DEBUGLEVEL >= 4)
1177 : {
1178 0 : if (!retries)
1179 0 : err_printf("Rho: searching small factor of %ld-bit integer\n", size);
1180 : else
1181 0 : err_printf("Rho: restarting for remaining rounds...\n");
1182 0 : err_printf("Rho: using X^2%+1ld for up to %ld rounds of 32 iterations\n",
1183 : delta, c >> 5);
1184 : }
1185 802 : x = gen_2; P = gen_1; g1 = NULL; k = 1; l = 1;
1186 802 : affui(2, y);
1187 802 : affui(2, x1);
1188 : for (;;) /* terminated under the control of c */
1189 : { /* use the polynomial x^2 + delta */
1190 13205185 : one_iter(&x, &P, x1, n, delta);
1191 :
1192 13206750 : if ((--c & 0x1f)==0)
1193 : { /* one round complete */
1194 412593 : g = gcdii(n, P); if (!is_pm1(g)) goto fin;
1195 412326 : if (c <= 0)
1196 : { /* getting bored */
1197 396 : if (DEBUGLEVEL >= 4)
1198 0 : err_printf("Rho: time = %6ld ms,\tPollard-Brent giving up.\n",
1199 : timer_delay(&T));
1200 396 : return NULL;
1201 : }
1202 411930 : P = gen_1;
1203 411930 : if (DEBUGLEVEL >= 4) rho_dbg(&T, c0-(c>>5), msg_mask);
1204 411930 : affii(x,y); x = y; set_avma(av);
1205 : }
1206 :
1207 13204676 : if (--k) continue; /* normal end of loop body */
1208 :
1209 8612 : if (c & 0x1f) /* otherwise, we already checked */
1210 : {
1211 4812 : g = gcdii(n, P); if (!is_pm1(g)) goto fin;
1212 4812 : P = gen_1;
1213 : }
1214 :
1215 : /* Fast forward phase, doing l inner iterations without computing gcds.
1216 : * Check first whether it would take us beyond the alloted time.
1217 : * Fast forward rounds count only half (although they're taking
1218 : * more like 2/3 the time of normal rounds). This to counteract the
1219 : * nuisance that all c0 between 4096 and 6144 would act exactly as
1220 : * 4096; with the halving trick only the range 4096..5120 collapses
1221 : * (similarly for all other powers of two) */
1222 8612 : if ((c -= (l>>1)) <= 0)
1223 : { /* got bored */
1224 179 : if (DEBUGLEVEL >= 4)
1225 0 : err_printf("Rho: time = %6ld ms,\tPollard-Brent giving up.\n",
1226 : timer_delay(&T));
1227 179 : return NULL;
1228 : }
1229 8433 : c &= ~0x1f; /* keep it on multiples of 32 */
1230 :
1231 : /* Fast forward loop */
1232 8433 : affii(x, x1); set_avma(av); x = x1;
1233 8434 : k = l; l <<= 1;
1234 : /* don't show this for the first several (short) fast forward phases. */
1235 8434 : if (DEBUGLEVEL >= 4 && (l>>7) > msg_mask)
1236 0 : err_printf("Rho: fast forward phase (%ld rounds of 64)...\n", l>>7);
1237 15080029 : for (k1=k; k1; k1--)
1238 : {
1239 15072194 : one_iter(&x, &P, x1, n, delta);
1240 15069711 : if ((k1 & 0x1f) == 0) gerepileall(av, 2, &x, &P);
1241 : }
1242 7835 : if (DEBUGLEVEL >= 4 && (l>>7) > msg_mask)
1243 0 : err_printf("Rho: time = %6ld ms,\t%3ld rounds, back to normal mode\n",
1244 0 : timer_delay(&T), c0-(c>>5));
1245 7835 : affii(x,y); P = gerepileuptoint(av, P); x = y;
1246 : } /* forever */
1247 :
1248 227 : fin:
1249 : /* An accumulated gcd was > 1 */
1250 227 : if (!equalii(g,n))
1251 : { /* if it isn't n, and looks prime, return it */
1252 227 : if (MR_Jaeschke(g))
1253 : {
1254 227 : if (DEBUGLEVEL >= 4)
1255 : {
1256 0 : rho_dbg(&T, c0-(c>>5), 0);
1257 0 : err_printf("\tfound factor = %Ps\n",g);
1258 : }
1259 227 : return g;
1260 : }
1261 0 : set_avma(av); g1 = icopy(g); /* known composite, keep it safe */
1262 0 : av = avma;
1263 : }
1264 0 : else g1 = n; /* and work modulo g1 for backtracking */
1265 :
1266 : /* Here g1 is known composite */
1267 0 : if (DEBUGLEVEL >= 4 && size > 192)
1268 0 : err_printf("Rho: hang on a second, we got something here...\n");
1269 0 : x = y;
1270 : for(;;)
1271 : { /* backtrack until period recovered. Must terminate */
1272 0 : x = addis(remii(sqri(x), g1), delta);
1273 0 : g = gcdii(subii(x1, x), g1); if (!is_pm1(g)) break;
1274 :
1275 0 : if (DEBUGLEVEL >= 4 && (--c & 0x1f) == 0) rho_dbg(&T, c0-(c>>5), msg_mask);
1276 : }
1277 :
1278 0 : if (g1 == n || equalii(g,g1))
1279 : {
1280 0 : if (g1 == n && equalii(g,g1))
1281 : { /* out of luck */
1282 0 : if (DEBUGLEVEL >= 4)
1283 : {
1284 0 : rho_dbg(&T, c0-(c>>5), 0);
1285 0 : err_printf("\tPollard-Brent failed.\n");
1286 : }
1287 0 : if (++retries >= 4) pari_err_BUG("");
1288 0 : goto PB_RETRY;
1289 : }
1290 : /* half lucky: we've split n, but g1 equals either g or n */
1291 0 : if (DEBUGLEVEL >= 4)
1292 : {
1293 0 : rho_dbg(&T, c0-(c>>5), 0);
1294 0 : err_printf("\tfound %sfactor = %Ps\n", (g1!=n ? "composite " : ""), g);
1295 : }
1296 0 : res = cgetg(7, t_VEC);
1297 : /* g^1: known composite when g1!=n */
1298 0 : INIT(res+1, g, gen_1, (g1!=n? gen_0: NULL));
1299 : /* cofactor^1: status unknown */
1300 0 : INIT(res+4, diviiexact(n,g), gen_1, NULL);
1301 0 : return res;
1302 : }
1303 : /* g < g1 < n : our lucky day -- we've split g1, too */
1304 0 : res = cgetg(10, t_VEC);
1305 : /* unknown status for all three factors */
1306 0 : INIT(res+1, g, gen_1, NULL);
1307 0 : INIT(res+4, diviiexact(g1,g), gen_1, NULL);
1308 0 : INIT(res+7, diviiexact(n,g1), gen_1, NULL);
1309 0 : if (DEBUGLEVEL >= 4)
1310 : {
1311 0 : rho_dbg(&T, c0-(c>>5), 0);
1312 0 : err_printf("\tfound factors = %Ps, %Ps,\n\tand %Ps\n",
1313 0 : gel(res,1), gel(res,4), gel(res,7));
1314 : }
1315 0 : return res;
1316 : }
1317 : /* decline if n < 2^96 */
1318 : static GEN
1319 3494 : pollardbrent(GEN n)
1320 : {
1321 3494 : const long tune = 14; /* FIXME: retune this */
1322 3494 : long c0, size, tf = lgefint(n);
1323 : #ifdef LONG_IS_64BIT
1324 3045 : if (tf < 4 || (tf == 4 && uel(n,2) < (1UL << 32))) return NULL;
1325 : #else /* 32 bits */
1326 449 : if (tf < 5) return NULL;
1327 : #endif
1328 802 : size = expi(n) + 1;
1329 : /* nonlinear increase in effort, kicking in around 80 bits */
1330 802 : if (size <= 301) /* 301 gives 48121 + tune */
1331 795 : c0 = tune + size-60 + ((size-73)>>1)*((size-70)>>3)*((size-56)>>4);
1332 : else
1333 7 : c0 = 49152; /* ECM is faster when it'd take longer */
1334 802 : return pollardbrent_i(n, size, c0, 0);
1335 : }
1336 : GEN
1337 0 : Z_pollardbrent(GEN n, long rounds, long seed)
1338 : {
1339 0 : pari_sp av = avma;
1340 0 : GEN v = pollardbrent_i(n, expi(n)+1, rounds, seed);
1341 0 : if (!v) return NULL;
1342 0 : if (typ(v) == t_INT) v = mkvec2(v, diviiexact(n,v));
1343 0 : else if (lg(v) == 7) v = mkvec2(gel(v,1), gel(v,4));
1344 0 : else v = mkvec3(gel(v,1), gel(v,4), gel(v,7));
1345 0 : return gerepilecopy(av, v);
1346 : }
1347 :
1348 : /***********************************************************************/
1349 : /** FACTORIZATION (Shanks' SQUFOF) --GN2000Sep30-Oct01 **/
1350 : /** squfof() returns a proper factor of n, or NULL (failure). Assume **/
1351 : /** n is composite, not a square, and has no small prime divisors. **/
1352 : /** Works on two discriminants at once: n and 5n or 3n and 4n **/
1353 : /** Present implementation is limited to input <2^59, and works most **/
1354 : /** of the time in signed arithmetic on integers <2^31 in absolute **/
1355 : /** size. (Cf. Algo 8.7.2 in ACiCNT) **/
1356 : /***********************************************************************/
1357 :
1358 : /* squfof_ambig walks back along the ambiguous cycle until we hit an ambiguous
1359 : * form and thus the desired factor, which it returns. Returs 0 on failure.
1360 : *
1361 : * Input: a form (A, B, -C) with A = a^2, where a isn't blacklisted and
1362 : * (a, B) = 1. We should now proceed reducing the form (a, -B, -aC), but the
1363 : * first reduction step always sends this to (-aC, B, a), and the next one,
1364 : * with q computed as usual from B and a (occupying the c position), gives a
1365 : * reduced form, whose third member is easiest to recover by going back to D.
1366 : * From this point onwards, we're once again working with single-word numbers.
1367 : * No need to track signs, just work with the abs values of the coefficients.
1368 : * HACK: if LONG_IS_64BIT, D is actually a typecast long */
1369 : static long
1370 2091 : squfof_ambig(long a, long B, long dd, GEN D)
1371 : {
1372 : long b, c, q, qa, a0, b0, b1;
1373 2091 : long cnt = 0; /* count reduction steps on the cycle */
1374 :
1375 2091 : q = (dd + (B>>1)) / a;
1376 2091 : qa = q * a;
1377 2091 : b = (qa - B) + qa; /* avoid overflow */
1378 : #ifdef LONG_IS_64BIT
1379 1620 : c = (((long)D - b*b) >> 2) / a;
1380 : #else
1381 : {
1382 471 : pari_sp av = avma;
1383 471 : c = itos(divis(shifti(subii(D, sqrs(b)), -2), a));
1384 471 : set_avma(av);
1385 : }
1386 : #endif
1387 2091 : a0 = a; b0 = b1 = b; /* end of loop detection and safeguard */
1388 : for (;;)
1389 957157 : { /* reduction step */
1390 959248 : long c0 = c, qc, qcb;
1391 959248 : if (c0 > dd)
1392 267802 : q = 1;
1393 : else
1394 691446 : q = (dd + (b>>1)) / c0;
1395 959248 : if (q == 1)
1396 : {
1397 396866 : qcb = c0 - b; b = c0 + qcb; c = a - qcb;
1398 : }
1399 : else
1400 : {
1401 562382 : qc = q*c0; qcb = qc - b; b = qc + qcb; c = a - q*qcb;
1402 : }
1403 959248 : a = c0;
1404 :
1405 959248 : cnt++; if (b == b1) break;
1406 :
1407 : /* safeguard against infinite loop: we walked the cycle in vain.
1408 : * (I don't think this can actually happen.) */
1409 957157 : if (b == b0 && a == a0) return 0;
1410 :
1411 957157 : b1 = b;
1412 : }
1413 2091 : q = a&1 ? a : a>>1;
1414 2091 : if (DEBUGLEVEL >= 4)
1415 : {
1416 0 : if (q > 1)
1417 0 : err_printf("SQUFOF: found factor %ld from ambiguous form\n"
1418 : "\tafter %ld steps on the ambiguous cycle\n",
1419 0 : q / ugcd(q,15), cnt);
1420 : else
1421 0 : err_printf("SQUFOF: ...found nothing on the ambiguous cycle\n"
1422 : "\tafter %ld steps there\n", cnt);
1423 0 : if (DEBUGLEVEL >= 6) err_printf("SQUFOF: squfof_ambig returned %ld\n", q);
1424 : }
1425 2091 : return q;
1426 : }
1427 :
1428 : #define SQUFOF_BLACKLIST_SZ 64
1429 :
1430 : /* assume (n,30) = 1 */
1431 : static GEN
1432 5125 : squfof(GEN n)
1433 : {
1434 : ulong d1, d2;
1435 : #ifdef LONG_IS_64BIT
1436 : ulong uD1, uD2;
1437 : #endif
1438 5125 : long tf = lgefint(n), nm4, cnt = 0;
1439 : long a1, b1, c1, dd1, L1, a2, b2, c2, dd2, L2, a, q, c, qc, qcb;
1440 : GEN D1, D2;
1441 5125 : pari_sp av = avma;
1442 : long blacklist1[SQUFOF_BLACKLIST_SZ], blacklist2[SQUFOF_BLACKLIST_SZ];
1443 5125 : long blp1 = 0, blp2 = 0;
1444 5125 : int act1 = 1, act2 = 1;
1445 :
1446 : #ifdef LONG_IS_64BIT
1447 4323 : if (tf > 3 || (tf == 3 && uel(n,2) >= (1UL << 46)))
1448 : #else /* 32 bits */
1449 802 : if (tf > 4 || (tf == 4 && (ulong)(*int_MSW(n)) >= (1UL << 17))) /* 49 */
1450 : #endif
1451 3494 : return NULL; /* n too large */
1452 :
1453 : /* now we have 5 < n < 2^59 */
1454 1631 : nm4 = mod4(n);
1455 : #ifdef LONG_IS_64BIT
1456 1278 : if (nm4 == 1)
1457 : { /* n = 1 (mod4): run one iteration on D1 = n, another on D2 = 5n */
1458 600 : uD1 = n[2];
1459 600 : uD2 = 5 * n[2]; d2 = usqrt(uD2); dd2 = (long)((d2>>1) + (d2&1));
1460 600 : b2 = (long)((d2-1) | 1); /* b1, b2 will always stay odd */
1461 : }
1462 : else
1463 : { /* n = 3 (mod4): run one iteration on D1 = 3n, another on D2 = 4n */
1464 678 : uD1 = 3 * n[2];
1465 678 : uD2 = 4 * n[2]; dd2 = usqrt(n[2]); d2 = dd2 << 1;
1466 678 : b2 = (long)(d2 & (~1UL)); /* largest even below d2, will stay even */
1467 : }
1468 1278 : D1 = (GEN)uD1;
1469 1278 : D2 = (GEN)uD2;
1470 1278 : d1 = usqrt(uD1);
1471 1278 : b1 = (long)((d1-1) | 1); /* largest odd number not exceeding d1 */
1472 : /* c1 != 0 else n or 3n would be a square */
1473 1278 : c1 = (uD1 - b1*b1) / 4;
1474 : /* c2 != 0 else 5n would be a square */
1475 1278 : c2 = (uD2 - b2*b2) / 4;
1476 : #else
1477 353 : if (nm4 == 1)
1478 : { /* n = 1 (mod4): run one iteration on D1 = n, another on D2 = 5n */
1479 179 : D1 = n;
1480 179 : D2 = mului(5,n); d2 = itou(sqrti(D2)); dd2 = (long)((d2>>1) + (d2&1));
1481 179 : b2 = (long)((d2-1) | 1); /* b1, b2 will always stay odd */
1482 : }
1483 : else
1484 : { /* n = 3 (mod4): run one iteration on D1 = 3n, another on D2 = 4n */
1485 174 : D1 = mului(3,n);
1486 174 : D2 = shifti(n,2); dd2 = itou(sqrti(n)); d2 = dd2 << 1;
1487 174 : b2 = (long)(d2 & (~1UL)); /* largest even below d2, will stay even */
1488 : }
1489 353 : d1 = itou(sqrti(D1));
1490 353 : b1 = (long)((d1-1) | 1); /* largest odd number not exceeding d1 */
1491 : /* c1 != 0 else n or 3n would be a square */
1492 353 : c1 = itos(shifti(subii(D1, sqru((ulong)b1)), -2));
1493 : /* c2 != 0 else 5n would be a square */
1494 353 : c2 = itos(shifti(subii(D2, sqru((ulong)b2)), -2));
1495 : #endif
1496 1631 : L1 = (long)usqrt(d1);
1497 1631 : L2 = (long)usqrt(d2);
1498 : /* dd1 used to compute floor((d1+b1)/2) as dd1+floor(b1/2), without
1499 : * overflowing the 31bit signed integer size limit. Same for dd2. */
1500 1631 : dd1 = (long) ((d1>>1) + (d1&1));
1501 1631 : a1 = a2 = 1;
1502 :
1503 : /* The two (identity) forms (a1,b1,-c1) and (a2,b2,-c2) are now set up.
1504 : *
1505 : * a1 and c1 represent the absolute values of the a,c coefficients; we keep
1506 : * track of the sign separately, via the iteration counter cnt: when cnt is
1507 : * even, c < 0 and a > 0, else c > 0 and a < 0.
1508 : *
1509 : * L1, L2 are the limits for blacklisting small leading coefficients
1510 : * on the principal cycle, to guarantee that when we find a square form,
1511 : * its square root will belong to an ambiguous cycle, i.e. won't be an
1512 : * earlier form on the principal cycle.
1513 : *
1514 : * When n = 3(mod 4), D2 = 12(mod 16), and b^2 is 0 or 4 mod 16.
1515 : * It follows that 4*a*c must be 4 or 8 mod 16, respectively, so at most
1516 : * one of a,c can be divisible by 2 at most to the first power. This fact
1517 : * is used a couple of times below.
1518 : *
1519 : * The flags act1, act2 remain true while the respective cycle is still
1520 : * active; we drop them to false when we return to the identity form with-
1521 : * out having found a square form (or when the blacklist overflows, which
1522 : * shouldn't happen). */
1523 1631 : if (DEBUGLEVEL >= 4)
1524 0 : err_printf("SQUFOF: entering main loop with forms\n"
1525 : "\t(1, %ld, %ld) and (1, %ld, %ld)\n", b1, -c1, b2, -c2);
1526 :
1527 : /* MAIN LOOP: walk around the principal cycle looking for a square form.
1528 : * Blacklist small leading coefficients.
1529 : *
1530 : * The reduction operator can be computed entirely in 32-bit arithmetic:
1531 : * Let q = floor(floor((d1+b1)/2)/c1) (when c1>dd1, q=1, which happens
1532 : * often enough to special-case it). Then the new b1 = (q*c1-b1) + q*c1,
1533 : * which does not overflow, and the new c1 = a1 - q*(q*c1-b1), which is
1534 : * bounded by d1 in abs size since both the old and the new a1 are positive
1535 : * and bounded by d1. */
1536 1384926 : while (act1 || act2)
1537 : {
1538 1384926 : if (act1)
1539 : { /* send first form through reduction operator if active */
1540 1384926 : c = c1;
1541 1384926 : q = (c > dd1)? 1: (dd1 + (b1>>1)) / c;
1542 1384926 : if (q == 1)
1543 574343 : { qcb = c - b1; b1 = c + qcb; c1 = a1 - qcb; }
1544 : else
1545 810583 : { qc = q*c; qcb = qc - b1; b1 = qc + qcb; c1 = a1 - q*qcb; }
1546 1384926 : a1 = c;
1547 :
1548 1384926 : if (a1 <= L1)
1549 : { /* blacklist this */
1550 1080 : if (blp1 >= SQUFOF_BLACKLIST_SZ) /* overflows: shouldn't happen */
1551 0 : act1 = 0;
1552 : else
1553 : {
1554 1080 : if (DEBUGLEVEL >= 6)
1555 0 : err_printf("SQUFOF: blacklisting a = %ld on first cycle\n", a1);
1556 1080 : blacklist1[blp1++] = a1;
1557 : }
1558 : }
1559 : }
1560 1384926 : if (act2)
1561 : { /* send second form through reduction operator if active */
1562 1384734 : c = c2;
1563 1384734 : q = (c > dd2)? 1: (dd2 + (b2>>1)) / c;
1564 1384734 : if (q == 1)
1565 574193 : { qcb = c - b2; b2 = c + qcb; c2 = a2 - qcb; }
1566 : else
1567 810541 : { qc = q*c; qcb = qc - b2; b2 = qc + qcb; c2 = a2 - q*qcb; }
1568 1384734 : a2 = c;
1569 :
1570 1384734 : if (a2 <= L2)
1571 : { /* blacklist this */
1572 1108 : if (blp2 >= SQUFOF_BLACKLIST_SZ) /* overflows: shouldn't happen */
1573 0 : act2 = 0;
1574 : else
1575 : {
1576 1108 : if (DEBUGLEVEL >= 6)
1577 0 : err_printf("SQUFOF: blacklisting a = %ld on second cycle\n", a2);
1578 1108 : blacklist2[blp2++] = a2;
1579 : }
1580 : }
1581 : }
1582 :
1583 1384926 : if (++cnt & 1) continue; /* odd iteration */
1584 : /* even iteration: the leading coefficients are positive */
1585 :
1586 : /* examine first form if active */
1587 692463 : if (act1 && a1 == 1) /* back to identity */
1588 : { /* drop this discriminant */
1589 0 : act1 = 0;
1590 0 : if (DEBUGLEVEL >= 4)
1591 0 : err_printf("SQUFOF: first cycle exhausted after %ld iterations,\n"
1592 : "\tdropping it\n", cnt);
1593 : }
1594 692463 : if (act1)
1595 : {
1596 692463 : if (uissquareall((ulong)a1, (ulong*)&a))
1597 : { /* square form */
1598 1309 : if (DEBUGLEVEL >= 4)
1599 0 : err_printf("SQUFOF: square form (%ld^2, %ld, %ld) on first cycle\n"
1600 : "\tafter %ld iterations\n", a, b1, -c1, cnt);
1601 1309 : if (a <= L1)
1602 : { /* blacklisted? */
1603 : long j;
1604 2436 : for (j = 0; j < blp1; j++)
1605 1613 : if (a == blacklist1[j]) { a = 0; break; }
1606 : }
1607 1309 : if (a > 0)
1608 : { /* not blacklisted */
1609 823 : q = ugcd(a, b1); /* imprimitive form? */
1610 823 : if (q > 1)
1611 : { /* q^2 divides D1 hence n [ assuming n % 3 != 0 ] */
1612 0 : set_avma(av);
1613 0 : if (DEBUGLEVEL >= 4) err_printf("SQUFOF: found factor %ld^2\n", q);
1614 0 : return mkvec3(utoipos(q), gen_2, NULL);/* q^2, unknown status */
1615 : }
1616 : /* chase the inverse root form back along the ambiguous cycle */
1617 823 : q = squfof_ambig(a, b1, dd1, D1);
1618 823 : if (q > 3)
1619 : {
1620 677 : if (nm4 == 3 && q % 3 == 0) q /= 3;
1621 677 : return gc_utoipos(av, q); /* SUCCESS! */
1622 : }
1623 : }
1624 486 : else if (DEBUGLEVEL >= 4) /* blacklisted */
1625 0 : err_printf("SQUFOF: ...but the root form seems to be on the "
1626 : "principal cycle\n");
1627 : }
1628 : }
1629 :
1630 : /* examine second form if active */
1631 691786 : if (act2 && a2 == 1) /* back to identity form */
1632 : { /* drop this discriminant */
1633 2 : act2 = 0;
1634 2 : if (DEBUGLEVEL >= 4)
1635 0 : err_printf("SQUFOF: second cycle exhausted after %ld iterations,\n"
1636 : "\tdropping it\n", cnt);
1637 : }
1638 691786 : if (act2)
1639 : {
1640 691690 : if (uissquareall((ulong)a2, (ulong*)&a))
1641 : { /* square form */
1642 1549 : if (DEBUGLEVEL >= 4)
1643 0 : err_printf("SQUFOF: square form (%ld^2, %ld, %ld) on second cycle\n"
1644 : "\tafter %ld iterations\n", a, b2, -c2, cnt);
1645 1549 : if (a <= L2)
1646 : { /* blacklisted? */
1647 : long j;
1648 2729 : for (j = 0; j < blp2; j++)
1649 1461 : if (a == blacklist2[j]) { a = 0; break; }
1650 : }
1651 1549 : if (a > 0)
1652 : { /* not blacklisted */
1653 1268 : q = ugcd(a, b2); /* imprimitive form? */
1654 : /* NB if b2 is even, a is odd, so the gcd is always odd */
1655 1268 : if (q > 1)
1656 : { /* q^2 divides D2 hence n [ assuming n % 5 != 0 ] */
1657 0 : set_avma(av);
1658 0 : if (DEBUGLEVEL >= 4) err_printf("SQUFOF: found factor %ld^2\n", q);
1659 0 : return mkvec3(utoipos(q), gen_2, NULL);/* q^2, unknown status */
1660 : }
1661 : /* chase the inverse root form along the ambiguous cycle */
1662 1268 : q = squfof_ambig(a, b2, dd2, D2);
1663 1268 : if (q > 5)
1664 : {
1665 954 : if (nm4 == 1 && q % 5 == 0) q /= 5;
1666 954 : return gc_utoipos(av, q); /* SUCCESS! */
1667 : }
1668 : }
1669 281 : else if (DEBUGLEVEL >= 4) /* blacklisted */
1670 0 : err_printf("SQUFOF: ...but the root form seems to be on the "
1671 : "principal cycle\n");
1672 : }
1673 : }
1674 : } /* end main loop */
1675 :
1676 0 : if (DEBUGLEVEL>=4) err_printf("SQUFOF: giving up\n");
1677 0 : return gc_NULL(av);
1678 : }
1679 :
1680 : /***********************************************************************/
1681 : /* DETECTING ODD POWERS --GN1998Jun28 */
1682 : /* Factoring engines like MPQS which ultimately rely on computing */
1683 : /* gcd(N, x^2-y^2) to find a nontrivial factor of N can't split */
1684 : /* N = p^k for an odd prime p, since (Z/p^k)^* is then cyclic. Here */
1685 : /* is an analogue of Z_issquareall() for 3rd, 5th and 7th powers. */
1686 : /* The general case is handled by is_kth_power */
1687 : /***********************************************************************/
1688 :
1689 : /* Multistage sieve. First stages work mod 211, 209, 61, 203 in this order
1690 : * (first reduce mod the product of these and then take the remainder apart).
1691 : * Second stages use 117, 31, 43, 71. Moduli which are no longer interesting
1692 : * are skipped. Everything is encoded in a table of 106 24-bit masks. We only
1693 : * need the first half of the residues. Three bits per modulus indicate which
1694 : * residues are 7th (bit 2), 5th (bit 1) or 3rd (bit 0) powers; the eight
1695 : * moduli above are assigned right-to-left. The table was generated using: */
1696 :
1697 : #if 0
1698 : L = [71, 43, 31, [O(3^2),O(13)], [O(7),O(29)], 61, [O(11),O(19)], 211];
1699 : ispow(x, N, k)=
1700 : {
1701 : if (type(N) == "t_INT", return (ispower(Mod(x,N), k)));
1702 : for (i = 1, #N, if (!ispower(x + N[i], k), return (0))); 1
1703 : }
1704 : check(r) =
1705 : {
1706 : print1(" 0");
1707 : for (i=1,#L,
1708 : N = 0;
1709 : if (ispow(r, L[i], 3), N += 1);
1710 : if (ispow(r, L[i], 5), N += 2);
1711 : if (ispow(r, L[i], 7), N += 4);
1712 : print1(N);
1713 : ); print("ul, /* ", r, " */")
1714 : }
1715 : for (r = 0, 105, check(r))
1716 : #endif
1717 : static ulong powersmod[106] = {
1718 : 077777777ul, /* 0 */
1719 : 077777777ul, /* 1 */
1720 : 013562440ul, /* 2 */
1721 : 012402540ul, /* 3 */
1722 : 013562440ul, /* 4 */
1723 : 052662441ul, /* 5 */
1724 : 016603440ul, /* 6 */
1725 : 016463450ul, /* 7 */
1726 : 013573551ul, /* 8 */
1727 : 012462540ul, /* 9 */
1728 : 012462464ul, /* 10 */
1729 : 013462771ul, /* 11 */
1730 : 012406473ul, /* 12 */
1731 : 012463641ul, /* 13 */
1732 : 052463646ul, /* 14 */
1733 : 012503446ul, /* 15 */
1734 : 013562440ul, /* 16 */
1735 : 052466440ul, /* 17 */
1736 : 012472451ul, /* 18 */
1737 : 012462454ul, /* 19 */
1738 : 032463550ul, /* 20 */
1739 : 013403664ul, /* 21 */
1740 : 013463460ul, /* 22 */
1741 : 032562565ul, /* 23 */
1742 : 012402540ul, /* 24 */
1743 : 052662441ul, /* 25 */
1744 : 032672452ul, /* 26 */
1745 : 013573551ul, /* 27 */
1746 : 012467541ul, /* 28 */
1747 : 012567640ul, /* 29 */
1748 : 032706450ul, /* 30 */
1749 : 012762452ul, /* 31 */
1750 : 033762662ul, /* 32 */
1751 : 012502562ul, /* 33 */
1752 : 032463562ul, /* 34 */
1753 : 013563440ul, /* 35 */
1754 : 016663440ul, /* 36 */
1755 : 036662550ul, /* 37 */
1756 : 012462552ul, /* 38 */
1757 : 033502450ul, /* 39 */
1758 : 012462643ul, /* 40 */
1759 : 033467540ul, /* 41 */
1760 : 017403441ul, /* 42 */
1761 : 017463462ul, /* 43 */
1762 : 017472460ul, /* 44 */
1763 : 033462470ul, /* 45 */
1764 : 052566450ul, /* 46 */
1765 : 013562640ul, /* 47 */
1766 : 032403640ul, /* 48 */
1767 : 016463450ul, /* 49 */
1768 : 016463752ul, /* 50 */
1769 : 033402440ul, /* 51 */
1770 : 012462540ul, /* 52 */
1771 : 012472540ul, /* 53 */
1772 : 053562462ul, /* 54 */
1773 : 012463465ul, /* 55 */
1774 : 012663470ul, /* 56 */
1775 : 052607450ul, /* 57 */
1776 : 012566553ul, /* 58 */
1777 : 013466440ul, /* 59 */
1778 : 012502741ul, /* 60 */
1779 : 012762744ul, /* 61 */
1780 : 012763740ul, /* 62 */
1781 : 012763443ul, /* 63 */
1782 : 013573551ul, /* 64 */
1783 : 013462471ul, /* 65 */
1784 : 052502460ul, /* 66 */
1785 : 012662463ul, /* 67 */
1786 : 012662451ul, /* 68 */
1787 : 012403550ul, /* 69 */
1788 : 073567540ul, /* 70 */
1789 : 072463445ul, /* 71 */
1790 : 072462740ul, /* 72 */
1791 : 012472442ul, /* 73 */
1792 : 012462644ul, /* 74 */
1793 : 013406650ul, /* 75 */
1794 : 052463471ul, /* 76 */
1795 : 012563474ul, /* 77 */
1796 : 013503460ul, /* 78 */
1797 : 016462441ul, /* 79 */
1798 : 016462440ul, /* 80 */
1799 : 012462540ul, /* 81 */
1800 : 013462641ul, /* 82 */
1801 : 012463454ul, /* 83 */
1802 : 013403550ul, /* 84 */
1803 : 057563540ul, /* 85 */
1804 : 017466441ul, /* 86 */
1805 : 017606471ul, /* 87 */
1806 : 053666573ul, /* 88 */
1807 : 012562561ul, /* 89 */
1808 : 013473641ul, /* 90 */
1809 : 032573440ul, /* 91 */
1810 : 016763440ul, /* 92 */
1811 : 016702640ul, /* 93 */
1812 : 033762552ul, /* 94 */
1813 : 012562550ul, /* 95 */
1814 : 052402451ul, /* 96 */
1815 : 033563441ul, /* 97 */
1816 : 012663561ul, /* 98 */
1817 : 012677560ul, /* 99 */
1818 : 012462464ul, /* 100 */
1819 : 032562642ul, /* 101 */
1820 : 013402551ul, /* 102 */
1821 : 032462450ul, /* 103 */
1822 : 012467445ul, /* 104 */
1823 : 032403440ul, /* 105 */
1824 : };
1825 :
1826 : static int
1827 3928662 : check_res(ulong x, ulong N, int shift, ulong *mask)
1828 : {
1829 3928662 : long r = x%N; if ((ulong)r> (N>>1)) r = N - r;
1830 3928662 : *mask &= (powersmod[r] >> shift);
1831 3928662 : return *mask;
1832 : }
1833 :
1834 : /* is x mod 211*209*61*203*117*31*43*71 a 3rd, 5th or 7th power ? */
1835 : int
1836 2422077 : uis_357_powermod(ulong x, ulong *mask)
1837 : {
1838 2422077 : if ( !check_res(x, 211UL, 0, mask)) return 0;
1839 978599 : if (*mask & 3 && !check_res(x, 209UL, 3, mask)) return 0;
1840 399592 : if (*mask & 3 && !check_res(x, 61UL, 6, mask)) return 0;
1841 222296 : if (*mask & 5 && !check_res(x, 203UL, 9, mask)) return 0;
1842 56862 : if (*mask & 1 && !check_res(x, 117UL,12, mask)) return 0;
1843 32886 : if (*mask & 3 && !check_res(x, 31UL,15, mask)) return 0;
1844 24697 : if (*mask & 5 && !check_res(x, 43UL,18, mask)) return 0;
1845 7452 : if (*mask & 6 && !check_res(x, 71UL,21, mask)) return 0;
1846 3786 : return 1;
1847 : }
1848 : /* asume x > 0 and pt != NULL */
1849 : int
1850 2366983 : uis_357_power(ulong x, ulong *pt, ulong *mask)
1851 : {
1852 : double logx;
1853 2366983 : if (!odd(x))
1854 : {
1855 9081 : long v = vals(x);
1856 9081 : if (v % 7) *mask &= ~4;
1857 9081 : if (v % 5) *mask &= ~2;
1858 9081 : if (v % 3) *mask &= ~1;
1859 9081 : if (!*mask) return 0;
1860 : }
1861 2359390 : if (!uis_357_powermod(x, mask)) return 0;
1862 2981 : logx = log((double)x);
1863 3853 : while (*mask)
1864 : {
1865 : long e, b;
1866 : ulong y, ye;
1867 2981 : if (*mask & 1) { b = 1; e = 3; }
1868 873 : else if (*mask & 2) { b = 2; e = 5; }
1869 355 : else { b = 4; e = 7; }
1870 2981 : y = (ulong)(exp(logx / e) + 0.5);
1871 2981 : ye = upowuu(y,e);
1872 2981 : if (ye == x) { *pt = y; return e; }
1873 : #ifdef LONG_IS_64BIT
1874 750 : if (ye > x) y--; else y++;
1875 750 : ye = upowuu(y,e);
1876 750 : if (ye == x) { *pt = y; return e; }
1877 : #endif
1878 872 : *mask &= ~b; /* turn the bit off */
1879 : }
1880 872 : return 0;
1881 : }
1882 :
1883 : #ifndef LONG_IS_64BIT
1884 : /* as above, split in two functions */
1885 : /* is x mod 211*209*61*203 a 3rd, 5th or 7th power ? */
1886 : static int
1887 14433 : uis_357_powermod_32bit_1(ulong x, ulong *mask)
1888 : {
1889 14433 : if ( !check_res(x, 211UL, 0, mask)) return 0;
1890 8008 : if (*mask & 3 && !check_res(x, 209UL, 3, mask)) return 0;
1891 4115 : if (*mask & 3 && !check_res(x, 61UL, 6, mask)) return 0;
1892 2914 : if (*mask & 5 && !check_res(x, 203UL, 9, mask)) return 0;
1893 828 : return 1;
1894 : }
1895 : /* is x mod 117*31*43*71 a 3rd, 5th or 7th power ? */
1896 : static int
1897 828 : uis_357_powermod_32bit_2(ulong x, ulong *mask)
1898 : {
1899 828 : if (*mask & 1 && !check_res(x, 117UL,12, mask)) return 0;
1900 656 : if (*mask & 3 && !check_res(x, 31UL,15, mask)) return 0;
1901 538 : if (*mask & 5 && !check_res(x, 43UL,18, mask)) return 0;
1902 266 : if (*mask & 6 && !check_res(x, 71UL,21, mask)) return 0;
1903 176 : return 1;
1904 : }
1905 : #endif
1906 :
1907 : /* Returns 3, 5, or 7 if x is a cube (but not a 5th or 7th power), a 5th
1908 : * power (but not a 7th), or a 7th power, and in this case creates the
1909 : * base on the stack and assigns its address to *pt. Otherwise returns 0.
1910 : * x must be of type t_INT and positive; this is not checked. The *mask
1911 : * argument tells us which things to check -- bit 0: 3rd, bit 1: 5th,
1912 : * bit 2: 7th pwr; set a bit to have the corresponding power examined --
1913 : * and is updated appropriately for a possible follow-up call */
1914 : int
1915 2801660 : is_357_power(GEN x, GEN *pt, ulong *mask)
1916 : {
1917 2801660 : long lx = lgefint(x);
1918 : ulong r;
1919 : pari_sp av;
1920 : GEN y;
1921 :
1922 2801660 : if (!*mask) return 0; /* useful when running in a loop */
1923 2430273 : if (DEBUGLEVEL>4)
1924 0 : err_printf("OddPwrs: examining %ld-bit integer\n", expi(x)+1);
1925 2430273 : if (lgefint(x) == 3) {
1926 : ulong t;
1927 2353154 : long e = uis_357_power(x[2], &t, mask);
1928 2353154 : if (e)
1929 : {
1930 2083 : if (pt) *pt = utoi(t);
1931 2083 : return e;
1932 : }
1933 2351071 : return 0;
1934 : }
1935 : #ifdef LONG_IS_64BIT
1936 62686 : r = (lx == 3)? uel(x,2): umodiu(x, 6046846918939827UL);
1937 62686 : if (!uis_357_powermod(r, mask)) return 0;
1938 : #else
1939 14433 : r = (lx == 3)? uel(x,2): umodiu(x, 211*209*61*203);
1940 14433 : if (!uis_357_powermod_32bit_1(r, mask)) return 0;
1941 828 : r = (lx == 3)? uel(x,2): umodiu(x, 117*31*43*71);
1942 828 : if (!uis_357_powermod_32bit_2(r, mask)) return 0;
1943 : #endif
1944 981 : av = avma;
1945 1602 : while (*mask)
1946 : {
1947 : long e, b;
1948 : /* priority to higher powers: if we have a 21st, it is easier to rediscover
1949 : * that its 7th root is a cube than that its cube root is a 7th power */
1950 1303 : if (*mask & 4) { b = 4; e = 7; }
1951 963 : else if (*mask & 2) { b = 2; e = 5; }
1952 366 : else { b = 1; e = 3; }
1953 1303 : y = mpround( sqrtnr(itor(x, nbits2prec(64 + bit_accuracy(lx) / e)), e) );
1954 1303 : if (equalii(powiu(y,e), x))
1955 : {
1956 682 : if (!pt) return gc_int(av,e);
1957 668 : set_avma((pari_sp)y); *pt = gerepileuptoint(av, y);
1958 668 : return e;
1959 : }
1960 621 : *mask &= ~b; /* turn the bit off */
1961 621 : set_avma(av);
1962 : }
1963 299 : return 0;
1964 : }
1965 :
1966 : /* Is x a n-th power ? If pt != NULL, it receives the n-th root */
1967 : ulong
1968 471038 : is_kth_power(GEN x, ulong n, GEN *pt)
1969 : {
1970 : forprime_t T;
1971 : long j;
1972 : ulong q, residue;
1973 : GEN y;
1974 471038 : pari_sp av = avma;
1975 :
1976 471038 : (void)u_forprime_arith_init(&T, odd(n)? 2*n+1: n+1, ULONG_MAX, 1,n);
1977 : /* we'll start at q, smallest prime >= n */
1978 :
1979 : /* Modular checks, use small primes q congruent 1 mod n */
1980 : /* A non n-th power nevertheless passes the test with proba n^(-#checks),
1981 : * We'd like this < 1e-6 but let j = floor(log(1e-6) / log(n)) which
1982 : * ensures much less. */
1983 470964 : if (n < 16)
1984 118068 : j = 5;
1985 352896 : else if (n < 32)
1986 154522 : j = 4;
1987 198374 : else if (n < 101)
1988 176423 : j = 3;
1989 21951 : else if (n < 1001)
1990 5390 : j = 2;
1991 16561 : else if (n < 17886697) /* smallest such that smallest suitable q is > 2^32 */
1992 16303 : j = 1;
1993 : else
1994 275 : j = 0;
1995 514422 : for (; j > 0; j--)
1996 : {
1997 511806 : if (!(q = u_forprime_next(&T))) break;
1998 : /* q a prime = 1 mod n */
1999 511899 : residue = umodiu(x, q);
2000 511893 : if (residue == 0)
2001 : {
2002 483 : if (Z_lval(x,q) % n) return gc_ulong(av,0);
2003 49 : continue;
2004 : }
2005 : /* n-th power mod q ? */
2006 511410 : if (Fl_powu(residue, (q-1)/n, q) != 1) return gc_ulong(av,0);
2007 : }
2008 2616 : set_avma(av);
2009 :
2010 2611 : if (DEBUGLEVEL>4) err_printf("\nOddPwrs: [%lu] passed modular checks\n",n);
2011 : /* go to the horse's mouth... */
2012 2611 : y = roundr( sqrtnr(itor(x, nbits2prec(expi(x)/n + 16)), n) );
2013 2611 : if (!equalii(powiu(y, n), x)) {
2014 1662 : if (DEBUGLEVEL>4) err_printf("\tBut it wasn't a pure power.\n");
2015 1662 : return gc_ulong(av,0);
2016 : }
2017 949 : if (!pt) set_avma(av); else { set_avma((pari_sp)y); *pt = gerepileuptoint(av,y); }
2018 949 : return 1;
2019 : }
2020 :
2021 : /* is x a p^i-th power, p >= 11 prime ? Similar to is_357_power(), but instead
2022 : * of the mask, we keep the current test exponent around. Cut off when
2023 : * log_2 x^(1/k) < cutoffbits since we would have found it by trial division.
2024 : * Everything needed here (primitive roots etc.) is computed from scratch on
2025 : * the fly; compared to the size of numbers under consideration, these
2026 : * word-sized computations take negligible time.
2027 : * Any cutoffbits > 0 is safe, but direct root extraction attempts are faster
2028 : * when trial division has been used to discover very small bases. We become
2029 : * competitive at cutoffbits ~ 10 */
2030 : int
2031 69277 : is_pth_power(GEN x, GEN *pt, forprime_t *T, ulong cutoffbits)
2032 : {
2033 69277 : long cnt=0, size = expi(x) /* not +1 */;
2034 : ulong p;
2035 69277 : pari_sp av = avma;
2036 486016 : while ((p = u_forprime_next(T)) && size/p >= cutoffbits) {
2037 416770 : long v = 1;
2038 416770 : if (DEBUGLEVEL>5 && cnt++==2000)
2039 0 : { cnt=0; err_printf("%lu%% ", 100*p*cutoffbits/size); }
2040 416828 : while (is_kth_power(x, p, pt)) {
2041 56 : v *= p; x = *pt;
2042 56 : size = expi(x);
2043 : }
2044 416781 : if (v > 1)
2045 : {
2046 42 : if (DEBUGLEVEL>5) err_printf("\nOddPwrs: is a %ld power\n",v);
2047 42 : return v;
2048 : }
2049 : }
2050 69224 : if (DEBUGLEVEL>5) err_printf("\nOddPwrs: not a power\n",p);
2051 69224 : return gc_int(av,0); /* give up */
2052 : }
2053 :
2054 : /***********************************************************************/
2055 : /** FACTORIZATION (master iteration) **/
2056 : /** Driver for the various methods of finding large factors **/
2057 : /** (after trial division has cast out the very small ones). **/
2058 : /** GN1998Jun24--30 **/
2059 : /***********************************************************************/
2060 :
2061 : /* Direct use:
2062 : * ifac_start_hint(n,moebius,hint) registers with the iterative factorizer
2063 : * - an integer n (without prime factors < tridiv_bound(n))
2064 : * - registers whether or not we should terminate early if we find a square
2065 : * factor,
2066 : * - a hint about which method(s) to use.
2067 : * This must always be called first. If input is not composite, oo loop.
2068 : * The routine decomposes n nontrivially into a product of two factors except
2069 : * in squarefreeness ('Moebius') mode.
2070 : *
2071 : * ifac_start(n,moebius) same using default hint.
2072 : *
2073 : * ifac_primary_factor() returns a prime divisor (not necessarily the
2074 : * smallest) and the corresponding exponent.
2075 : *
2076 : * Encapsulated user interface: Many arithmetic functions have a 'contributor'
2077 : * ifac_xxx, to be called on any large composite cofactor left over after trial
2078 : * division by small primes: xxx is one of moebius, issquarefree, totient, etc.
2079 : *
2080 : * We never test whether the input number is prime or composite, since
2081 : * presumably it will have come out of the small factors finder stage
2082 : * (which doesn't really exist yet but which will test the left-over
2083 : * cofactor for primality once it does). */
2084 :
2085 : /* The data structure in which we preserve whatever we know about our number N
2086 : * is kept on the PARI stack, and updated as needed.
2087 : * This makes the machinery re-entrant, and avoids memory leaks when a lengthy
2088 : * factorization is interrupted. We try to keep the whole affair connected,
2089 : * and the parent object is always older than its children. This may in
2090 : * rare cases lead to some extra copying around, and knowing what is garbage
2091 : * at any given time is not trivial. See below for examples how to do it right.
2092 : * (Connectedness is destroyed if callers of ifac_main() create stuff on the
2093 : * stack in between calls. This is harmless as long as ifac_realloc() is used
2094 : * to re-create a connected object at the head of the stack just before
2095 : * collecting garbage.)
2096 : * A t_INT may well have > 10^6 distinct prime factors larger than 2^16. Since
2097 : * we need not find factors in order of increasing size, we must be prepared to
2098 : * drag a very large amount of data around. We start with a small structure
2099 : * and extend it when necessary. */
2100 :
2101 : /* The idea of the algorithm is:
2102 : * Let N0 be whatever is currently left of N after dividing off all the
2103 : * prime powers we have already returned to the caller. Then we maintain
2104 : * N0 as a product
2105 : * (1) N0 = \prod_i P_i^{e_i} * \prod_j Q_j^{f_j} * \prod_k C_k^{g_k}
2106 : * where the P_i and Q_j are distinct primes, each C_k is known composite,
2107 : * none of the P_i divides any C_k, and we also know the total ordering
2108 : * of all the P_i, Q_j and C_k; in particular, we will never try to divide
2109 : * a C_k by a larger Q_j. Some of the C_k may have common factors.
2110 : *
2111 : * Caveat implementor: Taking gcds among C_k's is very likely to cost at
2112 : * least as much time as dividing off any primes as we find them, and book-
2113 : * keeping would be tough (since D=gcd(C_1,C_2) can still have common factors
2114 : * with both C_1/D and C_2/D, and so on...).
2115 : *
2116 : * At startup, we just initialize the structure to
2117 : * (2) N = C_1^1 (composite).
2118 : *
2119 : * Whenever ifac_primary_factor() or one of the arithmetic user interface
2120 : * routines needs a primary factor, and the smallest thing in our list is P_1,
2121 : * we return that and its exponent, and remove it from our list. (When nothing
2122 : * is left, we return a sentinel value -- gen_1. And in Moebius mode, when we
2123 : * see something with exponent > 1, whether prime or composite, we return gen_0
2124 : * or 0, depending on the function). In all other cases, ifac_main() iterates
2125 : * the following steps until we have a P_1 in the smallest position.
2126 : *
2127 : * When the smallest item is C_1, as it is initially:
2128 : * (3.1) Crack C_1 into a nontrivial product U_1 * U_2 by whatever method
2129 : * comes to mind for this size. (U for 'unknown'.) Cracking will detect
2130 : * perfect powers, so we may instead see a power of some U_1 here, or even
2131 : * something of the form U_1^k*U_2^k; of course the exponent already attached
2132 : * to C_1 is taken into account in the following.
2133 : * (3.2) If we have U_1*U_2, sort the two factors (distinct: squares are caught
2134 : * in stage 3.1). N.B. U_1 and U_2 are smaller than anything else in our list.
2135 : * (3.3) Check U_1 and U_2 for primality, and flag them accordingly.
2136 : * (3.4) Iterate.
2137 : *
2138 : * When the smallest item is Q_1:
2139 : * This is the unpleasant case. We go through the entire list and try to
2140 : * divide Q_1 off each of the current C_k's, which usually fails, but may
2141 : * succeed several times. When a division was successful, the corresponding
2142 : * C_k is removed from our list, and the cofactor becomes a U_l for the moment
2143 : * unless it is 1 (which happens when C_k was a power of Q_1). When we're
2144 : * through we upgrade Q_1 to P_1 status, then do a primality check on each U_l
2145 : * and sort it back into the list either as a Q_j or as a C_k. If during the
2146 : * insertion sort we discover that some U_l equals some P_i or Q_j or C_k we
2147 : * already have, we just add U_l's exponent to that of its twin. (The sorting
2148 : * therefore happens before the primality test). Since this may produce one or
2149 : * more elements smaller than the P_1 we just confirmed, we may have to repeat
2150 : * the iteration.
2151 : * A trick avoids some Q_1 instances: just after the sweep classifying
2152 : * all current unknowns as either composites or primes, we do another downward
2153 : * sweep beginning with the largest current factor and stopping just above the
2154 : * largest current composite. Every Q_j we pass is turned into a P_i.
2155 : * (Different primes are automatically coprime among each other, and primes do
2156 : * not divide smaller composites.)
2157 : * NB: We have no use for comparing the square of a prime to N0. Normally
2158 : * we will get called after casting out only the smallest primes, and
2159 : * since we cannot guarantee that we see the large prime factors in as-
2160 : * cending order, we cannot stop when we find one larger than sqrt(N0). */
2161 :
2162 : /* Data structure: We keep everything in a single t_VEC of t_INTs. The
2163 : * first 2 components are read-only:
2164 : * 1) the first records whether we're doing full (NULL) or Moebius (gen_1)
2165 : * factorization; in the latter case subroutines return a sentinel value as
2166 : * soon as they spot an exponent > 1.
2167 : * 2) the second records the hint from factorint()'s optional flag, for use by
2168 : * ifac_crack().
2169 : *
2170 : * The remaining components (initially 15) are used in groups of three:
2171 : * [ factor (t_INT), exponent (t_INT), factor class ], where factor class is
2172 : * NULL : unknown
2173 : * gen_0: known composite C_k
2174 : * gen_1: known prime Q_j awaiting trial division
2175 : * gen_2: finished prime P_i.
2176 : * When during the division stage we re-sort a C_k-turned-U_l to a lower
2177 : * position, we rotate any intervening material upward towards its old
2178 : * slot. When a C_k was divided down to 1, its slot is left empty at
2179 : * first; similarly when the re-sorting detects a repeated factor.
2180 : * After the sorting phase, we de-fragment the list and squeeze all the
2181 : * occupied slots together to the high end, so that ifac_crack() has room
2182 : * for new factors. When this doesn't suffice, we abandon the current vector
2183 : * and allocate a somewhat larger one, defragmenting again while copying.
2184 : *
2185 : * For internal use: note that all exponents will fit into C longs, given
2186 : * PARI's lgefint field size. When we work with them, we sometimes read
2187 : * out the GEN pointer, and sometimes do an itos, whatever is more con-
2188 : * venient for the task at hand. */
2189 :
2190 : /*** Overview ***/
2191 :
2192 : /* The '*where' argument in the following points into *partial at the first of
2193 : * the three fields of the first occupied slot. It's there because the caller
2194 : * would already know where 'here' is, so we don't want to search for it again.
2195 : * We do not preserve this from one user-interface call to the next. */
2196 :
2197 : /* In the most cases, control flows from the user interface to ifac_main() and
2198 : * then to a succession of ifac_crack()s and ifac_divide()s, with (typically)
2199 : * none of the latter finding anything. */
2200 :
2201 : #define LAST(x) x+lg(x)-3
2202 : #define FIRST(x) x+3
2203 :
2204 : #define MOEBIUS(x) gel(x,1)
2205 : #define HINT(x) gel(x,2)
2206 :
2207 : /* y <- x */
2208 : INLINE void
2209 0 : SHALLOWCOPY(GEN x, GEN y) {
2210 0 : VALUE(y) = VALUE(x);
2211 0 : EXPON(y) = EXPON(x);
2212 0 : CLASS(y) = CLASS(x);
2213 0 : }
2214 : /* y <- x */
2215 : INLINE void
2216 0 : COPY(GEN x, GEN y) {
2217 0 : icopyifstack(VALUE(x), VALUE(y));
2218 0 : icopyifstack(EXPON(x), EXPON(y));
2219 0 : CLASS(y) = CLASS(x);
2220 0 : }
2221 :
2222 : /* Diagnostics */
2223 : static void
2224 0 : ifac_factor_dbg(GEN x)
2225 : {
2226 0 : GEN c = CLASS(x), v = VALUE(x);
2227 0 : if (c == gen_2) err_printf("IFAC: factor %Ps\n\tis prime (finished)\n", v);
2228 0 : else if (c == gen_1) err_printf("IFAC: factor %Ps\n\tis prime\n", v);
2229 0 : else if (c == gen_0) err_printf("IFAC: factor %Ps\n\tis composite\n", v);
2230 0 : }
2231 : static void
2232 0 : ifac_check(GEN partial, GEN where)
2233 : {
2234 0 : if (!where || where < FIRST(partial) || where > LAST(partial))
2235 0 : pari_err_BUG("ifac_check ['where' out of bounds]");
2236 0 : }
2237 : static void
2238 0 : ifac_print(GEN part, GEN where)
2239 : {
2240 0 : long l = lg(part);
2241 : GEN p;
2242 :
2243 0 : err_printf("ifac partial factorization structure: %ld slots, ", (l-3)/3);
2244 0 : if (MOEBIUS(part)) err_printf("Moebius mode, ");
2245 0 : err_printf("hint = %ld\n", itos(HINT(part)));
2246 0 : ifac_check(part, where);
2247 0 : for (p = part+3; p < part + l; p += 3)
2248 : {
2249 0 : GEN v = VALUE(p), e = EXPON(p), c = CLASS(p);
2250 0 : const char *s = "";
2251 0 : if (!v) { err_printf("[empty slot]\n"); continue; }
2252 0 : if (c == NULL) s = "unknown";
2253 0 : else if (c == gen_0) s = "composite";
2254 0 : else if (c == gen_1) s = "unfinished prime";
2255 0 : else if (c == gen_2) s = "prime";
2256 0 : else pari_err_BUG("unknown factor class");
2257 0 : err_printf("[%Ps, %Ps, %s]\n", v, e, s);
2258 : }
2259 0 : err_printf("Done.\n");
2260 0 : }
2261 :
2262 : static const long decomp_default_hint = 0;
2263 : /* assume n > 0, which we can assign to */
2264 : /* return initial data structure, see ifac_crack() for the hint argument */
2265 : static GEN
2266 5940 : ifac_start_hint(GEN n, int moebius, long hint)
2267 : {
2268 5940 : const long ifac_initial_length = 3 + 7*3;
2269 : /* codeword, moebius, hint, 7 slots -- a 512-bit product of distinct 8-bit
2270 : * primes needs at most 7 slots at a time) */
2271 5940 : GEN here, part = cgetg(ifac_initial_length, t_VEC);
2272 :
2273 5940 : MOEBIUS(part) = moebius? gen_1 : NULL;
2274 5940 : HINT(part) = stoi(hint);
2275 : /* fill first slot at the top end */
2276 5940 : here = part + ifac_initial_length - 3; /* LAST(part) */
2277 5940 : INIT(here, n,gen_1,gen_0); /* n^1: composite */
2278 41580 : while ((here -= 3) > part) ifac_delete(here);
2279 5940 : return part;
2280 : }
2281 : GEN
2282 2509 : ifac_start(GEN n, int moebius)
2283 2509 : { return ifac_start_hint(n,moebius,decomp_default_hint); }
2284 :
2285 : /* Return next nonempty slot after 'here', NULL if none exist */
2286 : static GEN
2287 15293 : ifac_find(GEN partial)
2288 : {
2289 15293 : GEN scan, end = partial + lg(partial);
2290 :
2291 : #ifdef IFAC_DEBUG
2292 : ifac_check(partial, partial);
2293 : #endif
2294 111505 : for (scan = partial+3; scan < end; scan += 3)
2295 106781 : if (VALUE(scan)) return scan;
2296 4724 : return NULL;
2297 : }
2298 :
2299 : /* Defragment: squeeze out unoccupied slots above *where. Unoccupied slots
2300 : * arise when a composite factor dissolves completely whilst dividing off a
2301 : * prime, or when ifac_resort() spots a coincidence and merges two factors.
2302 : * Update *where */
2303 : static void
2304 14 : ifac_defrag(GEN *partial, GEN *where)
2305 : {
2306 14 : GEN scan_new = LAST(*partial), scan_old;
2307 :
2308 42 : for (scan_old = scan_new; scan_old >= *where; scan_old -= 3)
2309 : {
2310 28 : if (!VALUE(scan_old)) continue; /* empty slot */
2311 28 : if (scan_old < scan_new) SHALLOWCOPY(scan_old, scan_new);
2312 28 : scan_new -= 3; /* point at next slot to be written */
2313 : }
2314 14 : scan_new += 3; /* back up to last slot written */
2315 14 : *where = scan_new;
2316 84 : while ((scan_new -= 3) > *partial) ifac_delete(scan_new); /* erase junk */
2317 14 : }
2318 :
2319 : /* Move to a larger main vector, updating *where if it points into it, and
2320 : * *partial in any case. Can be used as a specialized gcopy before
2321 : * a gerepileupto() (pass 0 as the new length). Normally, one would pass
2322 : * new_lg=1 to let this function guess the new size. To be used sparingly.
2323 : * Complex version of ifac_defrag(), combined with reallocation. If new_lg
2324 : * is 0, use the old length, so this acts just like gcopy except that the
2325 : * 'where' pointer is carried along; if it is 1, we make an educated guess.
2326 : * Exception: If new_lg is 0, the vector is full to the brim, and the first
2327 : * entry is composite, we make it longer to avoid being called again a
2328 : * microsecond later. It is safe to call this with *where = NULL:
2329 : * if it doesn't point anywhere within the old structure, it is left alone */
2330 : static void
2331 0 : ifac_realloc(GEN *partial, GEN *where, long new_lg)
2332 : {
2333 0 : long old_lg = lg(*partial);
2334 : GEN newpart, scan_new, scan_old;
2335 :
2336 0 : if (new_lg == 1)
2337 0 : new_lg = 2*old_lg - 6; /* from 7 slots to 13 to 25... */
2338 0 : else if (new_lg <= old_lg) /* includes case new_lg == 0 */
2339 : {
2340 0 : GEN first = *partial + 3;
2341 0 : new_lg = old_lg;
2342 : /* structure full and first entry composite or unknown */
2343 0 : if (VALUE(first) && (CLASS(first) == gen_0 || CLASS(first)==NULL))
2344 0 : new_lg += 6; /* give it a little more breathing space */
2345 : }
2346 0 : newpart = cgetg(new_lg, t_VEC);
2347 0 : if (DEBUGMEM >= 3)
2348 0 : err_printf("IFAC: new partial factorization structure (%ld slots)\n",
2349 0 : (new_lg - 3)/3);
2350 0 : MOEBIUS(newpart) = MOEBIUS(*partial);
2351 0 : icopyifstack(HINT(*partial), HINT(newpart));
2352 : /* Downward sweep through the old *partial. Pick up 'where' and carry it
2353 : * over if we pass it. (Only useful if it pointed at a nonempty slot.)
2354 : * Factors are COPY'd so that we again have a nice object (parent older
2355 : * than children, connected), except the one factor that may still be living
2356 : * in a clone where n originally was; exponents are similarly copied if they
2357 : * aren't global constants; class-of-factor fields are global constants so we
2358 : * need only copy them as pointers. Caller may then do a gerepileupto() */
2359 0 : scan_new = newpart + new_lg - 3; /* LAST(newpart) */
2360 0 : scan_old = *partial + old_lg - 3; /* LAST(*partial) */
2361 0 : for (; scan_old > *partial + 2; scan_old -= 3)
2362 : {
2363 0 : if (*where == scan_old) *where = scan_new;
2364 0 : if (!VALUE(scan_old)) continue; /* skip empty slots */
2365 0 : COPY(scan_old, scan_new); scan_new -= 3;
2366 : }
2367 0 : scan_new += 3; /* back up to last slot written */
2368 0 : while ((scan_new -= 3) > newpart) ifac_delete(scan_new);
2369 0 : *partial = newpart;
2370 0 : }
2371 :
2372 : /* Re-sort one (typically unknown) entry from washere to a new position,
2373 : * rotating intervening entries upward to fill the vacant space. If the new
2374 : * position is the same as the old one, or the new value of the entry coincides
2375 : * with a value already occupying a lower slot, then we just add exponents (and
2376 : * use the 'more known' class, and return 1 immediately when in Moebius mode).
2377 : * Slots between *where and washere must be in sorted order, so a sweep using
2378 : * this to re-sort several unknowns must proceed upward, see ifac_resort().
2379 : * Bubble-sort-of-thing sort. Won't be exercised frequently, so this is ok */
2380 : static void
2381 7 : ifac_sort_one(GEN *where, GEN washere)
2382 : {
2383 7 : GEN old, scan = washere - 3;
2384 : GEN value, exponent, class0, class1;
2385 : long cmp_res;
2386 :
2387 7 : if (scan < *where) return; /* nothing to do, washere==*where */
2388 7 : value = VALUE(washere);
2389 7 : exponent = EXPON(washere);
2390 7 : class0 = CLASS(washere);
2391 7 : cmp_res = -1; /* sentinel */
2392 7 : while (scan >= *where) /* at least once */
2393 : {
2394 7 : if (VALUE(scan))
2395 : { /* current slot nonempty, check against where */
2396 7 : cmp_res = cmpii(value, VALUE(scan));
2397 7 : if (cmp_res >= 0) break; /* have found where to stop */
2398 : }
2399 : /* copy current slot upward by one position and move pointers down */
2400 0 : SHALLOWCOPY(scan, scan+3);
2401 0 : scan -= 3;
2402 : }
2403 7 : scan += 3;
2404 : /* At this point there are the following possibilities:
2405 : * 1) cmp_res == -1. Either value is less than that at *where, or *where was
2406 : * pointing at vacant slots and any factors we saw en route were larger than
2407 : * value. At any rate, scan == *where now, and scan is pointing at an empty
2408 : * slot, into which we'll stash our entry.
2409 : * 2) cmp_res == 0. The entry at scan-3 is the one, we compare class0
2410 : * fields and add exponents, and put it all into the vacated scan slot,
2411 : * NULLing the one at scan-3 (and possibly updating *where).
2412 : * 3) cmp_res == 1. The slot at scan is the one to store our entry into. */
2413 7 : if (cmp_res)
2414 : {
2415 7 : if (cmp_res < 0 && scan != *where)
2416 0 : pari_err_BUG("ifact_sort_one [misaligned partial]");
2417 7 : INIT(scan, value, exponent, class0); return;
2418 : }
2419 : /* case cmp_res == 0: repeated factor detected */
2420 0 : if (DEBUGLEVEL >= 4)
2421 0 : err_printf("IFAC: repeated factor %Ps\n\tin ifac_sort_one\n", value);
2422 0 : old = scan - 3;
2423 : /* if old class0 was composite and new is prime, or vice versa, complain
2424 : * (and if one class0 was unknown and the other wasn't, use the known one) */
2425 0 : class1 = CLASS(old);
2426 0 : if (class0) /* should never be used */
2427 : {
2428 0 : if (class1)
2429 : {
2430 0 : if (class0 == gen_0 && class1 != gen_0)
2431 0 : pari_err_BUG("ifac_sort_one (composite = prime)");
2432 0 : else if (class0 != gen_0 && class1 == gen_0)
2433 0 : pari_err_BUG("ifac_sort_one (prime = composite)");
2434 0 : else if (class0 == gen_2)
2435 0 : CLASS(scan) = class0;
2436 : }
2437 : else
2438 0 : CLASS(scan) = class0;
2439 : }
2440 : /* else stay with the existing known class0 */
2441 0 : CLASS(scan) = class1;
2442 : /* in any case, add exponents */
2443 0 : if (EXPON(old) == gen_1 && exponent == gen_1)
2444 0 : EXPON(scan) = gen_2;
2445 : else
2446 0 : EXPON(scan) = addii(EXPON(old), exponent);
2447 : /* move the value over and null out the vacated slot below */
2448 0 : old = scan - 3;
2449 0 : *scan = *old;
2450 0 : ifac_delete(old);
2451 : /* finally, see whether *where should be pulled in */
2452 0 : if (old == *where) *where += 3;
2453 : }
2454 :
2455 : /* Sort all current unknowns downward to where they belong. Sweeps in the
2456 : * upward direction. Not needed after ifac_crack(), only when ifac_divide()
2457 : * returned true. Update *where. */
2458 : static void
2459 7 : ifac_resort(GEN *partial, GEN *where)
2460 : {
2461 : GEN scan, end;
2462 7 : ifac_defrag(partial, where); end = LAST(*partial);
2463 21 : for (scan = *where; scan <= end; scan += 3)
2464 14 : if (VALUE(scan) && !CLASS(scan)) ifac_sort_one(where, scan); /*unknown*/
2465 7 : ifac_defrag(partial, where); /* remove newly created gaps */
2466 7 : }
2467 :
2468 : /* Let x be a t_INT known not to have small divisors (< 661, and < 2^14 for huge
2469 : * x > 2^512). Return 0 if x is a proven composite. Return 1 if we believe it
2470 : * to be prime (fully proven prime if factor_proven is set). */
2471 : int
2472 28148 : ifac_isprime(GEN x)
2473 : {
2474 28148 : if (!BPSW_psp_nosmalldiv(x)) return 0; /* composite */
2475 19818 : if (factor_proven && ! BPSW_isprime(x))
2476 : {
2477 0 : pari_warn(warner,
2478 : "IFAC: pseudo-prime %Ps\n\tis not prime. PLEASE REPORT!\n", x);
2479 0 : return 0;
2480 : }
2481 19818 : return 1;
2482 : }
2483 :
2484 : static int
2485 11291 : ifac_checkprime(GEN x)
2486 : {
2487 11291 : int res = ifac_isprime(VALUE(x));
2488 11291 : CLASS(x) = res? gen_1: gen_0;
2489 11291 : if (DEBUGLEVEL>2) ifac_factor_dbg(x);
2490 11291 : return res;
2491 : }
2492 :
2493 : /* Determine primality or compositeness of all current unknowns, and set
2494 : * class Q primes to finished (class P) if everything larger is already
2495 : * known to be prime. When after_crack >= 0, only look at the
2496 : * first after_crack things in the list (do nothing when it's 0) */
2497 : static void
2498 5794 : ifac_whoiswho(GEN *partial, GEN *where, long after_crack)
2499 : {
2500 5794 : GEN scan, scan_end = LAST(*partial);
2501 :
2502 : #ifdef IFAC_DEBUG
2503 : ifac_check(*partial, *where);
2504 : #endif
2505 5794 : if (after_crack == 0) return;
2506 5180 : if (after_crack > 0) /* check at most after_crack entries */
2507 5173 : scan = *where + 3*(after_crack - 1); /* assert(scan <= scan_end) */
2508 : else
2509 7 : for (scan = scan_end; scan >= *where; scan -= 3)
2510 : {
2511 7 : if (CLASS(scan))
2512 : { /* known class of factor */
2513 0 : if (CLASS(scan) == gen_0) break;
2514 0 : if (CLASS(scan) == gen_1)
2515 : {
2516 0 : if (DEBUGLEVEL>=3)
2517 : {
2518 0 : err_printf("IFAC: factor %Ps\n\tis prime (no larger composite)\n",
2519 0 : VALUE(*where));
2520 0 : err_printf("IFAC: prime %Ps\n\tappears with exponent = %ld\n",
2521 0 : VALUE(*where), itos(EXPON(*where)));
2522 : }
2523 0 : CLASS(scan) = gen_2;
2524 : }
2525 0 : continue;
2526 : }
2527 7 : if (!ifac_checkprime(scan)) break; /* must disable Q-to-P */
2528 0 : CLASS(scan) = gen_2; /* P_i, finished prime */
2529 0 : if (DEBUGLEVEL>2) ifac_factor_dbg(scan);
2530 : }
2531 : /* go on, Q-to-P trick now disabled */
2532 15809 : for (; scan >= *where; scan -= 3)
2533 : {
2534 10629 : if (CLASS(scan)) continue;
2535 10594 : (void)ifac_checkprime(scan); /* Qj | Ck */
2536 : }
2537 : }
2538 :
2539 : /* Divide all current composites by first (prime, class Q) entry, updating its
2540 : * exponent, and turning it into a finished prime (class P). Return 1 if any
2541 : * such divisions succeeded (in Moebius mode, the update may then not have
2542 : * been completed), or 0 if none of them succeeded. Doesn't modify *where.
2543 : * Here we normally do not check that the first entry is a not-finished
2544 : * prime. Stack management: we may allocate a new exponent */
2545 : static long
2546 9847 : ifac_divide(GEN *partial, GEN *where, long moebius_mode)
2547 : {
2548 9847 : GEN scan, scan_end = LAST(*partial);
2549 9847 : long res = 0, exponent, newexp, otherexp;
2550 :
2551 : #ifdef IFAC_DEBUG
2552 : ifac_check(*partial, *where);
2553 : if (CLASS(*where) != gen_1)
2554 : pari_err_BUG("ifac_divide [division by composite or finished prime]");
2555 : if (!VALUE(*where)) pari_err_BUG("ifac_divide [division by nothing]");
2556 : #endif
2557 9847 : newexp = exponent = itos(EXPON(*where));
2558 9847 : if (exponent > 1 && moebius_mode) return 1;
2559 : /* should've been caught by caller */
2560 :
2561 15559 : for (scan = *where+3; scan <= scan_end; scan += 3)
2562 : {
2563 5712 : if (CLASS(scan) != gen_0) continue; /* the other thing ain't composite */
2564 205 : otherexp = 0;
2565 : /* divide in place to keep stack clutter minimal */
2566 219 : while (dvdiiz(VALUE(scan), VALUE(*where), VALUE(scan)))
2567 : {
2568 14 : if (moebius_mode) return 1; /* immediately */
2569 14 : if (!otherexp) otherexp = itos(EXPON(scan));
2570 14 : newexp += otherexp;
2571 : }
2572 205 : if (newexp > exponent) /* did anything happen? */
2573 : {
2574 7 : EXPON(*where) = (newexp == 2 ? gen_2 : utoipos(newexp));
2575 7 : exponent = newexp;
2576 7 : if (is_pm1((GEN)*scan)) /* factor dissolved completely */
2577 : {
2578 0 : ifac_delete(scan);
2579 0 : if (DEBUGLEVEL >= 4)
2580 0 : err_printf("IFAC: a factor was a power of another prime factor\n");
2581 : } else {
2582 7 : CLASS(scan) = NULL; /* at any rate it's Unknown now */
2583 7 : if (DEBUGLEVEL >= 4)
2584 0 : err_printf("IFAC: a factor was divisible by another prime factor,\n"
2585 : "\tleaving a cofactor = %Ps\n", VALUE(scan));
2586 : }
2587 7 : res = 1;
2588 7 : if (DEBUGLEVEL >= 5)
2589 0 : err_printf("IFAC: prime %Ps\n\tappears at least to the power %ld\n",
2590 0 : VALUE(*where), newexp);
2591 : }
2592 : } /* for */
2593 9847 : CLASS(*where) = gen_2; /* make it a finished prime */
2594 9847 : if (DEBUGLEVEL >= 3)
2595 0 : err_printf("IFAC: prime %Ps\n\tappears with exponent = %ld\n",
2596 0 : VALUE(*where), newexp);
2597 9847 : return res;
2598 : }
2599 :
2600 : /* found out our integer was factor^exp. Update */
2601 : static void
2602 880 : update_pow(GEN where, GEN factor, long exp, pari_sp *av)
2603 : {
2604 880 : GEN ex = EXPON(where);
2605 880 : if (DEBUGLEVEL>3)
2606 0 : err_printf("IFAC: found %Ps =\n\t%Ps ^%ld\n", *where, factor, exp);
2607 880 : affii(factor, VALUE(where)); set_avma(*av);
2608 880 : if (ex == gen_1)
2609 677 : { EXPON(where) = exp == 2? gen_2: utoipos(exp); *av = avma; }
2610 203 : else if (ex == gen_2)
2611 182 : { EXPON(where) = utoipos(exp<<1); *av = avma; }
2612 : else
2613 21 : affsi(exp * itos(ex), EXPON(where));
2614 880 : }
2615 : /* hint = 0 : Use a default strategy
2616 : * hint & 1 : avoid MPQS
2617 : * hint & 2 : avoid first-stage ECM (may fall back to ECM if MPQS gives up)
2618 : * hint & 4 : avoid Pollard and SQUFOF stages.
2619 : * hint & 8 : avoid final ECM; may flag a composite as prime. */
2620 : #define get_hint(partial) (itos(HINT(*partial)) & 15)
2621 :
2622 : /* Complete ifac_crack's job when a factoring engine splits the current factor
2623 : * into a product of three or more new factors. Makes room for them if
2624 : * necessary, sorts them, gives them the right exponents and class. Returns the
2625 : * number of factors actually written, which may be less than #facvec if there
2626 : * are duplicates. Vectors of factors (cf pollardbrent() or mpqs()) contain
2627 : * 'slots' of 3 GENs per factor, interpreted as in our partial factorization
2628 : * data structure. Thus engines can tell us what they already know about
2629 : * factors being prime/composite or appearing to a power larger than thefirst.
2630 : * Don't collect garbage. No diagnostics: engine has printed what it found.
2631 : * facvec contains slots of three components per factor; repeated factors are
2632 : * allowed (their classes shouldn't contradict each other whereas their
2633 : * exponents will be added up) */
2634 : static long
2635 3284 : ifac_insert_multiplet(GEN *partial, GEN *where, GEN facvec, long moebius_mode)
2636 : {
2637 3284 : long j,k=1, lfv=lg(facvec)-1, nf=lfv/3, room=(long)(*where-*partial);
2638 : /* one of the factors will go into the *where slot, so room is now 3 times
2639 : * the number of slots we can use */
2640 3284 : long needroom = lfv - room;
2641 3284 : GEN newexp, cur, sorted, auxvec = cgetg(nf+1, t_VEC), factor;
2642 3284 : long E = itos(EXPON(*where)); /* the old exponent */
2643 :
2644 3284 : if (DEBUGLEVEL >= 5) /* squfof may return a single squared factor as a set */
2645 0 : err_printf("IFAC: incorporating set of %ld factor(s)\n", nf);
2646 3284 : if (needroom > 0) ifac_realloc(partial, where, lg(*partial) + needroom);
2647 :
2648 : /* create sort permutation from the values of the factors */
2649 10128 : for (j=nf; j; j--) gel(auxvec,j) = gel(facvec,3*j-2);
2650 3284 : sorted = ZV_indexsort(auxvec);
2651 : /* store factors, beginning at *where, and catching any duplicates */
2652 3284 : cur = facvec + 3*sorted[nf]-2; /* adjust for triple spacing */
2653 3284 : VALUE(*where) = VALUE(cur);
2654 3284 : newexp = EXPON(cur);
2655 : /* if new exponent is 1, the old exponent already in place will do */
2656 3284 : if (newexp != gen_1) EXPON(*where) = mului(E,newexp);
2657 3284 : CLASS(*where) = CLASS(cur);
2658 3284 : if (DEBUGLEVEL >= 6) err_printf("\tstored (largest) factor no. %ld...\n", nf);
2659 :
2660 6844 : for (j=nf-1; j; j--)
2661 : {
2662 : GEN e;
2663 3560 : cur = facvec + 3*sorted[j]-2;
2664 3560 : factor = VALUE(cur);
2665 3560 : if (equalii(factor, VALUE(*where)))
2666 : {
2667 0 : if (DEBUGLEVEL >= 6)
2668 0 : err_printf("\tfactor no. %ld is a duplicate%s\n", j, (j>1? "...": ""));
2669 : /* update exponent, ignore class which would already have been set,
2670 : * then forget current factor */
2671 0 : newexp = EXPON(cur);
2672 0 : if (newexp != gen_1) /* new exp > 1 */
2673 0 : e = addiu(EXPON(*where), E * itou(newexp));
2674 0 : else if (EXPON(*where) == gen_1 && E == 1)
2675 0 : e = gen_2;
2676 : else
2677 0 : e = addiu(EXPON(*where), E);
2678 0 : EXPON(*where) = e;
2679 :
2680 0 : if (moebius_mode) return 0; /* stop now, with exponent updated */
2681 0 : continue;
2682 : }
2683 :
2684 3560 : *where -= 3;
2685 3560 : CLASS(*where) = CLASS(cur); /* class as given */
2686 3560 : newexp = EXPON(cur);
2687 3560 : if (newexp != gen_1) /* new exp > 1 */
2688 99 : e = (E == 1 && newexp == gen_2)? gen_2: mului(E, newexp);
2689 : else /* inherit parent's exponent */
2690 3461 : e = (E == 1 ? gen_1 : (E == 2 ? gen_2 : utoipos(E)));
2691 3560 : EXPON(*where) = e;
2692 : /* keep components younger than *partial */
2693 3560 : icopyifstack(factor, VALUE(*where));
2694 3560 : k++;
2695 3560 : if (DEBUGLEVEL >= 6)
2696 0 : err_printf("\tfactor no. %ld was unique%s\n", j, j>1? " (so far)...": "");
2697 : }
2698 3284 : return k;
2699 : }
2700 :
2701 : /* x /= y; exact division */
2702 : static void
2703 1883 : diviiexact_inplace(GEN x, GEN y)
2704 1883 : { pari_sp av = avma; affii(diviiexact(x, y), x); set_avma(av); }
2705 :
2706 : /* Split the first (composite) entry. There must already be room for another
2707 : * factor below *where, and *where is updated. Two cases:
2708 : * - entry is a pure power: factor^k is inserted, leaving *where unchanged;
2709 : * - entry = factor * cofactor (not necessarily coprime): both factors are
2710 : * inserted in the correct order, updating *where
2711 : * The inserted factors class is set to unknown, they inherit the exponent
2712 : * (or a multiple thereof) of their ancestor.
2713 : *
2714 : * Returns number of factors written into the structure, usually 2; only 1
2715 : * if pure power, and > 2 if a factoring engine returned a vector of factors.
2716 : * Can reallocate the data structure in the rare > 2 case; may create one or
2717 : * more objects: new factors or exponents > 2 */
2718 : static long
2719 5789 : ifac_crack(GEN *partial, GEN *where, long moebius_mode)
2720 : {
2721 5789 : long hint = get_hint(partial);
2722 : GEN cofactor, factor, exponent;
2723 :
2724 : #ifdef IFAC_DEBUG
2725 : ifac_check(*partial, *where);
2726 : if (*where < *partial + 6)
2727 : pari_err_BUG("ifac_crack ['*where' out of bounds]");
2728 : if (!(VALUE(*where)) || typ(VALUE(*where)) != t_INT)
2729 : pari_err_BUG("ifac_crack [incorrect VALUE(*where)]");
2730 : if (CLASS(*where) != gen_0)
2731 : pari_err_BUG("ifac_crack [operand not known composite]");
2732 : #endif
2733 :
2734 5789 : if (DEBUGLEVEL>2) {
2735 0 : err_printf("IFAC: cracking composite\n\t%Ps\n", **where);
2736 0 : if (DEBUGLEVEL>3) err_printf("IFAC: checking for pure square\n");
2737 : }
2738 : /* MPQS cannot factor prime powers. Look for pure powers even if MPQS is
2739 : * blocked by hint: fast and useful in bounded factorization */
2740 : {
2741 : forprime_t T;
2742 5789 : ulong exp = 1, mask = 7;
2743 5789 : long good = 0;
2744 5789 : pari_sp av = avma;
2745 5789 : (void)u_forprime_init(&T, 11, ULONG_MAX);
2746 : /* crack squares */
2747 6626 : while (Z_issquareall(VALUE(*where), &factor))
2748 : {
2749 838 : good = 1; /* remember we succeeded once */
2750 838 : update_pow(*where, factor, 2, &av);
2751 1452 : if (moebius_mode) return 0; /* no need to carry on */
2752 : }
2753 5830 : while ( (exp = is_357_power(VALUE(*where), &factor, &mask)) )
2754 : {
2755 42 : good = 1; /* remember we succeeded once */
2756 42 : update_pow(*where, factor, exp, &av);
2757 42 : if (moebius_mode) return 0; /* no need to carry on */
2758 : }
2759 : /* cutoff at 14 bits: OK if tridiv_bound >= 2^14 OR if >= 661 for
2760 : * an integer < 701^11 (103 bits). */
2761 5788 : while ( (exp = is_pth_power(VALUE(*where), &factor, &T, 15)) )
2762 : {
2763 0 : good = 1; /* remember we succeeded once */
2764 0 : update_pow(*where, factor, exp, &av);
2765 0 : if (moebius_mode) return 0; /* no need to carry on */
2766 : }
2767 :
2768 5788 : if (good && hint != 15 && ifac_checkprime(*where))
2769 : { /* our composite was a prime power */
2770 614 : if (DEBUGLEVEL>3)
2771 0 : err_printf("IFAC: factor %Ps\n\tis prime\n", VALUE(*where));
2772 614 : return 0; /* bypass subsequent ifac_whoiswho() call */
2773 : }
2774 : } /* pure power stage */
2775 :
2776 5174 : factor = NULL;
2777 5174 : if (!(hint & 4))
2778 : { /* SQUFOF then Rho */
2779 5125 : if (DEBUGLEVEL >= 4)
2780 0 : err_printf("IFAC: trying Shanks' SQUFOF, will fail silently if input\n"
2781 : " is too large for it.\n");
2782 5125 : factor = squfof(VALUE(*where));
2783 5125 : if (!factor)
2784 : {
2785 3494 : if (DEBUGLEVEL >= 4) err_printf("IFAC: trying Pollard-Brent rho\n");
2786 3494 : factor = pollardbrent(VALUE(*where));
2787 : }
2788 : }
2789 5174 : if (!factor && !(hint & 2))
2790 : { /* First ECM stage */
2791 3267 : if (DEBUGLEVEL >= 4) err_printf("IFAC: trying Lenstra-Montgomery ECM\n");
2792 3267 : factor = ellfacteur(VALUE(*where), 0); /* do not insist */
2793 : }
2794 5174 : if (!factor && !(hint & 1))
2795 : { /* MPQS stage */
2796 3291 : if (DEBUGLEVEL >= 4) err_printf("IFAC: trying MPQS\n");
2797 3291 : factor = mpqs(VALUE(*where));
2798 : }
2799 5174 : if (!factor && !(hint & 8))
2800 : { /* Final ECM stage, guaranteed to succeed */
2801 0 : if (DEBUGLEVEL >= 4)
2802 0 : err_printf("IFAC: forcing ECM, may take some time\n");
2803 0 : factor = ellfacteur(VALUE(*where), 1);
2804 : }
2805 5174 : if (!factor)
2806 : { /* limited factorization */
2807 7 : if (DEBUGLEVEL >= 2)
2808 : {
2809 0 : pari_warn(warner, hint==15? "IFAC: untested integer declared prime"
2810 : : "IFAC: unfactored composite declared prime");
2811 : /* don't print it out at level 3 or above, where it would appear
2812 : * several times before and after this message already */
2813 0 : if (DEBUGLEVEL == 2) err_printf("\t%Ps\n", VALUE(*where));
2814 : }
2815 7 : CLASS(*where) = gen_1; /* might as well trial-divide by it... */
2816 7 : return 1;
2817 : }
2818 : /* At least two factors */
2819 5167 : if (typ(factor) == t_VEC)
2820 3284 : return ifac_insert_multiplet(partial, where, factor, moebius_mode);
2821 :
2822 : /* Single factor (t_INT): work out cofactor in place */
2823 1883 : diviiexact_inplace(VALUE(*where), factor);
2824 1883 : cofactor = VALUE(*where);
2825 : /* factoring engines reported factor; tell about the cofactor */
2826 1883 : if (DEBUGLEVEL >= 4) err_printf("IFAC: cofactor = %Ps\n", cofactor);
2827 1883 : CLASS(*where) = NULL; /* mark factor /cofactor 'unknown' */
2828 1883 : exponent = EXPON(*where); /* common exponent */
2829 :
2830 1883 : *where -= 3;
2831 1883 : CLASS(*where) = NULL; /* mark factor /cofactor 'unknown' */
2832 1883 : icopyifstack(exponent, EXPON(*where)); /* copy common exponent */
2833 :
2834 1883 : if (cmpii(factor, cofactor) < 0)
2835 1705 : VALUE(*where) = factor; /* common case */
2836 : else
2837 : { /* factor > cofactor, swap. Exponent is the same, so no need to swap. */
2838 178 : GEN old = *where + 3;
2839 178 : VALUE(*where) = VALUE(old); /* move cofactor pointer to lowest slot */
2840 178 : VALUE(old) = factor; /* save factor */
2841 : }
2842 1883 : return 2;
2843 : }
2844 :
2845 : /* main loop: iterate until smallest entry is a finished prime; returns
2846 : * a 'where' pointer, or NULL if nothing left, or gen_0 in Moebius mode if
2847 : * we aren't squarefree */
2848 : static GEN
2849 14213 : ifac_main(GEN *partial)
2850 : {
2851 14213 : const long moebius_mode = !!MOEBIUS(*partial);
2852 14213 : GEN here = ifac_find(*partial);
2853 : long nf;
2854 :
2855 14213 : if (!here) return NULL; /* nothing left */
2856 : /* loop until first entry is a finished prime. May involve reallocations,
2857 : * thus updates of *partial */
2858 25483 : while (CLASS(here) != gen_2)
2859 : {
2860 15636 : if (CLASS(here) == gen_0) /* composite: crack it */
2861 : { /* make sure there's room for another factor */
2862 5789 : if (here < *partial + 6)
2863 : {
2864 0 : ifac_defrag(partial, &here);
2865 0 : if (here < *partial + 6) ifac_realloc(partial, &here, 1); /* no luck */
2866 : }
2867 5789 : nf = ifac_crack(partial, &here, moebius_mode);
2868 5789 : if (moebius_mode && EXPON(here) != gen_1) /* that was a power */
2869 : {
2870 2 : if (DEBUGLEVEL >= 3)
2871 0 : err_printf("IFAC: main loop: repeated new factor\n\t%Ps\n", *here);
2872 2 : return gen_0;
2873 : }
2874 : /* deal with the new unknowns. No sort: ifac_crack did it */
2875 5787 : ifac_whoiswho(partial, &here, nf);
2876 5787 : continue;
2877 : }
2878 9847 : if (CLASS(here) == gen_1) /* prime but not yet finished: finish it */
2879 : {
2880 9847 : if (ifac_divide(partial, &here, moebius_mode))
2881 : {
2882 7 : if (moebius_mode)
2883 : {
2884 0 : if (DEBUGLEVEL >= 3)
2885 0 : err_printf("IFAC: main loop: another factor was divisible by\n"
2886 : "\t%Ps\n", *here);
2887 0 : return gen_0;
2888 : }
2889 7 : ifac_resort(partial, &here); /* sort new cofactors down */
2890 7 : ifac_whoiswho(partial, &here, -1);
2891 : }
2892 9847 : continue;
2893 : }
2894 0 : pari_err_BUG("ifac_main [nonexistent factor class]");
2895 : } /* while */
2896 9847 : if (moebius_mode && EXPON(here) != gen_1)
2897 : {
2898 0 : if (DEBUGLEVEL >= 3)
2899 0 : err_printf("IFAC: after main loop: repeated old factor\n\t%Ps\n", *here);
2900 0 : return gen_0;
2901 : }
2902 9847 : if (DEBUGLEVEL >= 4)
2903 : {
2904 0 : nf = (*partial + lg(*partial) - here - 3)/3;
2905 0 : if (nf)
2906 0 : err_printf("IFAC: main loop: %ld factor%s left\n", nf, (nf>1)? "s": "");
2907 : else
2908 0 : err_printf("IFAC: main loop: this was the last factor\n");
2909 : }
2910 9847 : if (factor_add_primes && !(get_hint(partial) & 8))
2911 : {
2912 0 : GEN p = VALUE(here);
2913 0 : if (lgefint(p)>3 || uel(p,2) > 0x1000000UL) (void)addprimes(p);
2914 : }
2915 9847 : return here;
2916 : }
2917 :
2918 : /* Encapsulated routines */
2919 :
2920 : /* prime/exponent pairs need to appear contiguously on the stack, but we also
2921 : * need our data structure somewhere, and we don't know in advance how many
2922 : * primes will turn up. The following discipline achieves this: When
2923 : * ifac_decomp() is called, n should point at an object older than the oldest
2924 : * small prime/exponent pair (ifactor() guarantees this).
2925 : * We allocate sufficient space to accommodate several pairs -- eleven pairs
2926 : * ought to fit in a space not much larger than n itself -- before calling
2927 : * ifac_start(). If we manage to complete the factorization before we run out
2928 : * of space, we free the data structure and cull the excess reserved space
2929 : * before returning. When we do run out, we have to leapfrog to generate more
2930 : * (guesstimating the requirements from what is left in the partial
2931 : * factorization structure); room for fresh pairs is allocated at the head of
2932 : * the stack, followed by an ifac_realloc() to reconnect the data structure and
2933 : * move it out of the way, followed by a few pointer tweaks to connect the new
2934 : * pairs space to the old one. This whole affair translates into a surprisingly
2935 : * compact routine. */
2936 :
2937 : /* find primary factors of n; destroy n */
2938 : static long
2939 2568 : ifac_decomp(GEN n, long hint)
2940 : {
2941 2568 : pari_sp av = avma;
2942 2568 : long nb = 0;
2943 2568 : GEN part, here, workspc, pairs = (GEN)av;
2944 :
2945 : /* workspc will be doled out in pairs of smaller t_INTs. For n = prod p^{e_p}
2946 : * (p not necessarily prime), need room to store all p and e_p [ cgeti(3) ],
2947 : * bounded by
2948 : * sum_{p | n} ( log_{2^BIL} (p) + 6 ) <= log_{2^BIL} n + 6 log_2 n */
2949 2568 : workspc = new_chunk((expi(n) + 1) * 7);
2950 2568 : part = ifac_start_hint(n, 0, hint);
2951 : for (;;)
2952 : {
2953 7562 : here = ifac_main(&part);
2954 7562 : if (!here) break;
2955 4994 : if (gc_needed(av,1))
2956 : {
2957 : long offset;
2958 0 : if(DEBUGMEM>1)
2959 : {
2960 0 : pari_warn(warnmem,"[2] ifac_decomp");
2961 0 : ifac_print(part, here);
2962 : }
2963 0 : ifac_realloc(&part, &here, 0);
2964 0 : offset = here - part;
2965 0 : part = gerepileupto((pari_sp)workspc, part);
2966 0 : here = part + offset;
2967 : }
2968 4994 : nb++;
2969 4994 : pairs = icopy_avma(VALUE(here), (pari_sp)pairs);
2970 4994 : pairs = icopy_avma(EXPON(here), (pari_sp)pairs);
2971 4994 : ifac_delete(here);
2972 : }
2973 2568 : set_avma((pari_sp)pairs);
2974 2568 : if (DEBUGLEVEL >= 3)
2975 0 : err_printf("IFAC: found %ld large prime (power) factor%s.\n",
2976 : nb, (nb>1? "s": ""));
2977 2568 : return nb;
2978 : }
2979 :
2980 : /***********************************************************************/
2981 : /** ARITHMETIC FUNCTIONS WITH EARLY-ABORT **/
2982 : /** needing direct access to the factoring machinery to avoid work: **/
2983 : /** e.g. if we find a square factor, moebius returns 0, core doesn't **/
2984 : /** need to factor it, etc. **/
2985 : /***********************************************************************/
2986 : /* memory management */
2987 : static void
2988 0 : ifac_GC(pari_sp av, GEN *part)
2989 : {
2990 0 : GEN here = NULL;
2991 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ifac_xxx");
2992 0 : ifac_realloc(part, &here, 0);
2993 0 : *part = gerepileupto(av, *part);
2994 0 : }
2995 :
2996 : /* destroys n */
2997 : static long
2998 236 : ifac_moebius(GEN n)
2999 : {
3000 236 : long mu = 1;
3001 236 : pari_sp av = avma;
3002 236 : GEN part = ifac_start(n, 1);
3003 : for(;;)
3004 468 : {
3005 : long v;
3006 : GEN p;
3007 704 : if (!ifac_next(&part,&p,&v)) return v? 0: mu;
3008 468 : mu = -mu;
3009 468 : if (gc_needed(av,1)) ifac_GC(av,&part);
3010 : }
3011 : }
3012 :
3013 : int
3014 760 : ifac_read(GEN part, GEN *p, long *e)
3015 : {
3016 760 : GEN here = ifac_find(part);
3017 760 : if (!here) return 0;
3018 400 : *p = VALUE(here);
3019 400 : *e = EXPON(here)[2];
3020 400 : return 1;
3021 : }
3022 : void
3023 320 : ifac_skip(GEN part)
3024 : {
3025 320 : GEN here = ifac_find(part);
3026 320 : if (here) ifac_delete(here);
3027 320 : }
3028 :
3029 : /* destroys n */
3030 : static int
3031 7 : ifac_ispowerful(GEN n)
3032 : {
3033 7 : pari_sp av = avma;
3034 7 : GEN part = ifac_start(n, 0);
3035 : for(;;)
3036 7 : {
3037 : long e;
3038 : GEN p;
3039 14 : if (!ifac_read(part,&p,&e)) return 1;
3040 : /* power: skip */
3041 7 : if (e != 1 || Z_isanypower(p,NULL)) { ifac_skip(part); continue; }
3042 0 : if (!ifac_next(&part,&p,&e)) return 1;
3043 0 : if (e == 1) return 0;
3044 0 : if (gc_needed(av,1)) ifac_GC(av,&part);
3045 : }
3046 : }
3047 : /* destroys n; assume n != 0 is composite */
3048 : static GEN
3049 353 : ifac_core(GEN n)
3050 : {
3051 353 : GEN m = gen_1, c = cgeti(lgefint(n));
3052 353 : pari_sp av = avma;
3053 353 : GEN part = ifac_start(n, 0);
3054 : for(;;)
3055 393 : {
3056 : long e;
3057 : GEN p;
3058 746 : if (!ifac_read(part,&p,&e)) return m;
3059 : /* square: skip */
3060 393 : if (!odd(e) || Z_issquare(p)) { ifac_skip(part); continue; }
3061 80 : if (!ifac_next(&part,&p,&e)) return m;
3062 80 : if (odd(e)) m = mulii(m, p);
3063 80 : if (gc_needed(av,1)) { affii(m,c); m=c; ifac_GC(av,&part); }
3064 : }
3065 : }
3066 :
3067 : /* must be >= 661 (various functions assume it in order to call uisprime_661
3068 : * instead of uisprime, and Z_isanypower_nosmalldiv instead of Z_isanypower) */
3069 : ulong
3070 4369152 : tridiv_boundu(ulong n)
3071 : {
3072 4369152 : long e = expu(n);
3073 4369017 : if(e<30) return 1UL<<12;
3074 : #ifdef LONG_IS_64BIT
3075 182925 : if(e<34) return 1UL<<13;
3076 100188 : if(e<37) return 1UL<<14;
3077 60054 : if(e<42) return 1UL<<15;
3078 26547 : if(e<47) return 1UL<<16;
3079 14882 : if(e<56) return 1UL<<17;
3080 4589 : if(e<56) return 1UL<<18;
3081 4589 : if(e<62) return 1UL<<19;
3082 1220 : return 1UL<<18;
3083 : #else
3084 6945 : return 1UL<<13;
3085 : #endif
3086 : }
3087 :
3088 : /* Where to stop trial dividing in factorization. Must be >= 661.
3089 : * If further n > 2^512, must be >= 2^14 */
3090 : ulong
3091 869536 : tridiv_bound(GEN n)
3092 : {
3093 869536 : if (lgefint(n)==3) return tridiv_boundu(n[2]);
3094 : else
3095 : {
3096 76318 : ulong l = (ulong)expi(n) + 1;
3097 76318 : if (l <= 512) return (l-16) << 10;
3098 1196 : return 1UL<<19; /* Rho is generally faster above this */
3099 : }
3100 : }
3101 :
3102 : /* destroys n */
3103 : static void
3104 863 : ifac_factoru(GEN n, long hint, GEN P, GEN E, long *pi)
3105 : {
3106 863 : GEN part = ifac_start_hint(n, 0, hint);
3107 863 : long i = *pi;
3108 : for(;;)
3109 1634 : {
3110 : long v;
3111 : GEN p;
3112 2497 : if (!ifac_next(&part,&p,&v)) { *pi = i; return; }
3113 1634 : P[i] = itou(p); E[i] = v; i++;
3114 : }
3115 : }
3116 : /* destroys n */
3117 : static long
3118 663 : ifac_moebiusu(GEN n)
3119 : {
3120 663 : GEN part = ifac_start(n, 1);
3121 663 : long s = 1;
3122 : for(;;)
3123 1326 : {
3124 : long v;
3125 : GEN p;
3126 1989 : if (!ifac_next(&part,&p,&v)) return v? 0: s;
3127 1326 : s = -s;
3128 : }
3129 : }
3130 :
3131 : INLINE ulong
3132 143807561 : u_forprime_next_fast(forprime_t *T)
3133 : {
3134 143807561 : if (++T->n <= pari_PRIMES[0])
3135 : {
3136 143808292 : T->p = pari_PRIMES[T->n];
3137 143808292 : return T->p > T->b ? 0: T->p;
3138 : }
3139 0 : return u_forprime_next(T);
3140 : }
3141 :
3142 : /* uisprime(n) knowing n has no prime divisor <= lim */
3143 : static int
3144 7866 : uisprime_nosmall(ulong n, ulong lim)
3145 7866 : { return (lim >= 661)? uisprime_661(n): uisprime(n); }
3146 :
3147 : static GEN factoru_sign(ulong n, ulong all, long hint, ulong *pU1, ulong *pU2);
3148 : static GEN ifactor_sign(GEN n, ulong all, long hint, long sn, GEN *pU);
3149 :
3150 : /* simplified version of factoru_sign, to be called on squarefree n whose
3151 : * prime divisors are in [minp, maxp]. In practice called with
3152 : * maxp <= maxprimelim() */
3153 : static GEN
3154 917212 : factoru_primes(ulong n, ulong minp, ulong maxp)
3155 : {
3156 : forprime_t S;
3157 : ulong p;
3158 : long i;
3159 : GEN P;
3160 :
3161 917212 : if (n < minp) return NULL;
3162 898337 : if (n <= maxp && PRIMES_search(n) > 0) return mkvecsmall(n);
3163 761722 : P = cgetg(16, t_VECSMALL); i = 1;
3164 761722 : u_forprime_init(&S, minp, maxp);
3165 39059687 : while ( (p = u_forprime_next_fast(&S)) )
3166 : {
3167 39059687 : ulong q = n / p;
3168 39059687 : if (n % p == 0)
3169 : {
3170 1246799 : P[i++] = p; n = q;
3171 1246799 : if (q <= p || (n <= maxp && PRIMES_search(n) > 0)) { P[i++] = n; break; }
3172 : }
3173 37812888 : else if (q <= p) { P[i++] = n; break; } /* n <= p^2: n is now prime */
3174 : }
3175 761722 : if (i == 1) return NULL;
3176 761722 : setlg(P, i); return P;
3177 : }
3178 : static GEN
3179 127536 : Z_factor_primes(GEN N, ulong minp, ulong maxp)
3180 : {
3181 : forprime_t S;
3182 127536 : ulong p, n = 0;
3183 : long i;
3184 : GEN P;
3185 127536 : if (lgefint(N) == 3) return factoru_primes(uel(N,2), minp, maxp);
3186 19086 : u_forprime_init(&S, minp, maxp);
3187 19086 : P = cgetg(expi(N) + 1, t_VECSMALL); i = 1;
3188 11437046 : while ( (p = u_forprime_next_fast(&S)) )
3189 : {
3190 : int stop;
3191 11437050 : long v = Z_lvalrem_stop(&N, p, &stop);
3192 11437046 : if (v) P[i++] = p;
3193 11437046 : if (stop)
3194 : {
3195 320 : if (!equali1(N)) P[i++] = uel(N,2);
3196 320 : goto END;
3197 : }
3198 11436726 : if (v && lgefint(N) == 3) { n = uel(N,2); break; }
3199 : }
3200 23972052 : if (n) while ( (p = u_forprime_next_fast(&S)) )
3201 : {
3202 23972052 : ulong q = n / p;
3203 23972052 : if (n % p == 0)
3204 : {
3205 64015 : P[i++] = p; n = q;
3206 64015 : if (q <= p || (n <= maxp && PRIMES_search(n) > 0)) { P[i++] = n; break; }
3207 : }
3208 23908037 : else if (q <= p) { P[i++] = n; break; } /* n <= p^2: n is now prime */
3209 : }
3210 0 : END:
3211 19086 : if (i == 1) return NULL;
3212 19086 : setlg(P, i); return P;
3213 : }
3214 :
3215 : /* N != 0. Product of odd prime divisors less than
3216 : * min(*pLIM, factorlimit) [WARNING!];
3217 : * with lim <= *pLIM < 2*lim and *pLIM prime
3218 : * Assume lim >= 128. Better for efficiency if N >= lim^2. */
3219 : static ulong
3220 809154 : u_oddprimedivisors_gcd(ulong N, ulong lim, ulong *pLIM)
3221 : {
3222 809154 : GEN PR = prodprimes(), LIM = prodprimeslim();
3223 809154 : long b = minss(lg(PR)-1, expu(lim)-6);
3224 : /* 2^{b+6} <= lim < 2^{b+7}, b >= 1 */
3225 809154 : *pLIM = LIM[b]; return ugcd(N, umodiu(gel(PR,b), N));
3226 : }
3227 : /* not GC-clean */
3228 : static GEN
3229 128213 : Z_oddprimedivisors_gcd(GEN N, ulong lim, ulong *pLIM)
3230 : {
3231 128213 : GEN PR = prodprimes(), LIM = prodprimeslim();
3232 128213 : long b = minss(lg(PR)-1, expu(lim)-6);
3233 128213 : *pLIM = LIM[b]; return gcdii(N, gel(PR,b));
3234 : }
3235 :
3236 : /* Assume lim >= 128 and N odd. */
3237 : static GEN
3238 126933 : Z_oddprimedivisors_fast(GEN N, ulong lim)
3239 : {
3240 126933 : pari_sp av = avma;
3241 126933 : GEN Nr = Z_oddprimedivisors_gcd(N, lim, &lim);
3242 126933 : GEN P = Z_factor_primes(Nr, 3, lim);
3243 126933 : return P? P: gc_NULL(av);
3244 : }
3245 : /* return mask with bit 0, 1, 2 set if respectively 3, 5, 7 divide n */
3246 : static int
3247 15807641 : u_357_divides(ulong n)
3248 : { /* vector (105, i, n = i-1; !(n%3) + 2 * !(n%5) + 4 * !(n%7)) */
3249 15807641 : const unsigned int tab[] = {
3250 : 7,0,0,1,0,2,1,4,0,1,2,0,1,0,4,3,0,0,1,0,2,5,0,0,1,2,0,1,4,0,3,0,0,1,0,6,1,0,
3251 : 0,1,2,0,5,0,0,3,0,0,1,4,2,1,0,0,1,2,4,1,0,0,3,0,0,5,0,2,1,0,0,1,6,0,1,0,0,3,
3252 : 0,4,1,0,2,1,0,0,5,2,0,1,0,0,3,4,0,1,0,2,1,0,4,1,2,0,1,0,0};
3253 15807641 : return tab[n % 105UL];
3254 : }
3255 :
3256 : static GEN
3257 17181630 : factoru_result(GEN P, GEN E, long i)
3258 : {
3259 17181630 : GEN P2, E2, f = cgetg(3,t_VEC);
3260 17181261 : gel(f,1) = P2 = cgetg(i, t_VECSMALL);
3261 17180670 : gel(f,2) = E2 = cgetg(i, t_VECSMALL);
3262 52490595 : while (--i >= 1) { P2[i] = P[i]; E2[i] = E[i]; }
3263 17180555 : return f;
3264 : }
3265 :
3266 : /* Factor n and output [p,e] where
3267 : * p, e are vecsmall with n = prod{p[i]^e[i]}. If all != 0:
3268 : * if pU1 is not NULL, set *pU1 and *pU2 so that unfactored composite is
3269 : * U1^U2 with U1 not a pure power; else include it in factorization */
3270 : static GEN
3271 17613518 : factoru_sign(ulong n, ulong all, long hint, ulong *pU1, ulong *pU2)
3272 : {
3273 17613518 : pari_sp av = avma;
3274 17613518 : ulong ALL, p, lim = 0;
3275 17613518 : long i, oldi = -1;
3276 : forprime_t S;
3277 : GEN E, P;
3278 :
3279 17613518 : if (pU1) *pU1 = *pU2 = 1;
3280 17613518 : if (n == 0) retmkvec2(mkvecsmall(0), mkvecsmall(1));
3281 17613518 : if (n == 1) retmkvec2(cgetg(1,t_VECSMALL), cgetg(1,t_VECSMALL));
3282 :
3283 : /* enough room to store <= 15 primes and exponents (OK if n < 2^64) */
3284 17181045 : (void)new_chunk(3 + 16*2);
3285 17181009 : P = cgetg(16, t_VECSMALL); i = 1;
3286 17180916 : E = cgetg(16, t_VECSMALL);
3287 17180551 : ALL = all? all: ULONG_MAX; /* (!all || all > n) == ALL > n */
3288 17180551 : if (ALL > 2)
3289 : {
3290 : ulong maxp;
3291 17180517 : long v = vals(n);
3292 17180477 : if (v)
3293 : {
3294 11638672 : P[1] = 2; E[1] = v; i = 2;
3295 11638672 : n >>= v; if (n == 1) goto END;
3296 : }
3297 14053473 : if (ALL > 3)
3298 : {
3299 14053932 : int mask = u_357_divides(n);
3300 14054408 : if (mask)
3301 : {
3302 9458704 : if (mask & 1)
3303 6001299 : { P[i] = 3; E[i] = 1 + u_lvalrem(n / 3, 3, &n); i++;
3304 6001355 : if (n == 1) goto END; }
3305 7144607 : if ((mask & 2) && ALL > 5)
3306 3468839 : { P[i] = 5; E[i] = 1 + u_lvalrem(n / 5, 5, &n); i++;
3307 3468855 : if (n == 1) goto END; }
3308 5578039 : if ((mask & 4) && ALL > 7)
3309 2079364 : { P[i] = 7; E[i] = 1 + u_lvalrem(n / 7, 7, &n); i++;
3310 2079371 : if (n == 1) goto END; }
3311 : }
3312 : }
3313 9013178 : maxp = maxprime();
3314 9014411 : if (n <= maxp && PRIMES_search(n) > 0) { P[i] = n; E[i] = 1; i++; goto END; }
3315 3254154 : lim = all? all-1: tridiv_boundu(n);
3316 3254152 : if (lim >= 128 && n >= 691 * 691) /* expu(lim) >= 7 */
3317 7777 : { /* fast trial division */
3318 804400 : ulong gcdlim, gcd, sqrtn = usqrt(n);
3319 : GEN Q;
3320 804400 : lim = minuu(lim, sqrtn);
3321 804400 : gcd = u_oddprimedivisors_gcd(n, lim, &gcdlim);
3322 804400 : Q = factoru_primes(gcd, 11, gcdlim);
3323 804400 : maxp = GP_DATA->factorlimit;
3324 804400 : if (Q)
3325 : {
3326 789362 : long j, l = lg(Q);
3327 789362 : int stop = 0;
3328 2578953 : for (j = 1; j < l; j++)
3329 : {
3330 1789597 : ulong p = uel(Q,j);
3331 1789597 : if (all && p >= all) { stop = 1; break; }
3332 1789591 : E[i] = u_lvalrem_stop(&n, p, &stop); /* > 0 */
3333 1789591 : P[i] = p; i++;
3334 : }
3335 825110 : if (n == 1) goto END;
3336 39113 : if (stop || (n <= maxp && PRIMES_search(n) > 0))
3337 35748 : { P[i] = n; E[i] = 1; i++; goto END; }
3338 : }
3339 15038 : else if (lim == sqrtn && lim <= maxp)
3340 10626 : { P[i] = n; E[i] = 1; i++; goto END; }
3341 : }
3342 : else
3343 : { /* naive trial division */
3344 2449752 : maxp = lim;
3345 2449752 : u_forprime_init(&S, 11, lim);
3346 17695257 : while ( (p = u_forprime_next_fast(&S)) )
3347 : {
3348 : int stop;
3349 : /* tiny integers without small factors are often primes */
3350 17694893 : if (p == 673)
3351 : {
3352 2449566 : if (uisprime_661(n)) { P[i] = n; E[i] = 1; i++; goto END; }
3353 190 : oldi = i;
3354 : }
3355 17694893 : v = u_lvalrem_stop(&n, p, &stop);
3356 17694886 : if (v) { P[i] = p; E[i] = v; i++; }
3357 17694886 : if (stop)
3358 : {
3359 2449376 : if (n != 1) { P[i] = n; E[i] = 1; i++; }
3360 2449376 : goto END;
3361 : }
3362 : }
3363 : }
3364 8168 : if (lim > maxp)
3365 : { /* second pass usually empty, outside fast trial division range */
3366 : long v;
3367 6 : u_forprime_init(&S, maxp+1, lim);
3368 5296866 : while ((p = u_forprime_next(&S)))
3369 : {
3370 : int stop;
3371 5296900 : v = u_lvalrem_stop(&n, p, &stop);
3372 5296866 : if (v) { P[i] = p; E[i] = v; i++; }
3373 5296866 : if (stop)
3374 : {
3375 6 : if (n != 1) { P[i] = n; E[i] = 1; i++; }
3376 6 : goto END;
3377 : }
3378 : }
3379 : }
3380 : }
3381 : /* if i > oldi (includes oldi = -1) we don't know that n is composite */
3382 8162 : if (all)
3383 : { /* smallfact: look for easy pure powers then stop */
3384 : #ifdef LONG_IS_64BIT
3385 1080 : ulong mask = all > 563 ? (all > 7129 ? 1: 3): 7;
3386 : #else
3387 15 : ulong mask = all > 22 ? (all > 83 ? 1: 3): 7;
3388 : #endif
3389 1095 : long k = 1, ex;
3390 1588 : while (uissquareall(n, &n)) k <<= 1;
3391 1108 : while ( (ex = uis_357_power(n, &n, &mask)) ) k *= ex;
3392 1095 : if (pU1 && (i == oldi || !uisprime_nosmall(n, lim)))
3393 266 : { *pU1 = n; *pU2 = (ulong)k; }
3394 : else
3395 829 : { P[i] = n; E[i] = k; i++; }
3396 1095 : goto END;
3397 : }
3398 : /* we don't known that n is composite ? */
3399 7067 : if (oldi != i && uisprime_nosmall(n, lim)) { P[i]=n; E[i]=1; i++; goto END; }
3400 :
3401 : {
3402 : GEN perm;
3403 863 : ifac_factoru(utoipos(n), hint, P, E, &i);
3404 863 : setlg(P, i);
3405 863 : perm = vecsmall_indexsort(P);
3406 863 : P = vecsmallpermute(P, perm);
3407 863 : E = vecsmallpermute(E, perm);
3408 : }
3409 17182456 : END:
3410 17182456 : set_avma(av); return factoru_result(P, E, i);
3411 : }
3412 : GEN
3413 3741503 : factoru(ulong n)
3414 3741503 : { return factoru_sign(n, 0, decomp_default_hint, NULL, NULL); }
3415 :
3416 : ulong
3417 0 : radicalu(ulong n)
3418 : {
3419 0 : pari_sp av = avma;
3420 0 : return gc_long(av, zv_prod(gel(factoru(n),1)));
3421 : }
3422 :
3423 : long
3424 54194 : moebiusu_fact(GEN f)
3425 : {
3426 54194 : GEN E = gel(f,2);
3427 54194 : long i, l = lg(E);
3428 93569 : for (i = 1; i < l; i++)
3429 57834 : if (E[i] > 1) return 0;
3430 35735 : return odd(l)? 1: -1;
3431 : }
3432 :
3433 : long
3434 2487299 : moebiusu(ulong n)
3435 : {
3436 : pari_sp av;
3437 : long s, v, test_prime;
3438 : ulong p, lim;
3439 : int mask;
3440 :
3441 2487299 : switch(n)
3442 : {
3443 0 : case 0: (void)check_arith_non0(gen_0,"moebius");/*error*/
3444 568237 : case 1: return 1;
3445 106667 : case 2: return -1;
3446 : }
3447 : /* n > 2 */
3448 1832112 : p = n & 3; if (!p) return 0;
3449 1743874 : if (p == 2) { n >>= 1; s = -1; } else s = 1;
3450 1743874 : mask = u_357_divides(n);
3451 1762802 : if (mask)
3452 : {
3453 677844 : if (mask & 1) { n /= 3; s = -s; if (n % 3 == 0) return 0; }
3454 613916 : if (mask & 2) { n /= 5; s = -s; if (n % 5 == 0) return 0; }
3455 572553 : if (mask & 4) { n /= 7; s = -s; if (n % 7 == 0) return 0; }
3456 : }
3457 1622366 : if (n <= maxprimelim() && PRIMES_search(n) > 0) return -s;
3458 661399 : av = avma; lim = tridiv_boundu(n);
3459 672117 : if (n >= 691 * 691)
3460 : {
3461 4507 : ulong gcdlim, gcd, sqrtn = usqrt(n);
3462 : GEN P;
3463 4507 : lim = minuu(sqrtn, lim);
3464 4507 : gcd = u_oddprimedivisors_gcd(n, lim, &gcdlim);
3465 4507 : if (gcd != 1)
3466 : {
3467 3523 : n /= gcd;
3468 3846 : if (ugcd(gcd, n) != 1) return 0;
3469 : }
3470 4362 : P = factoru_primes(gcd, 11, gcdlim);
3471 4362 : if (P && odd(lg(P) - 1)) s = -s;
3472 4362 : if (n == 1) return gc_long(av, s);
3473 4060 : if (lim == sqrtn && lim <= GP_DATA->factorlimit) return gc_long(av, -s);
3474 4039 : test_prime = 1;
3475 : }
3476 : else
3477 : {
3478 : forprime_t S;
3479 667610 : u_forprime_init(&S, 3, lim);
3480 666022 : test_prime = 0;
3481 7417388 : while ((p = u_forprime_next_fast(&S)))
3482 : {
3483 : int stop;
3484 : /* tiny integers without small factors are often primes */
3485 7416328 : if (p == 673)
3486 : {
3487 0 : test_prime = 0;
3488 665023 : if (uisprime_661(n)) return gc_long(av,-s);
3489 : }
3490 7416328 : v = u_lvalrem_stop(&n, p, &stop);
3491 7415431 : if (v) {
3492 635317 : if (v > 1) return gc_long(av,0);
3493 593329 : test_prime = 1;
3494 593329 : s = -s;
3495 : }
3496 7373443 : if (stop) return gc_long(av, n==1? s: -s);
3497 : }
3498 : }
3499 4039 : set_avma(av);
3500 4039 : if (test_prime && uisprime_661(n)) return -s;
3501 : else
3502 : {
3503 663 : long t = ifac_moebiusu(utoipos(n));
3504 663 : set_avma(av);
3505 663 : if (t == 0) return 0;
3506 663 : return (s == t)? 1: -1;
3507 : }
3508 : }
3509 :
3510 : long
3511 58815 : moebius(GEN n)
3512 : {
3513 58815 : pari_sp av = avma;
3514 : GEN F;
3515 : ulong p, lim, n357;
3516 : long i, l, s, v, copy;
3517 : int mask;
3518 :
3519 58815 : if ((F = check_arith_non0(n,"moebius")))
3520 : {
3521 : GEN E;
3522 728 : F = clean_Z_factor(F);
3523 728 : E = gel(F,2);
3524 728 : l = lg(E);
3525 1428 : for(i = 1; i < l; i++)
3526 980 : if (!equali1(gel(E,i))) return gc_long(av,0);
3527 448 : return gc_long(av, odd(l)? 1: -1);
3528 : }
3529 58050 : if (lgefint(n) == 3) return moebiusu(uel(n,2));
3530 1608 : p = mod4(n); if (!p) return 0;
3531 1408 : copy = s = 1;
3532 1408 : if (p == 2)
3533 : {
3534 358 : n = shifti(n, -1);
3535 358 : copy = 0; s = -1;
3536 : }
3537 1408 : n357 = umodiu(n, 9 * 25 * 49);
3538 1408 : mask = u_357_divides(n357);
3539 1408 : if (mask)
3540 : {
3541 764 : ulong m = 1;
3542 764 : if (mask & 1) { m = 3; s = -s; }
3543 764 : if (mask & 2) { m *= 5; s = -s; }
3544 764 : if (mask & 4) { m *= 7; s = -s; }
3545 764 : if (u_357_divides(n357 / m)) return gc_long(av, 0);
3546 549 : copy = 0; n = diviuexact(n, m);
3547 : }
3548 1193 : if (copy) n = icopy(n);
3549 701 : else if (lgefint(n) == 3) return gc_long(av, s * moebiusu(uel(n,2)));
3550 883 : setabssign(n); lim = tridiv_bound(n);
3551 883 : if (lim >= 128)
3552 : {
3553 : ulong gcdlim;
3554 883 : GEN gcd = Z_oddprimedivisors_gcd(n, lim, &gcdlim);
3555 883 : if (!equali1(gcd))
3556 : {
3557 : GEN P;
3558 622 : n = diviiexact(n, gcd);
3559 808 : if (!equali1(gcdii(gcd, n))) return gc_long(av, 0);
3560 603 : P = Z_factor_primes(gcd, 11, gcdlim);
3561 603 : if (P)
3562 : {
3563 603 : if (odd(lg(P) - 1)) s = -s;
3564 603 : if (is_pm1(n)) return gc_long(av, s);
3565 1158 : if (lim <= GP_DATA->factorlimit &&
3566 741 : cmpii(sqru(lim), n) >= 0) return gc_long(av, -s); /* n prime */
3567 : }
3568 : }
3569 : }
3570 : else
3571 : {
3572 : forprime_t S;
3573 0 : u_forprime_init(&S, 3, lim);
3574 0 : while ((p = u_forprime_next_fast(&S)))
3575 : {
3576 : int stop;
3577 0 : v = Z_lvalrem_stop(&n, p, &stop);
3578 0 : if (v)
3579 : {
3580 0 : if (v > 1) return gc_long(av,0);
3581 0 : s = -s;
3582 : }
3583 0 : if (stop) return gc_long(av, is_pm1(n)? s: -s);
3584 : }
3585 : }
3586 678 : l = lg(primetab);
3587 682 : for (i = 1; i < l; i++)
3588 : {
3589 7 : v = Z_pvalrem(n, gel(primetab,i), &n);
3590 7 : if (v)
3591 : {
3592 7 : if (v > 1) return gc_long(av,0);
3593 5 : s = -s;
3594 5 : if (is_pm1(n)) return gc_long(av,s);
3595 : }
3596 : }
3597 675 : if (ifac_isprime(n)) return gc_long(av,-s);
3598 : /* large composite without small factors */
3599 236 : v = ifac_moebius(n);
3600 236 : return gc_long(av, s < 0? -v: v); /* correct also if v==0 */
3601 : }
3602 :
3603 : long
3604 1708 : ispowerful(GEN n)
3605 : {
3606 1708 : pari_sp av = avma;
3607 : GEN F;
3608 : ulong p, lim, n357;
3609 : long i, l, v;
3610 1708 : int mask, copy = 1;
3611 :
3612 1708 : if ((F = check_arith_all(n, "ispowerful")))
3613 : {
3614 742 : GEN p, P = gel(F,1), E = gel(F,2);
3615 742 : if (lg(P) == 1) return 1; /* 1 */
3616 728 : p = gel(P,1);
3617 728 : if (!signe(p)) return 1; /* 0 */
3618 707 : i = is_pm1(p)? 2: 1; /* skip -1 */
3619 707 : l = lg(E);
3620 980 : for (; i < l; i++)
3621 847 : if (equali1(gel(E,i))) return 0;
3622 133 : return 1;
3623 : }
3624 966 : if (!signe(n)) return 1;
3625 952 : if (mod4(n) == 2) return 0;
3626 623 : n357 = umodiu(n, 9 * 25 * 49);
3627 623 : mask = u_357_divides(n357);
3628 623 : if (mask)
3629 : {
3630 315 : if ((mask & 1) && n357 % 9) return 0;
3631 203 : if ((mask & 2) && n357 % 25) return 0;
3632 126 : if ((mask & 4) && n357 % 49) return 0;
3633 98 : if (mask & 1) (void)Z_lvalrem(diviuexact(n,9), 3, &n);
3634 98 : if (mask & 2) (void)Z_lvalrem(diviuexact(n,25), 5, &n);
3635 98 : if (mask & 4) (void)Z_lvalrem(diviuexact(n,49), 7, &n);
3636 98 : copy = 0;
3637 : }
3638 406 : if (!mod2(n)) { n = shifti(n, -vali(n)); copy = 0; }
3639 406 : if (is_pm1(n)) return gc_long(av, 1);
3640 238 : if (copy) n = icopy(n);
3641 238 : setabssign(n); lim = tridiv_bound(n);
3642 238 : if (cmpiu(n, 691 * 691) >= 0)
3643 : {
3644 70 : ulong gcdlim, sqrtn = 0;
3645 : GEN gcd;
3646 70 : if (lgefint(n) == 3)
3647 : {
3648 6 : sqrtn = usqrt(n[2]);
3649 6 : lim = minuu(sqrtn, lim);
3650 : }
3651 70 : gcd = Z_oddprimedivisors_gcd(n, lim, &gcdlim);
3652 70 : if (!equali1(gcd))
3653 : {
3654 : GEN r;
3655 70 : n = diviiexact(n, gcd);
3656 70 : n = dvmdii(n, gcd, &r);
3657 70 : if (r != gen_0) return gc_long(av, 0);
3658 70 : n = Z_ppo(n, gcd);
3659 : }
3660 : /* prime divisors > gcdlim */
3661 70 : if (equali1(n)) return gc_long(av, 1);
3662 7 : if (sqrtn && gcdlim >= sqrtn) return gc_long(av, 0); /* prime */
3663 : }
3664 : else
3665 : {
3666 : forprime_t S;
3667 168 : u_forprime_init(&S, 3, lim);
3668 378 : while ((p = u_forprime_next_fast(&S)))
3669 : {
3670 : int stop;
3671 378 : v = Z_lvalrem_stop(&n, p, &stop);
3672 546 : if (v == 1) return gc_long(av,0);
3673 378 : if (stop) return gc_long(av, is_pm1(n)); /* n > 1 is now prime */
3674 : }
3675 : }
3676 7 : l = lg(primetab);
3677 7 : for (i = 1; i < l; i++)
3678 : {
3679 0 : v = Z_pvalrem(n, gel(primetab,i), &n);
3680 0 : if (v)
3681 : {
3682 0 : if (v == 1) return gc_long(av,0);
3683 0 : if (is_pm1(n)) return gc_long(av,1);
3684 : }
3685 : }
3686 : /* no need to factor: must be p^2 or not powerful */
3687 7 : if (cmpii(powuu(lim+1, 3), n) > 0) return gc_long(av, Z_issquare(n));
3688 :
3689 7 : if (ifac_isprime(n)) return gc_long(av,0);
3690 : /* large composite without small factors */
3691 7 : return gc_long(av, ifac_ispowerful(n));
3692 : }
3693 :
3694 : ulong
3695 0 : coreu_fact(GEN f)
3696 : {
3697 0 : GEN P = gel(f,1), E = gel(f,2);
3698 0 : long i, l = lg(P), m = 1;
3699 0 : for (i = 1; i < l; i++)
3700 : {
3701 0 : ulong p = P[i], e = E[i];
3702 0 : if (e & 1) m *= p;
3703 : }
3704 0 : return m;
3705 : }
3706 :
3707 : /* d = a squarefree divisor of n. Return n / (n, d^oo)
3708 : * and set *pcore = \prod_{p | (n,d), v_p(n) odd} p
3709 : * Simpified form of Z_cba. */
3710 : static GEN
3711 327 : core_init_from_squarefree(GEN n, GEN d, GEN *pcore)
3712 : {
3713 327 : GEN c = gen_1;
3714 : long v;
3715 :
3716 327 : if (equali1(d)) { *pcore = c; return n; }
3717 260 : v = Z_pvalrem(n, d, &n);
3718 59 : for (;; v++)
3719 59 : { /* d^v divides "original n" */
3720 319 : GEN newd = gcdii(n, d); /* newd^{v+1} divides original n */
3721 319 : if (!equalii(d, newd))
3722 : { /* new d loses primes dividing original n to exact power v */
3723 310 : if (odd(v)) c = mulii(c, diviiexact(d, newd)); /* lost primes */
3724 310 : d = newd; if (equali1(d)) break;
3725 : }
3726 63 : if (equalii(d, n))
3727 : {
3728 4 : if (odd(v + 1)) c = mulii(c, d);
3729 4 : *pcore = c; return gen_1;
3730 : }
3731 59 : n = diviiexact(n, d);
3732 : }
3733 256 : *pcore = c; return n;
3734 : }
3735 : static ulong
3736 247 : coreu_init_from_squarefree(ulong n, ulong d, ulong *pcore)
3737 : {
3738 247 : ulong c = 1;
3739 : long v;
3740 :
3741 247 : if (d == 1) { *pcore = c; return n; }
3742 205 : v = u_lvalrem(n, d, &n);
3743 18 : for (;; v++)
3744 18 : { /* d^v divides "original n" */
3745 223 : ulong newd = ugcd(n, d); /* newd^{v+1} divides original n */
3746 223 : if (d != newd)
3747 : { /* new d loses primes dividing original n to exact power v */
3748 211 : if (odd(v)) c *= d / newd; /* lost primes */
3749 211 : d = newd; if (d == 1) break;
3750 : }
3751 90 : if (d == n)
3752 : {
3753 72 : if (odd(v + 1)) c *= d;
3754 72 : *pcore = c; return 1;
3755 : }
3756 18 : n /= d;
3757 : }
3758 133 : *pcore = c; return n;
3759 : }
3760 :
3761 : ulong
3762 62112 : coreu(ulong n)
3763 : {
3764 : ulong p, lim, m;
3765 : long v;
3766 :
3767 62112 : if (!n) return 0;
3768 62112 : m = 1;
3769 62112 : v = vals(n);
3770 62112 : if (v)
3771 : {
3772 4133 : n >>= v;
3773 4133 : if (odd(v)) m = 2;
3774 : }
3775 62112 : v = u_lvalrem(n, 3, &n); if (odd(v)) m *= 3;
3776 62112 : v = u_lvalrem(n, 5, &n); if (odd(v)) m *= 5;
3777 62112 : v = u_lvalrem(n, 7, &n); if (odd(v)) m *= 7;
3778 62112 : if (n == 1) return m;
3779 3320 : if (n <= maxprimelim() && PRIMES_search(n) > 0) return m * n;
3780 408 : lim = tridiv_boundu(n);
3781 408 : if (n >= 691 * 691)
3782 : {
3783 247 : ulong mpart, gcd, gcdlim, sqrtn = usqrt(n);
3784 247 : lim = minuu(sqrtn, lim);
3785 247 : gcd = u_oddprimedivisors_gcd(n, lim, &gcdlim);
3786 247 : n = coreu_init_from_squarefree(n, gcd, &mpart);
3787 247 : m *= mpart;
3788 266 : if (n == 1) return m;
3789 : /* n has no prime divisor <= gcdlim */
3790 103 : if ((lim == sqrtn && lim <= GP_DATA->factorlimit)
3791 96 : || (gcdlim + 1) * (gcdlim + 1) > n)
3792 19 : return m * n; /* prime */
3793 : }
3794 : else
3795 : {
3796 : forprime_t S;
3797 161 : u_forprime_init(&S, 11, lim);
3798 3633 : while ((p = u_forprime_next_fast(&S)))
3799 : {
3800 : int stop;
3801 3633 : v = u_lvalrem_stop(&n, p, &stop);
3802 3633 : if (v & 1) m *= p;
3803 3633 : if (stop) return n == 1? m: m * n; /* n > 1 is now prime */
3804 : }
3805 : }
3806 84 : if (uisprime_661(n)) return m * n;
3807 : else
3808 : {
3809 84 : pari_sp av = avma;
3810 84 : m *= itou(ifac_core(utoipos(n)));
3811 84 : return gc_ulong(av, m);
3812 : }
3813 : }
3814 :
3815 : GEN
3816 708659 : core(GEN n)
3817 : {
3818 708659 : pari_sp av = avma;
3819 : GEN m, mpart, gcd, F;
3820 : ulong lim, gcdlim, mask, m0;
3821 : long i, l, v, s;
3822 708659 : int copy = 1;
3823 :
3824 708659 : if ((F = check_arith_all(n, "core")))
3825 : {
3826 646236 : GEN p, x, P = gel(F,1), E = gel(F,2);
3827 646236 : long j = 1;
3828 646236 : if (lg(P) == 1) return gen_1;
3829 646208 : p = gel(P,1);
3830 646208 : if (!signe(p)) return gen_0;
3831 646166 : l = lg(P); x = cgetg(l, t_VEC);
3832 2283030 : for (i = 1; i < l; i++)
3833 1636865 : if (mpodd(gel(E,i))) gel(x,j++) = gel(P,i);
3834 646165 : setlg(x, j); return ZV_prod(x);
3835 : }
3836 62410 : s = signe(n);
3837 62410 : if (!s) return gen_0;
3838 62382 : if (lgefint(n) == 3)
3839 : {
3840 61975 : ulong c = coreu(uel(n,2));
3841 61975 : return s < 0? utoineg(c): utoipos(c);
3842 : }
3843 407 : v = vali(n); m0 = 1;
3844 408 : if (v)
3845 : {
3846 123 : n = shifti(n, -v); if (odd(v)) m0 *= 2;
3847 123 : copy = 0;
3848 : }
3849 408 : if ((mask = u_357_divides(umodiu(n, 3 * 5 * 7))))
3850 : {
3851 275 : if (mask & 1) { v = Z_lvalrem(n, 3, &n); if (odd(v)) m0 *= 3; }
3852 275 : if (mask & 2) { v = Z_lvalrem(n, 5, &n); if (odd(v)) m0 *= 5; }
3853 275 : if (mask & 4) { v = Z_lvalrem(n, 7, &n); if (odd(v)) m0 *= 7; ; }
3854 275 : copy = 0;
3855 : }
3856 408 : if (copy) n = absi(n); /* ifac_core destroys n */
3857 289 : else if (lgefint(n) == 3)
3858 : {
3859 81 : ulong c = coreu(uel(n,2));
3860 81 : m = muluu(m0, c); if (s < 0) setsigne(m, -1);
3861 81 : return gerepileuptoint(av, m);
3862 : }
3863 327 : setabssign(n); lim = tridiv_bound(n);
3864 : /* n >= 691^2 */
3865 327 : gcd = Z_oddprimedivisors_gcd(n, lim, &gcdlim);
3866 327 : n = core_init_from_squarefree(n, gcd, &mpart);
3867 327 : m = mului(m0, mpart); if (s < 0) setsigne(m, -1);
3868 327 : if (equali1(n)) return gerepileuptoint(av, m);
3869 : /* n has no prime divisor <= gcdlim */
3870 276 : if (cmpii(sqru(gcdlim + 1), n) > 0)
3871 2 : return gerepileuptoint(av, mulii(m, n)); /* prime */
3872 274 : l = lg(primetab);
3873 750 : for (i = 1; i < l; i++)
3874 : {
3875 478 : GEN q = gel(primetab,i);
3876 478 : v = Z_pvalrem(n, q, &n);
3877 478 : if (v)
3878 : {
3879 8 : if (v & 1) m = mulii(m, q);
3880 8 : if (is_pm1(n)) return gerepileuptoint(av, m);
3881 : }
3882 : }
3883 272 : if (!ifac_isprime(n)) n = ifac_core(n); /* composite without small factors */
3884 272 : return gerepileuptoint(av, mulii(m, n));
3885 : }
3886 :
3887 : long
3888 0 : Z_issmooth(GEN m, ulong lim)
3889 : {
3890 0 : pari_sp av = avma;
3891 0 : ulong p = 2;
3892 : forprime_t S;
3893 0 : u_forprime_init(&S, 2, lim);
3894 0 : while ((p = u_forprime_next_fast(&S)))
3895 : {
3896 : int stop;
3897 0 : (void)Z_lvalrem_stop(&m, p, &stop);
3898 0 : if (stop) return gc_long(av, abscmpiu(m,lim) <= 0);
3899 : }
3900 0 : return gc_long(av,0);
3901 : }
3902 :
3903 : GEN
3904 178649 : Z_issmooth_fact(GEN m, ulong lim)
3905 : {
3906 178649 : pari_sp av = avma;
3907 : GEN F, P, E;
3908 : ulong p;
3909 178649 : long i = 1, l = expi(m)+1;
3910 : forprime_t S;
3911 178641 : P = cgetg(l, t_VECSMALL);
3912 178576 : E = cgetg(l, t_VECSMALL); F = mkmat2(P,E);
3913 178566 : if (l == 1) return F; /* m == 1 */
3914 178524 : u_forprime_init(&S, 2, lim);
3915 44193391 : while ((p = u_forprime_next_fast(&S)))
3916 : {
3917 : int stop;
3918 44151856 : long v = Z_lvalrem_stop(&m, p, &stop);
3919 44152639 : if (v) { P[i] = p; E[i] = v; i++; }
3920 44152639 : if (stop)
3921 : {
3922 137802 : if (abscmpiu(m,lim) > 0) break;
3923 111610 : if (m[2] > 1) { P[i] = m[2]; E[i] = 1; i++; }
3924 111610 : setlg(P, i);
3925 111692 : setlg(E, i); return gc_const((pari_sp)F, F);
3926 : }
3927 : }
3928 66057 : return gc_NULL(av);
3929 : }
3930 :
3931 : /* Is (a mod p^e) a K-th power ? p is prime and e > 0 */
3932 : static int
3933 798 : Zp_ispower(GEN a, GEN L, GEN K, GEN p, long e)
3934 : {
3935 798 : GEN t = gen_0;
3936 798 : long v = Z_pvalrem(a, p, &a), d = e - v;
3937 798 : if (d > 0)
3938 : { /* is a mod p^d a K-th power ? a a p-unit */
3939 : ulong r;
3940 763 : v = uabsdivui_rem(v, K, &r); if (r) return 0;
3941 651 : if (d == 1)
3942 : { /* mod p: faster */
3943 567 : if (!Fp_ispower(a, K, p)) return 0;
3944 560 : if (L) t = Fp_sqrtn(a, K, p, NULL);
3945 : }
3946 : else
3947 : { /* mod p^{2 +} */
3948 84 : if (!ispower(cvtop(a, p, d), K, L? &t: NULL)) return 0;
3949 63 : if (L) t = gtrunc(t);
3950 : }
3951 623 : if (L && v) t = mulii(t, powiu(p, v));
3952 : }
3953 686 : if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
3954 686 : return 1;
3955 : }
3956 : long
3957 756 : Zn_ispower(GEN a, GEN q, GEN K, GEN *pt)
3958 : {
3959 : GEN L, N;
3960 : pari_sp av;
3961 : long e, i, l;
3962 : ulong pp, lim;
3963 :
3964 756 : if (!signe(a))
3965 : {
3966 91 : if (pt) {
3967 91 : GEN t = cgetg(3, t_INTMOD);
3968 91 : gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
3969 : }
3970 91 : return 1;
3971 : }
3972 : /* a != 0 */
3973 665 : av = avma;
3974 :
3975 665 : if (typ(q) != t_INT) /* integer factorization */
3976 : {
3977 0 : GEN P = gel(q,1), E = gel(q,2);
3978 0 : l = lg(P);
3979 0 : L = pt? vectrunc_init(l): NULL;
3980 0 : for (i = 1; i < l; i++)
3981 : {
3982 0 : GEN p = gel(P,i);
3983 0 : long e = itos(gel(E,i));
3984 0 : if (!Zp_ispower(a, L, K, p, e)) return gc_long(av,0);
3985 : }
3986 0 : goto END;
3987 : }
3988 665 : e = vali(q); if (e) q = shifti(q, -e);
3989 665 : if (!mod2(K) && kronecker(a, q) == -1) return gc_long(av,0);
3990 658 : L = pt? vectrunc_init(expi(q)+2): NULL;
3991 658 : if (e)
3992 : {
3993 469 : if (!Zp_ispower(a, L, K, gen_2, e)) return gc_long(av,0);
3994 455 : a = modii(a, q);
3995 : }
3996 644 : lim = tridiv_bound(q);
3997 644 : if (cmpiu(q, 691 * 691) >= 0)
3998 : {
3999 161 : ulong sqrtq = lgefint(q) == 3? usqrt(q[2]): 0;
4000 : GEN P;
4001 161 : if (sqrtq) lim = minuu(sqrtq, lim);
4002 161 : P = Z_oddprimedivisors_fast(q, lim);
4003 161 : if (P)
4004 : {
4005 103 : long nP = lg(P) - 1;
4006 103 : int stop = 0;
4007 206 : for (i = 1; i <= nP; i++)
4008 : {
4009 151 : ulong pp = uel(P,i);
4010 151 : e = Z_lvalrem_stop(&q, pp, &stop);
4011 151 : if (!Zp_ispower(a, L, K, utoipos(pp), e)) return gc_long(av,0);
4012 103 : a = modii(a, q);
4013 : }
4014 55 : if (stop)
4015 : {
4016 48 : if (!is_pm1(q) && !Zp_ispower(a, L, K, q, 1)) return gc_long(av,0);
4017 48 : goto END;
4018 : }
4019 : }
4020 58 : else if (lim == sqrtq && lim <= GP_DATA->factorlimit)
4021 : {
4022 0 : if (!Zp_ispower(a, L, K, q, 1)) return gc_long(av,0);
4023 0 : goto END;
4024 : }
4025 : }
4026 : else
4027 : {
4028 : forprime_t S;
4029 483 : u_forprime_init(&S, 3, lim);
4030 237174 : while ((pp = u_forprime_next(&S)))
4031 : {
4032 : int stop;
4033 236754 : e = Z_lvalrem_stop(&q, pp, &stop);
4034 236754 : if (!e) continue;
4035 63 : if (!Zp_ispower(a, L, K, utoipos(pp), e)) return gc_long(av,0);
4036 56 : a = modii(a, q);
4037 56 : if (stop)
4038 : {
4039 56 : if (!is_pm1(q) && !Zp_ispower(a, L, K, q, 1)) return gc_long(av,0);
4040 56 : goto END;
4041 : }
4042 : }
4043 : }
4044 485 : l = lg(primetab);
4045 485 : for (i = 1; i < l; i++)
4046 : {
4047 0 : GEN p = gel(primetab,i);
4048 0 : e = Z_pvalrem(q, p, &q);
4049 0 : if (!e) continue;
4050 0 : if (!Zp_ispower(a, L, K, p, e)) return gc_long(av,0);
4051 0 : if (is_pm1(q)) goto END;
4052 0 : a = modii(a, q);
4053 : }
4054 485 : N = gcdii(a,q);
4055 485 : if (!is_pm1(N))
4056 : {
4057 52 : if (ifac_isprime(N))
4058 : {
4059 34 : e = Z_pvalrem(q, N, &q);
4060 34 : if (!Zp_ispower(a, L, K, N, e)) return gc_long(av,0);
4061 0 : a = modii(a, q);
4062 : }
4063 : else
4064 : {
4065 18 : GEN part = ifac_start(N, 0);
4066 : for(;;)
4067 18 : {
4068 : long e;
4069 : GEN p;
4070 36 : if (!ifac_next(&part, &p, &e)) break;
4071 18 : e = Z_pvalrem(q, p, &q);
4072 18 : if (!Zp_ispower(a, L, K, p, e)) return gc_long(av,0);
4073 18 : a = modii(a, q);
4074 : }
4075 : }
4076 : }
4077 451 : if (!is_pm1(q))
4078 : {
4079 31 : if (ifac_isprime(q))
4080 : {
4081 4 : if (!Zp_ispower(a, L, K, q, 1)) return gc_long(av,0);
4082 : }
4083 : else
4084 : { /* icopy needed: ifac_next would destroy q */
4085 27 : GEN part = ifac_start(icopy(q), 0);
4086 : for(;;)
4087 43 : {
4088 : long e;
4089 : GEN p;
4090 70 : if (!ifac_next(&part, &p, &e)) break;
4091 52 : if (!Zp_ispower(a, L, K, p, e)) return gc_long(av,0);
4092 43 : a = modii(a, q);
4093 : }
4094 : }
4095 : }
4096 420 : END:
4097 546 : if (!pt) return gc_long(av, 1);
4098 476 : *pt = gerepileupto(av, chinese1_coprime_Z(L)); return 1;
4099 : }
4100 :
4101 :
4102 : /***********************************************************************/
4103 : /** **/
4104 : /** COMPUTING THE MATRIX OF PRIME DIVISORS AND EXPONENTS **/
4105 : /** **/
4106 : /***********************************************************************/
4107 : static GEN
4108 128365 : aux_end(GEN M, GEN n, long nb)
4109 : {
4110 128365 : GEN P,E, z = (GEN)avma;
4111 : long i;
4112 :
4113 128365 : guncloneNULL(n);
4114 128365 : P = cgetg(nb+1,t_COL);
4115 128365 : E = cgetg(nb+1,t_COL);
4116 846197 : for (i=nb; i; i--)
4117 : { /* allow a stackdummy in the middle */
4118 792964 : while (typ(z) != t_INT) z += lg(z);
4119 717832 : gel(E,i) = z; z += lg(z);
4120 717832 : gel(P,i) = z; z += lg(z);
4121 : }
4122 128365 : gel(M,1) = P;
4123 128365 : gel(M,2) = E;
4124 128365 : return sort_factor(M, (void*)&abscmpii, cmp_nodata);
4125 : }
4126 :
4127 : static void
4128 712838 : STORE(long *nb, GEN x, long e) { (*nb)++; (void)x; (void)utoipos(e); }
4129 : static void
4130 681523 : STOREu(long *nb, ulong x, long e) { STORE(nb, utoipos(x), e); }
4131 : static void
4132 31163 : STOREi(long *nb, GEN x, long e) { STORE(nb, icopy(x), e); }
4133 : /* no prime less than p divides n; return 1 if factored completely */
4134 : static int
4135 39425 : special_primes(GEN n, ulong p, long *nb, GEN T)
4136 : {
4137 39425 : long i, l = lg(T);
4138 39425 : if (l > 1)
4139 : { /* pp = square of biggest p tried so far */
4140 540 : long pp[] = { evaltyp(t_INT)|_evallg(4), 0,0,0 };
4141 540 : pari_sp av = avma; affii(sqru(p), pp); set_avma(av);
4142 :
4143 1184 : for (i = 1; i < l; i++)
4144 777 : if (dvdiiz(n, gel(T,i), n))
4145 : {
4146 329 : long k = 1; while (dvdiiz(n, gel(T,i), n)) k++;
4147 231 : STOREi(nb, gel(T,i), k);
4148 231 : if (abscmpii(pp, n) > 0)
4149 : {
4150 133 : if (!is_pm1(n)) STOREi(nb, n, 1);
4151 133 : return 1;
4152 : }
4153 : }
4154 : }
4155 39292 : return 0;
4156 : }
4157 :
4158 : /* factor(sn*|n|), where sn = -1 or 1.
4159 : * all != 0 : only look for prime divisors < all. If pU is not NULL,
4160 : * set it to unfactored composite */
4161 : static GEN
4162 13999600 : ifactor_sign(GEN n, ulong all, long hint, long sn, GEN *pU)
4163 : {
4164 : GEN M, N;
4165 : pari_sp av;
4166 13999600 : long nb = 0, nb0 = -1, i;
4167 : ulong lim;
4168 : forprime_t T;
4169 :
4170 13999600 : if (lgefint(n) == 3)
4171 : { /* small integer */
4172 13871250 : GEN f, Pf, Ef, P, E, F = cgetg(3, t_MAT);
4173 : ulong U1, U2;
4174 : long l;
4175 13871752 : av = avma;
4176 : /* enough room to store <= 15 primes and exponents (OK if n < 2^64) */
4177 13871752 : (void)new_chunk((15*3 + 15 + 1) * 2);
4178 13871744 : f = factoru_sign(uel(n,2), all, hint, pU? &U1: NULL, pU? &U2: NULL);
4179 13871956 : set_avma(av);
4180 13871954 : Pf = gel(f,1);
4181 13871954 : Ef = gel(f,2);
4182 13871954 : l = lg(Pf);
4183 13871954 : if (sn < 0)
4184 : { /* add sign */
4185 6519 : long L = l+1;
4186 6519 : gel(F,1) = P = cgetg(L, t_COL);
4187 6519 : gel(F,2) = E = cgetg(L, t_COL);
4188 6519 : gel(P,1) = gen_m1; P++;
4189 6519 : gel(E,1) = gen_1; E++;
4190 : }
4191 : else
4192 : {
4193 13865435 : gel(F,1) = P = cgetg(l, t_COL);
4194 13865455 : gel(F,2) = E = cgetg(l, t_COL);
4195 : }
4196 42749118 : for (i = 1; i < l; i++)
4197 : {
4198 28877210 : gel(P,i) = utoipos(Pf[i]);
4199 28877152 : gel(E,i) = utoipos(Ef[i]);
4200 : }
4201 13871908 : if (pU) *pU = U1 == 1? NULL: mkvec2(utoipos(U1), utoipos(U2));
4202 13871908 : return F;
4203 : }
4204 128350 : if (pU) *pU = NULL;
4205 128350 : M = cgetg(3,t_MAT);
4206 128365 : if (sn < 0) STORE(&nb, utoineg(1), 1);
4207 128365 : if (is_pm1(n)) return aux_end(M,NULL,nb);
4208 :
4209 128365 : n = N = gclone(n); setabssign(n);
4210 : /* trial division bound; look for primes <= lim; nb is the number of
4211 : * distinct prime factors so far; if nb0 >= 0, it records the value of nb
4212 : * for which we made a successful compositeness test: if later nb = nb0,
4213 : * we know that n is composite */
4214 128365 : lim = 1;
4215 128365 : if (!all || all > 2)
4216 : { /* trial divide p < all if all != 0, else p <= tridiv_bound() */
4217 : ulong maxp, p;
4218 : pari_sp av2;
4219 128351 : i = vali(n);
4220 128351 : if (i)
4221 : {
4222 76854 : STOREu(&nb, 2, i);
4223 76854 : av = avma; affii(shifti(n,-i), n); set_avma(av);
4224 : }
4225 128351 : if (is_pm1(n)) return aux_end(M,n,nb);
4226 128235 : lim = all? all-1: tridiv_bound(n);
4227 128235 : av = avma;
4228 128235 : if (lim >= 128)
4229 : { /* fast trial division */
4230 126772 : GEN Q = Z_oddprimedivisors_fast(n, lim);
4231 126772 : av2 = avma;
4232 126772 : if (Q)
4233 : {
4234 123977 : long l = lg(Q);
4235 724350 : for (i = 1; i < l; i++)
4236 : {
4237 601361 : pari_sp av3 = avma;
4238 601361 : ulong p = uel(Q, i);
4239 : long k;
4240 601361 : if (all && p >= all) break;
4241 600373 : k = Z_lvalrem(n, p, &n); /* > 0 */
4242 600373 : affii(n, N); n = N; set_avma(av3);
4243 600373 : STOREu(&nb, p, k);
4244 : }
4245 123977 : if (is_pm1(n))
4246 : {
4247 88502 : stackdummy(av, av2);
4248 88502 : return aux_end(M,n,nb);
4249 : }
4250 : }
4251 38270 : maxp = GP_DATA->factorlimit;
4252 : }
4253 : else
4254 : { /* naive trial division */
4255 1463 : maxp = maxprime();
4256 1463 : u_forprime_init(&T, 3, minuu(lim, maxp)); av2 = avma;
4257 : /* first pass: known to fit in private prime table */
4258 31133 : while ((p = u_forprime_next_fast(&T)))
4259 : {
4260 29991 : pari_sp av3 = avma;
4261 : int stop;
4262 29991 : long k = Z_lvalrem_stop(&n, p, &stop);
4263 29991 : if (k)
4264 : {
4265 4295 : affii(n, N); n = N; set_avma(av3);
4266 4295 : STOREu(&nb, p, k);
4267 : }
4268 : /* prodeuler(p=2,16381,1-1/p) ~ 0.0578; if probability of being prime
4269 : * knowing P^-(n) > 16381 is at least 10%, try BPSW */
4270 29991 : if (!stop && p == 16381)
4271 : {
4272 0 : if (bit_accuracy_mul(lgefint(n), 0.0578 * M_LN2) < 10)
4273 0 : { nb0 = nb; stop = ifac_isprime(n); }
4274 : }
4275 29991 : if (stop)
4276 : {
4277 321 : if (!is_pm1(n)) STOREi(&nb, n, 1);
4278 321 : stackdummy(av, av2);
4279 321 : return aux_end(M,n,nb);
4280 : }
4281 : }
4282 : }
4283 39412 : stackdummy(av, av2);
4284 39412 : if (lim > maxp)
4285 : { /* second pass usually empty, outside fast trial division range */
4286 1 : av = avma; u_forprime_init(&T, maxp+1, lim); av2 = avma;
4287 882811 : while ((p = u_forprime_next(&T)))
4288 : {
4289 882811 : pari_sp av3 = avma;
4290 : int stop;
4291 882811 : long k = Z_lvalrem_stop(&n, p, &stop);
4292 882811 : if (k)
4293 : {
4294 1 : affii(n, N); n = N; set_avma(av3);
4295 1 : STOREu(&nb, p, k);
4296 : }
4297 882811 : if (stop)
4298 : {
4299 1 : if (!is_pm1(n)) STOREi(&nb, n, 1);
4300 1 : stackdummy(av, av2);
4301 1 : return aux_end(M,n,nb);
4302 : }
4303 : }
4304 0 : stackdummy(av, av2);
4305 : }
4306 : }
4307 39425 : if (special_primes(n, lim, &nb, primetab)) return aux_end(M,n, nb);
4308 : /* if nb > nb0 (includes nb0 = -1) we don't know that n is composite */
4309 39292 : if (all)
4310 : { /* smallfact: look for easy pure powers then stop. Cf Z_isanypower */
4311 : GEN x;
4312 31290 : long k, e = expu(lim);
4313 31290 : av = avma;
4314 29976 : k = e >= 10? Z_isanypower_nosmalldiv(n, e, &x)
4315 31290 : : Z_isanypower(n, &x);
4316 31290 : if (k > 1) { affii(x, n); nb0 = -1; } else if (k < 1) k = 1;
4317 31290 : if (pU)
4318 : {
4319 : GEN F;
4320 13156 : if (abscmpiu(n, lim) <= 0
4321 13156 : || cmpii(n, sqru(lim)) <= 0
4322 8833 : || ((e >= 14) &&
4323 7880 : (nb>nb0 && bit_accuracy(lgefint(n))<2048 && ifac_isprime(n))))
4324 13156 : { set_avma(av); STOREi(&nb, n, k); return aux_end(M,n, nb); }
4325 5946 : set_avma(av); F = aux_end(M, NULL, nb); /* don't destroy n */
4326 5946 : *pU = mkvec2(icopy(n), utoipos(k)); /* composite cofactor */
4327 5946 : gunclone(n); return F;
4328 : }
4329 18134 : set_avma(av); STOREi(&nb, n, k);
4330 18134 : if (DEBUGLEVEL >= 2) {
4331 0 : pari_warn(warner,
4332 0 : "IFAC: untested %ld-bit integer declared prime", expi(n)+1);
4333 0 : if (expi(n) <= 256) err_printf("\t%Ps\n", n);
4334 : }
4335 : }
4336 8002 : else if (nb > nb0 && ifac_isprime(n)) STOREi(&nb, n, 1);
4337 2568 : else nb += ifac_decomp(n, hint);
4338 26136 : return aux_end(M,n, nb);
4339 : }
4340 :
4341 : static GEN
4342 9541902 : ifactor(GEN n, ulong all, long hint)
4343 : {
4344 9541902 : long s = signe(n);
4345 9541902 : if (!s) retmkmat2(mkcol(gen_0), mkcol(gen_1));
4346 9541853 : return ifactor_sign(n, all, hint, s, NULL);
4347 : }
4348 :
4349 : int
4350 6651 : ifac_next(GEN *part, GEN *p, long *e)
4351 : {
4352 6651 : GEN here = ifac_main(part);
4353 6651 : if (here == gen_0) { *p = NULL; *e = 1; return 0; }
4354 6649 : if (!here) { *p = NULL; *e = 0; return 0; }
4355 4853 : *p = VALUE(here);
4356 4853 : *e = EXPON(here)[2];
4357 4853 : ifac_delete(here); return 1;
4358 : }
4359 :
4360 : /* see before ifac_crack for current semantics of 'hint' (factorint's 'flag') */
4361 : GEN
4362 10290 : factorint(GEN n, long flag)
4363 : {
4364 : GEN F;
4365 10290 : if ((F = check_arith_all(n,"factorint"))) return gcopy(F);
4366 10276 : return ifactor(n,0,flag);
4367 : }
4368 :
4369 : GEN
4370 49797 : Z_factor_limit(GEN n, ulong all)
4371 : {
4372 49797 : if (!all) all = GP_DATA->factorlimit + 1;
4373 49797 : return ifactor(n, all, decomp_default_hint);
4374 : }
4375 : GEN
4376 886941 : absZ_factor_limit_strict(GEN n, ulong all, GEN *pU)
4377 : {
4378 : GEN F, U;
4379 886941 : if (!signe(n))
4380 : {
4381 0 : if (pU) *pU = NULL;
4382 0 : retmkmat2(mkcol(gen_0), mkcol(gen_1));
4383 : }
4384 886941 : if (!all) all = GP_DATA->factorlimit + 1;
4385 886941 : F = ifactor_sign(n, all, decomp_default_hint, 1, &U);
4386 886972 : if (pU) *pU = U;
4387 886972 : return F;
4388 : }
4389 : GEN
4390 290811 : absZ_factor_limit(GEN n, ulong all)
4391 : {
4392 290811 : if (!signe(n)) retmkmat2(mkcol(gen_0), mkcol(gen_1));
4393 290811 : if (!all) all = GP_DATA->factorlimit + 1;
4394 290811 : return ifactor_sign(n, all, decomp_default_hint, 1, NULL);
4395 : }
4396 : GEN
4397 9481780 : Z_factor(GEN n)
4398 9481780 : { return ifactor(n,0,decomp_default_hint); }
4399 : GEN
4400 3277039 : absZ_factor(GEN n)
4401 : {
4402 3277039 : if (!signe(n)) retmkmat2(mkcol(gen_0), mkcol(gen_1));
4403 3277025 : return ifactor_sign(n, 0, decomp_default_hint, 1, NULL);
4404 : }
4405 : /* Factor until the unfactored part is smaller than limit. Return the
4406 : * factored part. Hence factorback(output) may be smaller than n */
4407 : GEN
4408 3045 : Z_factor_until(GEN n, GEN limit)
4409 : {
4410 3045 : pari_sp av = avma;
4411 3045 : long s = signe(n), eq;
4412 : GEN q, F, U;
4413 :
4414 3045 : if (!s) retmkmat2(mkcol(gen_0), mkcol(gen_1));
4415 3045 : F = ifactor_sign(n, tridiv_bound(n), decomp_default_hint, s, &U);
4416 3045 : if (!U) return F;
4417 1155 : q = gel(U,1); /* composite, q^eq = unfactored part */
4418 1155 : eq = itou(gel(U,2));
4419 1155 : if (cmpii(eq == 1? q: powiu(q, eq), limit) > 0)
4420 : { /* factor further */
4421 1022 : long l2 = expi(q)+1;
4422 : GEN P2, E2, F2, part;
4423 1022 : if (eq > 1) limit = sqrtnint(limit, eq);
4424 1022 : P2 = coltrunc_init(l2);
4425 1022 : E2 = coltrunc_init(l2); F2 = mkmat2(P2,E2);
4426 1022 : part = ifac_start(icopy(q), 0); /* ifac_next would destroy q */
4427 : for(;;)
4428 70 : {
4429 : long e;
4430 : GEN p;
4431 1092 : if (!ifac_next(&part,&p,&e)) break;
4432 1092 : vectrunc_append(P2, p);
4433 1092 : vectrunc_append(E2, utoipos(e * eq));
4434 1092 : q = diviiexact(q, powiu(p, e));
4435 1092 : if (cmpii(q, limit) <= 0) break;
4436 : }
4437 1022 : F2 = sort_factor(F2, (void*)&abscmpii, cmp_nodata);
4438 1022 : F = merge_factor(F, F2, (void*)&abscmpii, cmp_nodata);
4439 : }
4440 1155 : return gerepilecopy(av, F);
4441 : }
4442 :
4443 : static void
4444 96569823 : matsmalltrunc_append(GEN m, ulong p, ulong e)
4445 : {
4446 96569823 : GEN P = gel(m,1), E = gel(m,2);
4447 96569823 : long l = lg(P);
4448 96569823 : P[l] = p; lg_increase(P);
4449 96573968 : E[l] = e; lg_increase(E);
4450 96554374 : }
4451 : static GEN
4452 37767785 : matsmalltrunc_init(long l)
4453 : {
4454 37767785 : GEN P = vecsmalltrunc_init(l);
4455 37730662 : GEN E = vecsmalltrunc_init(l); return mkvec2(P,E);
4456 : }
4457 :
4458 : /* return optimal N s.t. omega(b) <= N for all b <= x */
4459 : long
4460 71606 : maxomegau(ulong x)
4461 : { /* P=primes(15); for(i=1,15, print([i, vecprod(P[1..i])])) */
4462 71606 : if (x < 30030UL)/* rare trivial cases */
4463 : {
4464 37455 : if (x < 2UL) return 0;
4465 19199 : if (x < 6UL) return 1;
4466 13599 : if (x < 30UL) return 2;
4467 12892 : if (x < 210UL) return 3;
4468 12633 : if (x < 2310UL) return 4;
4469 11604 : return 5;
4470 : }
4471 34151 : if (x < 510510UL) return 6; /* most frequent case */
4472 18753 : if (x < 9699690UL) return 7;
4473 7 : if (x < 223092870UL) return 8;
4474 : #ifdef LONG_IS_64BIT
4475 6 : if (x < 6469693230UL) return 9;
4476 0 : if (x < 200560490130UL) return 10;
4477 0 : if (x < 7420738134810UL) return 11;
4478 0 : if (x < 304250263527210UL) return 12;
4479 0 : if (x < 13082761331670030UL) return 13;
4480 0 : if (x < 614889782588491410UL) return 14;
4481 0 : return 15;
4482 : #else
4483 1 : return 9;
4484 : #endif
4485 : }
4486 : /* return optimal N s.t. omega(b) <= N for all odd b <= x */
4487 : long
4488 2229 : maxomegaoddu(ulong x)
4489 : { /* P=primes(15+1); for(i=1,15, print([i, vecprod(P[2..i+1])])) */
4490 2229 : if (x < 255255UL)/* rare trivial cases */
4491 : {
4492 1355 : if (x < 3UL) return 0;
4493 1355 : if (x < 15UL) return 1;
4494 1355 : if (x < 105UL) return 2;
4495 1355 : if (x < 1155UL) return 3;
4496 1327 : if (x < 15015UL) return 4;
4497 1327 : return 5;
4498 : }
4499 874 : if (x < 4849845UL) return 6; /* most frequent case */
4500 0 : if (x < 111546435UL) return 7;
4501 0 : if (x < 3234846615UL) return 8;
4502 : #ifdef LONG_IS_64BIT
4503 0 : if (x < 100280245065UL) return 9;
4504 0 : if (x < 3710369067405UL) return 10;
4505 0 : if (x < 152125131763605UL) return 11;
4506 0 : if (x < 6541380665835015UL) return 12;
4507 0 : if (x < 307444891294245705UL) return 13;
4508 0 : if (x < 16294579238595022365UL) return 14;
4509 0 : return 15;
4510 : #else
4511 0 : return 9;
4512 : #endif
4513 : }
4514 :
4515 : /* If a <= c <= b , factoru(c) = L[c-a+1] */
4516 : GEN
4517 31925 : vecfactoru_i(ulong a, ulong b)
4518 : {
4519 31925 : ulong k, p, n = b-a+1, N = maxomegau(b) + 1;
4520 31925 : GEN v = const_vecsmall(n, 1);
4521 31925 : GEN L = cgetg(n+1, t_VEC);
4522 : forprime_t T;
4523 20783074 : for (k = 1; k <= n; k++) gel(L,k) = matsmalltrunc_init(N);
4524 31925 : u_forprime_init(&T, 2, usqrt(b));
4525 883441 : while ((p = u_forprime_next(&T)))
4526 : { /* p <= sqrt(b) */
4527 853170 : ulong pk = p, K = ulogint(b, p);
4528 2967666 : for (k = 1; k <= K; k++)
4529 : {
4530 2116150 : ulong j, t = a / pk, ap = t * pk;
4531 2116150 : if (ap < a) { ap += pk; t++; }
4532 : /* t = (j+a-1) \ pk */
4533 2116150 : t %= p;
4534 60953760 : for (j = ap-a+1; j <= n; j += pk)
4535 : {
4536 58839268 : if (t) { v[j] *= pk; matsmalltrunc_append(gel(L,j), p,k); }
4537 58837610 : if (++t == p) t = 0;
4538 : }
4539 2114492 : pk *= p;
4540 : }
4541 : }
4542 : /* complete factorisation of non-sqrt(b)-smooth numbers */
4543 21150542 : for (k = 1, N = a; k <= n; k++, N++)
4544 21118326 : if (uel(v,k) != N) matsmalltrunc_append(gel(L,k), N/uel(v,k),1UL);
4545 32216 : return L;
4546 : }
4547 : GEN
4548 0 : vecfactoru(ulong a, ulong b)
4549 : {
4550 0 : pari_sp av = avma;
4551 0 : return gerepilecopy(av, vecfactoru_i(a,b));
4552 : }
4553 :
4554 : /* Assume a and b odd, return L s.t. L[k] = factoru(a + 2*(k-1))
4555 : * If a <= c <= b odd, factoru(c) = L[(c-a)>>1 + 1] */
4556 : GEN
4557 2229 : vecfactoroddu_i(ulong a, ulong b)
4558 : {
4559 2229 : ulong k, p, n = ((b-a)>>1) + 1, N = maxomegaoddu(b) + 1;
4560 2229 : GEN v = const_vecsmall(n, 1);
4561 2229 : GEN L = cgetg(n+1, t_VEC);
4562 : forprime_t T;
4563 :
4564 17321271 : for (k = 1; k <= n; k++) gel(L,k) = matsmalltrunc_init(N);
4565 2229 : u_forprime_init(&T, 3, usqrt(b));
4566 185277 : while ((p = u_forprime_next(&T)))
4567 : { /* p <= sqrt(b) */
4568 183064 : ulong pk = p, K = ulogint(b, p);
4569 623601 : for (k = 1; k <= K; k++)
4570 : {
4571 440553 : ulong j, t = (a / pk) | 1UL, ap = t * pk;
4572 : /* t and ap are odd, ap multiple of pk = p^k */
4573 440553 : if (ap < a) { ap += pk<<1; t+=2; }
4574 : /* c=t*p^k by steps of 2*p^k; factorization of c*=p^k if (t,p)=1 */
4575 440553 : t %= p;
4576 32504310 : for (j = ((ap-a)>>1)+1; j <= n; j += pk)
4577 : {
4578 32063881 : if (t) { v[j] *= pk; matsmalltrunc_append(gel(L,j), p,k); }
4579 32063757 : t += 2; if (t >= p) t -= p;
4580 : }
4581 440429 : pk *= p;
4582 : }
4583 : }
4584 : /* complete factorisation of non-sqrt(b)-smooth numbers */
4585 17688492 : for (k = 1, N = a; k <= n; k++, N+=2)
4586 17686437 : if (uel(v,k) != N) matsmalltrunc_append(gel(L,k), N/uel(v,k),1UL);
4587 2055 : return L;
4588 : }
4589 : GEN
4590 0 : vecfactoroddu(ulong a, ulong b)
4591 : {
4592 0 : pari_sp av = avma;
4593 0 : return gerepilecopy(av, vecfactoroddu_i(a,b));
4594 : }
4595 :
4596 : /* If 0 <= a <= c <= b; L[c-a+1] = factoru(c)[,1] if c squarefree, else NULL */
4597 : GEN
4598 7014 : vecfactorsquarefreeu(ulong a, ulong b)
4599 : {
4600 7014 : ulong k, p, n = b-a+1, N = maxomegau(b) + 1;
4601 7014 : GEN v = const_vecsmall(n, 1);
4602 7014 : GEN L = cgetg(n+1, t_VEC);
4603 : forprime_t T;
4604 14007238 : for (k = 1; k <= n; k++) gel(L,k) = vecsmalltrunc_init(N);
4605 7014 : u_forprime_init(&T, 2, usqrt(b));
4606 838334 : while ((p = u_forprime_next(&T)))
4607 : { /* p <= sqrt(b), kill nonsquarefree */
4608 831320 : ulong j, pk = p*p, t = a / pk, ap = t * pk;
4609 831320 : if (ap < a) ap += pk;
4610 7160090 : for (j = ap-a+1; j <= n; j += pk) gel(L,j) = NULL;
4611 :
4612 831320 : t = a / p; ap = t * p;
4613 831320 : if (ap < a) { ap += p; t++; }
4614 30551556 : for (j = ap-a+1; j <= n; j += p, t++)
4615 29720236 : if (gel(L,j)) { v[j] *= p; vecsmalltrunc_append(gel(L,j), p); }
4616 : }
4617 : /* complete factorisation of non-sqrt(b)-smooth numbers */
4618 14007238 : for (k = 1, N = a; k <= n; k++, N++)
4619 14000224 : if (gel(L,k) && uel(v,k) != N) vecsmalltrunc_append(gel(L,k), N/uel(v,k));
4620 7014 : return L;
4621 : }
4622 : /* If 0 <= a <= c <= b; L[c-a+1] = factoru(c)[,1] if c squarefree and coprime
4623 : * to all the primes in sorted zv P, else NULL */
4624 : GEN
4625 32590 : vecfactorsquarefreeu_coprime(ulong a, ulong b, GEN P)
4626 : {
4627 32590 : ulong k, p, n = b-a+1, sqb = usqrt(b), N = maxomegau(b) + 1;
4628 32590 : GEN v = const_vecsmall(n, 1);
4629 32591 : GEN L = cgetg(n+1, t_VEC);
4630 : forprime_t T;
4631 90637851 : for (k = 1; k <= n; k++) gel(L,k) = vecsmalltrunc_init(N);
4632 32592 : u_forprime_init(&T, 2, sqb);
4633 3680347 : while ((p = u_forprime_next(&T)))
4634 : { /* p <= sqrt(b), kill nonsquarefree */
4635 3648305 : ulong j, t, ap, bad = zv_search(P, p), pk = bad ? p: p * p;
4636 3648553 : t = a / pk; ap = t * pk; if (ap < a) ap += pk;
4637 80847095 : for (j = ap-a+1; j <= n; j += pk) gel(L,j) = NULL;
4638 3648553 : if (bad) continue;
4639 :
4640 3586400 : t = a / p; ap = t * p;
4641 3586400 : if (ap < a) { ap += p; t++; }
4642 116461758 : for (j = ap-a+1; j <= n; j += p, t++)
4643 112876156 : if (gel(L,j)) { v[j] *= p; vecsmalltrunc_append(gel(L,j), p); }
4644 : }
4645 32592 : if (uel(P,lg(P)-1) <= sqb) P = NULL;
4646 : /* complete factorisation of non-sqrt(b)-smooth numbers */
4647 90824890 : for (k = 1, N = a; k <= n; k++, N++)
4648 90792409 : if (gel(L,k) && uel(v,k) != N)
4649 : {
4650 25662448 : ulong q = N / uel(v,k);
4651 25662448 : if (!P || !zv_search(P, q)) vecsmalltrunc_append(gel(L,k), q);
4652 : }
4653 32481 : return L;
4654 : }
4655 :
4656 : GEN
4657 49 : vecsquarefreeu(ulong a, ulong b)
4658 : {
4659 49 : ulong j, k, p, n = b-a+1;
4660 49 : GEN L = const_vecsmall(n, 1);
4661 : forprime_t T;
4662 49 : u_forprime_init(&T, 2, usqrt(b));
4663 462 : while ((p = u_forprime_next(&T)))
4664 : { /* p <= sqrt(b), kill nonsquarefree */
4665 413 : ulong pk = p*p, t = a / pk, ap = t * pk;
4666 413 : if (ap < a) { ap += pk; t++; }
4667 : /* t = (j+a-1) \ pk */
4668 21777 : for (j = ap-a+1; j <= n; j += pk, t++) L[j] = 0;
4669 : }
4670 48258 : for (k = j = 1; k <= n; k++)
4671 48209 : if (L[k]) L[j++] = a+k-1;
4672 49 : setlg(L,j); return L;
4673 : }
|
