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