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