Annotation of OpenXM_contrib/pari/src/basemath/ifactor1.c, Revision 1.1.1.1
1.1 maekawa 1: /********************************************************************/
2: /** **/
3: /** INTEGER FACTORIZATION **/
4: /** **/
5: /********************************************************************/
6: /* $Id: ifactor1.c,v 1.1.1.1 1999/09/16 13:47:33 karim Exp $ */
7: #include "pari.h"
8:
9: /*********************************************************************/
10: /** **/
11: /** PSEUDO PRIMALITY **/
12: /** **/
13: /*********************************************************************/
14: static GEN sqrt1, sqrt2, t1, t;
15: static long r1;
16:
17: /* The following two internal routines don't restore avma -- the caller
18: must do so at the end. */
19: static GEN
20: init_miller(GEN n)
21: {
22: if (signe(n) < 0) n = absi(n);
23: t=addsi(-1,n); r1=vali(t); t1 = shifti(t,-r1);
24: sqrt1=cgeti(lg(t)); sqrt1[1]=evalsigne(0)|evallgefint(2);
25: sqrt2=cgeti(lg(t)); sqrt2[1]=evalsigne(0)|evallgefint(2);
26: return n;
27: }
28:
29: /* is n strong pseudo-prime for base a ? `End matching' (check for square
30: * roots of -1) added by GN */
31: /* TODO: If ends do mismatch, then we have factored n, and this information
32: should somehow be made available to the factoring machinery. --GN */
33: static int
34: bad_for_base(GEN n, GEN a)
35: {
36: long r, av=avma, lim=stack_lim(av,1);
37: GEN c2, c = powmodulo(a,t1,n);
38:
39: if (!is_pm1(c) && !egalii(t,c)) /* go fishing for -1, not for 1 */
40: {
41: for (r=r1-1; r; r--) /* (this saves one squaring/reduction) */
42: {
43: c2=c; c=resii(sqri(c),n);
44: if (egalii(t,c)) break;
45: if (low_stack(lim, stack_lim(av,1)))
46: {
47: GEN *gsav[2]; gsav[0]=&c; gsav[1]=&c2;
48: if(DEBUGMEM>1) err(warnmem,"miller(rabin)");
49: gerepilemany(av, gsav, 2);
50: }
51: }
52: if (!r) return 1;
53: /* sqrt(-1) seen, compare or remember */
54: if (signe(sqrt1)) /* we saw one earlier: compare */
55: {
56: /* check if too many sqrt(-1)s mod n */
57: if (!egalii(c2,sqrt1) && !egalii(c2,sqrt2)) return 1;
58: }
59: else { affii(c2,sqrt1); affii(subii(n,c2),sqrt2); } /* remember */
60: }
61: return 0;
62: }
63:
64: /* Miller-Rabin test for k random bases */
65: long
66: millerrabin(GEN n, long k)
67: {
68: long r,i,av2, av = avma;
69:
70: if (!signe(n)) return 0;
71: /* If |n| <= 3, check if n = +- 1 */
72: if (lgefint(n)==3 && (ulong)(n[2])<=3) return (n[2] != 1);
73:
74: if (!mod2(n)) return 0;
75: n = init_miller(n); av2=avma;
76: for (i=1; i<=k; i++)
77: {
78: do r = smodsi(mymyrand(),n); while (!r);
79: if (DEBUGLEVEL > 4)
80: fprintferr("Miller-Rabin: testing base %ld\n",
81: r);
82: if (bad_for_base(n, stoi(r))) { avma=av; return 0; }
83: avma=av2;
84: }
85: avma=av; return 1;
86: }
87:
88: /* As above for k bases taken in pr (i.e not random).
89: * We must have |n|>2 and 1<=k<=11 (not checked) or k in {16,17} to select
90: * some special sets of bases.
91: *
92: * By computations of Gerhard Jaeschke, `On strong pseudoprimes to several
93: * bases', Math.Comp. 61 (1993), 915--926 (see also Chris Caldwell's Prime
94: * Number Pages at http://www.utm.edu/research/primes/prove2.html), we have:
95: *
96: * k == 4 (bases 2,3,5,7) correctly detects all composites
97: * n < 118 670 087 467 == 172243 * 688969 with the single exception of
98: * n == 3 215 031 751 == 151 * 751 * 28351,
99: *
100: * k == 5 (bases 2,3,5,7,11) correctly detects all composites
101: * n < 2 152 302 898 747 == 6763 * 10627 * 29947,
102: *
103: * k == 6 (bases 2,3,...,13) correctly detects all composites
104: * n < 3 474 749 660 383 == 1303 * 16927 * 157543,
105: *
106: * k == 7 (bases 2,3,...,17) correctly detects all composites
107: * n < 341 550 071 728 321 == 10670053 * 32010157,
108: * and even this limiting value is caught by an end mismatch between bases
109: * 2 and 5 (or 5 and 17).
110: *
111: * Moreover, the four bases chosen at
112: *
113: * k == 16 (2,13,23,1662803) will correctly detect all composites up
114: * to at least 10^12, and the combination at
115: *
116: * k == 17 (31,73) detects most odd composites without prime factors > 100
117: * in the range n < 2^36 (with less than 250 exceptions, indeed with fewer
118: * than 1400 exceptions up to 2^42). --GN
119: * (DATA TO BE COMPLETED)
120: */
121: int /* no longer static -- needed in mpqs.c */
122: miller(GEN n, long k)
123: {
124: long r,i,av2, av = avma;
125: static long pr[] =
126: { 0, 2,3,5,7,11,13,17,19,23,29, 31,73, 2,13,23,1662803UL, };
127: long *p;
128:
129: if (!mod2(n)) return 0;
130: if (k==16)
131: { /* use smaller (faster) bases if possible */
132: if (lgefint(n)==3 && (ulong)(n[2]) < 3215031751UL) p = pr; /* 2,3,5,7 */
133: else p = pr+13; /* 2,13,23,1662803 */
134: k=4;
135: }
136: else if (k==17)
137: {
138: if (lgefint(n)==3 && (ulong)(n[2]) < 1373653) p = pr; /* 2,3 */
139: else p = pr+11; /* 31,73 */
140: k=2;
141: }
142: else p = pr; /* 2,3,5,... */
143: n = init_miller(n); av2=avma;
144: for (i=1; i<=k; i++)
145: {
146: r = smodsi(p[i],n); if (!r) break;
147: if (bad_for_base(n, stoi(r))) { avma = av; return 0; }
148: avma=av2;
149: }
150: avma=av; return 1;
151: }
152:
153: /***********************************************************************/
154: /** **/
155: /** PRIMES IN SUCCESSION **/
156: /** (abstracted by GN 1998Aug21 mainly for use in ellfacteur() below) **/
157: /** **/
158: /***********************************************************************/
159:
160: /* map from prime residue classes mod 210 to their numbers in {0...47}.
161: Subscripts into this array take the form ((k-1)%210)/2, ranging from
162: 0 to 104. Unused entries are 128 */
163: #define NPRC 128
164:
165: static
166: unsigned char prc210_no[] =
167: {
168: 0, NPRC, NPRC, NPRC, NPRC, 1, 2, NPRC, 3, 4, NPRC, /* 21 */
169: 5, NPRC, NPRC, 6, 7, NPRC, NPRC, 8, NPRC, 9, /* 41 */
170: 10, NPRC, 11, NPRC, NPRC, 12, NPRC, NPRC, 13, 14, NPRC, /* 63 */
171: NPRC, 15, NPRC, 16, 17, NPRC, NPRC, 18, NPRC, 19, /* 83 */
172: NPRC, NPRC, 20, NPRC, NPRC, NPRC, 21, NPRC, 22, 23, NPRC, /* 105 */
173: 24, 25, NPRC, 26, NPRC, NPRC, NPRC, 27, NPRC, NPRC, /* 125 */
174: 28, NPRC, 29, NPRC, NPRC, 30, 31, NPRC, 32, NPRC, NPRC, /* 147 */
175: 33, 34, NPRC, NPRC, 35, NPRC, NPRC, 36, NPRC, 37, /* 167 */
176: 38, NPRC, 39, NPRC, NPRC, 40, 41, NPRC, NPRC, 42, NPRC, /* 189 */
177: 43, 44, NPRC, 45, 46, NPRC, NPRC, NPRC, NPRC, 47, /* 209 */
178: };
179:
180: /* map from prime residue classes mod 210 (by number) to their smallest
181: positive representatives */
182: static
183: unsigned char prc210_rp[] =
184: {
185: 1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79,
186: 83, 89, 97, 101, 103, 107, 109, 113, 121, 127, 131, 137, 139, 143, 149,
187: 151, 157, 163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199, 209,
188: };
189:
190: /* first differences of the preceding */
191: static
192: unsigned char prc210_d1[] =
193: {
194: 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6,
195: 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6,
196: 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2,
197: };
198:
199: GEN
200: nextprime(GEN n)
201: {
202: long rc,rc0,rcd,rcn,av1,av2, av = avma;
203:
204: if (typ(n) != t_INT) n=gceil(n); /* accept arguments in R --GN */
205: if (typ(n) != t_INT) err(arither1);
206: if (signe(n) <= 0) { avma=av; return gdeux; }
207: if (lgefint(n) <= 3)
208: { /* check if n <= 7 */
209: ulong k = n[2];
210: if (k <= 2) { avma=av; return gdeux; }
211: if (k == 3) { avma = av; return stoi(3); }
212: if (k <= 5) { avma = av; return stoi(5); }
213: if (k <= 7) { avma = av; return stoi(7); }
214: }
215: /* here n > 7 */
216: if (!(mod2(n))) n = addsi(1,n);
217: rc = rc0 = smodis(n, 210);
218: rcn = (long)(prc210_no[rc0>>1]);
219: /* find next prime residue class mod 210 */
220: while (rcn == NPRC)
221: {
222: rc += 2; /* cannot wrap since 209 is coprime */
223: rcn = (long)(prc210_no[rc>>1]);
224: }
225: if (rc > rc0) n = addsi(rc - rc0, n);
226: /* now find an actual prime */
227: av2 = av1 = avma;
228: for(;;)
229: {
230: if (miller(n,10)) break;
231: av1 = avma;
232: rcd = prc210_d1[rcn];
233: if (++rcn > 47) rcn = 0;
234: n = addsi(rcd,n);
235: }
236: if (av1!=av2) return gerepile(av,av1,n);
237: return (av1==av)? icopy(n): n;
238: }
239:
240: GEN
241: precprime(GEN n)
242: {
243: long rc,rc0,rcd,rcn,av1,av2, av = avma;
244:
245: if (typ(n) != t_INT) n=gfloor(n); /* accept arguments in R --GN */
246: if (typ(n) != t_INT) err(arither1);
247: if (signe(n)<=0) { avma=av; return gzero; }
248: if (lgefint(n) <= 3)
249: { /* check if n <= 10 */
250: ulong k = n[2];
251: if (k <= 1) { avma=av; return gzero; }
252: if (k == 2) { avma=av; return gdeux; }
253: if (k <= 4) { avma=av; return stoi(3); }
254: if (k <= 6) { avma=av; return stoi(5); }
255: if (k <= 10) { avma=av; return stoi(7); }
256: }
257: /* here n >= 11 */
258: if (!(mod2(n))) n = addsi(-1,n);
259: rc = rc0 = smodis(n, 210);
260: rcn = (long)(prc210_no[rc0>>1]);
261: /* find last prime residue class mod 210 */
262: while (rcn == NPRC)
263: {
264: rc -= 2; /* cannot wrap since 1 is coprime */
265: rcn = (long)(prc210_no[rc>>1]);
266: }
267: if (rc < rc0) n = addsi(rc - rc0, n);
268: /* now find an actual prime */
269: av2 = av1 = avma;
270: for(;;)
271: {
272: if (miller(n,10)) break;
273: av1 = avma;
274: if (rcn == 0)
275: { rcd = 2; rcn = 47; }
276: else
277: rcd = prc210_d1[--rcn];
278: n = addsi(-rcd,n);
279: }
280: if (av1!=av2) return gerepile(av,av1,n);
281: return (av1==av)? icopy(n): n;
282: }
283:
284: /* find next single-word prime strictly larger than p. If **d is non-NULL,
285: this will be p + *(*d)++, using the diffptr table. Otherwise imitate
286: nextprime(). Apart from *d, caller must supply a long variable to which
287: rcn points, initialized either to NPRC or to the correct residue class
288: number for the current p; we'll use this to track the current prime
289: residue class mod 210 once we're out of range of the diffptr table, and
290: we'll update it before that if it isn't NPRC. *q is incremented when-
291: ever q!=NULL and we wrap from 209 mod 210 to 1 mod 210; this make sense
292: only when *rcn already held the correct value. Caller must also supply
293: the second argument for miller(). --GN1998Aug22 */
294: ulong
295: snextpr(ulong p, byteptr *d, long *rcn, long *q, long k)
296: {
297: static ulong pp[] =
298: { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0 };
299: static ulong *pp2 = pp + 2;
300: static GEN gp = (GEN)pp;
301: long d1 = **d, rcn0;
302:
303: if (d1)
304: {
305: if (*rcn != NPRC)
306: {
307: rcn0 = *rcn;
308: while (d1 > 0)
309: {
310: d1 -= prc210_d1[*rcn];
311: if (++*rcn > 47) { *rcn = 0; if (q) (*q)++; }
312: }
313: if (d1 < 0)
314: {
315: fprintferr("snextpr: prime %lu wasn\'t %lu mod 210\n",
316: p, (ulong)prc210_rp[rcn0]);
317: err(bugparier, "[caller of] snextpr");
318: }
319: }
320: return p + *(*d)++;
321: }
322: /* we are beyond the diffptr table */
323: if (*rcn == NPRC) /* we need to initialize this now */
324: {
325: *rcn = prc210_no[(p % 210) >> 1];
326: if (*rcn == NPRC)
327: {
328: fprintferr("snextpr: %lu should have been prime but isn\'t\n", p);
329: err(bugparier, "[caller of] snextpr");
330: }
331: }
332: /* look for the next one */
333: *pp2 = p;
334: *pp2 += prc210_d1[*rcn];
335: if (++*rcn > 47) *rcn = 0;
336: while (!miller(gp, k))
337: {
338: *pp2 += prc210_d1[*rcn];
339: if (++*rcn > 47) { *rcn = 0; if (q) (*q)++; }
340: if (*pp2 <= 11) /* wraparound mod 2^BITS_IN_LONG */
341: {
342: fprintferr("snextpr: integer wraparound after prime %lu\n", p);
343: err(bugparier, "[caller of] snextpr");
344: }
345: }
346: return *pp2;
347: }
348:
349:
350: /***********************************************************************/
351: /** **/
352: /** FACTORIZATION (ECM) **/
353: /** Integer factorization using the elliptic curves method (ECM). **/
354: /** ellfacteur() returns a non trivial factor of N, assuming N>0, **/
355: /** is composite, and has no prime divisor below 2^14 or so. **/
356: /** Extensively modified by GN Jul-Aug 1998, with much helpful **/
357: /** advice by Paul Zimmermann. Thanks also to Guillaume Hanrot **/
358: /** and Igor Schein for providing many CPU cycles whilst testing. **/
359: /** **/
360: /***********************************************************************/
361:
362: static GEN N, gl, *XAUX;
363: #define nbcmax 64 /* max number of simultaneous curves */
364: #define bstpmax 1024 /* max number of baby step table entries */
365:
366: /* addition/doubling/multiplication of a point on an `elliptic curve'
367: mod N may result in one of three things: a new bona fide point,
368: a point at infinity (betraying itself by a denominator divisible
369: by N), or a point which is at infinity mod some nontrivial factor
370: of N but finite mod some other factor (betraying itself by a denom-
371: inator which has nontrivial gcd with N, and this is of course what
372: we want). */
373: /* (In the second case, addition/doubling will simply abort, copying one
374: of the summands to the destination array of points unless they coincide.
375: Multiplication will stop at some unpredictable intermediate stage: The
376: destination will contain _some_ multiple of the input point, but not
377: necessarily the desired one, which doesn't matter. As long as we're
378: multiplying (B1 phase) we simply carry on with the next multiplier.
379: During the B2 phase, the only additions are the giant steps, and the
380: worst that can happen here is that we lose one residue class mod 210
381: of prime multipliers on 4 of the curves, so again, we ignore the problem
382: and just carry on.) */
383: /* The idea is: Select a handful of curves mod N and one point P on each of
384: them. Try to compute, for each such point, the multiple [M]P = Q where
385: M is the product of all powers <= B2 of primes <= nextprime(B1), for some
386: suitably chosen B1 and B2. Then check whether multiplying Q by one of the
387: primes < nextprime(B2) would betray a factor. This second stage proceeds
388: by looking separately at the primes in each residue class mod 210, four
389: curves at a time, and stepping additively to ever larger multipliers,
390: by comparing X coordinates of points which we would need to add in order
391: to reach another prime multiplier in the same residue class. `Comparing'
392: means that we accumulate a product of differences of X coordinates, and
393: from time to time take a gcd of this product with N. */
394: /* Montgomery's trick of hiding the cost of computing inverses mod N at a
395: price of three extra multiplications mod N, by working on up to 64 or
396: even 128 points in parallel, is used heavily. --GN */
397:
398: /* *** auxiliary functions for ellfacteur: *** */
399:
400: /* Parallel addition on nbc curves, assigning the result to locations at and
401: following *X3, *Y3. Safe to be called with X3,Y3 equal to X2,Y2 (_not_
402: to X1,Y1). It is also safe to overwrite Y2 with X3. (If Y coords of
403: result not desired, set Y3=NULL.) If nbc1 < nbc, the first summand is
404: assumed to hold only nbc1 distinct points, which are repeated as often
405: as we need them (useful for adding one point on each of a few curves
406: to several other points on the same curves).
407: Return 0 when successful, 1 when we hit a denominator divisible by N,
408: and 2 when gcd(denominator, N) is a nontrivial factor of N, which will
409: be preserved in gl.
410: We use more stack space than the old code did, and thus run a bit of a
411: risk of overflowing it, but it's still bounded by a constant multiple
412: of lgefint(N)*nbc, as it was in the old version --GN1998Jul02,Aug12 */
413: /* (Lessee: Second phase creates 12 items on the stack, per iteration,
414: of which four are twice as long and one is thrice as long as N --
415: makes 18 units per iteration. First phase creates 4 units. Total
416: can be as large as about 4*nbcmax+18*8 units. And elladd2() is just
417: as bad, and elldouble() comes to about 3*nbcmax+29*8 units. A few
418: strategic garbage collections every 8 iterations should help when nbc
419: is large...) --GN1998Aug23 */
420:
421: static int
422: elladd0(long nbc, long nbc1,
423: GEN *X1, GEN *Y1, GEN *X2, GEN *Y2, GEN *X3, GEN *Y3)
424: {
425: GEN lambda;
426: GEN W[2*nbcmax], *A=W+nbc; /* W[0],A[0] never used */
427: long i, av=avma, tetpil;
428: ulong mask = ~0UL;
429:
430: /* actually, this is only ever called with nbc1==nbc or nbc1==4, so: */
431: if (nbc1 == 4) mask = 3;
432: else if (nbc1 < nbc) err(bugparier, "[caller of] elladd0");
433:
434: /* W[0] = gun; */
435: W[1] = /* A[0] =*/ subii(X1[0], X2[0]);
436: for (i=1; i<nbc; i++)
437: {
438: A[i] = subii(X1[i&mask], X2[i]); /* don't waste time reducing mod N here */
439: W[i+1] = modii(mulii(A[i], W[i]), N);
440: }
441: tetpil = avma;
442:
443: /* if gl != N we have a factor */
444: if (!invmod(W[nbc], N, &gl))
445: {
446: if (!egalii(N,gl)) { gl = gerepile(av,tetpil,gl); return 2; }
447: if (X2 != X3)
448: {
449: long k;
450: /* cannot add on one of the curves mod N: make sure X3 contains
451: something useful before letting the caller proceed */
452: for (k = 2*nbc; k--; ) affii(X2[k],X3[k]);
453: }
454: avma = av; return 1;
455: }
456:
457: while (i--) /* nbc times, actually */
458: {
459: lambda = modii(mulii(subii(Y1[i&mask], Y2[i]),
460: i?mulii(gl, W[i]):gl), N);
461: modiiz(subii(sqri(lambda), addii(X2[i], X1[i&mask])), N, X3[i]);
462: if (Y3)
463: modiiz(subii(mulii(lambda, subii(X1[i&mask], X3[i])),
464: Y1[i&mask]),
465: N, Y3[i]);
466: if (!i) break;
467: gl = modii(mulii(gl, A[i]), N);
468: if (!(i&7)) gl = gerepileupto(tetpil, gl);
469: }
470: avma=av; return 0;
471: }
472:
473: /* Shortcut variant, for use in cases where Y coordinates follow their
474: corresponding X coordinates, and the first summand doesn't need to be
475: repeated */
476: static int
477: elladd(long nbc, GEN *X1, GEN *X2, GEN *X3)
478: {
479: return elladd0(nbc, nbc, X1, X1+nbc, X2, X2+nbc, X3, X3+nbc);
480: }
481: /* this could perhaps become a macro --GN */
482:
483: /* The next one is exactly the same except it does twice as many additions
484: (and thus hides even more of the cost of the modular inverse); the net
485: effect is the same as elladd(nbc,X1,X2,X3) followed by elladd(nbc,X4,X5,X6).
486: Safe to have X2==X3 and/or X5==X6, and of course safe to have X1 or X2
487: coincide with X4 or X5, in any order. */
488:
489: static int
490: elladd2(long nbc, GEN *X1, GEN *X2, GEN *X3, GEN *X4, GEN *X5, GEN *X6)
491: {
492: GEN lambda, *Y1 = X1+nbc, *Y2 = X2+nbc, *Y3 = X3+nbc;
493: GEN *Y4 = X4+nbc, *Y5 = X5+nbc, *Y6 = X6+nbc;
494: GEN W[4*nbcmax], *A=W+2*nbc; /* W[0],A[0] never used */
495: long i,j, av=avma, tetpil;
496: /* W[0] = gun; */
497: W[1] = /* A[0] =*/ subii(X1[0], X2[0]);
498: for (i=1; i<nbc; i++)
499: {
500: A[i] = subii(X1[i], X2[i]); /* don't waste time reducing mod N here */
501: W[i+1] = modii(mulii(A[i], W[i]), N);
502: }
503: for (j=0; j<nbc; i++,j++)
504: {
505: A[i] = subii(X4[j], X5[j]);
506: W[i+1] = modii(mulii(A[i], W[i]), N);
507: }
508: tetpil = avma;
509:
510: /* if gl != N we have a factor */
511: if (!invmod(W[2*nbc], N, &gl))
512: {
513: if (!egalii(N,gl)) { gl = gerepile(av,tetpil,gl); return 2; }
514: if (X2 != X3)
515: {
516: long k;
517: /* cannot add on one of the curves mod N: make sure X3 contains
518: something useful before letting the caller proceed */
519: for (k = 2*nbc; k--; ) affii(X2[k],X3[k]);
520: }
521: if (X5 != X6)
522: {
523: long k;
524: /* same for X6 */
525: for (k = 2*nbc; k--; ) affii(X5[k],X6[k]);
526: }
527: avma = av; return 1;
528: }
529:
530: while (j--) /* nbc times, actually */
531: {
532: i--;
533: lambda = modii(mulii(subii(Y4[j], Y5[j]),
534: mulii(gl, W[i])), N);
535: modiiz(subii(sqri(lambda), addii(X5[j], X4[j])), N, X6[j]);
536: modiiz(subii(mulii(lambda, subii(X4[j], X6[j])), Y4[j]), N, Y6[j]);
537: gl = modii(mulii(gl, A[i]), N);
538: if (!(i&7)) gl = gerepileupto(tetpil, gl);
539: }
540: while (i--) /* nbc times */
541: {
542: lambda = modii(mulii(subii(Y1[i], Y2[i]),
543: i?mulii(gl, W[i]):gl), N);
544: modiiz(subii(sqri(lambda), addii(X2[i], X1[i])), N, X3[i]);
545: modiiz(subii(mulii(lambda, subii(X1[i], X3[i])), Y1[i]), N, Y3[i]);
546: if (!i) break;
547: gl = modii(mulii(gl, A[i]), N);
548: if (!(i&7)) gl = gerepileupto(tetpil, gl);
549: }
550: avma=av; return 0;
551: }
552:
553: /* Parallel doubling on nbc curves, assigning the result to locations at
554: and following *X2. Safe to be called with X2 equal to X1. Return
555: value as for elladd() above. If we find a point at infinity mod N,
556: and if X1 != X2, we copy the points at X1 to X2.
557: Use fewer assignments than the old code. Strangely, whereas this gains
558: about 3% on my P133 with elladd(), it makes hardly any difference here
559: with elldouble() --GN */
560: static int
561: elldouble(long nbc, GEN *X1, GEN *X2)
562: {
563: GEN lambda,v, *Y1 = X1+nbc, *Y2 = X2+nbc;
564: GEN W[nbcmax+1]; /* W[0] never used */
565: long i, av=avma, tetpil;
566: /*W[0] = gun;*/ W[1] = Y1[0];
567: for (i=1; i<nbc; i++)
568: W[i+1] = modii(mulii(Y1[i], W[i]), N);
569: tetpil = avma;
570:
571: if (!invmod(W[nbc], N, &gl))
572: {
573: if (!egalii(N,gl)) { gl = gerepile(av,tetpil,gl); return 2; }
574: if (X1 != X2)
575: {
576: long k;
577: /* cannot double on one of the curves mod N: make sure X2 contains
578: something useful before letting the caller proceed */
579: for (k = 2*nbc; k--; ) affii(X1[k],X2[k]);
580: }
581: avma = av; return 1;
582: }
583:
584: while (i--) /* nbc times, actually */
585: {
586: lambda = modii(mulii(addsi(1, mulsi(3, sqri(X1[i]))),
587: i?mulii(gl,W[i]):gl), N);
588: if (signe(lambda)) /* half of zero is still zero */
589: lambda = shifti(mod2(lambda)? addii(lambda, N): lambda, -1);
590: v = modii(subii(sqri(lambda), shifti(X1[i],1)), N);
591: if (i) gl = modii(mulii(gl, Y1[i]), N);
592: modiiz(subii(mulii(lambda, subii(X1[i], v)), Y1[i]), N, Y2[i]);
593: affii(v, X2[i]);
594: if (!(i&7) && i) gl = gerepileupto(tetpil, gl);
595: }
596: avma = av; return 0;
597: }
598:
599: /* Parallel multiplication by an odd prime k on nbc curves, storing the
600: result to locations at and following *X2. Safe to be called with X2
601: equal to X1. Return values as for elladd() and elldouble().
602: Uses (a simplified variant of) Peter Montgomery's PRAC (PRactical Addition
603: Chain) algorithm; see ftp://ftp.cwi.nl/pub/pmontgom/Lucas.ps.gz .
604: With thanks to Paul Zimmermann for the reference. --GN1998Aug13 */
605:
606: /* We use an array of GENs pointed at by XAUX as a scratchpad; this will
607: have been set up by ellfacteur() (so we don't need to reinitialize it
608: each time). */
609:
610: static int
611: ellmult(long nbc, ulong k, GEN *X1, GEN *X2) /* k>2 prime, not checked */
612: {
613: long i,d,e,e1,r,av=avma,tetpil;
614: int res;
615: GEN *A=X2, *B=XAUX, *S, *T=XAUX+2*nbc;
616:
617: for (i = 2*nbc; i--; ) { affii(X1[i],XAUX[i]); }
618: tetpil = avma;
619:
620: /* first doubling picks up X1; after this we'll be working in XAUX and
621: X2 only, mostly via A and B and T */
622: if ((res = elldouble(nbc, X1, X2)) != 0)
623: {
624: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
625: return res;
626: }
627:
628: /* split the work at the golden ratio */
629: r = (long)(k*0.61803398875 + .5);
630: d = k - r; e = r - d; /* NB d+e == r, so no danger of ofl below */
631:
632: while (d != e)
633: {
634:
635: /* apply one of the nine transformations from PM's Table 4. We first
636: figure out which, and then go into an eight-way switch, because
637: some of the transformations are similar enough to share code. */
638:
639: if (d <= e + (e>>2)) /* floor(1.25*e) */
640: {
641: if ((d+e)%3 == 0)
642: { i = 0; goto apply; } /* Table 4, rule 1 */
643: else if ((d-e)%6 == 0)
644: { i = 1; goto apply; } /* rule 2 */
645: /* else fall through */
646: }
647: if ((d+3)>>2 <= e) /* equiv to d <= 4*e but cannot ofl */
648: { i = 2; goto apply; } /* rule 3, the most common case */
649: if ((d&1)==(e&1))
650: { i = 1; goto apply; } /* rule 4, which does the same as rule 2 */
651: if (!(d&1))
652: { i = 3; goto apply; } /* rule 5 */
653: if (d%3 == 0)
654: { i = 4; goto apply; } /* rule 6 */
655: if ((d+e)%3 == 0)
656: { i = 5; goto apply; } /* rule 7 */
657: if ((d-e)%3 == 0)
658: { i = 6; goto apply; } /* rule 8 */
659: /* when we get here, e must be even, for otherwise one of rules 4,5
660: would have applied */
661: i = 7; /* rule 9 */
662:
663: apply:
664: switch(i) /* i takes values in {0,...,7} here */
665: {
666: case 0: /* rule 1 */
667: e1 = d - e; d = (d + e1)/3; e = (e - e1)/3;
668: if ((res = elladd(nbc, A, B, T)) != 0)
669: {
670: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
671: return res;
672: }
673: if ((res = elladd2(nbc, T, A, A, T, B, B)) != 0)
674: {
675: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
676: return res;
677: }
678: break; /* end of rule 1 */
679: case 1: /* rules 2 and 4, part 1 */
680: d -= e;
681: if ((res = elladd(nbc, A, B, B)) != 0)
682: {
683: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
684: return res;
685: }
686: /* FALL THROUGH */
687: case 3: /* rule 5, and 2nd part of rules 2 and 4 */
688: d >>= 1;
689: if ((res = elldouble(nbc, A, A)) != 0)
690: {
691: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
692: return res;
693: }
694: break; /* end of rules 2, 4, and 5 */
695: case 4: /* rule 6 */
696: d /= 3;
697: if ((res = elldouble(nbc, A, T)) != 0)
698: {
699: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
700: return res;
701: }
702: if ((res = elladd(nbc, T, A, A)) != 0)
703: {
704: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
705: return res;
706: }
707: /* FALL THROUGH */
708: case 2: /* rule 3, and 2nd part of rule 6 */
709: d -= e;
710: if ((res = elladd(nbc, A, B, B)) != 0)
711: {
712: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
713: return res;
714: }
715: break; /* end of rules 3 and 6 */
716: case 5: /* rule 7 */
717: d = (d - e - e)/3;
718: if ((res = elldouble(nbc, A, T)) != 0)
719: {
720: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
721: return res;
722: }
723: if ((res = elladd2(nbc, T, A, A, T, B, B)) != 0)
724: {
725: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
726: return res;
727: }
728: break; /* end of rule 7 */
729: case 6: /* rule 8 */
730: d = (d - e)/3;
731: if ((res = elladd(nbc, A, B, B)) != 0)
732: {
733: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
734: return res;
735: }
736: if ((res = elldouble(nbc, A, T)) != 0)
737: {
738: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
739: return res;
740: }
741: if ((res = elladd(nbc, T, A, A)) != 0)
742: {
743: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
744: return res;
745: }
746: break; /* end of rule 8 */
747: case 7: /* rule 9 */
748: e >>= 1;
749: if ((res = elldouble(nbc, B, B)) != 0)
750: {
751: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
752: return res;
753: }
754: break; /* end of rule 9 */
755: default: /* never reached */
756: break;
757: }
758: /* end of Table 4 processing */
759:
760: /* swap d <-> e and A <-> B if necessary */
761: if (d < e) { r = d; d = e; e = r; S = A; A = B; B = S; }
762: } /* while */
763: if ((res = elladd(nbc, XAUX, X2, X2)) != 0)
764: {
765: if (res > 1) { gl = gerepile(av,tetpil,gl); } else avma = av;
766: return res;
767: }
768: avma = av; return 0;
769: }
770:
771: /* PRAC implementation notes - main changes against the paper version:
772: (1) The general function [m+n]P = f([m]P,[n]P,[m-n]P) collapses (for
773: m!=n) to an elladd() which does not depend on the third argument; and
774: thus all references to the third variable (C in the paper) can be elimi-
775: nated. (2) Since our multipliers are prime, the outer loop of the paper
776: version executes only once, and thus is invisible above. (3) The first
777: step in the inner loop of the paper version will always be rule 3, but
778: the addition requested by this rule amounts to a doubling, and it will
779: always be followed by a swap, so we have unrolled this first iteration.
780: (4) Some simplifications in rules 6 and 7 are possible given the above,
781: and we can save one addition in each of the two cases. NB one can show
782: that none of the other elladd()s in the loop can ever turn out to de-
783: generate into an elldouble. (5) I tried to optimize for rule 3, which
784: is used far more frequently than all others together, but it didn't
785: improve things, so I removed the nested tight loop again. --GN */
786:
787: /* The main loop body of ellfacteur() runs slightly _slower_ under PRAC than
788: under a straightforward left-shift binary multiplication algorithm when
789: N has <30 digits and B1 is small; PRAC wins when N and B1 get larger.
790: Weird. --GN */
791:
792:
793: /* memory layout in ellfacteur(): We'll have a large-ish array of GEN
794: pointers, and one huge chunk of memory containing all the actual GEN
795: (t_INT) objects.
796: nbc will be held constant throughout the invocation. */
797: /* The B1 stage of each iteration through the main loop needs little
798: space: enough for the X and Y coordinates of the current points,
799: and twice as much again as scratchpad for ellmult(). */
800: /* The B2 stage, starting from some current set of points Q, needs, in
801: succession:
802: - space for [2]Q, [4]Q, ..., [10]Q, and [p]Q for building the helix;
803: - space for 48*nbc X and Y coordinates to hold the helix. Now this
804: could re-use [2]Q,...,[8]Q, but only with difficulty, since we don't
805: know in advance which residue class mod 210 our p is going to be in.
806: It can and should re-use [p]Q, though;
807: - space for (temporarily [30]Q and then) [210]Q, [420]Q, and several
808: further doublings until the giant step multiplier is reached. This
809: _can_ re-use the remaining cells from above. The computation of [210]Q
810: will have been the last call to ellmult() within this iteration of the
811: main loop, so the scratchpad is now also free to be re-used. We also
812: compute [630]Q by a parallel addition; we'll need it later to get the
813: baby-step table bootstrapped a little faster. */
814: /* Finally, for no more than 4 curves at a time, room for up to 1024 X
815: coordinates only (the Y coordinates needed whilst setting up this baby
816: step table are temporarily stored in the upper half, and overwritten
817: during the last series of additions). */
818:
819: /* Graphically: after end of B1 stage (X,Y are the coords of Q):
820: +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
821: | X Y | scratch | [2]Q| [4]Q| [6]Q| [8]Q|[10]Q| ... | ...
822: +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
823: *X *XAUX *XT *XD *XB
824:
825: [30]Q is computed from [10]Q. [210]Q can go into XY, etc:
826: +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
827: |[210]|[420]|[630]|[840]|[1680,3360,6720,...,2048*210] |bstp table...
828: +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
829: *X *XAUX *XT *XD [*XG, somewhere here] *XB .... *XH
830:
831: So we need (13 + 48) * 2 * nbc slots here, and another 4096 slots for
832: the baby step table (not all of which will be used when we start with a
833: small B1, but it's better to allocate and initialize ahead of time all
834: the slots that might be needed later). */
835:
836: /* Note on memory locality: During the B2 phase, accesses to the helix
837: (once it has been set up) will be clustered by curves (4 out of nbc at
838: a time). Accesses to the baby steps table will wander from one end of
839: the array to the other and back, one such cycle per giant step, and
840: during a full cycle we would expect on the order of 2E4 accesses when
841: using the largest giant step size. Thus we shouldn't be doing too bad
842: with respect to thrashing a (512KBy) L2 cache. However, we don't want
843: the baby step table to grow larger than this, even if it would reduce
844: the number of E.C. operations by a few more per cent for very large B2,
845: lest cache thrashing slow down everything disproportionally. --GN */
846:
847: /* parameters for miller() via snextpr(), for use by ellfacteur() */
848: #define miller_k1 16 /* B1 phase, foolproof below 10^12 */
849: #define miller_k2 1 /* B2 phase, not foolproof, much faster */
850: /* (miller_k2 will let thousands of composites slip through, which doesn't
851: harm ECM, but ellmult() during the B1 phase should only be fed primes
852: which really are prime) */
853:
854: /* ellfacteur() has been re-tuned to be useful as a first stage before
855: MPQS, especially for _large_ arguments, when insist is false, and now
856: also for the case when insist is true, vaguely following suggestions
857: by Paul Zimmermann (see http://www.loria.fr/~zimmerma/ and especially
858: http://www.loria.fr/~zimmerma/records/ecmnet.html) of INRIA/LORIA.
859: --GN 1998Jul,Aug */
860: GEN
861: ellfacteur(GEN n, int insist)
862: {
863: static ulong TB1[] =
864: {
865: /* table revised, cf. below 1998Aug15 --GN */
866: 142,172,208,252,305,370,450,545,661,801,972,1180,1430,
867: 1735,2100,2550,3090,3745,4540,5505,6675,8090,9810,11900,
868: 14420,17490,21200,25700,31160,37780UL,45810UL,55550UL,67350UL,
869: 81660UL,99010UL,120050UL,145550UL,176475UL,213970UL,259430UL,
870: 314550UL,381380UL,462415UL,560660UL,679780UL,824220UL,999340UL,
871: 1211670UL,1469110UL,1781250UL,2159700UL,2618600UL,3175000UL,
872: 3849600UL,4667500UL,5659200UL,6861600UL,8319500UL,10087100UL,
873: 12230300UL,14828900UL,17979600UL,21799700UL,26431500UL,
874: 32047300UL,38856400UL, /* 110 times that still fits into 32bits */
875: #ifdef LONG_IS_64BITS
876: 47112200UL,57122100UL,69258800UL,83974200UL,101816200UL,
877: 123449000UL,149678200UL,181480300UL,220039400UL,266791100UL,
878: 323476100UL,392204900UL,475536500UL,576573500UL,699077800UL,
879: 847610500UL,1027701900UL,1246057200UL,1510806400UL,1831806700UL,
880: 2221009800UL,2692906700UL,3265067200UL,3958794400UL,4799917500UL,
881: /* the only reason to stop here is that I got bored (and that users
882: will get bored watching their 64bit machines churning on such large
883: numbers for month after month). Someone can extend this table when
884: the hardware has gotten 100 times faster than now --GN */
885: #endif
886: };
887: static ulong TB1_for_stage[] =
888: {
889: /* table revised 1998Aug11 --GN. The idea is to start a little below
890: the optimal B1 for finding factors which would just have been missed
891: by pollardbrent(), and escalate gradually, changing curves suf-
892: ficiently frequently to give good coverage of the small factor
893: ranges. The table entries grow a bit faster than what Paul says
894: would be optimal, but having a single table instead of a 2D array
895: keeps the code simple */
896: 500,520,560,620,700,800,900,1000,1150,1300,1450,1600,1800,2000,
897: 2200,2450,2700,2950,3250,3600,4000,4400,4850,5300,5800,6400,
898: 7100,7850,8700,9600,10600,11700,12900,14200,15700,17300,
899: 19000,21000,23200,25500,28000,31000,34500UL,38500UL,43000UL,
900: 48000UL,53800UL,60400UL,67750UL,76000UL,85300UL,95700UL,
901: 107400UL,120500UL,135400UL,152000UL,170800UL,191800UL,215400UL,
902: 241800UL,271400UL,304500UL,341500UL,383100UL,429700UL,481900UL,
903: 540400UL,606000UL,679500UL,761800UL,854100UL,957500UL,1073500UL,
904: };
905: long nbc,nbc2,dsn,dsnmax,rep,spc,gse,gss,rcn,rcn0,bstp,bstp0;
906: long a,i,j,k, av,av1,lim, size = expi(n) + 1, tf = lgefint(n);
907: ulong B1,B2,B2_p,B2_rt,m,p,p0,p2,dp;
908: GEN w,w0,x,*X,*XT,*XD,*XG,*YG,*XH,*XB,*XB2,*Xh,*Yh,*Xb, res = cgeti(tf);
909: int rflag, use_clones = 0;
910: byteptr d, d0;
911:
912: av = avma; /* taking res into account */
913: N = n; /* make n known to auxiliary functions */
914: /* determine where we'll start, how long we'll persist, and how many
915: curves we'll use in parallel */
916: if (insist)
917: {
918: dsnmax = (size >> 2) - 10;
919: if (dsnmax < 0) dsnmax = 0;
920: #ifdef LONG_IS_64BITS
921: else if (dsnmax > 90) dsnmax = 90;
922: #else
923: else if (dsnmax > 65) dsnmax = 65;
924: #endif
925: dsn = (size >> 3) - 5;
926: if (dsn < 0) dsn = 0;
927: else if (dsn > 47) dsn = 47;
928: /* pick up the torch where non-insistent stage would have given up */
929: nbc = dsn + (dsn >> 2) + 9; /* 8 or more curves in parallel */
930: nbc &= ~3; /* nbc is always a multiple of 4 */
931: if (nbc > nbcmax) nbc = nbcmax;
932: a = 1 + (nbcmax<<7); /* seed for choice of curves */
933: }
934: else
935: {
936: dsn = (size - 140) >> 3;
937: if (dsn > 12) dsn = 12;
938: dsnmax = 72;
939: if (dsn < 0) /* < 140 bits: decline the task */
940: {
941: #ifdef __EMX__
942: /* MPQS's disk access under DOS/EMX would be abysmally slow, so... */
943: dsn = 0;
944: rep = 20;
945: nbc = 8;
946: #else
947: if (DEBUGLEVEL >= 4)
948: {
949: fprintferr("ECM: number too small to justify this stage\n");
950: flusherr();
951: }
952: avma = av; return NULL;
953: #endif
954: }
955: else
956: {
957: rep = (size <= 248 ?
958: (size <= 176 ? (size - 124) >> 4 : (size - 148) >> 3) :
959: (size - 224) >> 1);
960: nbc = ((size >> 3) << 2) - 80;
961: if (nbc < 8) nbc = 8;
962: else if (nbc > nbcmax) nbc = nbcmax;
963: #ifdef __EMX__
964: rep += 20;
965: #endif
966: }
967:
968: /* it may be convenient to use disjoint sets of curves for the non-insist
969: and insist phases; moreover, repeated non-insistent calls acting on
970: factors of the same original number should try to use fresh curves.
971: The following achieves this */
972: a = 1 + (nbcmax<<3)*(size & 0xf);
973: }
974: if (dsn > dsnmax) dsn = dsnmax;
975:
976: if (DEBUGLEVEL >= 4)
977: {
978: (void) timer2(); /* clear timer */
979: fprintferr("ECM: working on %ld curves at a time; initializing", nbc);
980: if (!insist)
981: {
982: if (rep == 1)
983: fprintferr(" for one round");
984: else
985: fprintferr(" for up to %ld rounds", rep);
986: }
987: fprintferr("...\n");
988: }
989:
990: /* The auxiliary routines above need < (3*nbc+240)*tf words on the PARI
991: stack, in addition to the spc*(tf+1) words occupied by our main table.
992: If stack space is already tight, try the heap, using newbloc() and
993: killbloc() */
994: nbc2 = nbc << 1;
995: spc = (13 + 48) * nbc2 + bstpmax * 4;
996: if ((long)((GEN)avma - (GEN)bot) < spc + 385 + (spc + 3*nbc + 240)*tf)
997: {
998: if (DEBUGLEVEL >= 5)
999: {
1000: fprintferr("ECM: stack tight, using clone space on the heap\n");
1001: }
1002: use_clones = 1;
1003: x = newbloc(spc + 385);
1004: }
1005: else
1006: x = new_chunk(spc + 385);
1007: X = 1 + (GEN*)x; /* B1 phase: current point */
1008: XAUX = X + nbc2; /* scratchpad for ellmult() */
1009: XT = XAUX + nbc2; /* ditto, will later hold [3*210]Q */
1010: XD = XT + nbc2; /* room for various multiples */
1011: XB = XD + 20*nbc; /* start of baby steps table */
1012: XB2 = XB + 2 * bstpmax; /* middle of baby steps table */
1013: XH = XB2 + 2 * bstpmax; /* end of bstps table, start of helix */
1014: Xh = XH + 96*nbc; /* little helix, X coords */
1015: Yh = XH + 192; /* ditto, Y coords */
1016: /* XG will be set later, inside the main loop, since it depends on B2 */
1017:
1018: {
1019: long tw = evallg(tf) | evaltyp(t_INT);
1020:
1021: if (use_clones)
1022: w = newbloc(spc*tf);
1023: else
1024: w = new_chunk(spc*tf);
1025: w0 = w; /* remember this for later... */
1026: for (i = spc; i--; )
1027: {
1028: *w = tw; X[i] = w; w += tf; /* hack for: w = cgeti(tf) */
1029: }
1030: /* Xh range of 384 pointers not set; these will later duplicate the
1031: pointers in the XH range, 4 curves at a time. Some of the cells
1032: reserved here for the XB range will never be used, instead, we'll
1033: warp the pointers to connect to (read-only) GENs in the X/XD range;
1034: it would be complicated to skip them here to conserve merely a few
1035: KBy of stack or heap space. --GN */
1036: }
1037:
1038: /* *** ECM MAIN LOOP *** */
1039: for(;;)
1040: {
1041: d = diffptr; rcn = NPRC; /* multipliers begin at the beginning */
1042:
1043: /* pick curves */
1044: for (i = nbc2; i--; ) affsi(a++, X[i]);
1045: /* pick bounds */
1046: B1 = insist ? TB1[dsn] : TB1_for_stage[dsn];
1047: B2 = 110*B1;
1048: B2_rt = (ulong)(sqrt((double)B2));
1049: /* pick giant step exponent and size.
1050: With 32 baby steps, a giant step corresponds to 32*420 = 13440, appro-
1051: priate for the smallest B2s. With 1024, a giant step will be 430080;
1052: this will be appropriate for B1 >~ 42000, where 512 baby steps would
1053: imply roughly the same number of E.C. additions. */
1054: gse = (B1 < 656 ?
1055: (B1 < 200 ? 5 : 6) :
1056: (B1 < 10500 ?
1057: (B1 < 2625 ? 7 : 8) :
1058: (B1 < 42000 ? 9 : 10)
1059: )
1060: );
1061: gss = 1UL << gse;
1062: XG = XT + gse*nbc2; /* will later hold [2^(gse+1)*210]Q */
1063: YG = XG + nbc;
1064:
1065: if (DEBUGLEVEL >= 4)
1066: {
1067: fprintferr("ECM: time = %6ld ms\nECM: dsn = %2ld,\tB1 = %4lu,",
1068: timer2(), dsn, B1);
1069: fprintferr("\tB2 = %6lu,\tgss = %4ld*420\n", B2, gss);
1070: flusherr();
1071: }
1072: p = *d++;
1073:
1074: /* ---B1 PHASE--- */
1075: /* treat p=2 separately */
1076: B2_p = B2 >> 1;
1077: for (m=1; m<=B2_p; m<<=1)
1078: {
1079: if ((rflag = elldouble(nbc, X, X)) > 1) goto fin;
1080: else if (rflag) break;
1081: }
1082:
1083: /* p=3,...,nextprime(B1) */
1084: while (p < B1 && p <= B2_rt)
1085: {
1086: p = snextpr(p, &d, &rcn, NULL, miller_k1);
1087: B2_p = B2/p; /* beware integer overflow on 32-bit CPUs */
1088: for (m=1; m<=B2_p; m*=p)
1089: {
1090: if ((rflag = ellmult(nbc, p, X, X)) > 1) goto fin;
1091: else if (rflag) break;
1092: }
1093: }
1094: /* primes p larger than sqrt(B2) can appear only to the 1st power */
1095: while (p < B1)
1096: {
1097: p = snextpr(p, &d, &rcn, NULL, miller_k1);
1098: if (ellmult(nbc, p, X, X) > 1) goto fin; /* p^2 > B2: no loop */
1099: }
1100:
1101: if (DEBUGLEVEL >= 4)
1102: {
1103: fprintferr("ECM: time = %6ld ms, B1 phase done, ", timer2());
1104: fprintferr("p = %lu, setting up for B2\n", p);
1105: }
1106:
1107: /* ---B2 PHASE--- */
1108: /* compute [2]Q,...,[10]Q, which we need to build the helix */
1109: if (elldouble(nbc, X, XD) > 1)
1110: goto fin; /* [2]Q */
1111: if (elldouble(nbc, XD, XD + nbc2) > 1)
1112: goto fin; /* [4]Q */
1113: if (elladd(nbc, XD, XD + nbc2, XD + (nbc<<2)) > 1)
1114: goto fin; /* [6]Q */
1115: if (elladd2(nbc,
1116: XD, XD + (nbc<<2), XT + (nbc<<3),
1117: XD + nbc2, XD + (nbc<<2), XD + (nbc<<3)) > 1)
1118: goto fin; /* [8]Q and [10]Q */
1119: if (DEBUGLEVEL >= 7)
1120: fprintferr("\t(got [2]Q...[10]Q)\n");
1121:
1122: /* get next prime (still using the foolproof test) */
1123: p = snextpr(p, &d, &rcn, NULL, miller_k1);
1124: /* make sure we have the residue class number (mod 210) */
1125: if (rcn == NPRC)
1126: {
1127: rcn = prc210_no[(p % 210) >> 1];
1128: if (rcn == NPRC)
1129: {
1130: fprintferr("ECM: %lu should have been prime but isn\'t\n", p);
1131: err(bugparier, "ellfacteur");
1132: }
1133: }
1134:
1135: /* compute [p]Q and put it into its place in the helix */
1136: if (ellmult(nbc, p, X, XH + rcn*nbc2) > 1) goto fin;
1137: if (DEBUGLEVEL >= 7)
1138: fprintferr("\t(got [p]Q, p = %lu = %lu mod 210)\n",
1139: p, (ulong)(prc210_rp[rcn]));
1140:
1141: /* save current p, d, and rcn; we'll need them more than once below */
1142: p0 = p;
1143: d0 = d;
1144: rcn0 = rcn; /* remember where the helix wraps */
1145: bstp0 = 0; /* p is at baby-step offset 0 from itself */
1146:
1147: /* fill up the helix, stepping forward through the prime residue classes
1148: mod 210 until we're back at the r'class of p0. Keep updating p so
1149: that we can print meaningful diagnostics if a factor shows up; but
1150: don't bother checking which of these p's are in fact prime */
1151: for (i = 47; i; i--) /* 47 iterations */
1152: {
1153: p += (dp = (ulong)prc210_d1[rcn]);
1154: if (rcn == 47)
1155: { /* wrap mod 210 */
1156: if (elladd(nbc, XT + dp*nbc, XH + rcn*nbc2, XH) > 1)
1157: goto fin;
1158: rcn = 0;
1159: continue;
1160: }
1161: if (elladd(nbc, XT + dp*nbc, XH + rcn*nbc2, XH + rcn*nbc2 + nbc2) > 1)
1162: goto fin;
1163: rcn++;
1164: }
1165: if (DEBUGLEVEL >= 7)
1166: fprintferr("\t(got initial helix)\n");
1167:
1168: /* compute [210]Q etc, which will be needed for the baby step table */
1169: if (ellmult(nbc, 3, XD + (nbc<<3), X) > 1) goto fin;
1170: if (ellmult(nbc, 7, X, X) > 1) goto fin; /* [210]Q */
1171: /* this was the last call to ellmult() in the main loop body; may now
1172: overwrite XAUX and slots XD and following */
1173: if (elldouble(nbc, X, XAUX) > 1) goto fin; /* [420]Q */
1174: if (elladd(nbc, X, XAUX, XT) > 1) goto fin; /* [630]Q */
1175: if (elladd(nbc, X, XT, XD) > 1) goto fin; /* [840]Q */
1176: for (i=1; i <= gse; i++) /* gse successive doublings */
1177: {
1178: if (elldouble(nbc, XT + i*nbc2, XD + i*nbc2) > 1) goto fin;
1179: }
1180: /* (the last iteration has initialized XG to [210*2^(gse+1)]Q) */
1181:
1182: if (DEBUGLEVEL >= 4)
1183: {
1184: fprintferr("ECM: time = %6ld ms, entering B2 phase, p = %lu\n",
1185: timer2(), p);
1186: }
1187:
1188: /* inner loop over small sets of 4 curves at a time */
1189: for (i = nbc - 4; i >= 0; i -= 4)
1190: {
1191: if (DEBUGLEVEL >= 6)
1192: fprintferr("ECM: finishing curves %ld...%ld\n", i, i+3);
1193: /* copy relevant pointers from XH to Xh. Recall memory layout in XH
1194: is: nbc X coordinates followed by nbc Y coordinates for residue
1195: class 1 mod 210, then the same for r.c. 11 mod 210, etc. Memory
1196: layout for Xh is: four X coords for 1 mod 210, four for 11 mod 210,
1197: etc, four for 209 mod 210, and then the corresponding Y coordinates
1198: in the same order. This will allow us to do a giant step on Xh
1199: using just three calls to elladd0() each acting on 64 points in
1200: parallel */
1201: for (j = 48; j--; )
1202: {
1203: k = nbc2*j + i;
1204: m = j << 2; /* X coordinates */
1205: Xh[m] = XH[k]; Xh[m+1] = XH[k+1];
1206: Xh[m+2] = XH[k+2]; Xh[m+3] = XH[k+3];
1207: k += nbc; /* Y coordinates */
1208: Yh[m] = XH[k]; Yh[m+1] = XH[k+1];
1209: Yh[m+2] = XH[k+2]; Yh[m+3] = XH[k+3];
1210: }
1211: /* build baby step table of X coords of multiples of [210]Q. XB[4*j]
1212: will point at X coords on four curves from [(j+1)*210]Q. Until
1213: we're done, we need some Y coords as well, which we keep in the
1214: second half of the table, overwriting them at the end when gse==10.
1215: Those multiples which we already have (by 1,2,3,4,8,16,...,2^gse)
1216: are entered simply by copying the pointers, ignoring the small
1217: number of slots in w that were initially reserved for them.
1218: Here are the initial entries... */
1219: for (Xb=XB,k=2,j=i; k--; Xb=XB2,j+=nbc) /* do first X, then Y coords */
1220: {
1221: Xb[0] = X[j]; Xb[1] = X[j+1]; /* [210]Q */
1222: Xb[2] = X[j+2]; Xb[3] = X[j+3];
1223: Xb[4] = XAUX[j]; Xb[5] = XAUX[j+1]; /* [420]Q */
1224: Xb[6] = XAUX[j+2]; Xb[7] = XAUX[j+3];
1225: Xb[8] = XT[j]; Xb[9] = XT[j+1]; /* [630]Q */
1226: Xb[10] = XT[j+2]; Xb[11] = XT[j+3];
1227: Xb += 4; /* this points at [420]Q */
1228: /* ... entries at powers of 2 times 210 .... */
1229: for (m = 2; m < gse+k; m++) /* omit Y coords of [2^gse*210]Q */
1230: {
1231: long m2 = m*nbc2 + j;
1232: Xb += (2UL<<m); /* points now at [2^m*210]Q */
1233: Xb[0] = XAUX[m2]; Xb[1] = XAUX[m2+1];
1234: Xb[2] = XAUX[m2+2]; Xb[3] = XAUX[m2+3];
1235: }
1236: }
1237: if (DEBUGLEVEL >= 7)
1238: fprintferr("\t(extracted precomputed helix / baby step entries)\n");
1239: /* ... glue in between, up to 16*210 ... */
1240: if (elladd0(12, 4, /* 12 pts + (4 pts replicated thrice) */
1241: XB + 12, XB2 + 12,
1242: XB, XB2,
1243: XB + 16, XB2 + 16)
1244: > 1) goto fin; /* 4 + {1,2,3} = {5,6,7} */
1245: if (elladd0(28, 4, /* 28 pts + (4 pts replicated 7fold) */
1246: XB + 28, XB2 + 28,
1247: XB, XB2,
1248: XB + 32, XB2 + 32)
1249: > 1) goto fin; /* 8 + {1,...,7} = {9,...,15} */
1250: /* ... and the remainder of the lot */
1251: for (m = 5; m <= gse; m++)
1252: {
1253: /* fill in from 2^(m-1)+1 to 2^m-1 in chunks of 64 and 60 points */
1254: ulong m2 = 2UL << m; /* will point at 2^(m-1)+1 */
1255: for (j = 0; j < m2-64; j+=64) /* executed 0 times when m == 5 */
1256: {
1257: if (elladd0(64, 4,
1258: XB + m2 - 4, XB2 + m2 - 4,
1259: XB + j, XB2 + j,
1260: XB + m2 + j,
1261: (m<gse ? XB2 + m2 + j : NULL))
1262: > 1) goto fin;
1263: } /* j == m2-64 here, 60 points left */
1264: if (elladd0(60, 4,
1265: XB + m2 - 4, XB2 + m2 - 4,
1266: XB + j, XB2 + j,
1267: XB + m2 + j,
1268: (m<gse ? XB2 + m2 + j : NULL))
1269: > 1) goto fin;
1270: /* (when m==gse, drop Y coords of result, and when both equal 1024,
1271: overwrite Y coords of second argument with X coords of result) */
1272: }
1273: if (DEBUGLEVEL >= 7)
1274: fprintferr("\t(baby step table complete)\n");
1275: /* initialize a few other things */
1276: bstp = bstp0;
1277: p = p0; d = d0; rcn = rcn0;
1278: gl = gun;
1279: av1 = avma;
1280: lim=stack_lim(av1,1);
1281: /* the correct entry in XB to use depends on bstp and on where we are
1282: on the helix. As we skip from prime to prime, bstp will be incre-
1283: mented by snextpr() each time we wrap around through residue class
1284: number 0 (1 mod 210), but the baby step should not be taken until
1285: rcn>=rcn0 (i.e. until we pass again the residue class of p0).
1286: The correct signed multiplier is thus k = bstp - (rcn < rcn0),
1287: and the offset from XB is four times (|k| - 1). When k==0, we may
1288: ignore the current prime (if it had led to a factorization, this
1289: would have been noted during the last giant step, or -- when we
1290: first get here -- whilst initializing the helix). When k > gss,
1291: we must do a giant step and bump bstp back by -2*gss.
1292: The gcd of the product of X coord differences against N is taken just
1293: before we do a giant step. */
1294:
1295: /* loop over probable primes p0 < p <= nextprime(B2),
1296: inserting giant steps as necessary */
1297: while (p < B2)
1298: {
1299: /* save current p for diagnostics */
1300: p2 = p;
1301: /* get next probable prime */
1302: p = snextpr(p, &d, &rcn, &bstp, miller_k2);
1303: /* work out the corresponding baby-step multiplier */
1304: k = bstp - (rcn < rcn0 ? 1 : 0);
1305: /* check whether it's giant-step time */
1306: if (k > gss)
1307: {
1308: /* take gcd */
1309: gl = mppgcd(gl, n);
1310: if (!is_pm1(gl) && !egalii(gl, n)) { p = p2; goto fin; }
1311: gl = gun;
1312: avma = av1;
1313: while (k > gss) /* hm, just how large are those prime gaps? */
1314: {
1315: /* giant step */
1316: if (DEBUGLEVEL >= 7)
1317: fprintferr("\t(giant step at p = %lu)\n", p);
1318: if (elladd0(64, 4,
1319: XG + i, YG + i,
1320: Xh, Yh, Xh, Yh) > 1) goto fin;
1321: if (elladd0(64, 4,
1322: XG + i, YG + i,
1323: Xh + 64, Yh + 64, Xh + 64, Yh + 64) > 1) goto fin;
1324: if (elladd0(64, 4,
1325: XG + i, YG + i,
1326: Xh + 128, Yh + 128, Xh + 128, Yh + 128)
1327: > 1) goto fin;
1328: bstp -= (gss << 1);
1329: /* recompute multiplier */
1330: k = bstp - (rcn < rcn0 ? 1 : 0);
1331: }
1332: }
1333: if (!k) continue; /* point of interest is already in Xh */
1334: if (k < 0) k = -k;
1335: m = ((ulong)k - 1) << 2;
1336: /* accumulate product of differences of X coordinates */
1337: j = rcn<<2;
1338: gl = modii(mulii(gl, subii(XB[m], Xh[j])), n);
1339: gl = modii(mulii(gl, subii(XB[m+1], Xh[j+1])), n);
1340: gl = modii(mulii(gl, subii(XB[m+2], Xh[j+2])), n);
1341: gl = modii(mulii(gl, subii(XB[m+3], Xh[j+3])), n);
1342: if (low_stack(lim, stack_lim(av1,1)))
1343: {
1344: if(DEBUGMEM>1) err(warnmem,"ellfacteur");
1345: gl = gerepileupto(av1, gl);
1346: }
1347: } /* loop over p */
1348: avma = av1;
1349: } /* for i (loop over sets of 4 curves) */
1350:
1351: /* continuation part of main loop */
1352:
1353: if (dsn < dsnmax)
1354: {
1355: dsn += insist ? 1 : 2;
1356: if (dsn > dsnmax) dsn = dsnmax;
1357: }
1358:
1359: if (!insist && !--rep)
1360: {
1361: if (DEBUGLEVEL >= 4)
1362: {
1363: fprintferr("ECM: time = %6ld ms,\tellfacteur giving up.\n",
1364: timer2());
1365: flusherr();
1366: }
1367: avma = av;
1368: if (use_clones) { gunclone(w0); gunclone(x); }
1369: return NULL;
1370: }
1371: }
1372: /* *** END OF ECM MAIN LOOP *** */
1373: fin:
1374: affii(gl, res);
1375:
1376: if (DEBUGLEVEL >= 4)
1377: {
1378: fprintferr("ECM: time = %6ld ms,\tp <= %6lu,\n\tfound factor = %Z\n",
1379: timer2(), p, res);
1380: flusherr();
1381: }
1382: avma=av;
1383: if (use_clones) { gunclone(w0); gunclone(x); }
1384: return res;
1385: }
1386:
1387: /***********************************************************************/
1388: /** **/
1389: /** FACTORIZATION (Pollard-Brent rho) **/
1390: /** pollardbrent() returns a non trivial factor of n, assuming n is **/
1391: /** composite and has no small prime divisor, or NULL if going on **/
1392: /** would take more time than we want to spend. GN1998Jun18-26 **/
1393: /** (Cf. Algorithm 8.5.2 in ACiCNT) **/
1394: /** **/
1395: /***********************************************************************/
1396: static void
1397: rho_dbg(long c, long msg_mask)
1398: {
1399: if (c & msg_mask) return;
1400: fprintferr("Rho: time = %6ld ms,\t%3ld round%s\n",
1401: timer2(), c, (c==1?"":"s"));
1402: flusherr();
1403: }
1404:
1405: /* Tuning parameter: for input up to 64 bits long, we must not spend more
1406: * than a very short time, for fear of slowing things down on average.
1407: * With the current tuning formula, increase our efforts somewhat at 49 bit
1408: * input (an extra round for each bit at first), and go up more and more
1409: * rapidly after we pass 80 bits. */
1410:
1411: #define tune_pb_min 14 /* even 15 seems too much */
1412:
1413: /* We return NULL when we run out of time, or a single t_INT containing a
1414: nontrivial factor of n, or a vector of t_INTs, each triple of successive
1415: entries containing a factor, an exponent (equal to un), and a factor
1416: class (NULL for unknown or zero for known composite), matching the
1417: internal representation used by the ifac_*() routines below. Repeated
1418: factors can arise and are legal; the caller will be sorting the factors
1419: anyway. */
1420: GEN
1421: pollardbrent(GEN n)
1422: {
1423: long tf = lgefint(n), size = 0, delta, retries = 0, msg_mask;
1424: long c0, c, k, k1, l, avP, avx, GGG, av = avma;
1425: GEN x, x1, y, P, g, g1, res;
1426:
1427: if (DEBUGLEVEL > 3) (void)timer2(); /* clear timer */
1428:
1429: if (tf >= 4)
1430: size = expi(n) + 1;
1431: else if (tf == 3) /* try to keep purify happy... */
1432: size = BITS_IN_LONG - bfffo(n[2]);
1433:
1434: if (size <= 32)
1435: c0 = 32; /* amounts very nearly to `insist' */
1436: else if (size <= 48)
1437: c0 = tune_pb_min;
1438: else if (size <= 72)
1439: c0 = tune_pb_min + size - 24;
1440: else if (size <= 301)
1441: /* nonlinear increase in effort, kicking in around 80 bits */
1442: /* 301 gives 48121 + tune_pb_min */
1443: c0 = tune_pb_min + size - 60 +
1444: ((size-73)>>1)*((size-70)>>3)*((size-56)>>4);
1445: else
1446: c0 = 49152; /* ECM is faster when it'd take longer */
1447:
1448: c = c0 << 5; /* 32 iterations per round */
1449: msg_mask = (size >= 448? 0x1fff:
1450: (size >= 192? (256L<<((size-128)>>6))-1: 0xff));
1451: PB_RETRY:
1452: /* trick to make a `random' choice determined by n. Don't use x^2+0 or
1453: * x^2-2, ever. Don't use x^2-3 or x^2-7 with a starting value of 2.
1454: * x^2+4, x^2+9 are affine conjugate to x^2+1, so don't use them either.
1455: *
1456: * (the point being that when we get called again on a composite cofactor
1457: * of something we've already seen, we had better avoid the same delta) */
1458: switch ((size + retries) & 7)
1459: {
1460: case 0: delta= 1; break;
1461: case 1: delta= -1; break;
1462: case 2: delta= 3; break;
1463: case 3: delta= 5; break;
1464: case 4: delta= -5; break;
1465: case 5: delta= 7; break;
1466: case 6: delta= 11; break;
1467: case 7: delta=-11; break;
1468: }
1469: if (DEBUGLEVEL > 3)
1470: {
1471: if (!retries)
1472: {
1473: if (size < 1536)
1474: fprintferr("Rho: searching small factor of %ld-bit integer\n", size);
1475: else
1476: fprintferr("Rho: searching small factor of %ld-word integer\n", tf-2);
1477: }
1478: else
1479: fprintferr("Rho: restarting for remaining rounds...\n");
1480: fprintferr("Rho: using X^2%+1ld for up to %ld rounds of 32 iterations\n",
1481: delta, c >> 5);
1482: flusherr();
1483: }
1484: x=gdeux; P=gun; g1 = NULL; k = 1; l = 1;
1485: (void)new_chunk(10 + 6 * tf); /* enough for cgetg(10) + 3 divii */
1486: y = cgeti(tf); affsi(2, y);
1487: x1= cgeti(tf); affsi(2, x1);
1488: avx = avma;
1489: avP = (long)new_chunk(2 * tf); /* enough for x = addsi(tf+1) */
1490: GGG = (long)new_chunk(4 * tf); /* enough for P = modii(2tf+1, tf) */
1491:
1492: for (;;) /* terminated under the control of c */
1493: {
1494: /* use the polynomial x^2 + delta */
1495: #define one_iter() {\
1496: avma = GGG; x = resii(sqri(x), n); /* to garbage zone */\
1497: avma = avx; x = addsi(delta,x); /* erase garbage */\
1498: avma = GGG; P = mulii(P, subii(x1, x));\
1499: avma = avP; P = modii(P,n); }
1500:
1501: one_iter();
1502:
1503: if ((--c & 0x1f)==0) /* one round complete */
1504: {
1505: g = mppgcd(n, P);
1506: if (!is_pm1(g)) goto fin; /* caught something */
1507: if (c <= 0)
1508: { /* getting bored */
1509: if (DEBUGLEVEL > 3)
1510: {
1511: fprintferr("Rho: time = %6ld ms,\tPollard-Brent giving up.\n",
1512: timer2());
1513: flusherr();
1514: }
1515: avma=av; return NULL;
1516: }
1517: P = gun; /* not necessary, but saves 1 mulii/round */
1518: if (DEBUGLEVEL > 3) rho_dbg(c0-(c>>5), msg_mask);
1519: affii(x,y);
1520: }
1521:
1522: if (--k) continue; /* normal end of loop body */
1523:
1524: if (c & 0x1f) /* otherwise, we already checked */
1525: {
1526: g = mppgcd(n, P);
1527: if (!is_pm1(g)) goto fin;
1528: P = gun;
1529: }
1530:
1531: /* Fast forward phase, doing l inner iterations without computing gcds.
1532: * Check first whether it would take us beyond the alloted time.
1533: * Fast forward rounds count only half (although they're taking
1534: * more like 2/3 the time of normal rounds). This to counteract the
1535: * nuisance that all c0 between 4096 and 6144 would act exactly as
1536: * 4096; with the halving trick only the range 4096..5120 collapses
1537: * (similarly for all other powers of two) */
1538: if ((c-=(l>>1)) <= 0)
1539: { /* got bored */
1540: if (DEBUGLEVEL > 3)
1541: {
1542: fprintferr("Rho: time = %6ld ms,\tPollard-Brent giving up.\n",
1543: timer2());
1544: flusherr();
1545: }
1546: avma=av; return NULL;
1547: }
1548: c &= ~0x1f; /* keep it on multiples of 32 */
1549:
1550: /* Fast forward loop */
1551: affii(x, x1); k = l; l <<= 1;
1552: /* don't show this for the first several (short) fast forward phases. */
1553: if (DEBUGLEVEL > 3 && (l>>7) > msg_mask)
1554: {
1555: fprintferr("Rho: fast forward phase (%ld rounds of 64)...\n", l>>7);
1556: flusherr();
1557: }
1558: for (k1=k; k1; k1--) one_iter();
1559: if (DEBUGLEVEL > 3 && (l>>7) > msg_mask)
1560: {
1561: fprintferr("Rho: time = %6ld ms,\t%3ld rounds, back to normal mode\n",
1562: timer2(), c0-(c>>5));
1563: flusherr();
1564: }
1565:
1566: affii(x,y);
1567: } /* forever */
1568:
1569: fin:
1570: /* An accumulated gcd was > 1 */
1571: /* if it isn't n, and looks prime, return it */
1572: if (!egalii(g,n))
1573: {
1574: if (miller(g,17))
1575: {
1576: if (DEBUGLEVEL > 3)
1577: {
1578: rho_dbg(c0-(c>>5), 0);
1579: fprintferr("\tfound factor = %Z\n",g);
1580: flusherr();
1581: }
1582: avma=av; return icopy(g);
1583: }
1584: avma = avx; g1 = icopy(g); /* known composite, keep it safe */
1585: avx = avma;
1586: }
1587: else g1 = n; /* and work modulo g1 for backtracking */
1588:
1589: /* Here g1 is known composite */
1590: if (DEBUGLEVEL > 3 && size > 192)
1591: {
1592: fprintferr("Rho: hang on a second, we got something here...\n");
1593: flusherr();
1594: }
1595: for(;;) /* backtrack until period recovered. Must terminate */
1596: {
1597: avma = GGG; y = resii(sqri(y), g1);
1598: avma = avx; y = addsi(delta,y);
1599: g = mppgcd(subii(x1, y), g1);
1600: if (!is_pm1(g)) break;
1601:
1602: if (DEBUGLEVEL > 3 && (--c & 0x1f) == 0) rho_dbg(c0-(c>>5), msg_mask);
1603: }
1604:
1605: avma = av; /* safe */
1606: if (g1 == n || egalii(g,g1))
1607: {
1608: if (g1 == n && egalii(g,g1))
1609: { /* out of luck */
1610: if (DEBUGLEVEL > 3)
1611: {
1612: rho_dbg(c0-(c>>5), 0);
1613: fprintferr("\tPollard-Brent failed.\n"); flusherr();
1614: }
1615: if (++retries >= 4) return NULL;
1616: goto PB_RETRY;
1617: }
1618: /* half lucky: we've split n, but g1 equals either g or n */
1619: if (DEBUGLEVEL > 3)
1620: {
1621: rho_dbg(c0-(c>>5), 0);
1622: fprintferr("\tfound %sfactor = %Z\n",
1623: (g1!=n ? "composite " : ""), g);
1624: flusherr();
1625: }
1626: res = cgetg(7, t_VEC);
1627: res[1] = licopy(g); /* factor */
1628: res[2] = un; /* exponent 1 */
1629: res[3] = (g1!=n? zero: (long)NULL); /* known composite when g1!=n */
1630:
1631: res[4] = ldivii(n,g); /* cofactor */
1632: res[5] = un; /* exponent 1 */
1633: res[6] = (long)NULL; /* unknown */
1634: return res;
1635: }
1636: /* g < g1 < n : our lucky day -- we've split g1, too */
1637: res = cgetg(10, t_VEC);
1638: /* unknown status for all three factors */
1639: res[1] = licopy(g); res[2] = un; res[3] = (long)NULL;
1640: res[4] = ldivii(g1,g); res[5] = un; res[6] = (long)NULL;
1641: res[7] = ldivii(n,g1); res[8] = un; res[9] = (long)NULL;
1642: if (DEBUGLEVEL > 3)
1643: {
1644: rho_dbg(c0-(c>>5), 0);
1645: fprintferr("\tfound factors = %Z, %Z,\n\tand %Z\n",
1646: res[1], res[4], res[7]);
1647: flusherr();
1648: }
1649: return res;
1650: }
1651:
1652: /***********************************************************************/
1653: /** **/
1654: /** DETECTING ODD POWERS **/
1655: /** Factoring engines like MPQS which ultimately rely on computing **/
1656: /** gcd(N, x^2-y^2) to find a nontrivial factor of N are fundamen- **/
1657: /** tally incapable of splitting a proper power of an odd prime, **/
1658: /** because of the cyclicity of the prime residue class group. We **/
1659: /** already have a square-detection function carrecomplet(), which **/
1660: /** also returns the square root if appropriate. Here's an analogue **/
1661: /** for cubes, fifth and 7th powers. 11th powers are a non-issue so **/
1662: /** long as mpqs() gives up beyond 100 decimal digits (since ECM **/
1663: /** easily find a 10-digit prime factor of a 100-digit number). **/
1664: /** GN1998Jun28 **/
1665: /** **/
1666: /***********************************************************************/
1667:
1668: /* Use a multistage sieve. First stages work mod 211, 209, 61, 203;
1669: if the argument is larger than a word, we first reduce mod the product
1670: of these and then take the remainder apart. Second stages use 117,
1671: 31, 43, 71 in this order. Moduli which are no longer interesting are
1672: skipped. Everything is encoded in a single table of 106 24-bit masks.
1673: We only need the first half of the residues. Three bits per modulus
1674: indicate which residues are 7th (bit 2), 5th (bit 1) powers or cubes
1675: (bit 0); the eight moduli above are assigned right-to-left. The table
1676: will err on the side of safety if one of the moduli divides the number
1677: to be tested, but as this leads to inefficiency it should still be
1678: avoided. */
1679:
1680: static ulong powersmod[106] = {
1681: 077777777ul, /* 0 */
1682: 077777777ul, /* 1 */
1683: 013562440ul, /* 2 */
1684: 012462540ul, /* 3 */
1685: 013562440ul, /* 4 */
1686: 052662441ul, /* 5 */
1687: 016663440ul, /* 6 */
1688: 016463450ul, /* 7 */
1689: 013573551ul, /* 8 */
1690: 012462540ul, /* 9 */
1691: 012462464ul, /* 10 */
1692: 013462771ul, /* 11 */
1693: 012466473ul, /* 12 */
1694: 012463641ul, /* 13 */
1695: 052463646ul, /* 14 */
1696: 012563446ul, /* 15 */
1697: 013762440ul, /* 16 */
1698: 052766440ul, /* 17 */
1699: 012772451ul, /* 18 */
1700: 012762454ul, /* 19 */
1701: 032763550ul, /* 20 */
1702: 013763664ul, /* 21 */
1703: 017763460ul, /* 22 */
1704: 037762565ul, /* 23 */
1705: 017762540ul, /* 24 */
1706: 057762441ul, /* 25 */
1707: 037772452ul, /* 26 */
1708: 017773551ul, /* 27 */
1709: 017767541ul, /* 28 */
1710: 017767640ul, /* 29 */
1711: 037766450ul, /* 30 */
1712: 017762752ul, /* 31 */
1713: 037762762ul, /* 32 */
1714: 017762742ul, /* 33 */
1715: 037763762ul, /* 34 */
1716: 017763740ul, /* 35 */
1717: 077763740ul, /* 36 */
1718: 077762750ul, /* 37 */
1719: 077762752ul, /* 38 */
1720: 077762750ul, /* 39 */
1721: 077762743ul, /* 40 */
1722: 077767740ul, /* 41 */
1723: 077763741ul, /* 42 */
1724: 077763762ul, /* 43 */
1725: 077772760ul, /* 44 */
1726: 077762770ul, /* 45 */
1727: 077766750ul, /* 46 */
1728: 077762740ul, /* 47 */
1729: 077763740ul, /* 48 */
1730: 077763750ul, /* 49 */
1731: 077763752ul, /* 50 */
1732: 077762740ul, /* 51 */
1733: 077762740ul, /* 52 */
1734: 077772740ul, /* 53 */
1735: 077762762ul, /* 54 */
1736: 077763765ul, /* 55 */
1737: 077763770ul, /* 56 */
1738: 077767750ul, /* 57 */
1739: 077766753ul, /* 58 */
1740: 077776740ul, /* 59 */
1741: 077772741ul, /* 60 */
1742: 077772744ul, /* 61 */
1743: 077773740ul, /* 62 */
1744: 077773743ul, /* 63 */
1745: 077773751ul, /* 64 */
1746: 077772771ul, /* 65 */
1747: 077772760ul, /* 66 */
1748: 077772763ul, /* 67 */
1749: 077772751ul, /* 68 */
1750: 077773750ul, /* 69 */
1751: 077777740ul, /* 70 */
1752: 077773745ul, /* 71 */
1753: 077772740ul, /* 72 */
1754: 077772742ul, /* 73 */
1755: 077772744ul, /* 74 */
1756: 077776750ul, /* 75 */
1757: 077773771ul, /* 76 */
1758: 077773774ul, /* 77 */
1759: 077773760ul, /* 78 */
1760: 077772741ul, /* 79 */
1761: 077772740ul, /* 80 */
1762: 077772740ul, /* 81 */
1763: 077772741ul, /* 82 */
1764: 077773754ul, /* 83 */
1765: 077773750ul, /* 84 */
1766: 077773740ul, /* 85 */
1767: 077776741ul, /* 86 */
1768: 077776771ul, /* 87 */
1769: 077776773ul, /* 88 */
1770: 077772761ul, /* 89 */
1771: 077773741ul, /* 90 */
1772: 077773740ul, /* 91 */
1773: 077773740ul, /* 92 */
1774: 077772740ul, /* 93 */
1775: 077772752ul, /* 94 */
1776: 077772750ul, /* 95 */
1777: 077772751ul, /* 96 */
1778: 077773741ul, /* 97 */
1779: 077773761ul, /* 98 */
1780: 077777760ul, /* 99 */
1781: 077772765ul, /* 100 */
1782: 077772742ul, /* 101 */
1783: 077777751ul, /* 102 */
1784: 077777750ul, /* 103 */
1785: 077777745ul, /* 104 */
1786: 077777770ul /* 105 */
1787: };
1788:
1789: /* Returns 3, 5, or 7 if x is a cube (but not a 5th or 7th power), a 5th
1790: power (but not a 7th), or a 7th power, and in this case creates the
1791: base on the stack and assigns its address to *pt. Otherwise returns 0.
1792: x must be of type t_INT and nonzero; this is not checked. The *mask
1793: argument tells us which things to check -- bit 0: 3rd, bit 1: 5th,
1794: bit 2: 7th pwr; set a bit to have the corresponding power examined --
1795: and is updated appropriately for a possible follow-up call */
1796:
1797: long /* no longer static -- used in mpqs.c */
1798: is_odd_power(GEN x, GEN *pt, long *mask)
1799: {
1800: long av=avma, tetpil, lgx=lgefint(x), exponent=0, residue, resbyte;
1801: GEN y;
1802:
1803: *mask &= 7; /* paranoia */
1804: if (!*mask) return 0; /* useful when running in a loop */
1805: if (signe(x) < 0) x=absi(x);
1806:
1807: if (DEBUGLEVEL >= 5)
1808: {
1809: fprintferr("OddPwrs: is %Z\n\t...a", x);
1810: if (*mask&1) fprintferr(" 3rd%s",
1811: (*mask==7?",":(*mask!=1?" or":"")));
1812: if (*mask&2) fprintferr(" 5th%s",
1813: (*mask==7?", or":(*mask&4?" or":"")));
1814: if (*mask&4) fprintferr(" 7th");
1815: fprintferr(" power?\n");
1816: }
1817: if (lgx > 3) residue = smodis(x, 211*209*61*203);
1818: else residue = x[2];
1819:
1820: resbyte=residue%211; if (resbyte > 105) resbyte = 211 - resbyte;
1821: *mask &= powersmod[resbyte];
1822: if (DEBUGLEVEL >= 5)
1823: {
1824: fprintferr("\tmodulo: resid. (remaining possibilities)\n");
1825: fprintferr("\t 211: %3ld (3rd %ld, 5th %ld, 7th %ld)\n",
1826: resbyte, *mask&1, (*mask>>1)&1, (*mask>>2)&1);
1827: }
1828: if (!*mask) { avma=av; return 0; }
1829:
1830: if (*mask & 3)
1831: {
1832: resbyte=residue%209; if (resbyte > 104) resbyte = 209 - resbyte;
1833: *mask &= (powersmod[resbyte] >> 3);
1834: if (DEBUGLEVEL >= 5)
1835: fprintferr("\t 209: %3ld (3rd %ld, 5th %ld, 7th %ld)\n",
1836: resbyte, *mask&1, (*mask>>1)&1, (*mask>>2)&1);
1837: if (!*mask) { avma=av; return 0; }
1838: }
1839: if (*mask & 3)
1840: {
1841: resbyte=residue%61; if (resbyte > 30) resbyte = 61 - resbyte;
1842: *mask &= (powersmod[resbyte] >> 6);
1843: if (DEBUGLEVEL >= 5)
1844: fprintferr("\t 61: %3ld (3rd %ld, 5th %ld, 7th %ld)\n",
1845: resbyte, *mask&1, (*mask>>1)&1, (*mask>>2)&1);
1846: if (!*mask) { avma=av; return 0; }
1847: }
1848: if (*mask & 5)
1849: {
1850: resbyte=residue%203; if (resbyte > 101) resbyte = 203 - resbyte;
1851: *mask &= (powersmod[resbyte] >> 9);
1852: if (DEBUGLEVEL >= 5)
1853: fprintferr("\t 203: %3ld (3rd %ld, 5th %ld, 7th %ld)\n",
1854: resbyte, *mask&1, (*mask>>1)&1, (*mask>>2)&1);
1855: if (!*mask) { avma=av; return 0; }
1856: }
1857:
1858: if (lgx > 3) residue = smodis(x, 117*31*43*71);
1859: else residue = x[2];
1860:
1861: if (*mask & 1)
1862: {
1863: resbyte=residue%117; if (resbyte > 58) resbyte = 117 - resbyte;
1864: *mask &= (powersmod[resbyte] >> 12);
1865: if (DEBUGLEVEL >= 5)
1866: fprintferr("\t 117: %3ld (3rd %ld, 5th %ld, 7th %ld)\n",
1867: resbyte, *mask&1, (*mask>>1)&1, (*mask>>2)&1);
1868: if (!*mask) { avma=av; return 0; }
1869: }
1870: if (*mask & 3)
1871: {
1872: resbyte=residue%31; if (resbyte > 15) resbyte = 31 - resbyte;
1873: *mask &= (powersmod[resbyte] >> 15);
1874: if (DEBUGLEVEL >= 5)
1875: fprintferr("\t 31: %3ld (3rd %ld, 5th %ld, 7th %ld)\n",
1876: resbyte, *mask&1, (*mask>>1)&1, (*mask>>2)&1);
1877: if (!*mask) { avma=av; return 0; }
1878: }
1879: if (*mask & 5)
1880: {
1881: resbyte=residue%43; if (resbyte > 21) resbyte = 43 - resbyte;
1882: *mask &= (powersmod[resbyte] >> 18);
1883: if (DEBUGLEVEL >= 5)
1884: fprintferr("\t 43: %3ld (3rd %ld, 5th %ld, 7th %ld)\n",
1885: resbyte, *mask&1, (*mask>>1)&1, (*mask>>2)&1);
1886: if (!*mask) { avma=av; return 0; }
1887: }
1888: if (*mask & 6)
1889: {
1890: resbyte=residue%71; if (resbyte > 35) resbyte = 71 - resbyte;
1891: *mask &= (powersmod[resbyte] >> 21);
1892: if (DEBUGLEVEL >= 5)
1893: fprintferr("\t 71: %3ld (3rd %ld, 5th %ld, 7th %ld)\n",
1894: resbyte, *mask&1, (*mask>>1)&1, (*mask>>2)&1);
1895: if (!*mask) { avma=av; return 0; }
1896: }
1897:
1898: /* priority to higher powers -- if we have a 21st, it'll be easier to
1899: rediscover that its 7th root is a cube than that its cube root is
1900: a 7th power */
1901: if ((resbyte = *mask & 4)) /* assignment */
1902: exponent = 7;
1903: else if ((resbyte = *mask & 2))
1904: exponent = 5;
1905: else
1906: { resbyte = 1; exponent = 3; }
1907: /* leave that mask bit on for the moment, we might need it for a
1908: subsequent call */
1909:
1910: /* precision in the following is one extra significant word (overkill) */
1911: y=ground(gpow(x, ginv(stoi(exponent)), lgx));
1912: if (!egalii(gpowgs(y, exponent), x))
1913: {
1914: if (DEBUGLEVEL >= 5)
1915: {
1916: if (exponent == 3)
1917: fprintferr("\tBut it nevertheless wasn't a cube.\n");
1918: else
1919: fprintferr("\tBut it nevertheless wasn't a %ldth power.\n",
1920: exponent);
1921: }
1922: *mask &= ~resbyte; /* _now_ turn the bit off */
1923: avma=av; return 0;
1924: }
1925: /* caller (ifac_crack() below) will report the final result if it was
1926: a pure power, so no further diagnostics here */
1927:
1928: tetpil=avma;
1929: if (!pt) { avma=av; return exponent; } /* this branch not used */
1930: *pt=gerepile(av,tetpil,icopy(y));
1931: return exponent;
1932: }
1933:
1934: /***********************************************************************/
1935: /** **/
1936: /** FACTORIZATION (master iteration) **/
1937: /** Driver for the various methods of finding large factors **/
1938: /** (after trial division has cast out the very small ones). **/
1939: /** GN1998Jun24--30 **/
1940: /** **/
1941: /***********************************************************************/
1942:
1943: /** Direct use:
1944: ** ifac_start() registers a number (without prime factors < 100)
1945: ** with the iterative factorizer, and also registers whether or
1946: ** not we should terminate early if we find that the number is
1947: ** not squarefree, and a hint about which method(s) to use. This
1948: ** must always be called first. The input _must_ have been checked
1949: ** to be composite by the caller. The routine immediately tries
1950: ** to decompose it nontrivially into a product of two factors,
1951: ** except in squarefreeness (`Moebius') mode.
1952: ** ifac_primary_factor() returns a prime divisor (not necessarily
1953: ** the smallest) and the corresponding exponent. */
1954:
1955: /** Encapsulated user interface:
1956: ** ifac_decomp() does the right thing for auxdecomp() (put a succession
1957: ** of prime divisor / exponent pairs onto the stack, not necessarily
1958: ** sorted, although in practice they will tend not to be too far from
1959: ** the correct order).
1960: **
1961: ** For each of the additive/multiplicative arithmetic functions, there is
1962: ** a `contributor' below, to be called on any large composite cofactor
1963: ** left over after trial division by small primes, whose result can then
1964: ** be added to or multiplied with whatever we already have:
1965: ** ifac_moebius() ifac_issquarefree() ifac_totient() ifac_omega()
1966: ** ifac_bigomega() ifac_numdiv() ifac_sumdiv() ifac_sumdivk() */
1967:
1968: /* We never test whether the input number is prime or composite, since
1969: presumably it will have come out of the small factors finder stage
1970: (which doesn't really exist yet but which will test the left-over
1971: cofactor for primality once it does). */
1972:
1973: /* The data structure in which we preserve whatever we know at any given
1974: time about our number N is kept on the PARI stack, and updated as needed.
1975: This makes the machinery re-entrant (you can have more than one fac-
1976: torization using ifac_start()/ifac_primary_factor() in progress simul-
1977: taneously so long as you preserve the GEN across garbage collections),
1978: and which avoids memory leaks when a lengthy factorization is interrupted.
1979: We also make an effort to keep the whole affair connected, and the parent
1980: object will always be older than its children. This may in rare cases
1981: lead to some extra copying around, and knowing what is garbage at any
1982: given time is not entirely trivial. See below for examples how to do
1983: it right. (Connectedness can be destroyed if callers of ifac_main()
1984: create other stuff on the stack in between calls. This is harmless
1985: as long as ifac_realloc() is used to re-create a connected object at
1986: the head of the stack just before collecting garbage.) */
1987:
1988: /* Note that a PARI integer can have hundreds of millions of distinct prime
1989: factors larger than 2^16, given enough memory. And since there's no
1990: guarantee that we will find factors in order of increasing size, we must
1991: be prepared to drag a very large amount of data around (although this
1992: will _very_ rarely happen for random input!). So we start with a small
1993: structure and extend it when necessary. */
1994:
1995: /* The idea of data structure and algorithm is:
1996: Let N0 be whatever is currently left of N after dividing off all the
1997: prime powers we have already returned to the caller. Then we maintain
1998: N0 as a product
1999: (1) N0 = \prod_i P_i^{e_i} * \prod_j Q_j^{f_j} * \prod_k C_k^{g_k}
2000: where the P_i and Q_j are distinct primes, each C_k is known composite,
2001: none of the P_i divides any C_k, and we also know the total ordering
2002: of all the P_i, Q_j and C_k (in particular, we will never try to divide
2003: a C_k by a larger Q_j). Some of the C_k may have common factors, although
2004: this will not often be the case. */
2005:
2006: /* Caveat implementor: Taking gcds among C_k's is very likely to cost at
2007: least as much time as dividing off any primes as we find them, and book-
2008: keeping would be a nightmare (since D=gcd(C_1,C_2) can still have common
2009: factors with both C_1/D and C_2/D, and so on...). */
2010:
2011: /* At startup, we just initialize the structure to
2012: (2) N = C_1^1 (composite). */
2013:
2014: /* Whenever ifac_primary_factor() or ifac_decomp() (or, mutatis mutandis,
2015: one of the three arithmetic user interface routines) needs a primary
2016: factor, and the smallest thing in our list is P_1, we return that and
2017: its exponent, and remove it from our list.
2018: (When nothing is left, we return a sentinel value -- gun. And in Moebius
2019: mode, when we see something with exponent > 1, whether prime or composite,
2020: we yell at our caller by returning gzero or 0, depending on the function).
2021: In all other cases, ifac_main() iterates the following steps until we have
2022: a P_1 in the smallest position. */
2023:
2024: /* When the smallest item is C_1 (as it is initially):
2025: (3.1) Crack C_1 into a nontrivial product U_1 * U_2 by whatever method
2026: comes to mind for this size. (U for `unknown'.) Cracking will detect
2027: squares (and biquadrates etc), and it may detect odd powers, so we
2028: might instead see a power of some U_1 here, or even something of the form
2029: U_1^k*U_2^k. (Of course the exponent already attached to C_1 is taken
2030: into account in the following.)
2031: (3.2) If we have U_1*U_2, sort the two factors; convert to U_1^2 if they
2032: happen to be equal (which they shouldn't -- squares should have been
2033: caught at the preceding stage). Note that U_1 and (if it exists) U_2
2034: are automatically smaller than anything else in our list.
2035: (3.3) Check U_1 (and U_2) for primality, and flag them accordingly.
2036: (3.4) Iterate. */
2037:
2038: /* When the smallest item is Q_1:
2039: This is the potentially unpleasant case. The idea is to go through the
2040: entire list and try to divide Q_1 off each of the current C_k's, which
2041: will usually fail, but may succeed several times. When a division was
2042: successful, the corresponding C_k is removed from our list, and the co-
2043: factor becomes a U_l for the moment unless it is 1 (which happens when
2044: C_k was a power of Q_1). When we're through we upgrade Q_1 to P_1 status,
2045: and then do a primality check on each U_l and sort it back into the list
2046: either as a Q_j or as a C_k. If during the insertion sort we discover
2047: that some U_l equals some P_i or Q_j or C_k we already have, we just add
2048: U_l's exponent to that of its twin. (The sorting should therefore happen
2049: before the primality test).
2050: Note that this may produce one or more elements smaller than the P_1
2051: we just confirmed, so we may have to repeat the iteration. */
2052:
2053: /* There's a little trick that avoids some Q_1 instances. Just after we do
2054: a sweep to classify all current unknowns as either composites or primes,
2055: we do another downward sweep beginning with the largest current factor
2056: and stopping just above the largest current composite. Every Q_j we
2057: pass is turned into a P_i. (Different primes are automatically coprime
2058: among each other, and primes tend not to divide smaller composites.) */
2059:
2060: /* (We have no use for comparing the square of a prime to N0. Normally
2061: we will get called after casting out only the smallest primes, and
2062: since we cannot guarantee that we see the large prime factors in as-
2063: cending order, we cannot stop when we find one larger than sqrt(N0).) */
2064:
2065: /* Data structure: We keep everything in a single t_VEC of t_INTs. The
2066: first component records whether we're doing full (NULL) or Moebius (un)
2067: factorization; in the latter case many subroutines return a sentinel
2068: value as soon as they spot an exponent > 1. The second component records
2069: the hint from factorint()'s optional flag, for use by ifac_crack().
2070: The remaining components (initially 15) are used in groups of three:
2071: a GEN pointer at the t_INT value of the factor, a pointer at the t_INT
2072: exponent (usually gun or gdeux so we don't clutter up the stack too
2073: much), and another t_INT GEN pointer to record the class of the factor:
2074: NULL for unknown, zero for known composite C_k, un for known prime Q_j
2075: awaiting trial division, and deux for finished prime P_i. */
2076:
2077: /* When during the division stage we re-sort a C_k-turned-U_l to a lower
2078: position, we rotate any intervening material upward towards its old
2079: slot. When a C_k was divided down to 1, its slot is left empty at
2080: first; similarly when the re-sorting detects a repeated factor.
2081: After the sorting phase, we de-fragment the list and squeeze all the
2082: occupied slots together to the high end, so that ifac_crack() has room
2083: for new factors. When this doesn't suffice, we abandon the current
2084: vector and allocate a somewhat larger one, defragmenting again during
2085: copying. */
2086:
2087: /* (For internal use, note that all exponents will fit into C longs, given
2088: PARI's lgefint field size. When we work with them, we sometimes read
2089: out the GEN pointer, and sometimes do an itos, whatever is more con-
2090: venient for the task at hand.) */
2091:
2092:
2093: /*** Overview and forward declarations: ***/
2094:
2095: /* The `*where' argument in the following points into *partial at the
2096: first of the three fields of the first occupied slot. It's there
2097: because the caller would already know where `here' is, so we don't
2098: want to search for it again, although it wouldn't take much time.
2099: On the other hand, we do not preserve this from one user-interface
2100: call to the next. */
2101:
2102: static GEN
2103: ifac_find(GEN *partial, GEN *where);
2104: /* Return GEN pointing at the first nonempty slot strictly behind the
2105: current *where, or NULL if such doesn't exist. Can be used to skip
2106: a range of vacant slots, or to initialize *where in the first place
2107: (pass partial in both args). Does not modify its argument pointers. */
2108:
2109: void
2110: ifac_realloc(GEN *partial, GEN *where, long new_lg);
2111: /* Move to a larger main vector, updating *where if it points into it.
2112: Certainly updates *partial. Can be used as a specialized gcopy before
2113: a gerepileupto()/gerepilemanysp() (pass 0 as the new length).
2114: Normally, one would pass new_lg=1 to let this function guess the
2115: new size. To be used sparingly. */
2116:
2117: static long
2118: ifac_crack(GEN *partial, GEN *where);
2119: /* Split the first (composite) entry. There _must_ already be room for
2120: another factor below *where, and *where will be updated. Factor and
2121: cofactor will be inserted in the correct order, updating *where, or
2122: factor^k will be inserted if such should be the case (leaving *where
2123: unchanged). The factor or factors will be set to unknown, and inherit
2124: the exponent (or a multiple thereof) of its/their ancestor. Returns
2125: number of factors written into the structure (normally 2, but 1 if a
2126: factor equalled its cofactor, and may be more than 1 if a factoring
2127: engine returned a vector of factors instead of a single factor). Can
2128: reallocate the data structure in the vector-of-factors case (but not
2129: in the more common single-factor case) */
2130:
2131: static long
2132: ifac_insert_multiplet(GEN *partial, GEN *where, GEN facvec);
2133: /* Gets called to complete ifac_crack()'s job when a factoring engine
2134: splits the current factor into a product of three or more new factors.
2135: Makes room for them if necessary, sorts them, gives them the right
2136: exponents and class etc. Also returns the number of factors actually
2137: written, which may be less than the number of components in facvec
2138: if there are duplicates.--- Vectors of factors (cf pollardbrent()
2139: above) actually contain `slots' of three GENs per factor with the
2140: three fields being interpreted exactly as in our partial factorization
2141: data structure. Thus `engines' can tell us what they already happen to
2142: know about factors being prime or composite and/or appearing to a power
2143: larger than the first */
2144:
2145: static long
2146: ifac_divide(GEN *partial, GEN *where);
2147: /* Divide all current composites by first (prime, class Q) entry, updating
2148: its exponent, and turning it into a finished prime (class P). Return 1
2149: if any such divisions succeeded (in Moebius mode, the update may then
2150: not have been completed), or 0 if none of them succeeded. Doesn't
2151: modify *where. */
2152:
2153: static long
2154: ifac_sort_one(GEN *partial, GEN *where, GEN washere);
2155: /* re-sort one (typically unknown) entry from washere to a new position,
2156: rotating intervening entries upward to fill the vacant space. It may
2157: happen (rarely) that the new position is the same as the old one, or
2158: that the new value of the entry coincides with a value already occupying
2159: a lower slot, in which latter case we just add exponents (and use the
2160: `more known' class, and return 1 immediately when in Moebius mode).
2161: The slots between *where and washere must be in sorted order, so a
2162: sweep using this to re-sort several unknowns must proceed upward (see
2163: ifac_resort() below). Return 1 if we see an exponent > 1 (in Moebius
2164: mode without completing the update), 0 otherwise. */
2165:
2166: static long
2167: ifac_resort(GEN *partial, GEN *where);
2168: /* sort all current unknowns downward to where they belong. Sweeps
2169: in the upward direction. Not needed after ifac_crack(), only when
2170: ifac_divide() returned true. May update *where. Returns 1 when an
2171: ifac_sort_one() call does so to indicate a repeated factor, or 0 if
2172: any and all such calls returned 0 */
2173:
2174: static void
2175: ifac_defrag(GEN *partial, GEN *where);
2176: /* defragment: collect and squeeze out any unoccupied slots above *where
2177: during a downward sweep. Unoccupied slots arise when a composite factor
2178: dissolves completely whilst dividing off a prime, or when ifac_resort()
2179: spots a coincidence and merges two factors. *where will be updated */
2180:
2181: static void
2182: ifac_whoiswho(GEN *partial, GEN *where, long after_crack);
2183: /* determine primality or compositeness of all current unknowns, and set
2184: class Q primes to finished (class P) if everything larger is already
2185: known to be prime. When after_crack is nonnegative, only look at the
2186: first after_crack things in the list (do nothing when it's zero) */
2187:
2188: static GEN
2189: ifac_main(GEN *partial);
2190: /* main loop: iterate until smallest entry is a finished prime; returns
2191: a `where' pointer, or gun if nothing left, or gzero in Moebius mode if
2192: we aren't squarefree */
2193:
2194: /* NB In the most common cases, control flows from the user interface to
2195: ifac_main() and then to a succession of ifac_crack()s and ifac_divide()s,
2196: with (typically) none of the latter finding anything. */
2197:
2198: /** user interface: **/
2199: /* return initial data structure, see ifac_crack() below for semantics
2200: of the hint argument */
2201: GEN
2202: ifac_start(GEN n, long moebius, long hint);
2203:
2204: /* run main loop until primary factor is found, return the prime and
2205: assign the exponent. If nothing left, return gun and set exponent
2206: to 0; if in Moebius mode and a square factor is discovered, return
2207: gzero and set exponent to 0 */
2208: GEN
2209: ifac_primary_factor(GEN *partial, long *exponent);
2210:
2211: /* call ifac_start() and run main loop until factorization is complete,
2212: accumulating prime / exponent pairs on the PARI stack to be picked up
2213: by aux_end(). Return number of distinct primes found */
2214: long
2215: ifac_decomp(GEN n, long hint);
2216:
2217: /* completely encapsulated functions; these call ifac_start() themselves,
2218: and ensure proper stack housekeeping etc. Call them on any large
2219: composite left over after trial division, and multiply/add the result
2220: onto whatever you already have from the small factors. Don't call
2221: them on large primes; they will run into trouble */
2222: long
2223: ifac_moebius(GEN n, long hint);
2224:
2225: long
2226: ifac_issquarefree(GEN n, long hint);
2227:
2228: long
2229: ifac_omega(GEN n, long hint);
2230:
2231: long
2232: ifac_bigomega(GEN n, long hint);
2233:
2234: GEN
2235: ifac_totient(GEN n, long hint); /* for gp's eulerphi() */
2236:
2237: GEN
2238: ifac_numdiv(GEN n, long hint);
2239:
2240: GEN
2241: ifac_sumdiv(GEN n, long hint);
2242:
2243: GEN
2244: ifac_sumdivk(GEN n, long k, long hint);
2245:
2246: /*** implementation ***/
2247:
2248: #define ifac_initial_length 24 /* codeword, moebius flag, hint, 7 slots */
2249: /* (more than enough in most cases -- a 512-bit product of distinct 8-bit
2250: primes needs at most 7 slots at a time) */
2251:
2252: GEN
2253: ifac_start(GEN n, long moebius, long hint)
2254: {
2255: GEN part, here;
2256:
2257: if (typ(n) != t_INT) err(typeer, "ifac_start");
2258: if (signe(n) == 0)
2259: err(talker, "factoring 0 in ifac_start");
2260:
2261: part = cgetg(ifac_initial_length, t_VEC);
2262: here = part + ifac_initial_length;
2263: part[1] = moebius? un : (long)NULL;
2264: switch(hint)
2265: {
2266: case 0:
2267: part[2] = zero; break;
2268: case 1:
2269: part[2] = un; break;
2270: case 2:
2271: part[2] = deux; break;
2272: default:
2273: part[2] = (long)stoi(hint);
2274: }
2275: if (isonstack(n))
2276: n = absi(n);
2277: /* make copy, because we'll later want to mpdivis() into it in place.
2278: If it's not on stack, then we assume it is a clone made for us by
2279: auxdecomp0(), and we assume the sign has already been set positive */
2280: /* fill first slot at the top end */
2281: *--here = zero; /* initially composite */
2282: *--here = un; /* initial exponent 1 */
2283: *--here = (long) n;
2284: /* and NULL out the remaining slots */
2285: while (here > part + 3) *--here = (long)NULL;
2286: return part;
2287: }
2288:
2289: static GEN
2290: ifac_find(GEN *partial, GEN *where)
2291: {
2292: long lgp = lg(*partial);
2293: GEN end = *partial + lgp;
2294: GEN scan = *where + 3;
2295:
2296: if (DEBUGLEVEL >= 5)
2297: {
2298: if (!*partial || typ(*partial) != t_VEC)
2299: err(typeer, "ifac_find");
2300: if (lg(*partial) < ifac_initial_length)
2301: err(talker, "partial impossibly short in ifac_find");
2302: if (!(*where) ||
2303: *where > *partial + lgp - 3 ||
2304: *where < *partial) /* sic */
2305: err(talker, "`*where\' out of bounds in ifac_find");
2306: }
2307: while (scan < end && !*scan) scan += 3;
2308: /* paranoia -- check completely NULLed ? nope -- we never inspect the
2309: exponent field for deciding whether a slot is empty or occupied */
2310: if (scan < end)
2311: {
2312: if (DEBUGLEVEL >= 5)
2313: {
2314: if (!scan[1])
2315: err(talker, "factor has NULL exponent in ifac_find");
2316: }
2317: return scan;
2318: }
2319: return NULL;
2320: }
2321:
2322: /* simple defragmenter */
2323: static void
2324: ifac_defrag(GEN *partial, GEN *where)
2325: {
2326: long lgp = lg(*partial);
2327: GEN scan_new = *partial + lgp - 3, scan_old = scan_new;
2328:
2329: while (scan_old >= *where)
2330: {
2331: if (*scan_old) /* slot occupied? */
2332: {
2333: if (scan_old < scan_new)
2334: {
2335: scan_new[2] = scan_old[2];
2336: scan_new[1] = scan_old[1];
2337: *scan_new = *scan_old;
2338: }
2339: scan_new -= 3; /* point at next slot to be written */
2340: }
2341: scan_old -= 3;
2342: }
2343: scan_new += 3; /* back up to last slot written */
2344: *where = scan_new;
2345: while (scan_new > *partial + 3)
2346: *--scan_new = (long)NULL; /* erase junk */
2347: }
2348:
2349: /* and complex version combined with reallocation. If new_lg is 0, we
2350: use the old length, so this acts just like gcopy except that the where
2351: pointer is carried along; if it is 1, we make an educated guess.
2352: Exception: If new_lg is 0, the vector is full to the brim, and the
2353: first entry is composite, we make it longer to avoid being called again
2354: a microsecond later (at significant cost).
2355: It is safe to call this with NULL for the where argument; if it doesn't
2356: point anywhere within the old structure, it will be left alone */
2357: void
2358: ifac_realloc(GEN *partial, GEN *where, long new_lg)
2359: {
2360: long old_lg = lg(*partial);
2361: GEN newpart, scan_new, scan_old;
2362:
2363: if (DEBUGLEVEL >= 5) /* none of these should ever happen */
2364: {
2365: if (!*partial || typ(*partial) != t_VEC)
2366: err(typeer, "ifac_realloc");
2367: if (lg(*partial) < ifac_initial_length)
2368: err(talker, "partial impossibly short in ifac_realloc");
2369: }
2370:
2371: if (new_lg == 1)
2372: new_lg = 2*old_lg - 6; /* from 7 slots to 13 to 25... */
2373: else if (new_lg <= old_lg) /* includes case new_lg == 0 */
2374: {
2375: new_lg = old_lg;
2376: if ((*partial)[3] && /* structure full */
2377: ((*partial)[5]==zero || (*partial)[5]==(long)NULL))
2378: /* and first entry composite or unknown */
2379: new_lg += 6; /* give it a little more breathing space */
2380: }
2381: newpart = cgetg(new_lg, t_VEC);
2382: if (DEBUGMEM >= 3)
2383: {
2384: fprintferr("IFAC: new partial factorization structure (%ld slots)\n",
2385: (new_lg - 3)/3);
2386: flusherr();
2387: }
2388: newpart[1] = (*partial)[1]; /* moebius */
2389: newpart[2] = (*partial)[2]; /* hint */
2390: /* downward sweep through the old *partial, picking up where1 and carry-
2391: ing it over if and when we pass it. (This will only be useful if
2392: it pointed at a non-empty slot.) Factors are licopy()d so that we
2393: again have a nice object (parent older than children, connected),
2394: except the one factor that may still be living in a clone where n
2395: originally was; exponents are similarly copied if they aren't global
2396: constants; class-of-factor fields are always global constants so we
2397: need only copy them as pointers. Caller may then do a gerepileupto()
2398: or a gerepilemanysp() */
2399: scan_new = newpart + new_lg - 3;
2400: scan_old = *partial + old_lg - 3;
2401: for (; scan_old > *partial + 2; scan_old -= 3)
2402: {
2403: if (*where == scan_old) *where = scan_new;
2404: if (!*scan_old) continue; /* skip empty slots */
2405:
2406: *scan_new =
2407: isonstack((GEN)(*scan_old)) ?
2408: licopy((GEN)(*scan_old)) : *scan_old;
2409: scan_new[1] =
2410: isonstack((GEN)(scan_old[1])) ?
2411: licopy((GEN)(scan_old[1])) : scan_old[1];
2412: scan_new[2] = scan_old[2];
2413: scan_new -= 3;
2414: }
2415: scan_new += 3; /* back up to last slot written */
2416: while (scan_new > newpart + 3)
2417: *--scan_new = (long)NULL;
2418: *partial = newpart;
2419: }
2420:
2421: #define moebius_mode ((*partial)[1])
2422:
2423: /* Bubble-sort-of-thing sort. Won't be exercised frequently,
2424: so this is ok. */
2425: static long
2426: ifac_sort_one(GEN *partial, GEN *where, GEN washere)
2427: {
2428: GEN scan = washere - 3;
2429: GEN value, exponent, class0, class1;
2430: long cmp_res;
2431:
2432: if (DEBUGLEVEL >= 5) /* none of these should ever happen */
2433: {
2434: long lgp;
2435: if (!*partial || typ(*partial) != t_VEC)
2436: err(typeer, "ifac_sort_one");
2437: if ((lgp = lg(*partial)) < ifac_initial_length)
2438: err(talker, "partial impossibly short in ifac_sort_one");
2439: if (!(*where) ||
2440: *where < *partial + 3 ||
2441: *where > *partial + lgp - 3)
2442: err(talker, "`*where\' out of bounds in ifac_sort_one");
2443: if (!washere ||
2444: washere < *where ||
2445: washere > *partial + lgp - 3)
2446: err(talker, "`washere\' out of bounds in ifac_sort_one");
2447: }
2448: value = (GEN)(*washere);
2449: exponent = (GEN)(washere[1]);
2450: if (exponent != gun && moebius_mode && cmpsi(1,exponent) < 0)
2451: return 1; /* should have been detected by caller */
2452: class0 = (GEN)(washere[2]);
2453:
2454: if (scan < *where) return 0; /* nothing to do, washere==*where */
2455:
2456: cmp_res = -1; /* sentinel */
2457: while (scan >= *where) /* therefore at least once */
2458: {
2459: if (*scan) /* current slot nonempty */
2460: {
2461: /* check against where */
2462: cmp_res = cmpii(value, (GEN)(*scan));
2463: if (cmp_res >= 0) break; /* have found where to stop */
2464: }
2465: /* copy current slot upward by one position and move pointers down */
2466: scan[5] = scan[2];
2467: scan[4] = scan[1];
2468: scan[3] = *scan;
2469: scan -= 3;
2470: }
2471: scan += 3;
2472: /* at this point there are the following possibilities:
2473: (*) cmp_res == -1. Either value is less than that at *where, or for
2474: some reason *where was pointing at one or more vacant slots and any
2475: factors we saw en route were larger than value. At any rate,
2476: scan == *where now, and scan is pointing at an empty slot, into
2477: which we'll stash our entry.
2478: (*) cmp_res == 0. The entry at scan-3 is the one, we compare class0
2479: fields and add exponents, and put it all into the vacated scan slot,
2480: NULLing the one at scan-3 (and possibly updating *where).
2481: (*) cmp_res == 1. The slot at scan is the one to store our entry
2482: into. */
2483: if (cmp_res != 0)
2484: {
2485: if (cmp_res < 0 && scan != *where)
2486: err(talker, "misaligned partial detected in ifac_sort_one");
2487: *scan = (long)value;
2488: scan[1] = (long)exponent;
2489: scan[2] = (long)class0;
2490: return 0;
2491: }
2492: /* case cmp_res == 0: repeated factor detected */
2493: if (DEBUGLEVEL >= 4)
2494: {
2495: fprintferr("IFAC: repeated factor %Z\n\tdetected in ifac_sort_one\n",
2496: value);
2497: flusherr();
2498: }
2499: if (moebius_mode) return 1; /* not squarefree */
2500: /* if old class0 was composite and new is prime, or vice versa,
2501: complain (and if one class0 was unknown and the other wasn't,
2502: use the known one) */
2503: class1 = (GEN)(scan[-1]);
2504: if (class0) /* should never be used */
2505: {
2506: if(class1)
2507: {
2508: if (class0 == gzero && class1 != gzero)
2509: err(talker, "composite equals prime in ifac_sort_one");
2510: else if (class0 != gzero && class1 == gzero)
2511: err(talker, "prime equals composite in ifac_sort_one");
2512: else if (class0 == gdeux) /* should happen even less */
2513: scan[2] = (long)class0; /* use it */
2514: }
2515: else /* shouldn't happen either */
2516: scan[2] = (long)class0; /* use it */
2517: }
2518: /* else stay with the existing known class0 */
2519: scan[2] = (long)class1;
2520: /* in any case, add exponents */
2521: if (scan[-2] == un && exponent == gun)
2522: scan[1] = deux;
2523: else
2524: scan[1] = laddii((GEN)(scan[-2]), exponent);
2525: /* move the value over */
2526: *scan = scan[-3];
2527: /* null out the vacated slot below */
2528: *--scan = (long)NULL;
2529: *--scan = (long)NULL;
2530: *--scan = (long)NULL;
2531: /* finally, see whether *where should be pulled in */
2532: if (scan == *where) *where += 3;
2533: return 0;
2534: }
2535:
2536: /* the following loop around the former doesn't need to check moebius_mode
2537: because ifac_sort_one() never returns 1 in normal mode */
2538: static long
2539: ifac_resort(GEN *partial, GEN *where)
2540: {
2541: long lgp = lg(*partial), res = 0;
2542: GEN scan = *where;
2543:
2544: for (; scan < *partial + lgp; scan += 3)
2545: {
2546: if (*scan && /* slot occupied */
2547: !scan[2]) /* with an unknown */
2548: {
2549: res |= ifac_sort_one(partial, where, scan);
2550: if (res) return res; /* early exit */
2551: }
2552: }
2553: return res;
2554: }
2555:
2556: /* sweep downward so we can with luck turn some Qs into Ps */
2557: static void
2558: ifac_whoiswho(GEN *partial, GEN *where, long after_crack)
2559: {
2560: long lgp = lg(*partial), larger_compos = 0;
2561: GEN scan, scan_end = *partial + lgp - 3;
2562:
2563: if (DEBUGLEVEL >= 5)
2564: {
2565: if (!*partial || typ(*partial) != t_VEC)
2566: err(typeer, "ifac_whoiswho");
2567: if (lg(*partial) < ifac_initial_length)
2568: err(talker, "partial impossibly short in ifac_whoiswho");
2569: if (!(*where) ||
2570: *where > scan_end ||
2571: *where < *partial + 3)
2572: err(talker, "`*where\' out of bounds in ifac_whoiswho");
2573: }
2574:
2575: if (after_crack == 0) return;
2576: if (after_crack > 0)
2577: {
2578: larger_compos = 1; /* disable Q-to-P trick */
2579: scan = *where + 3*(after_crack - 1);
2580: /* check at most after_crack entries */
2581: if (scan > scan_end) /* ooops... */
2582: {
2583: err(warner, "avoiding nonexistent factors in ifac_whoiswho");
2584: scan = scan_end;
2585: }
2586: }
2587: else { larger_compos = 0; scan = scan_end; }
2588:
2589: for (; scan >= *where; scan -= 3)
2590: {
2591: if (scan[2]) /* known class of factor */
2592: {
2593: if (scan[2] == zero) larger_compos = 1;
2594: else if (!larger_compos && scan[2] == un)
2595: {
2596: if (DEBUGLEVEL >= 3)
2597: {
2598: fprintferr("IFAC: factor %Z\n\tis prime (no larger composite)\n",
2599: **where);
2600: fprintferr("IFAC: prime %Z\n\tappears with exponent = %ld\n",
2601: **where, itos((GEN)(*where)[1]));
2602: }
2603: scan[2] = deux;
2604: } /* no else case */
2605: continue;
2606: }
2607: scan[2] =
2608: (isprime((GEN)(*scan)) ?
2609: (larger_compos ? un : deux) : /* un- or finished prime */
2610: zero); /* composite */
2611:
2612: if (scan[2] == zero) larger_compos = 1;
2613: if (DEBUGLEVEL >= 3)
2614: {
2615: fprintferr("IFAC: factor %Z\n\tis %s\n", *scan,
2616: (scan[2] == zero ? "composite" : "prime"));
2617: }
2618: }
2619: }
2620:
2621: /* Here we normally do not check that the first entry is a not-finished
2622: prime. Stack management: we may allocate a new exponent */
2623: static long
2624: ifac_divide(GEN *partial, GEN *where)
2625: {
2626: long lgp = lg(*partial);
2627: GEN scan = *where + 3;
2628: long res = 0, exponent, newexp, otherexp;
2629:
2630: if (DEBUGLEVEL >= 5) /* none of these should ever happen */
2631: {
2632: if (!*partial || typ(*partial) != t_VEC)
2633: err(typeer, "ifac_divide");
2634: if (lg(*partial) < ifac_initial_length)
2635: err(talker, "partial impossibly short in ifac_divide");
2636: if (!(*where) ||
2637: *where > *partial + lgp - 3 ||
2638: *where < *partial + 3)
2639: err(talker, "`*where\' out of bounds in ifac_divide");
2640: if ((*where)[2] != un)
2641: err(talker, "division by composite or finished prime in ifac_divide");
2642: }
2643: if (!(**where)) /* always test just this one */
2644: err(talker, "division by nothing in ifac_divide");
2645:
2646: newexp = exponent = itos((GEN)((*where)[1]));
2647: if (exponent > 1 && moebius_mode) return 1;
2648: /* should've been caught by caller already */
2649:
2650: /* go for it */
2651: for (; scan < *partial + lgp; scan += 3)
2652: {
2653: if (scan[2] != zero) continue; /* the other thing ain't composite */
2654: otherexp = 0;
2655: /* let mpdivis divide the other factor in place to keep stack clutter
2656: minimal */
2657: while (mpdivis((GEN)(*scan), (GEN)(**where), (GEN)(*scan)))
2658: {
2659: if (moebius_mode) return 1; /* immediately */
2660: if (!otherexp) otherexp = itos((GEN)(scan[1]));
2661: newexp += otherexp;
2662: }
2663: if (newexp > exponent) /* did anything happen? */
2664: {
2665: (*where)[1] = (newexp == 2 ? deux : (long)(stoi(newexp)));
2666: exponent = newexp;
2667: if (is_pm1((GEN)(*scan))) /* factor dissolved completely */
2668: {
2669: *scan = scan[1] = (long)NULL;
2670: if (DEBUGLEVEL >= 4)
2671: fprintferr("IFAC: a factor was a power of another prime factor\n");
2672: }
2673: else if (DEBUGLEVEL >= 4)
2674: {
2675: fprintferr("IFAC: a factor was divisible by another prime factor,\n");
2676: fprintferr("\tleaving a cofactor = %Z\n", *scan);
2677: }
2678: scan[2] = (long)NULL; /* at any rate it's Unknown now */
2679: res = 1;
2680: if (DEBUGLEVEL >= 5)
2681: {
2682: fprintferr("IFAC: prime %Z\n\tappears at least to the power %ld\n",
2683: **where, newexp);
2684: }
2685: }
2686: } /* for */
2687: (*where)[2] = deux; /* make it a finished prime */
2688: if (DEBUGLEVEL >= 3)
2689: {
2690: fprintferr("IFAC: prime %Z\n\tappears with exponent = %ld\n",
2691: **where, newexp);
2692: }
2693: return res;
2694: }
2695:
2696:
2697: GEN mpqs(GEN N); /* in src/modules/mpqs.c, maybe a dummy,
2698: returns a factor, or a vector of factors,
2699: or NULL */
2700:
2701: /* The following takes the place of 2.0.9.alpha's find_factor(). */
2702:
2703: /* The meaning of the hint changes against 2.0.9.alpha to:
2704: hint == 0 : Use our own strategy, and this should be the default
2705: hint & 1 : Avoid mpqs(), use ellfacteur() after pollardbrent()
2706: hint & 2 : Avoid first-stage ellfacteur() in favour of mpqs()
2707: (which may still fall back to ellfacteur() if mpqs() is not installed
2708: or gives up)
2709: hint & 4 : Avoid even the pollardbrent() stage
2710: hint & 8 : Avoid final ellfacteur(); this may `declare' a composite
2711: to be prime. */
2712:
2713: /* stack housekeeping: this routine may create one or more objects (a new
2714: factor, or possibly several, and perhaps one or more new exponents > 2) */
2715: static long
2716: ifac_crack(GEN *partial, GEN *where)
2717: {
2718: long hint, cmp_res, exp1 = 1, exp2 = 1, av;
2719: GEN factor = NULL, exponent;
2720:
2721: if (DEBUGLEVEL >= 5) /* none of these should ever happen */
2722: {
2723: long lgp;
2724: if (!*partial || typ(*partial) != t_VEC)
2725: err(typeer, "ifac_crack");
2726: if ((lgp = lg(*partial)) < ifac_initial_length)
2727: err(talker, "partial impossibly short in ifac_crack");
2728: if (!(*where) ||
2729: *where < *partial + 6 || /* sic -- caller must realloc first */
2730: *where > *partial + lgp - 3)
2731: err(talker, "`*where\' out of bounds in ifac_crack");
2732: if (!(**where) || typ((GEN)(**where)) != t_INT)
2733: err(typeer, "ifac_crack");
2734: if ((*where)[2] != zero)
2735: err(talker, "operand not known composite in ifac_crack");
2736: }
2737: hint = itos((GEN)((*partial)[2])) & 15;
2738: exponent = (GEN)((*where)[1]);
2739:
2740: if (DEBUGLEVEL >= 3)
2741: fprintferr("IFAC: cracking composite\n\t%Z\n", **where);
2742:
2743: /* crack squares. Quite fast due to the initial square residue test */
2744: if (DEBUGLEVEL >= 4)
2745: fprintferr("IFAC: checking for pure square\n");
2746: av = avma;
2747: while (carrecomplet((GEN)(**where), &factor))
2748: {
2749: if (DEBUGLEVEL >= 4)
2750: fprintferr("IFAC: found %Z =\n\t%Z ^2\n", **where, factor);
2751: affii(factor, (GEN)(**where)); avma = av; factor = NULL;
2752: if (exponent == gun)
2753: (*where)[1] = deux;
2754: else if (exponent == gdeux)
2755: { (*where)[1] = (long)stoi(4); av = avma; }
2756: else
2757: { affii(shifti(exponent, 1), (GEN)((*where)[1])); avma = av; }
2758: exponent = (GEN)((*where)[1]);
2759: if (moebius_mode) return 0; /* no need to carry on... */
2760: exp1 = 2;
2761: } /* while carrecomplet */
2762:
2763: /* check whether our composite hasn't become prime */
2764: if (exp1 > 1 && isprime((GEN)(**where)))
2765: {
2766: (*where)[2] = un;
2767: if (DEBUGLEVEL >= 4)
2768: {
2769: fprintferr("IFAC: factor %Z\n\tis prime\n",**where);
2770: flusherr();
2771: }
2772: return 0; /* bypass subsequent ifac_whoiswho() call */
2773: }
2774: /* still composite -- carry on */
2775:
2776: /* MPQS cannot factor prime powers; check for cubes/5th/7th powers.
2777: Do this even if MPQS is blocked by hint -- it still serves a useful
2778: purpose in bounded factorization */
2779: {
2780: long mask = 7;
2781: if (DEBUGLEVEL == 4)
2782: fprintferr("IFAC: checking for odd power\n");
2783: /* (At debug levels > 4, is_odd_power() itself prints something more
2784: informative) */
2785: av = avma;
2786: while ((exp1 = /* assignment */
2787: is_odd_power((GEN)(**where), &factor, &mask)))
2788: {
2789: if (exp2 == 1) exp2 = exp1; /* remember this after the loop */
2790: if (DEBUGLEVEL >= 4)
2791: fprintferr("IFAC: found %Z =\n\t%Z ^%ld\n", **where, factor, exp1);
2792: affii(factor, (GEN)(**where)); avma = av; factor = NULL;
2793: if (exponent == gun)
2794: { (*where)[1] = (long)stoi(exp1); av = avma; }
2795: else if (exponent == gdeux)
2796: { (*where)[1] = (long)stoi(exp1<<1); av = avma; }
2797: else
2798: { affii(mulsi(exp1, exponent), (GEN)((*where)[1])); avma = av; }
2799: exponent = (GEN)((*where)[1]);
2800: if (moebius_mode) return 0; /* no need to carry on... */
2801: } /* while is_odd_power */
2802:
2803: if (exp2 > 1)
2804: { /* Something nice has happened */
2805: /* check whether our composite hasn't become prime */
2806: if (isprime((GEN)(**where)))
2807: {
2808: (*where)[2] = un;
2809: if (DEBUGLEVEL >= 4)
2810: {
2811: fprintferr("IFAC: factor %Z\n\tis prime\n", **where);
2812: flusherr();
2813: }
2814: return 0; /* bypass subsequent ifac_whoiswho() call */
2815: }
2816: /* base of power is still composite (an exceedingly rare case),
2817: fall through */
2818: }
2819: } /* odd power stage */
2820:
2821: /* pollardbrent() Rho usually gets a first chance */
2822: if (!(hint & 4))
2823: {
2824: if (DEBUGLEVEL >= 4)
2825: fprintferr("IFAC: trying Pollard-Brent rho method first\n");
2826: factor = pollardbrent((GEN)(**where));
2827: } /* Rho stage */
2828:
2829: /* if this didn't work, try one of our high-power beasties */
2830: if (!factor && !(hint & 2))
2831: {
2832: if (DEBUGLEVEL >= 4)
2833: fprintferr("IFAC: trying Lenstra-Montgomery ECM\n");
2834: factor = ellfacteur((GEN)(**where), 0); /* do not insist */
2835: } /* First ECM stage */
2836:
2837: if (!factor && !(hint & 1))
2838: {
2839: if (DEBUGLEVEL >= 4)
2840: fprintferr("IFAC: trying Multi-Polynomial Quadratic Sieve\n");
2841: factor = mpqs((GEN)(**where));
2842: } /* MPQS stage */
2843:
2844: if (!factor)
2845: {
2846: if (!(hint & 8)) /* still no luck? force it */
2847: {
2848: if (DEBUGLEVEL >= 4)
2849: fprintferr("IFAC: forcing ECM, may take some time\n");
2850: factor = ellfacteur((GEN)(**where), 1);
2851: } /* final ECM stage, guaranteed to succeed */
2852: else /* limited factorization */
2853: {
2854: if (DEBUGLEVEL >= 2)
2855: {
2856: err(warner, "IFAC: unfactored composite declared prime");
2857: /* don't print it out at level 3 or above, where it would appear
2858: several times before and after this message already */
2859: if (DEBUGLEVEL == 2)
2860: {
2861: fprintferr("\t%Z\n",**where);
2862: flusherr();
2863: }
2864: }
2865: (*where)[2] = un; /* might as well trial-divide by it... */
2866: return 1;
2867: }
2868: } /* Final ECM stage */
2869:
2870: if (DEBUGLEVEL >= 1)
2871: {
2872: if (!factor) /* never reached */
2873: err(talker, "all available factoring methods failed in ifac_crack");
2874: }
2875: if (typ(factor) == t_VEC) /* delegate this case */
2876: return ifac_insert_multiplet(partial, where, factor);
2877:
2878: else if (typ(factor) != t_INT)
2879: {
2880: fprintferr("IFAC: factorizer returned strange object to ifac_crack\n");
2881: outerr(factor);
2882: err(bugparier, "factoring");
2883: }
2884:
2885: /* got single integer back: work out the cofactor (in place) */
2886: if (!mpdivis((GEN)(**where), factor, (GEN)(**where)))
2887: {
2888: fprintferr("IFAC: factoring %Z\n", **where);
2889: fprintferr("\tyielded `factor\' %Z\n\twhich isn't!\n", factor);
2890: err(bugparier, "factoring");
2891: }
2892:
2893: /* the factoring engines report the factor found when DEBUGLEVEL is
2894: large enough; let's tell about the cofactor */
2895: if (DEBUGLEVEL >= 4)
2896: fprintferr("IFAC: cofactor = %Z\n", **where);
2897:
2898: /* ok, now `factor' is one factor and **where is the other, find out
2899: which is larger */
2900: cmp_res = cmpii(factor, (GEN)(**where));
2901: if (cmp_res < 0) /* common case */
2902: {
2903: (*where)[2] = (long)NULL; /* mark cofactor `unknown' */
2904: (*where)[-1] = (long)NULL; /* mark factor `unknown' */
2905: (*where)[-2] =
2906: isonstack(exponent) ? licopy(exponent) : (long)exponent;
2907: *where -= 3;
2908: **where = (long)factor;
2909: return 2;
2910: }
2911: else if (cmp_res == 0) /* hep, split a square in the middle */
2912: {
2913: err(warner,
2914: "square not found by carrecomplet, ifac_crack recovering");
2915: cgiv(factor);
2916: (*where)[2] = (long)NULL; /* mark the sqrt `unknown' */
2917: if (exponent == gun) /* double the exponent */
2918: (*where)[1] = deux;
2919: else if (exponent == gdeux)
2920: (*where)[1] = (long)stoi(4); /* make a new one */
2921: else /* overwrite old exponent */
2922: {
2923: av = avma;
2924: affii(shifti(exponent, 1), (GEN)((*where)[1]));
2925: avma = av;
2926: /* leave *where unchanged */
2927: }
2928: if (moebius_mode) return 0;
2929: return 1;
2930: }
2931: else /* factor > cofactor, rearrange */
2932: {
2933: (*where)[2] = (long)NULL; /* mark factor `unknown' */
2934: (*where)[-1] = (long)NULL; /* mark cofactor `unknown' */
2935: (*where)[-2] =
2936: isonstack(exponent) ? licopy(exponent) : (long)exponent;
2937: *where -= 3;
2938: **where = (*where)[3]; /* move cofactor pointer to lowest slot */
2939: (*where)[3] = (long)factor; /* save factor */
2940: return 2;
2941: }
2942: }
2943:
2944: /* the following doesn't collect garbage; caller's caller should do it
2945: (which means ifac_main()). No diagnostics either, the factoring engine
2946: should have printed what it found when DEBUGLEVEL>=4 or so. Note facvec
2947: contains slots of three components per factor; repeated factors are
2948: expressly allowed (and their classes shouldn't contradict each other
2949: whereas their exponents will be added up) */
2950: static long
2951: ifac_insert_multiplet(GEN *partial, GEN *where, GEN facvec)
2952: {
2953: long j,k=1, lfv=lg(facvec)-1, nf=lfv/3, room=(long)(*where-*partial);
2954: /* one of the factors will go into the *where slot, so room is now
2955: 3 times the number of slots we can use */
2956: long needroom = lfv - room;
2957: GEN sorted, auxvec = cgetg(nf+1, t_VEC), factor;
2958: long exponent = itos((GEN)((*where)[1])); /* the old exponent */
2959: GEN newexp;
2960:
2961: if (DEBUGLEVEL >= 5)
2962: fprintferr("IFAC: incorporating set of %ld factors%s\n",
2963: nf, (DEBUGLEVEL >=6 ? "..." : ""));
2964: if (needroom > 0)
2965: ifac_realloc(partial, where, lg(*partial) + needroom + 3);
2966: /* one extra slot for paranoia, errm, future use */
2967:
2968: /* create sort permutation from the values of the factors */
2969: for (j=nf; j; j--) auxvec[j] = facvec[3*j-2]; /* just the pointers */
2970: sorted = sindexsort(auxvec);
2971: /* and readjust the result for the triple spacing */
2972: for (j=nf; j; j--) sorted[j] = 3*sorted[j]-2;
2973: if (DEBUGLEVEL >= 6)
2974: fprintferr("\tsorted them...\n");
2975:
2976: /* store factors, beginning at *where, and catching any duplicates */
2977: **where = facvec[sorted[nf]];
2978: if ((newexp = (GEN)(facvec[sorted[nf]+1])) != gun) /* new exponent > 1 */
2979: {
2980: if (exponent == 1)
2981: (*where)[1] = isonstack(newexp) ? licopy(newexp) : (long)newexp;
2982: else
2983: (*where)[1] = lmulsi(exponent, newexp);
2984: } /* if new exponent is 1, the old exponent already in place will do */
2985: (*where)[2] = facvec[sorted[nf]+2]; /* copy class */
2986: if (DEBUGLEVEL >= 6)
2987: fprintferr("\tstored (largest) factor no. %ld...\n", nf);
2988:
2989: for (j=nf-1; j; j--)
2990: {
2991: factor = (GEN)(facvec[sorted[j]]);
2992: if (egalii(factor, (GEN)(**where)))
2993: {
2994: if (DEBUGLEVEL >= 6)
2995: fprintferr("\tfactor no. %ld is a duplicate%s\n",
2996: j, (j>1 ? "..." : ""));
2997: /* update exponent, ignore class which would already have been set,
2998: and then forget current factor */
2999: if ((newexp = (GEN)(facvec[sorted[j]+1])) != gun) /* new exp > 1 */
3000: { /* now we have at least 3 */
3001: (*where)[1] = laddii((GEN)((*where)[1]),
3002: mulsi(exponent, newexp));
3003: }
3004: else
3005: {
3006: if ((*where)[1] == un && exponent == 1)
3007: (*where)[1] = deux;
3008: else
3009: (*where)[1] = laddsi(exponent, (GEN)((*where)[1]));
3010: /* not safe to add 1 in place -- that might overwrite gdeux,
3011: with `interesting' consequences */
3012: }
3013: if (moebius_mode) return 0; /* stop now, but with exponent updated */
3014: continue;
3015: }
3016: (*where)[-1] = facvec[sorted[j]+2]; /* class as given */
3017: if ((newexp = (GEN)(facvec[sorted[j]+1])) != gun) /* new exp > 1 */
3018: {
3019: if (exponent == 1 && newexp == gdeux)
3020: (*where)[-2] = deux;
3021: else /* exponent*newexp > 2 */
3022: (*where)[-2] = lmulsi(exponent, newexp);
3023: }
3024: else
3025: {
3026: (*where)[-2] = (exponent == 1 ? un :
3027: (exponent == 2 ? deux :
3028: (long)stoi(exponent))); /* inherit parent's exponent */
3029: }
3030: (*where)[-3] = isonstack(factor) ? licopy(factor) : (long)factor;
3031: /* keep components younger than *partial */
3032: *where -= 3;
3033: k++;
3034: if (DEBUGLEVEL >= 6)
3035: fprintferr("\tfactor no. %ld was unique%s\n",
3036: j, (j>1 ? " (so far)..." : ""));
3037: }
3038: /* make the `sorted' object safe for garbage collection (probably not a
3039: problem, since it should be in the garbage zone from everybody's
3040: perspective, but it's easy to do it) */
3041: *sorted = evaltyp(t_INT) | evallg(nf+1);
3042: return k;
3043: }
3044:
3045: static GEN
3046: ifac_main(GEN *partial)
3047: {
3048: /* leave the basic error checking to ifac_find() */
3049: GEN here = ifac_find(partial, partial);
3050: long res, nf;
3051:
3052: /* if nothing left, return gun */
3053: if (!here) return gun;
3054:
3055: /* if we are in Moebius mode and have already detected a repeated factor,
3056: stop right here. Shouldn't really happen */
3057: if (moebius_mode && here[1] != un)
3058: {
3059: if (DEBUGLEVEL >= 3)
3060: {
3061: fprintferr("IFAC: main loop: repeated old factor\n\t%Z\n", *here);
3062: flusherr();
3063: }
3064: return gzero;
3065: }
3066:
3067: /* loop until first entry is a finished prime. May involve reallocations
3068: and thus updates of *partial */
3069: while (here[2] != deux)
3070: {
3071: /* if it's unknown, something has gone wrong; try to recover */
3072: if (!(here[2]))
3073: {
3074: err(warner, "IFAC: unknown factor seen in main loop");
3075: res = ifac_resort(partial, &here);
3076: if (res) return gzero; /* can only happen in Moebius mode */
3077: ifac_whoiswho(partial, &here, -1);
3078: /* defrag for good measure */
3079: ifac_defrag(partial, &here);
3080: continue;
3081: }
3082: /* if it's composite, crack it */
3083: if (here[2] == zero)
3084: {
3085: /* make sure there's room for another factor */
3086: if (here < *partial + 6)
3087: { /* try defrag first */
3088: ifac_defrag(partial, &here);
3089: if (here < *partial + 6) /* no luck */
3090: {
3091: ifac_realloc(partial, &here, 1); /* guaranteed to work */
3092: /* Unfortunately, we can't do a garbage collection here since we
3093: know too little about where in the stack the old components
3094: were. */
3095: }
3096: }
3097: nf = ifac_crack(partial, &here);
3098: if (moebius_mode && here[1] != un) /* that was a power */
3099: {
3100: if (DEBUGLEVEL >= 3)
3101: {
3102: fprintferr("IFAC: main loop: repeated new factor\n\t%Z\n", *here);
3103: flusherr();
3104: }
3105: return gzero;
3106: }
3107: /* deal with the new unknowns. No resort, since ifac_crack will
3108: already have sorted them */
3109: ifac_whoiswho(partial, &here, nf);
3110: continue;
3111: }
3112: /* if it's prime but not yet finished, finish it */
3113: if (here[2] == un)
3114: {
3115: res = ifac_divide(partial, &here);
3116: if (res)
3117: {
3118: if (moebius_mode)
3119: {
3120: if (DEBUGLEVEL >= 3)
3121: {
3122: fprintferr("IFAC: main loop: another factor was divisible by\n");
3123: fprintferr("\t%Z\n", *here); flusherr();
3124: }
3125: return gzero;
3126: }
3127: ifac_defrag(partial, &here);
3128: (void)(ifac_resort(partial, &here)); /* sort new cofactors down */
3129: /* it doesn't matter right now whether this finds a repeated factor,
3130: since we never get to this point in Moebius mode */
3131: ifac_defrag(partial, &here); /* resort may have created new gaps */
3132: ifac_whoiswho(partial, &here, -1);
3133: }
3134: continue;
3135: }
3136: /* there are no other cases, never reached */
3137: err(talker, "non-existent factor class in ifac_main");
3138: } /* while */
3139: if (moebius_mode && here[1] != un)
3140: {
3141: if (DEBUGLEVEL >= 3)
3142: {
3143: fprintferr("IFAC: after main loop: repeated old factor\n\t%Z\n", *here);
3144: flusherr();
3145: }
3146: return gzero; /* just a safety net */
3147: }
3148: if (DEBUGLEVEL >= 4)
3149: {
3150: long nf = (*partial + lg(*partial) - here - 3)/3;
3151: if (nf)
3152: fprintferr("IFAC: main loop: %ld factor%s left\n",
3153: nf, (nf>1 ? "s" : ""));
3154: else
3155: fprintferr("IFAC: main loop: this was the last factor\n");
3156: flusherr();
3157: }
3158: return here;
3159: }
3160:
3161: /* Caller of the following should worry about stack management, it makes
3162: a rather shameless mess :^) */
3163: GEN
3164: ifac_primary_factor(GEN *partial, long *exponent)
3165: {
3166: GEN here = ifac_main(partial);
3167: GEN res;
3168:
3169: if (here == gun) { *exponent = 0; return gun; }
3170: else if (here == gzero) { *exponent = 0; return gzero; }
3171:
3172: res = icopy((GEN)(*here));
3173: *exponent = itos((GEN)(here[1]));
3174: here[2] = here[1] = *here = (long)NULL;
3175: return res;
3176: }
3177:
3178: /* encapsulated routines */
3179:
3180: /* prime/exponent pairs need to appear contiguously on the stack, but we
3181: also need to have our data structure somewhere, and we don't know in
3182: advance how many primes will turn up. The following discipline achieves
3183: this: When ifac_decomp() is called, n should point at an object older
3184: than the oldest small prime/exponent pair (auxdecomp0() guarantees
3185: this easily since it mpdivis()es any divisors it discovers off its own
3186: copy of the original N). We allocate sufficient space to accommodate
3187: several pairs -- eleven pairs ought to fit in a space not much larger
3188: than n itself -- before calling ifac_start(). If we manage to complete
3189: the factorization before we run out of space, we free the data structure
3190: and cull the excess reserved space before returning. When we do run out,
3191: we have to leapfrog to generate more (guesstimating the requirements
3192: from what is left in the partial factorization structure); room for
3193: fresh pairs is allocated at the head of the stack, followed by an
3194: ifac_realloc() to reconnect the data structure and move it out of the
3195: way, followed by a few pointer tweaks to connect the new pairs space
3196: to the old one.-- This whole affair translates into a surprisingly
3197: compact little routine. */
3198:
3199: #define ifac_overshoot 64 /* lgefint(n)+64 words reserved */
3200:
3201: long
3202: ifac_decomp(GEN n, long hint)
3203: {
3204: long tf=lgefint(n), av=avma, lim=stack_lim(av,1);
3205: long nb=0;
3206: GEN part, here, workspc = new_chunk(tf + ifac_overshoot), pairs = (GEN)av;
3207: /* workspc will be doled out by us in pairs of smaller t_INTs */
3208: long tetpil = avma; /* remember head of workspc zone */
3209:
3210: if (!n || typ(n) != t_INT) err(typeer, "ifac_decomp");
3211: if (!signe(n) || tf < 3) err(talker, "factoring 0 in ifac_decomp");
3212:
3213: part = ifac_start(n, 0, hint);
3214: here = ifac_main(&part);
3215:
3216: while (here != gun)
3217: {
3218: long lf=lgefint((GEN)(*here));
3219: if (pairs - workspc < lf + 3) /* out of room, leapfrog */
3220: {
3221: /* the ifac_realloc() below will clear tetpil - avma words
3222: on the stack, which should be about enough for the extra
3223: primes we're going to see, and we'll want several more to
3224: accommodate further exponents. In most cases, the lf + 3
3225: below is pure paranoia, but the factor we're about to copy
3226: might be the one sitting off the stack in the original n,
3227: so let's play safe */
3228: workspc = new_chunk(lf + 3 + ifac_overshoot);
3229: ifac_realloc(&part, &here, 0);
3230: here = ifac_find(&part, &part);
3231: tetpil = (long)workspc;
3232: }
3233: /* room enough now */
3234: nb++;
3235: pairs -= lf;
3236: *pairs = evaltyp(t_INT) | evallg(lf);
3237: affii((GEN)(*here), pairs);
3238: pairs -= 3;
3239: *pairs = evaltyp(t_INT) | evallg(3);
3240: affii((GEN)(here[1]), pairs);
3241: here[2] = here[1] = *here = (long)NULL;
3242: here = ifac_main(&part);
3243: if (low_stack(lim, stack_lim(av,1)))
3244: {
3245: if(DEBUGMEM>1) err(warnmem,"[2] ifac_decomp");
3246: ifac_realloc(&part, &here, 0);
3247: part = gerepileupto(tetpil, part);
3248: }
3249: }
3250: avma = (long)pairs;
3251: if (DEBUGLEVEL >= 3)
3252: {
3253: fprintferr("IFAC: found %ld large prime (power) factor%s.\n",
3254: nb, (nb>1? "s": ""));
3255: flusherr();
3256: }
3257: return nb;
3258: }
3259:
3260: long
3261: ifac_moebius(GEN n, long hint)
3262: {
3263: long mu=1, av=avma, lim=stack_lim(av,1);
3264: GEN part = ifac_start(n, 1, hint);
3265: GEN here = ifac_main(&part);
3266:
3267: while (here != gun && here != gzero)
3268: {
3269: if (itos((GEN)(here[1])) > 1)
3270: { here = gzero; break; } /* shouldn't happen */
3271: mu = -mu;
3272: here[2] = here[1] = *here = (long)NULL;
3273: here = ifac_main(&part);
3274: if (low_stack(lim, stack_lim(av,1)))
3275: {
3276: if(DEBUGMEM>1) err(warnmem,"ifac_moebius");
3277: ifac_realloc(&part, &here, 0);
3278: part = gerepileupto(av, part);
3279: }
3280: }
3281: avma = av;
3282: return (here == gun ? mu : 0);
3283: }
3284:
3285: long
3286: ifac_issquarefree(GEN n, long hint)
3287: {
3288: long av=avma, lim=stack_lim(av,1);
3289: GEN part = ifac_start(n, 1, hint);
3290: GEN here = ifac_main(&part);
3291:
3292: while (here != gun && here != gzero)
3293: {
3294: if (itos((GEN)(here[1])) > 1)
3295: { here = gzero; break; } /* shouldn't happen */
3296: here[2] = here[1] = *here = (long)NULL;
3297: here = ifac_main(&part);
3298: if (low_stack(lim, stack_lim(av,1)))
3299: {
3300: if(DEBUGMEM>1) err(warnmem,"ifac_issquarefree");
3301: ifac_realloc(&part, &here, 0);
3302: part = gerepileupto(av, part);
3303: }
3304: }
3305: avma = av;
3306: return (here == gun ? 1 : 0);
3307: }
3308:
3309: long
3310: ifac_omega(GEN n, long hint)
3311: {
3312: long omega=0, av=avma, lim=stack_lim(av,1);
3313: GEN part = ifac_start(n, 0, hint);
3314: GEN here = ifac_main(&part);
3315:
3316: while (here != gun)
3317: {
3318: omega++;
3319: here[2] = here[1] = *here = (long)NULL;
3320: here = ifac_main(&part);
3321: if (low_stack(lim, stack_lim(av,1)))
3322: {
3323: if(DEBUGMEM>1) err(warnmem,"ifac_omega");
3324: ifac_realloc(&part, &here, 0);
3325: part = gerepileupto(av, part);
3326: }
3327: }
3328: avma = av;
3329: return omega;
3330: }
3331:
3332: long
3333: ifac_bigomega(GEN n, long hint)
3334: {
3335: long Omega=0, av=avma, lim=stack_lim(av,1);
3336: GEN part = ifac_start(n, 0, hint);
3337: GEN here = ifac_main(&part);
3338:
3339: while (here != gun)
3340: {
3341: Omega += itos((GEN)(here[1]));
3342: here[2] = here[1] = *here = (long)NULL;
3343: here = ifac_main(&part);
3344: if (low_stack(lim, stack_lim(av,1)))
3345: {
3346: if(DEBUGMEM>1) err(warnmem,"ifac_bigomega");
3347: ifac_realloc(&part, &here, 0);
3348: part = gerepileupto(av, part);
3349: }
3350: }
3351: avma = av;
3352: return Omega;
3353: }
3354:
3355: GEN
3356: ifac_totient(GEN n, long hint)
3357: {
3358: GEN res = cgeti(lgefint(n));
3359: long exponent, av=avma, tetpil, lim=stack_lim(av,1);
3360: GEN phi = gun;
3361: GEN part = ifac_start(n, 0, hint);
3362: GEN here = ifac_main(&part);
3363:
3364: while (here != gun)
3365: {
3366: phi = mulii(phi, addsi(-1, (GEN)(*here)));
3367: if (here[1] != un)
3368: {
3369: if (here[1] == deux)
3370: {
3371: phi = mulii(phi, (GEN)(*here));
3372: }
3373: else
3374: {
3375: exponent = itos((GEN)(here[1]));
3376: phi = mulii(phi, gpowgs((GEN)(*here), exponent-1));
3377: }
3378: }
3379: here[2] = here[1] = *here = (long)NULL;
3380: here = ifac_main(&part);
3381: if (low_stack(lim, stack_lim(av,1)))
3382: {
3383: GEN *gsav[2];
3384: if(DEBUGMEM>1) err(warnmem,"ifac_totient");
3385: tetpil = avma;
3386: ifac_realloc(&part, &here, 0);
3387: phi = icopy(phi);
3388: gsav[0] = φ gsav[1] = ∂
3389: gerepilemanysp(av, tetpil, gsav, 2);
3390: /* don't try to preserve here, safer to pick it up again
3391: (and ifac_find does a lot of sanity checking at high
3392: debuglevels) */
3393: here = ifac_find(&part, &part);
3394: }
3395: }
3396: affii(phi, res);
3397: avma = av;
3398: return res;
3399: }
3400:
3401: GEN
3402: ifac_numdiv(GEN n, long hint)
3403: {
3404: /* we don't preallocate since it's too hard to guess the right
3405: size here */
3406: GEN res;
3407: long av=avma, tetpil, lim=stack_lim(av,1);
3408: GEN exponent, tau = gun;
3409: GEN part = ifac_start(n, 0, hint);
3410: GEN here = ifac_main(&part);
3411:
3412: while (here != gun)
3413: {
3414: exponent = (GEN)(here[1]);
3415: tau = mulii(tau, addsi(1, exponent));
3416: here[2] = here[1] = *here = (long)NULL;
3417: here = ifac_main(&part);
3418: if (low_stack(lim, stack_lim(av,1)))
3419: {
3420: GEN *gsav[2];
3421: if(DEBUGMEM>1) err(warnmem,"ifac_numdiv");
3422: tetpil = avma;
3423: ifac_realloc(&part, &here, 0);
3424: tau = icopy(tau);
3425: gsav[0] = τ gsav[1] = ∂
3426: gerepilemanysp(av, tetpil, gsav, 2);
3427: /* (see ifac_totient()) */
3428: here = ifac_find(&part, &part);
3429: }
3430: }
3431: tetpil = avma;
3432: res = icopy(tau);
3433: return gerepile(av, tetpil, res);
3434: }
3435:
3436: GEN
3437: ifac_sumdiv(GEN n, long hint)
3438: {
3439: /* don't preallocate */
3440: GEN res;
3441: long exponent, av=avma, tetpil, lim=stack_lim(av,1);
3442: GEN contrib, sigma = gun;
3443: GEN part = ifac_start(n, 0, hint);
3444: GEN here = ifac_main(&part);
3445:
3446: while (here != gun)
3447: {
3448: exponent = itos((GEN)(here[1]));
3449: contrib = addsi(1, (GEN)(*here));
3450: for (; exponent > 1; exponent--)
3451: contrib = addsi(1, mulii((GEN)(*here), contrib));
3452: sigma = mulii(sigma, contrib);
3453: here[2] = here[1] = *here = (long)NULL;
3454: here = ifac_main(&part);
3455: if (low_stack(lim, stack_lim(av,1)))
3456: {
3457: GEN *gsav[2];
3458: if(DEBUGMEM>1) err(warnmem,"ifac_sumdiv");
3459: tetpil = avma;
3460: ifac_realloc(&part, &here, 0);
3461: sigma = icopy(sigma);
3462: gsav[0] = σ gsav[1] = ∂
3463: gerepilemanysp(av, tetpil, gsav, 2);
3464: /* (see ifac_totient()) */
3465: here = ifac_find(&part, &part);
3466: }
3467: }
3468: tetpil = avma;
3469: res = icopy(sigma);
3470: return gerepile(av, tetpil, res);
3471: }
3472:
3473: /* k should be positive, and indeed it had better be > 1 (not checked).
3474: The calling function knows what to do with the other cases. */
3475: GEN
3476: ifac_sumdivk(GEN n, long k, long hint)
3477: {
3478: /* don't preallocate */
3479: GEN res;
3480: long exponent, av=avma, tetpil, lim=stack_lim(av,1);
3481: GEN contrib, q, sigma = gun;
3482: GEN part = ifac_start(n, 0, hint);
3483: GEN here = ifac_main(&part);
3484:
3485: while (here != gun)
3486: {
3487: exponent = itos((GEN)(here[1]));
3488: q = gpowgs((GEN)(*here), k);
3489: contrib = addsi(1, q);
3490: for (; exponent > 1; exponent--)
3491: contrib = addsi(1, mulii(q, contrib));
3492: sigma = mulii(sigma, contrib);
3493: here[2] = here[1] = *here = (long)NULL;
3494: here = ifac_main(&part);
3495: if (low_stack(lim, stack_lim(av,1)))
3496: {
3497: GEN *gsav[2];
3498: if(DEBUGMEM>1) err(warnmem,"ifac_sumdivk");
3499: tetpil = avma;
3500: ifac_realloc(&part, &here, 0);
3501: sigma = icopy(sigma);
3502: gsav[0] = σ gsav[1] = ∂
3503: gerepilemanysp(av, tetpil, gsav, 2);
3504: /* (see ifac_totient()) */
3505: here = ifac_find(&part, &part);
3506: }
3507: }
3508: tetpil = avma;
3509: res = icopy(sigma);
3510: return gerepile(av, tetpil, res);
3511: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ifactor1.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / basemath