version 1.1, 2001/10/02 11:17:04 |
version 1.2, 2002/09/11 07:26:51 |
Line 19 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
Line 19 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
/** **/ |
/** **/ |
/********************************************************************/ |
/********************************************************************/ |
#include "pari.h" |
#include "pari.h" |
|
extern GEN decomp_limit(GEN n, GEN limit); |
|
extern int BSW_isprime(GEN x); |
|
extern int BSW_isprime_small(GEN x); |
|
extern ulong ucarrecomplet(ulong A); |
|
extern GEN mpqs(GEN N);/* in src/modules/mpqs.c, |
|
* returns a factor, a vector of factors, or NULL */ |
|
|
/*********************************************************************/ |
/*********************************************************************/ |
/** **/ |
/** **/ |
Line 47 init_miller(GEN n) |
|
Line 53 init_miller(GEN n) |
|
static int |
static int |
bad_for_base(GEN n, GEN a) |
bad_for_base(GEN n, GEN a) |
{ |
{ |
long r, av=avma, lim=stack_lim(av,1); |
long r; |
|
gpmem_t av=avma, lim=stack_lim(av, 1); |
GEN c2, c = powmodulo(a,t1,n); |
GEN c2, c = powmodulo(a,t1,n); |
|
|
if (!is_pm1(c) && !egalii(t,c)) /* go fishing for -1, not for 1 */ |
if (!is_pm1(c) && !egalii(t,c)) /* go fishing for -1, not for 1 */ |
Line 79 bad_for_base(GEN n, GEN a) |
|
Line 86 bad_for_base(GEN n, GEN a) |
|
long |
long |
millerrabin(GEN n, long k) |
millerrabin(GEN n, long k) |
{ |
{ |
long r,i,av2, av = avma; |
long r, i; |
|
gpmem_t av2, av = avma; |
|
|
if (!signe(n)) return 0; |
if (!signe(n)) return 0; |
/* If |n| <= 3, check if n = +- 1 */ |
/* If |n| <= 3, check if n = +- 1 */ |
Line 135 millerrabin(GEN n, long k) |
|
Line 143 millerrabin(GEN n, long k) |
|
int /* no longer static -- needed in mpqs.c */ |
int /* no longer static -- needed in mpqs.c */ |
miller(GEN n, long k) |
miller(GEN n, long k) |
{ |
{ |
long r,i,av2, av = avma; |
long r, i; |
|
gpmem_t av2, av = avma; |
static long pr[] = |
static long pr[] = |
{ 0, 2,3,5,7,11,13,17,19,23,29, 31,73, 2,13,23,1662803UL, }; |
{ 0, 2,3,5,7,11,13,17,19,23,29, 31,73, 2,13,23,1662803UL, }; |
long *p; |
long *p; |
Line 163 miller(GEN n, long k) |
|
Line 172 miller(GEN n, long k) |
|
} |
} |
avma=av; return 1; |
avma=av; return 1; |
} |
} |
|
|
|
/* compute n-th term of Lucas sequence modulo N. |
|
* v_{k+2} = P v_{k+1} - v_k, v_0 = 2, v_1 = P. |
|
* Assume n > 0 */ |
|
static GEN |
|
LucasMod(GEN n, long P, GEN N) |
|
{ |
|
gpmem_t av = avma, lim = stack_lim(av, 1); |
|
GEN nd = n+2; |
|
long i, m = *nd, j = 1+bfffo((ulong)m); |
|
GEN v = stoi(P), v1 = stoi(P*P - 2); |
|
|
|
m <<= j; j = BITS_IN_LONG - j; |
|
for (i=lgefint(n)-2;;) /* cf. leftright_pow */ |
|
{ |
|
for (; j; m<<=1,j--) |
|
{ /* v = v_k, v1 = v_{k+1} */ |
|
if (m < 0) |
|
{ /* set v = v_{2k+1}, v1 = v_{2k+2} */ |
|
v = subis(mulii(v,v1), P); |
|
v1= subis(sqri(v1), 2); |
|
} |
|
else |
|
{/* set v = v_{2k}, v1 = v_{2k+1} */ |
|
v1= subis(mulii(v,v1), P); |
|
v = subis(sqri(v), 2); |
|
} |
|
v = modii(v, N); |
|
v1= modii(v1,N); |
|
if (low_stack(lim,stack_lim(av,1))) |
|
{ |
|
GEN *gptr[2]; gptr[0]=&v; gptr[1]=&v1; |
|
if(DEBUGMEM>1) err(warnmem,"LucasMod"); |
|
gerepilemany(av,gptr,2); |
|
} |
|
} |
|
if (--i == 0) return v; |
|
j = BITS_IN_LONG; |
|
m = *++nd; |
|
} |
|
} |
|
|
|
/* check that N not a square first (taken care of here, but inefficient) |
|
* assume N > 3 */ |
|
static int |
|
IsLucasPsP0(GEN N) |
|
{ |
|
long b, i, v; |
|
GEN N_2, m, z; |
|
|
|
for (b=3, i=0;; b+=2, i++) |
|
{ |
|
if (i == 64 && carreparfait(N)) return 0; /* avoid oo loop if N = m^2 */ |
|
if (krosg(b*b - 4, N) < 0) break; |
|
} |
|
|
|
m = addis(N,1); v = vali(m); m = shifti(m,-v); |
|
z = LucasMod(m, b, N); |
|
if (egalii(z, gdeux)) return 1; |
|
N_2 = subis(N,2); |
|
if (egalii(z, N_2)) return 1; |
|
for (i=1; i<v; i++) |
|
{ |
|
if (!signe(z)) return 1; |
|
z = modii(subis(sqri(z), 2), N); |
|
if (egalii(z, gdeux)) return 0; |
|
} |
|
return 0; |
|
} |
|
|
|
long |
|
BSW_psp(GEN N) |
|
{ |
|
gpmem_t av = avma; |
|
int k; |
|
GEN T; |
|
|
|
if (typ(N) != t_INT) err(arither1); |
|
if (signe(N) <= 0) return 0; |
|
if (!is_bigint(N)) |
|
switch(itos(N)) |
|
{ |
|
case 1: return 0; |
|
case 2: |
|
case 3: return 1; |
|
} |
|
if (!mod2(N)) return 0; |
|
|
|
T = init_miller(N); |
|
k = bad_for_base(T,gdeux); |
|
avma = av; |
|
k = (!k && IsLucasPsP0(N)); |
|
avma = av; return k; |
|
} |
|
|
/***********************************************************************/ |
/***********************************************************************/ |
/** **/ |
/** **/ |
/** Pocklington-Lehmer **/ |
/** Pocklington-Lehmer **/ |
Line 173 miller(GEN n, long k) |
|
Line 277 miller(GEN n, long k) |
|
/*assume n>=2*/ |
/*assume n>=2*/ |
static long pl831(GEN N, GEN p) |
static long pl831(GEN N, GEN p) |
{ |
{ |
ulong ltop=avma,av; |
gpmem_t ltop=avma, av; |
long a; |
long a; |
GEN Nmun,Nmunp; |
GEN Nmun,Nmunp; |
Nmun=addis(N,-1); |
Nmun=addis(N,-1); |
Line 187 static long pl831(GEN N, GEN p) |
|
Line 291 static long pl831(GEN N, GEN p) |
|
{ |
{ |
GEN g; |
GEN g; |
g=mppgcd(addis(b,-1),N); |
g=mppgcd(addis(b,-1),N); |
if (gcmp1(g)) |
if (gcmp1(g)) { avma=ltop; return a; } |
{ |
if (!gegal(g,N)) { avma=ltop; return 0; } |
avma=ltop; |
|
return a; |
|
} |
|
if (!gegal(g,N)) |
|
{ |
|
avma=ltop; |
|
return 0; |
|
} |
|
} |
} |
else |
else { avma=ltop; return 0; } |
{ |
|
avma=ltop; |
|
return 0; |
|
} |
|
avma=av; |
avma=av; |
} |
} |
} |
} |
/* |
/* Assume N is a strong BSW pseudoprime |
|
* |
* flag 0: return gun (prime), gzero (composite) |
* flag 0: return gun (prime), gzero (composite) |
* flag 1: return gzero (composite), gun (small prime), matrix (large prime) |
* flag 1: return gzero (composite), gun (small prime), matrix (large prime) |
* |
* |
* The matrix has 3 columns, [a,b,c] with |
* The matrix has 3 columns, [a,b,c] with |
* a[i] prime factor of N-1, |
* a[i] prime factor of N-1, |
* b[i] witness for a[i] as in pl831 |
* b[i] witness for a[i] as in pl831 |
* c[i] plisprime(a[i]) |
* c[i] plisprime(a[i]) |
*/ |
*/ |
extern GEN decomp_limit(GEN n, GEN limit); |
|
GEN |
GEN |
plisprime(GEN N, long flag) |
plisprime(GEN N, long flag) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop = avma; |
long i; |
long i, l, t = typ(N); |
int eps; |
int eps; |
GEN C,F; |
GEN C, F = NULL; |
if ( typ(N) != t_INT ) err(arither1); |
|
|
if (t == t_VEC) |
|
{ /* [ N, [p1,...,pk] ], pi list of pseudoprime divisors of N */ |
|
F = (GEN)N[2]; |
|
N = (GEN)N[1]; t = typ(N); |
|
} |
|
if (t != t_INT) err(arither1); |
|
if (DEBUGLEVEL>3) fprintferr("PL: proving primality of N = %Z\n", N); |
|
|
eps = absi_cmp(N,gdeux); |
eps = absi_cmp(N,gdeux); |
if (eps<=0) return eps? gzero: gun; |
if (eps<=0) return eps? gzero: gun; |
|
|
N = absi(N); |
N = absi(N); |
/* Use Jaeschke results. See above */ |
if (!F) |
if (miller(N,7)) |
{ |
{ /* compare to 341550071728321 */ |
F = (GEN)decomp_limit(addis(N,-1), racine(N))[1]; |
if (cmpii(N, u2toi(0x136a3, 0x52b2c8c1)) < 0) { avma=ltop; return gun; } |
if (DEBUGLEVEL>3) fprintferr("PL: N-1 factored!\n"); |
} |
} |
else { avma=ltop; return gzero; } |
|
F=(GEN)decomp_limit(addis(N,-1),racine(N))[1]; |
C = cgetg(4,t_MAT); l = lg(F); |
if (DEBUGLEVEL>=3) fprintferr("P.L.:factor O.K.\n"); |
C[1] = lgetg(l,t_COL); |
C=cgetg(4,t_MAT); |
C[2] = lgetg(l,t_COL); |
C[1]=lgetg(lg(F),t_COL); |
C[3] = lgetg(l,t_COL); |
C[2]=lgetg(lg(F),t_COL); |
for(i=1; i<l; i++) |
C[3]=lgetg(lg(F),t_COL); |
|
for(i=1;i<lg(F);i++) |
|
{ |
{ |
long witness; |
GEN p = (GEN)F[i], r; |
GEN p; |
long witness = pl831(N,p); |
p=(GEN)F[i]; |
|
witness=pl831(N,p); |
if (!witness) { avma = ltop; return gzero; } |
if (!witness) |
mael(C,1,i) = licopy(p); |
|
mael(C,2,i) = lstoi(witness); |
|
if (!flag) r = BSW_isprime(p)? gun: gzero; |
|
else |
{ |
{ |
avma=ltop; |
if (BSW_isprime_small(p)) r = gun; |
return gzero; |
else if (expi(p) > 250) r = isprimeAPRCL(p)? gdeux: gzero; |
|
else r = plisprime(p,flag); |
} |
} |
mael(C,1,i)=lcopy(p); |
mael(C,3,i) = (long)r; |
mael(C,2,i)=lstoi(witness); |
if (r == gzero) err(talker,"False prime number %Z in plisprime", p); |
mael(C,3,i)=(long)plisprime(p,flag); |
|
if (gmael(C,3,i)==gzero) |
|
err(talker,"Sorry false prime number %Z in plisprime",p); |
|
} |
} |
if (!flag) { avma=ltop; return gun; } |
if (!flag) { avma = ltop; return gun; } |
return gerepileupto(ltop,C); |
return gerepileupto(ltop,C); |
} |
} |
|
|
Line 311 unsigned char prc210_d1[] = |
|
Line 411 unsigned char prc210_d1[] = |
|
GEN |
GEN |
nextprime(GEN n) |
nextprime(GEN n) |
{ |
{ |
long rc,rc0,rcd,rcn,av1,av2, av = avma; |
long rc, rc0, rcd, rcn; |
|
gpmem_t av1, av2, av = avma; |
|
|
if (typ(n) != t_INT) n=gceil(n); /* accept arguments in R --GN */ |
if (typ(n) != t_INT) n=gceil(n); /* accept arguments in R --GN */ |
if (typ(n) != t_INT) err(arither1); |
if (typ(n) != t_INT) err(arither1); |
Line 339 nextprime(GEN n) |
|
Line 440 nextprime(GEN n) |
|
av2 = av1 = avma; |
av2 = av1 = avma; |
for(;;) |
for(;;) |
{ |
{ |
if (miller(n,10)) break; |
if (BSW_psp(n)) break; |
av1 = avma; |
av1 = avma; |
rcd = prc210_d1[rcn]; |
rcd = prc210_d1[rcn]; |
if (++rcn > 47) rcn = 0; |
if (++rcn > 47) rcn = 0; |
Line 352 nextprime(GEN n) |
|
Line 453 nextprime(GEN n) |
|
GEN |
GEN |
precprime(GEN n) |
precprime(GEN n) |
{ |
{ |
long rc,rc0,rcd,rcn,av1,av2, av = avma; |
long rc, rc0, rcd, rcn; |
|
gpmem_t av1, av2, av = avma; |
|
|
if (typ(n) != t_INT) n=gfloor(n); /* accept arguments in R --GN */ |
if (typ(n) != t_INT) n=gfloor(n); /* accept arguments in R --GN */ |
if (typ(n) != t_INT) err(arither1); |
if (typ(n) != t_INT) err(arither1); |
Line 381 precprime(GEN n) |
|
Line 483 precprime(GEN n) |
|
av2 = av1 = avma; |
av2 = av1 = avma; |
for(;;) |
for(;;) |
{ |
{ |
if (miller(n,10)) break; |
if (BSW_psp(n)) break; |
av1 = avma; |
av1 = avma; |
if (rcn == 0) |
if (rcn == 0) |
{ rcd = 2; rcn = 47; } |
{ rcd = 2; rcn = 47; } |
|
|
snextpr(ulong p, byteptr *d, long *rcn, long *q, long k) |
snextpr(ulong p, byteptr *d, long *rcn, long *q, long k) |
{ |
{ |
static ulong pp[] = |
static ulong pp[] = |
{ evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0 }; |
{ evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0 }; |
static ulong *pp2 = pp + 2; |
static ulong *pp2 = pp + 2; |
static GEN gp = (GEN)pp; |
static GEN gp = (GEN)pp; |
long d1 = **d, rcn0; |
long rcn0; |
|
|
if (d1) |
if (**d) |
{ |
{ |
|
byteptr dd = *d; |
|
long d1 = 0; |
|
|
|
NEXT_PRIME_VIADIFF(d1,dd); |
if (*rcn != NPRC) |
if (*rcn != NPRC) |
{ |
{ |
rcn0 = *rcn; |
rcn0 = *rcn; |
Line 430 snextpr(ulong p, byteptr *d, long *rcn, long *q, long |
|
Line 536 snextpr(ulong p, byteptr *d, long *rcn, long *q, long |
|
err(bugparier, "[caller of] snextpr"); |
err(bugparier, "[caller of] snextpr"); |
} |
} |
} |
} |
return p + *(*d)++; |
NEXT_PRIME_VIADIFF(p,*d); |
|
return p; |
} |
} |
/* we are beyond the diffptr table */ |
/* we are beyond the diffptr table */ |
if (*rcn == NPRC) /* we need to initialize this now */ |
if (*rcn == NPRC) /* we need to initialize this now */ |
Line 542 elladd0(long nbc, long nbc1, |
|
Line 649 elladd0(long nbc, long nbc1, |
|
{ |
{ |
GEN lambda; |
GEN lambda; |
GEN W[2*nbcmax], *A=W+nbc; /* W[0],A[0] never used */ |
GEN W[2*nbcmax], *A=W+nbc; /* W[0],A[0] never used */ |
long i, av=avma, tetpil; |
long i; |
|
gpmem_t av=avma, tetpil; |
ulong mask = ~0UL; |
ulong mask = ~0UL; |
|
|
/* actually, this is only ever called with nbc1==nbc or nbc1==4, so: */ |
/* actually, this is only ever called with nbc1==nbc or nbc1==4, so: */ |
Line 613 elladd2(long nbc, GEN *X1, GEN *X2, GEN *X3, GEN *X4, |
|
Line 721 elladd2(long nbc, GEN *X1, GEN *X2, GEN *X3, GEN *X4, |
|
GEN lambda, *Y1 = X1+nbc, *Y2 = X2+nbc, *Y3 = X3+nbc; |
GEN lambda, *Y1 = X1+nbc, *Y2 = X2+nbc, *Y3 = X3+nbc; |
GEN *Y4 = X4+nbc, *Y5 = X5+nbc, *Y6 = X6+nbc; |
GEN *Y4 = X4+nbc, *Y5 = X5+nbc, *Y6 = X6+nbc; |
GEN W[4*nbcmax], *A=W+2*nbc; /* W[0],A[0] never used */ |
GEN W[4*nbcmax], *A=W+2*nbc; /* W[0],A[0] never used */ |
long i,j, av=avma, tetpil; |
long i, j; |
|
gpmem_t av=avma, tetpil; |
/* W[0] = gun; */ |
/* W[0] = gun; */ |
W[1] = /* A[0] =*/ subii(X1[0], X2[0]); |
W[1] = /* A[0] =*/ subii(X1[0], X2[0]); |
for (i=1; i<nbc; i++) |
for (i=1; i<nbc; i++) |
Line 685 elldouble(long nbc, GEN *X1, GEN *X2) |
|
Line 794 elldouble(long nbc, GEN *X1, GEN *X2) |
|
{ |
{ |
GEN lambda,v, *Y1 = X1+nbc, *Y2 = X2+nbc; |
GEN lambda,v, *Y1 = X1+nbc, *Y2 = X2+nbc; |
GEN W[nbcmax+1]; /* W[0] never used */ |
GEN W[nbcmax+1]; /* W[0] never used */ |
long i, av=avma, tetpil; |
long i; |
|
gpmem_t av=avma, tetpil; |
/*W[0] = gun;*/ W[1] = Y1[0]; |
/*W[0] = gun;*/ W[1] = Y1[0]; |
for (i=1; i<nbc; i++) |
for (i=1; i<nbc; i++) |
W[i+1] = modii(mulii(Y1[i], W[i]), N); |
W[i+1] = modii(mulii(Y1[i], W[i]), N); |
Line 736 elldouble(long nbc, GEN *X1, GEN *X2) |
|
Line 846 elldouble(long nbc, GEN *X1, GEN *X2) |
|
static int |
static int |
ellmult(long nbc, ulong k, GEN *X1, GEN *X2) /* k>2 prime, not checked */ |
ellmult(long nbc, ulong k, GEN *X1, GEN *X2) /* k>2 prime, not checked */ |
{ |
{ |
long i,d,e,e1,r,av=avma,tetpil; |
long i, d, e, e1, r; |
|
gpmem_t av=avma, tetpil; |
int res; |
int res; |
GEN *A=X2, *B=XAUX, *S, *T=XAUX+2*nbc; |
GEN *A=X2, *B=XAUX, *S, *T=XAUX+2*nbc; |
|
|
Line 1039 ellfacteur(GEN n, int insist) |
|
Line 1150 ellfacteur(GEN n, int insist) |
|
540400UL,606000UL,679500UL,761800UL,854100UL,957500UL,1073500UL, |
540400UL,606000UL,679500UL,761800UL,854100UL,957500UL,1073500UL, |
}; |
}; |
long nbc,nbc2,dsn,dsnmax,rep,spc,gse,gss,rcn,rcn0,bstp,bstp0; |
long nbc,nbc2,dsn,dsnmax,rep,spc,gse,gss,rcn,rcn0,bstp,bstp0; |
long a,i,j,k, av,av1,avtmp, size = expi(n) + 1, tf = lgefint(n); |
long a, i, j, k, size = expi(n) + 1, tf = lgefint(n); |
|
gpmem_t av, av1, avtmp; |
ulong B1,B2,B2_p,B2_rt,m,p,p0,p2,dp; |
ulong B1,B2,B2_p,B2_rt,m,p,p0,p2,dp; |
GEN w,w0,x,*X,*XT,*XD,*XG,*YG,*XH,*XB,*XB2,*Xh,*Yh,*Xb, res = cgeti(tf); |
GEN w,w0,x,*X,*XT,*XD,*XG,*YG,*XH,*XB,*XB2,*Xh,*Yh,*Xb, res = cgeti(tf); |
int rflag, use_clones = 0; |
int rflag, use_clones = 0; |
Line 1211 ellfacteur(GEN n, int insist) |
|
Line 1323 ellfacteur(GEN n, int insist) |
|
fprintferr("\tB2 = %6lu,\tgss = %4ld*420\n", B2, gss); |
fprintferr("\tB2 = %6lu,\tgss = %4ld*420\n", B2, gss); |
flusherr(); |
flusherr(); |
} |
} |
p = *d++; |
p = 0; |
|
NEXT_PRIME_VIADIFF(p,d); |
|
|
/* ---B1 PHASE--- */ |
/* ---B1 PHASE--- */ |
/* treat p=2 separately */ |
/* treat p=2 separately */ |
Line 1425 ellfacteur(GEN n, int insist) |
|
Line 1538 ellfacteur(GEN n, int insist) |
|
gl = gun; |
gl = gun; |
av1 = avma; |
av1 = avma; |
/* scratchspace for prod (x_i-x_j) */ |
/* scratchspace for prod (x_i-x_j) */ |
avtmp = (long)new_chunk(8 * lgefint(n)); |
avtmp = (gpmem_t)new_chunk(8 * lgefint(n)); |
/* the correct entry in XB to use depends on bstp and on where we are |
/* the correct entry in XB to use depends on bstp and on where we are |
* on the helix. As we skip from prime to prime, bstp will be incre- |
* on the helix. As we skip from prime to prime, bstp will be incre- |
* mented by snextpr() each time we wrap around through residue class |
* mented by snextpr() each time we wrap around through residue class |
|
|
pollardbrent(GEN n) |
pollardbrent(GEN n) |
{ |
{ |
long tf = lgefint(n), size = 0, delta, retries = 0, msg_mask; |
long tf = lgefint(n), size = 0, delta, retries = 0, msg_mask; |
long c0, c, k, k1, l, avP, avx, GGG, av = avma; |
long c0, c, k, k1, l; |
|
gpmem_t GGG, avP, avx, av = avma; |
GEN x, x1, y, P, g, g1, res; |
GEN x, x1, y, P, g, g1, res; |
|
|
if (DEBUGLEVEL >= 4) (void)timer2(); /* clear timer */ |
if (DEBUGLEVEL >= 4) (void)timer2(); /* clear timer */ |
Line 1582 pollardbrent(GEN n) |
|
Line 1696 pollardbrent(GEN n) |
|
if (tf >= 4) |
if (tf >= 4) |
size = expi(n) + 1; |
size = expi(n) + 1; |
else if (tf == 3) /* try to keep purify happy... */ |
else if (tf == 3) /* try to keep purify happy... */ |
size = BITS_IN_LONG - bfffo(n[2]); |
size = BITS_IN_LONG - bfffo((ulong)n[2]); |
|
|
if (size <= 28) |
if (size <= 28) |
c0 = 32; /* amounts very nearly to `insist'. |
c0 = 32; /* amounts very nearly to `insist'. |
Line 1600 pollardbrent(GEN n) |
|
Line 1714 pollardbrent(GEN n) |
|
/* 301 gives 48121 + tune_pb_min */ |
/* 301 gives 48121 + tune_pb_min */ |
c0 = tune_pb_min + size - 60 + |
c0 = tune_pb_min + size - 60 + |
((size-73)>>1)*((size-70)>>3)*((size-56)>>4); |
((size-73)>>1)*((size-70)>>3)*((size-56)>>4); |
else |
else |
c0 = 49152; /* ECM is faster when it'd take longer */ |
c0 = 49152; /* ECM is faster when it'd take longer */ |
|
|
c = c0 << 5; /* 32 iterations per round */ |
c = c0 << 5; /* 32 iterations per round */ |
|
|
y = cgeti(tf); affsi(2, y); |
y = cgeti(tf); affsi(2, y); |
x1= cgeti(tf); affsi(2, x1); |
x1= cgeti(tf); affsi(2, x1); |
avx = avma; |
avx = avma; |
avP = (long)new_chunk(2 * tf); /* enough for x = addsi(tf+1) */ |
avP = (gpmem_t)new_chunk(2 * tf); /* enough for x = addsi(tf+1) */ |
GGG = (long)new_chunk(4 * tf); /* enough for P = modii(2tf+1, tf) */ |
GGG = (gpmem_t)new_chunk(4 * tf); /* enough for P = modii(2tf+1, tf) */ |
|
|
for (;;) /* terminated under the control of c */ |
for (;;) /* terminated under the control of c */ |
{ |
{ |
|
|
/** (Cf. Algorithm 8.7.2 in ACiCNT) **/ |
/** (Cf. Algorithm 8.7.2 in ACiCNT) **/ |
/** **/ |
/** **/ |
/***********************************************************************/ |
/***********************************************************************/ |
extern ulong ucarrecomplet(ulong A); |
|
static long squfof_ambig(long a, long B, long dd, GEN D, long *cntamb); |
static long squfof_ambig(long a, long B, long dd, GEN D, long *cntamb); |
|
|
#define SQUFOF_BLACKLIST_SZ 64 |
#define SQUFOF_BLACKLIST_SZ 64 |
Line 1837 squfof(GEN n, long quiet) |
|
Line 1950 squfof(GEN n, long quiet) |
|
long tf = lgefint(n), nm4, cnt = 0, cntamb; |
long tf = lgefint(n), nm4, cnt = 0, cntamb; |
long a1, b1, c1, d1, dd1, L1, a2, b2, c2, d2, dd2, L2, a, q, c, qc, qcb; |
long a1, b1, c1, d1, dd1, L1, a2, b2, c2, d2, dd2, L2, a, q, c, qc, qcb; |
GEN D1, D2, Q, res; |
GEN D1, D2, Q, res; |
long av = avma; |
gpmem_t av = avma; |
static long blacklist1[SQUFOF_BLACKLIST_SZ], blacklist2[SQUFOF_BLACKLIST_SZ]; |
static long blacklist1[SQUFOF_BLACKLIST_SZ], blacklist2[SQUFOF_BLACKLIST_SZ]; |
long blp1 = 0, blp2 = 0; |
long blp1 = 0, blp2 = 0; |
long mydebug = DEBUGLEVEL - quiet; |
long mydebug = DEBUGLEVEL - quiet; |
Line 1846 squfof(GEN n, long quiet) |
|
Line 1959 squfof(GEN n, long quiet) |
|
if (cmpis(n,5) <= 0) return NULL; /* input n <= 5 */ |
if (cmpis(n,5) <= 0) return NULL; /* input n <= 5 */ |
|
|
#ifdef LONG_IS_64BIT |
#ifdef LONG_IS_64BIT |
if (tf > 3 || (tf == 3 && bfffo(n[2]) < 5)) /* n too large */ |
if (tf > 3 || (tf == 3 && bfffo((ulong)n[2]) < 5)) /* n too large */ |
return NULL; |
return NULL; |
#else /* 32 bits */ |
#else /* 32 bits */ |
if (tf > 4 || (tf == 4 && bfffo(n[2]) < 5)) /* n too large */ |
if (tf > 4 || (tf == 4 && bfffo((ulong)n[2]) < 5)) /* n too large */ |
return NULL; |
return NULL; |
#endif |
#endif |
/* now we have 5 < n < 2^59 */ |
/* now we have 5 < n < 2^59 */ |
|
|
long |
long |
squfof_ambig(long a, long B, long dd, GEN D, long *cntamb) |
squfof_ambig(long a, long B, long dd, GEN D, long *cntamb) |
{ |
{ |
long b, c, q, qc, qcb, av = avma; |
long b, c, q, qc, qcb; |
|
gpmem_t av = avma; |
long a0, b0, b1, c0; |
long a0, b0, b1, c0; |
|
|
q = (dd + (B>>1)) / a; qc = q*a; qcb = qc - B; |
q = (dd + (B>>1)) / a; qc = q*a; qcb = qc - B; |
Line 2490 static ulong powersmod[106] = { |
|
Line 2604 static ulong powersmod[106] = { |
|
long /* no longer static -- used in mpqs.c */ |
long /* no longer static -- used in mpqs.c */ |
is_odd_power(GEN x, GEN *pt, long *mask) |
is_odd_power(GEN x, GEN *pt, long *mask) |
{ |
{ |
long av=avma, tetpil, lgx=lgefint(x), exponent=0, residue, resbyte; |
long lgx=lgefint(x), exponent=0, residue, resbyte; |
|
gpmem_t av=avma, tetpil; |
GEN y; |
GEN y; |
|
|
*mask &= 7; /* paranoia */ |
*mask &= 7; /* paranoia */ |
Line 3325 ifac_whoiswho(GEN *partial, GEN *where, long after_cra |
|
Line 3440 ifac_whoiswho(GEN *partial, GEN *where, long after_cra |
|
continue; |
continue; |
} |
} |
scan[2] = |
scan[2] = |
(isprime((GEN)(*scan)) ? |
(BSW_psp((GEN)(*scan)) ? |
(larger_compos ? un : deux) : /* un- or finished prime */ |
(larger_compos ? un : deux) : /* un- or finished prime */ |
zero); /* composite */ |
zero); /* composite */ |
|
|
Line 3415 ifac_divide(GEN *partial, GEN *where) |
|
Line 3530 ifac_divide(GEN *partial, GEN *where) |
|
} |
} |
|
|
|
|
GEN mpqs(GEN N); /* in src/modules/mpqs.c, maybe a dummy, |
|
* returns a factor, or a vector of factors, |
|
* or NULL |
|
*/ |
|
|
|
/* The following takes the place of 2.0.9.alpha's find_factor(). */ |
/* The following takes the place of 2.0.9.alpha's find_factor(). */ |
|
|
/* The meaning of the hint changes against 2.0.9.alpha to: |
/* The meaning of the hint changes against 2.0.9.alpha to: |
Line 3440 GEN mpqs(GEN N); /* in src/modules/mpqs.c, maybe a du |
|
Line 3550 GEN mpqs(GEN N); /* in src/modules/mpqs.c, maybe a du |
|
static long |
static long |
ifac_crack(GEN *partial, GEN *where) |
ifac_crack(GEN *partial, GEN *where) |
{ |
{ |
long hint, cmp_res, exp1 = 1, exp2 = 1, av; |
long hint, cmp_res, exp1 = 1, exp2 = 1; |
|
gpmem_t av; |
GEN factor = NULL, exponent; |
GEN factor = NULL, exponent; |
|
|
if (DEBUGLEVEL >= 5) /* none of these should ever happen */ |
if (DEBUGLEVEL >= 5) /* none of these should ever happen */ |
Line 3486 ifac_crack(GEN *partial, GEN *where) |
|
Line 3597 ifac_crack(GEN *partial, GEN *where) |
|
} /* while carrecomplet */ |
} /* while carrecomplet */ |
|
|
/* check whether our composite hasn't become prime */ |
/* check whether our composite hasn't become prime */ |
if (exp1 > 1 && hint != 15 && isprime((GEN)(**where))) |
if (exp1 > 1 && hint != 15 && BSW_psp((GEN)(**where))) |
{ |
{ |
(*where)[2] = un; |
(*where)[2] = un; |
if (DEBUGLEVEL >= 4) |
if (DEBUGLEVEL >= 4) |
Line 3527 ifac_crack(GEN *partial, GEN *where) |
|
Line 3638 ifac_crack(GEN *partial, GEN *where) |
|
if (moebius_mode) return 0; /* no need to carry on... */ |
if (moebius_mode) return 0; /* no need to carry on... */ |
} /* while is_odd_power */ |
} /* while is_odd_power */ |
|
|
if (exp2 > 1 && hint != 15 && isprime((GEN)(**where))) |
if (exp2 > 1 && hint != 15 && BSW_psp((GEN)(**where))) |
{ /* Something nice has happened and our composite has become prime */ |
{ /* Something nice has happened and our composite has become prime */ |
(*where)[2] = un; |
(*where)[2] = un; |
if (DEBUGLEVEL >= 4) |
if (DEBUGLEVEL >= 4) |
Line 3949 ifac_primary_factor(GEN *partial, long *exponent) |
|
Line 4060 ifac_primary_factor(GEN *partial, long *exponent) |
|
*/ |
*/ |
|
|
#define ifac_overshoot 64 /* lgefint(n)+64 words reserved */ |
#define ifac_overshoot 64 /* lgefint(n)+64 words reserved */ |
/* ifac_decomp_break: |
/* ifac_decomp_break: |
* |
* |
* Find primary factors of n until ifac_break return true, or n is |
* Find primary factors of n until ifac_break return true, or n is |
* factored if ifac_break is NULL. |
* factored if ifac_break is NULL. |
|
|
ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN pairs,GEN here,GEN state), |
ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN pairs,GEN here,GEN state), |
GEN state, long hint) |
GEN state, long hint) |
{ |
{ |
long tf=lgefint(n), av=avma, lim=stack_lim(av,1); |
long tf=lgefint(n); |
|
gpmem_t av=avma, lim=stack_lim(av, 1); |
long nb=0; |
long nb=0; |
GEN part, here, workspc = new_chunk(tf + ifac_overshoot), pairs = (GEN)av; |
GEN part, here, workspc = new_chunk(tf + ifac_overshoot), pairs = (GEN)av; |
/* workspc will be doled out by us in pairs of smaller t_INTs */ |
/* workspc will be doled out by us in pairs of smaller t_INTs */ |
long tetpil = avma; /* remember head of workspc zone */ |
gpmem_t tetpil = avma; /* remember head of workspc zone */; |
|
|
if (!n || typ(n) != t_INT) err(typeer, "ifac_decomp"); |
if (!n || typ(n) != t_INT) err(typeer, "ifac_decomp"); |
if (!signe(n) || tf < 3) err(talker, "factoring 0 in ifac_decomp"); |
if (!signe(n) || tf < 3) err(talker, "factoring 0 in ifac_decomp"); |
Line 3996 ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN |
|
Line 4108 ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN |
|
workspc = new_chunk(lf + 3 + ifac_overshoot); |
workspc = new_chunk(lf + 3 + ifac_overshoot); |
ifac_realloc(&part, &here, 0); |
ifac_realloc(&part, &here, 0); |
here = ifac_find(&part, &part); |
here = ifac_find(&part, &part); |
tetpil = (long)workspc; |
tetpil = (gpmem_t)workspc; |
} |
} |
/* room enough now */ |
/* room enough now */ |
nb++; |
nb++; |
Line 4021 ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN |
|
Line 4133 ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN |
|
part = gerepileupto(tetpil, part); |
part = gerepileupto(tetpil, part); |
} |
} |
} |
} |
avma = (long)pairs; |
avma = (gpmem_t)pairs; |
if (DEBUGLEVEL >= 3) |
if (DEBUGLEVEL >= 3) |
{ |
{ |
fprintferr("IFAC: found %ld large prime (power) factor%s.\n", |
fprintferr("IFAC: found %ld large prime (power) factor%s.\n", |
Line 4040 ifac_decomp(GEN n, long hint) |
|
Line 4152 ifac_decomp(GEN n, long hint) |
|
long |
long |
ifac_moebius(GEN n, long hint) |
ifac_moebius(GEN n, long hint) |
{ |
{ |
long mu=1, av=avma, lim=stack_lim(av,1); |
long mu=1; |
|
gpmem_t av=avma, lim=stack_lim(av, 1); |
GEN part = ifac_start(n, 1, hint); |
GEN part = ifac_start(n, 1, hint); |
GEN here = ifac_main(&part); |
GEN here = ifac_main(&part); |
|
|
Line 4065 ifac_moebius(GEN n, long hint) |
|
Line 4178 ifac_moebius(GEN n, long hint) |
|
long |
long |
ifac_issquarefree(GEN n, long hint) |
ifac_issquarefree(GEN n, long hint) |
{ |
{ |
long av=avma, lim=stack_lim(av,1); |
gpmem_t av=avma, lim=stack_lim(av, 1); |
GEN part = ifac_start(n, 1, hint); |
GEN part = ifac_start(n, 1, hint); |
GEN here = ifac_main(&part); |
GEN here = ifac_main(&part); |
|
|
Line 4089 ifac_issquarefree(GEN n, long hint) |
|
Line 4202 ifac_issquarefree(GEN n, long hint) |
|
long |
long |
ifac_omega(GEN n, long hint) |
ifac_omega(GEN n, long hint) |
{ |
{ |
long omega=0, av=avma, lim=stack_lim(av,1); |
long omega=0; |
|
gpmem_t av=avma, lim=stack_lim(av, 1); |
GEN part = ifac_start(n, 0, hint); |
GEN part = ifac_start(n, 0, hint); |
GEN here = ifac_main(&part); |
GEN here = ifac_main(&part); |
|
|
Line 4112 ifac_omega(GEN n, long hint) |
|
Line 4226 ifac_omega(GEN n, long hint) |
|
long |
long |
ifac_bigomega(GEN n, long hint) |
ifac_bigomega(GEN n, long hint) |
{ |
{ |
long Omega=0, av=avma, lim=stack_lim(av,1); |
long Omega=0; |
|
gpmem_t av=avma, lim=stack_lim(av, 1); |
GEN part = ifac_start(n, 0, hint); |
GEN part = ifac_start(n, 0, hint); |
GEN here = ifac_main(&part); |
GEN here = ifac_main(&part); |
|
|
|
|
ifac_totient(GEN n, long hint) |
ifac_totient(GEN n, long hint) |
{ |
{ |
GEN res = cgeti(lgefint(n)); |
GEN res = cgeti(lgefint(n)); |
long exponent, av=avma, tetpil, lim=stack_lim(av,1); |
long exponent; |
|
gpmem_t av=avma, tetpil, lim=stack_lim(av, 1); |
GEN phi = gun; |
GEN phi = gun; |
GEN part = ifac_start(n, 0, hint); |
GEN part = ifac_start(n, 0, hint); |
GEN here = ifac_main(&part); |
GEN here = ifac_main(&part); |
Line 4186 ifac_numdiv(GEN n, long hint) |
|
Line 4302 ifac_numdiv(GEN n, long hint) |
|
* size here |
* size here |
*/ |
*/ |
GEN res; |
GEN res; |
long av=avma, tetpil, lim=stack_lim(av,1); |
gpmem_t av=avma, tetpil, lim=stack_lim(av, 1); |
GEN exponent, tau = gun; |
GEN exponent, tau = gun; |
GEN part = ifac_start(n, 0, hint); |
GEN part = ifac_start(n, 0, hint); |
GEN here = ifac_main(&part); |
GEN here = ifac_main(&part); |
Line 4220 ifac_sumdiv(GEN n, long hint) |
|
Line 4336 ifac_sumdiv(GEN n, long hint) |
|
{ |
{ |
/* don't preallocate */ |
/* don't preallocate */ |
GEN res; |
GEN res; |
long exponent, av=avma, tetpil, lim=stack_lim(av,1); |
long exponent; |
|
gpmem_t av=avma, tetpil, lim=stack_lim(av, 1); |
GEN contrib, sigma = gun; |
GEN contrib, sigma = gun; |
GEN part = ifac_start(n, 0, hint); |
GEN part = ifac_start(n, 0, hint); |
GEN here = ifac_main(&part); |
GEN here = ifac_main(&part); |
Line 4260 ifac_sumdivk(GEN n, long k, long hint) |
|
Line 4377 ifac_sumdivk(GEN n, long k, long hint) |
|
{ |
{ |
/* don't preallocate */ |
/* don't preallocate */ |
GEN res; |
GEN res; |
long exponent, av=avma, tetpil, lim=stack_lim(av,1); |
long exponent; |
|
gpmem_t av=avma, tetpil, lim=stack_lim(av, 1); |
GEN contrib, q, sigma = gun; |
GEN contrib, q, sigma = gun; |
GEN part = ifac_start(n, 0, hint); |
GEN part = ifac_start(n, 0, hint); |
GEN here = ifac_main(&part); |
GEN here = ifac_main(&part); |