File: [local] / OpenXM_contrib / pari-2.2 / src / modules / Attic / aprcl.c (download)
Revision 1.1, Wed Sep 11 07:27:05 2002 UTC (21 years, 9 months ago) by noro
Branch: MAIN
CVS Tags: RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2
Upgraded pari-2.2 to pari-2.2.4.
|
/* $Id: aprcl.c,v 1.23 2002/07/15 13:26:20 karim Exp $
Copyright (C) 2000 The PARI group.
This file is part of the PARI/GP package.
PARI/GP is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software
Foundation. It is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY WHATSOEVER.
Check the License for details. You should have received a copy of it, along
with the package; see the file 'COPYING'. If not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
/*******************************************************************/
/* */
/* THE APRCL PRIMALITY TEST */
/* */
/*******************************************************************/
#include "pari.h"
static ulong kglob,ctglob,NBITSN,sgtaut,sgt6,sgtjac;
static int isinstep5,ishack;
static GEN *tabaall,*tabtall,*tabcyc,tabE,tabTH,tabeta,*tabmatvite,*tabmatinvvite,globfa,tabfaq,tabj,sgt,ctsgt,tabavite,tabpkvite,errfill;
#define pkfalse (ishack ? 1 : pk)
#define dotime (DEBUGLEVEL>=3)
typedef struct red_t {
long n;
GEN C; /* polcyclo(n) */
GEN N; /* prime we are certifying */
GEN (*red)(GEN x, struct red_t *);
} red_t;
/* powering t_POLMOD mod N, flexible window */
static GEN
makepoldeg1(GEN c, GEN d)
{
GEN res;
if (signe(c)) {
res = cgetg(4,t_POL); res[1] = evalsigne(1)|evallgef(4);
res[2] = (long)d; res[3] = (long)c;
} else if (signe(d)) {
res = cgetg(3,t_POL); res[1] = evalsigne(1)|evallgef(3);
res[2] = (long)d;
} else {
res = cgetg(2,t_POL); res[1] = evalsigne(0)|evallgef(2);
}
return res;
}
/* T mod polcyclo(p), assume deg(T) < 2p */
static GEN
red_mod_cyclo(GEN T, long p)
{
long i, d;
GEN y, *z;
d = degpol(T) - p; /* < p */
if (d <= -2) return T;
/* reduce mod (x^p - 1) */
y = dummycopy(T); z = (GEN*)(y+2);
for (i = 0; i<=d; i++) z[i] = addii(z[i], z[i+p]);
/* reduce mod x^(p-1) + ... + 1 */
d = p-1;
if (signe(z[d]))
for (i=0; i<d; i++) z[i] = subii(z[i], z[d]);
return normalizepol_i(y, d+2);
}
/* x t_VECSMALL, n > 0. Return x mod polcyclo(2^n) = (x^(2^(n-1)) + 1). DESTROY x */
static GEN
u_red_mod_cyclo2n(GEN x, int n)
{
long i, pow2 = 1<<(n-1);
GEN z;
for (i = lg(x)-1; i>pow2; i--) x[i-pow2] -= x[i];
for (; i>0; i--)
if (x[i]) break;
i += 2;
z = cgetg(i, t_POL); z[1] = evalsigne(1)|evallgef(i);
for (i--; i>=2; i--) z[i] = lstoi(x[i-1]);
return z;
}
/* x t_POL, n > 0. Return x mod polcyclo(2^n) = (x^(2^(n-1)) + 1). DESTROY x */
static GEN
red_mod_cyclo2n(GEN x, int n)
{
long i, pow2 = 1<<(n-1);
for (i = lgef(x)-1; i>pow2+1; i--)
if (signe(x[i])) x[i-pow2] = lsubii((GEN)x[i-pow2], (GEN)x[i]);
for (; i>1; i--)
if (signe(x[i])) break;
i += 1; setlgef(x,i);
return x;
}
/* x a non-zero VECSMALL */
static GEN
smallpolrev(GEN x)
{
long i,j, lx = lg(x);
GEN y;
while (lx-- && x[lx]==0) /* empty */;
i = lx+2; y = cgetg(i,t_POL);
y[1] = evallgef(i) | evalsigne(1);
for (j=2; j<i; j++) y[j] = lstoi(x[j-1]);
return y;
}
/* x polynomial in t_VECSMALL form, T t_POL return x mod T */
static GEN
u_red(GEN x, GEN T)
{
return gres(smallpolrev(x), T);
}
/* special case R->C = polcyclo(p) */
static GEN
_red2(GEN x, red_t *R) {
return FpX_red(red_mod_cyclo(x, R->n), R->N);
}
static GEN
_red(GEN x, red_t *R) {
return FpX_red(gres(x, R->C), R->N);
}
static GEN
_redsimple(GEN x, red_t *R) {
return modii(x, R->N);
}
static GEN
sqrmod(GEN x, red_t *R) {
return R->red(gsqr(x), R);
}
static GEN
sqrconst(GEN pol, red_t *R)
{
GEN p1 = cgetg(3,t_POL);
p1[2] = (long)modii(sqri((GEN)pol[2]), R->N);
p1[1] = pol[1]; return p1;
}
/* pol^2 mod (x^2+x+1, N) */
static GEN
sqrmod3(GEN pol, red_t *R)
{
GEN a,b,bma,A,B;
long lv=lgef(pol);
if (lv==2) return pol;
if (lv==3) return sqrconst(pol, R);
a = (GEN)pol[3];
b = (GEN)pol[2]; bma = subii(b,a);
A = modii(mulii(a,addii(b,bma)), R->N);
B = modii(mulii(bma,addii(a,b)), R->N);
return makepoldeg1(A,B);
}
/* pol^2 mod (x^2+1,N) */
static GEN
sqrmod4(GEN pol, red_t *R)
{
GEN a,b,A,B;
long lv=lgef(pol);
if (lv==2) return pol;
if (lv==3) return sqrconst(pol, R);
a = (GEN)pol[3];
b = (GEN)pol[2];
A = modii(mulii(a, shifti(b,1)), R->N);
B = modii(mulii(subii(b,a),addii(b,a)), R->N);
return makepoldeg1(A,B);
}
/* pol^2 mod (polcyclo(5),N) */
static GEN
sqrmod5(GEN pol, red_t *R)
{
GEN c2,b,c,d,A,B,C,D;
long lv=lgef(pol);
if (lv==2) return pol;
if (lv==3) return sqrconst(pol, R);
b = (GEN)pol[4];
c = (GEN)pol[3]; c2 = shifti(c,1);
d = (GEN)pol[2];
if (lv==4)
{
A = mulii(b, subii(c2,b));
B = addii(sqri(c), mulii(b, subii(shifti(d,1),b)));
C = subii(mulii(c2,d), sqri(b));
D = mulii(subii(d,b), addii(d,b));
}
else
{
GEN a = (GEN)pol[5], a2 = shifti(a,1);
/* 2a(d - c) + b(2c - b) */
A = addii(mulii(a2, subii(d,c)), mulii(b, subii(c2,b)));
/* c(c - 2a) + b(2d - b) */
B = addii(mulii(c, subii(c,a2)), mulii(b, subii(shifti(d,1),b)));
/* (a-b)(a+b) + 2c(d - a) */
C = addii(mulii(subii(a,b),addii(a,b)), mulii(c2,subii(d,a)));
/* 2a(b - c) + (d-b)(d+b) */
D = addii(mulii(a2, subii(b,c)), mulii(subii(d,b), addii(d,b)));
}
A = modii(A, R->N);
B = modii(B, R->N);
C = modii(C, R->N);
D = modii(D, R->N);
return coefs_to_pol(4,A,B,C,D);
}
static GEN
_mul(GEN x, GEN y, red_t *R) {
return R->red(gmul(x,y), R);
}
static GEN
_powpolmod(int pk, GEN jac, red_t *R, GEN (*_sqr)(GEN, red_t *))
{
const GEN taba = tabaall[pkfalse];
const GEN tabt = tabtall[pkfalse];
const int efin = lg(taba)-1;
GEN res,pol2, *vz;
int lv,tf,f,i;
gpmem_t av;
lv = 1 << (kglob-1);
vz = (GEN*)cgetg(lv+1,t_VEC);
res = lift(jac); pol2 = _sqr(res, R);
vz[1] = res;
for (i=2; i<=lv; i++) vz[i] = _mul(vz[i-1],pol2, R);
av = avma;
for (f=efin; f>=1; f--)
{
res = f==efin ? vz[taba[f]]
: _mul(vz[taba[f]],res, R);
tf = tabt[f];
for (i=1; i<=tf; i++) res = _sqr(res, R);
if ((f&7) == 0) res = gerepilecopy(av, res);
}
return res;
}
static GEN
_powpolmodsimple(red_t *R, int pk, GEN jac)
{
GEN vres,w,wpow,wnew,ma;
long lvres;
int ph,j;
vres = gtovec(lift(jac)); settyp(vres,t_COL);
lvres = lg(vres);
ma = tabmatvite[pkfalse]; ph = lg(ma);
w = cgetg(lvres,t_MAT);
for (j=1; j<lvres; j++) w[j] = ma[j+ph-lvres];
w = gmul(w,vres);
R->red = &_redsimple;
wpow = cgetg(ph,t_COL);
for (j=1; j<ph; j++) wpow[j] = (long)_powpolmod(pk,(GEN)w[j],R,&sqrmod);
wnew = lift(gmul(tabmatinvvite[pkfalse],wpow));
return gtopoly(wnew,0);
}
GEN
powpolmod(red_t *R, int k, int pk, GEN jac)
{
GEN (*_sqr)(GEN, red_t *);
if (!gcmp0(tabmatvite[pkfalse])) return _powpolmodsimple(R, pk, jac);
if (pk == 3) {
R->red = &_red; _sqr = &sqrmod3;
} else if (pk == 4) {
R->red = &_red; _sqr = &sqrmod4;
} else if (pk == 5) {
R->red = &_red; _sqr = &sqrmod5;
} else if (k == 1 && pk >= 7) {
R->n = pk;
R->red = &_red2; _sqr = &sqrmod;
} else {
R->red = &_red; _sqr = &sqrmod;
}
return _powpolmod(pk, jac, R, _sqr);
}
static GEN
e(ulong t)
{
GEN fa,pr,ex,s;
ulong nbd,m,k,d;
int lfa,i,j;
fa = decomp(stoi(t));
pr = (GEN)fa[1]; settyp(pr, t_VECSMALL);
ex = (GEN)fa[2]; settyp(ex, t_VECSMALL); lfa = lg(pr);
nbd = 1;
for (i=1; i<lfa; i++)
{
ex[i] = itos((GEN)ex[i]) + 1;
pr[i] = itos((GEN)pr[i]);
nbd *= ex[i];
}
s = gdeux; /* nbd = number of divisors */
for (k=0; k<nbd; k++)
{
m = k; d = 1;
for (j=1; m; j++)
{
d *= u_pow(pr[j], m % ex[j]);
m /= ex[j];
}
/* d runs through the divisors of t */
if (isprime(stoi(++d))) s = muliu(s, u_pow(d, 1 + u_val(t,d)));
}
return s;
}
static ulong
compt(GEN N)
{
gpmem_t av0 = avma;
ulong Bint,t;
GEN B;
B = mulsr(100, divrr(glog(N,DEFAULTPREC), dbltor(log(10.))));
Bint = itos(gceil(B));
avma = av0;
/* Bint < [200*log_10 e(t)] ==> return t. For e(t) < 10^529, N < 10^1058 */
if (Bint < 540) return 6;
if (Bint < 963) return 12;
if (Bint < 1023) return 24;
if (Bint < 1330) return 48;
if (Bint < 1628) return 36;
if (Bint < 1967) return 60;
if (Bint < 2349) return 120;
if (Bint < 3083) return 180;
if (Bint < 3132) return 240;
if (Bint < 3270) return 504;
if (Bint < 3838) return 360;
if (Bint < 4115) return 420;
if (Bint < 4621) return 720;
if (Bint < 4987) return 840;
if (Bint < 5079) return 1440;
if (Bint < 6212) return 1260;
if (Bint < 6686) return 1680;
if (Bint < 8137) return 2520;
if (Bint < 8415) return 3360;
if (Bint < 10437) return 5040;
if (Bint < 11643) return 13860;
if (Bint < 12826) return 10080;
if (Bint < 11643) return 13860;
if (Bint < 12826) return 10080;
if (Bint < 13369) return 16380;
if (Bint < 13540) return 21840;
if (Bint < 15060) return 18480;
if (Bint < 15934) return 27720;
if (Bint < 17695) return 32760;
if (Bint < 18816) return 36960;
if (Bint < 21338) return 55440;
if (Bint < 23179) return 65520;
if (Bint < 23484) return 98280;
if (Bint < 27465) return 110880;
if (Bint < 30380) return 131040;
if (Bint < 31369) return 166320;
if (Bint < 33866) return 196560;
if (Bint < 34530) return 262080;
if (Bint < 36195) return 277200;
if (Bint < 37095) return 360360;
if (Bint < 38179) return 480480;
if (Bint < 41396) return 332640;
if (Bint < 43301) return 554400;
if (Bint < 47483) return 720720;
if (Bint < 47742) return 665280;
if (Bint < 50202) return 831600;
if (Bint < 52502) return 1113840;
if (Bint < 60245) return 1441440;
if (Bint < 63112) return 1663200;
if (Bint < 65395) return 2227680;
if (Bint < 69895) return 2162160;
if (Bint < 71567) return 2827440;
if (Bint < 75708) return 3326400;
if (Bint < 79377) return 3603600;
if (Bint < 82703) return 6126120;
if (Bint < 91180) return 4324320;
if (Bint < 93978) return 6683040;
if (Bint < 98840) return 7207200;
if (Bint < 99282) return 11138400;
if (Bint < 105811) return 8648640;
B = racine(N);
for (t = 8648640+840;; t+=840)
{
gpmem_t av = avma;
if (cmpii(e(t),B) > 0) break;
avma = av;
}
avma = av0; return t;
}
extern ulong u_gener(ulong p);
/* tabdl[i] = discrete log of i+1 in (Z/q)^*, q odd prime */
static GEN
computetabdl(ulong q)
{
GEN v = cgetg(q-1,t_VECSMALL), w = v-1; /* w[i] = dl(i) */
ulong g,qm3s2,qm1s2,a,i;
g = u_gener(q);
qm3s2 = (q-3)>>1;
qm1s2 = qm3s2+1;
w[q-1] = qm1s2; a = 1;
for (i=1; i<=qm3s2; i++)
{
a = mulssmod(g,a,q);
w[a] = i;
w[q-a] = i+qm1s2;
}
return v;
}
static void
compute_fg(ulong q, int half, GEN *tabf, GEN *tabg)
{
const long qm3s2 = (q-3)>>1;
const long l = half? qm3s2: q-2;
ulong x;
*tabf = computetabdl(q);
*tabg = cgetg(l+1,t_VECSMALL);
for (x=1; x<=l; x++) (*tabg)[x] = (*tabf)[x] + (*tabf)[q-x-1];
}
/* p odd prime */
static GEN
get_jac(ulong q, ulong p, int k, GEN tabf, GEN tabg)
{
GEN ze, vpk;
long i, qm3s2;
ulong x, pk;
pk = u_pow(p,k);
vpk = cgetg(pk+1,t_VECSMALL);
for (i=1; i<=pk; i++) vpk[i] = 0;
ze = tabcyc[pkfalse];
qm3s2 = (q-3)>>1;
for (x=1; x<=qm3s2; x++) vpk[ tabg[x]%pk + 1 ] += 2;
vpk[ (2*tabf[qm3s2+1])%pk + 1 ]++;
return u_red(vpk,ze);
}
/* p = 2 */
static GEN
get_jac2(GEN N, ulong q, int k, GEN *j2q, GEN *j3q)
{
GEN jpq, vpk, tabf, tabg;
long i, qm3s2;
ulong x, pk;
if (k == 1) return NULL;
compute_fg(q,0, &tabf,&tabg);
pk = u_pow(2,k);
vpk = cgetg(pk+1,t_VECSMALL);
for (i=1; i<=pk; i++) vpk[i] = 0;
qm3s2 = (q-3)>>1;
for (x=1; x<=qm3s2; x++) vpk[ tabg[x]%pk + 1 ] += 2;
vpk[ (2*tabf[qm3s2+1])%pk + 1 ]++;
jpq = u_red_mod_cyclo2n(vpk, k);
if (k == 2) return jpq;
if (mod8(N) >= 5)
{
GEN v8 = cgetg(9,t_VECSMALL);
for (x=1; x<=8; x++) v8[x] = 0;
for (x=1; x<=q-2; x++) v8[ ((2*tabf[x]+tabg[x])&7) + 1 ]++;
*j2q = polinflate(gsqr(u_red_mod_cyclo2n(v8,3)), pk>>3);
*j2q = red_mod_cyclo2n(*j2q, k);
}
for (i=1; i<=pk; i++) vpk[i] = 0;
for (x=1; x<=q-2; x++) vpk[ (tabf[x]+tabg[x])%pk + 1 ]++;
*j3q = gmul(jpq, u_red_mod_cyclo2n(vpk,k));
*j3q = red_mod_cyclo2n(*j3q, k);
return jpq;
}
static void
calcjac(GEN et)
{
gpmem_t av;
ulong q, l;
int lfaq,p,k,i,j;
GEN J,tabf,tabg,faq,faqpr,faqex;
globfa = (GEN)decomp(shifti(et,-vali(et)))[1];
l = lg(globfa);
tabfaq= cgetg(l,t_VEC);
tabj = cgetg(l,t_VEC);
for (i=1; i<l; i++)
{
q = itos((GEN)globfa[i]); /* odd prime */
faq = decomp(stoi(q-1)); tabfaq[i] = (long)faq;
faqpr = (GEN)faq[1];
faqex = (GEN)faq[2]; lfaq = lg(faqpr);
av = avma;
compute_fg(q, 1, &tabf, &tabg);
J = cgetg(lfaq,t_VEC);
J[1] = lgetg(1,t_STR); /* dummy */
for (j=2; j<lfaq; j++) /* skip p = faqpr[1] = 2 */
{
p = itos((GEN)faqpr[j]);
k = itos((GEN)faqex[j]);
J[j] = (long)get_jac(q,p,k,tabf,tabg);
}
tabj[i] = (long)gerepilecopy(av, J);
}
}
static void
inittabs(int lv)
{
int i;
tabaall = (GEN*)cgetg(lv,t_VEC);
tabtall = (GEN*)cgetg(lv,t_VEC);
tabcyc = (GEN*)cgetg(lv,t_VEC);
for (i=1; i<lv; i++) tabcyc[i] = gzero;
tabE = cgetg(lv,t_VEC);
tabTH= cgetg(lv,t_VEC);
tabeta=cgetg(lv,t_VEC);
tabmatvite = (GEN*)cgetg(lv,t_VEC);
tabmatinvvite = (GEN*)cgetg(lv,t_VEC);
tabavite = cgetg(lv,t_VEC);
tabpkvite = cgetg(lv,t_VECSMALL);
for (i=1; i<lv; i++)
{
tabE[i] = tabTH[i] = tabeta[i] = tabavite[i] = zero;
tabmatvite[i] = tabmatinvvite[i] = gzero;
}
sgt = cgetg(lv,t_VECSMALL);
ctsgt= cgetg(lv,t_VECSMALL);
for (i=1; i<lv; i++) tabpkvite[i] = sgt[i] = ctsgt[i] = 0;
}
static GEN
finda(GEN N, int pk, int p)
{
int ph,k;
ulong u;
GEN a,b,N1;
if (!isinstep5 && tabpkvite[p])
return gpowgs((GEN)tabavite[p],tabpkvite[p]/pk);
else
{
k = pvaluation(addis(N,-1),stoi(p),&N1);
ph = u_pow(p,k-1); tabpkvite[p] = p*ph;
u = 2;
if (p>2)
{
for(;;u++)
{
a = powgi(gmodulcp(stoi(u),N),N1);
b = gpowgs(a,ph);
if (!gcmp1(b)) break;
}
}
else
{
while(krosg(u,N)!=-1) u++;
a = powgi(gmodulcp(stoi(u),N),N1);
b = gpowgs(a,ph);
}
if (!isinstep5) tabavite[p] = (long)a;
/* Checking b^p = 1 mod N is done economically in filltabs */
b = mppgcd(addis(lift(b),-1),N);
if (!gcmp1(b)) {errfill = b; return NULL;}
return gpowgs(a,(p*ph)/pk);
}
}
static void
filltabs(GEN N, int p, int k, ulong ltab)
{
const ulong mask = (1<<kglob)-1;
gpmem_t av;
int pk,pk2,i,j,LE=0;
ulong s,e;
GEN tabt,taba,m,E=gzero,p1,a=gzero,a2=gzero;
pk = u_pow(p,k);
ishack = isinstep5 && (pk >= lg((GEN)tabcyc) || !signe(tabcyc[pk]));
pk2 = pkfalse;
tabcyc[pk2] = cyclo(pk,0);
p1 = cgetg(pk+1,t_VEC);
for (i=1; i<=pk; i++)
p1[i] = (long)FpX_res(gpowgs(polx[0],i-1), tabcyc[pk2], N);
tabeta[pk2] = (long)p1;
if (p > 2)
{
LE = pk - pk/p + 1; E = cgetg(LE,t_VECSMALL);
for (i=1,j=0; i<pk; i++)
if (i%p) E[++j] = i;
}
else if (k >= 3)
{
LE = (pk>>2) + 1; E = cgetg(LE,t_VECSMALL);
for (i=1,j=0; i<pk; i++)
if ((i%8)==1 || (i%8)==3) E[++j] = i;
}
if (p>2 || k>=3)
{
tabE[pk2] = (long)E;
p1 = cgetg(LE,t_VEC);
for (i=1; i<LE; i++)
{
GEN p2 = cgetg(3,t_VECSMALL);
p2[1] = p2[2] = E[i];
p1[i] = (long)p2;
}
tabTH[pk2] = (long)p1;
}
if (pk > 2 && smodis(N,pk) == 1)
{
int ph,jj;
GEN vpa,p1,p2,p3;
if (DEBUGLEVEL&& !isinstep5) fprintferr("%ld ",pk);
a = finda(N,pk,p); if (!a) return;
ph = pk - pk/p;
vpa = cgetg(ph+1,t_COL); vpa[1] = (long)a;
if (k>1) a2 = gsqr(a);
jj = 1;
for (i=2; i<pk; i++) /* vpa[i] = a^i */
if (i%p) { jj++; vpa[jj] = lmul((i%p==1) ? a2 : a, (GEN)vpa[jj-1]); }
if (jj!=ph) err(bugparier,"filltabs1");
if (!gcmp1(gmul(a,(GEN)vpa[ph])))
{
if (signe(errfill)<=0) errfill = stoi(-1);
return;
}
p1 = cgetg(ph+1,t_MAT);
p2 = cgetg(ph+1,t_COL); p1[ph] = (long)p2;
for (i=1; i<=ph; i++) p2[i] = un;
p1[ph-1] = (long)vpa; p3 = vpa;
for (j=3; j<=ph; j++)
{
p2 = cgetg(ph+1,t_COL); p1[ph+1-j] = (long)p2;
for (i=1; i<=ph; i++) p2[i] = lmul((GEN)vpa[i],(GEN)p3[i]);
p3 = p2;
}
tabmatvite[pk2] = p1;
tabmatinvvite[pk2] = invmat(p1);
}
tabt = cgetg(ltab+1,t_VECSMALL);
taba = cgetg(ltab+1,t_VECSMALL);
av = avma;
m = divis(N,pk);
for (e=1; e<=ltab && signe(m); e++)
{
s = vali(m); m = shifti(m,-s);
tabt[e] = e==1? s: s+kglob;
taba[e] = signe(m)? ((modBIL(m) & mask)+1)>>1: 0;
m = shifti(m, -kglob);
}
avma = av;
if (e > ltab) err(bugparier,"filltabs");
setlg(taba, e); tabaall[pk2] = taba;
setlg(tabt, e); tabtall[pk2] = tabt;
}
static GEN
calcglobs(GEN N, ulong t)
{
GEN fat,fapr,faex;
int lfa,pk,p,ex,i,k;
ulong ltab, ltaball;
NBITSN = bit_accuracy(lgefint(N)) - 1;
while (bittest(N,NBITSN)==0) NBITSN--;
NBITSN++;
kglob=3; while (((kglob+1)*(kglob+2) << (kglob-1)) < NBITSN) kglob++;
ltab = (NBITSN/kglob) + 2;
fat = decomp(stoi(t));
fapr = (GEN)fat[1];
faex = (GEN)fat[2];
lfa = lg(fapr); ltaball = 1;
for (i=1; i<lfa; i++)
{
pk = u_pow(itos((GEN)fapr[i]), itos((GEN)faex[i]));
if (pk > ltaball) ltaball = pk;
}
inittabs(ltaball+1);
if (DEBUGLEVEL) fprintferr("Fast pk-values: ");
for (i=1; i<lfa; i++)
{
p = itos((GEN)fapr[i]);
ex= itos((GEN)faex[i]);
for (k=1; k<=ex; k++)
{
filltabs(N,p,k,ltab);
if (signe(errfill)) break;
}
}
if (DEBUGLEVEL) fprintferr("\n");
return fapr;
}
/* Calculer sig_a^{-1}(z) pour z dans Q(ze) et sig_a: ze -> ze^a */
static GEN
aut(int pk, GEN z,int a)
{
GEN v = cgetg(pk+1,t_VEC);
int i;
for (i=1; i<=pk; i++)
v[i] = (long)polcoeff0(z, (a*(i-1))%pk, 0);
return gmodulcp(gtopolyrev(v,0), tabcyc[pkfalse]);
}
/* calculer z^v pour v dans Z[G] represente par des couples [sig_x^{-1},y] */
static GEN
autvec(int pk, GEN z, GEN v)
{
int i,lv = lg(v);
GEN s = gun;
for (i=1; i<lv; i++)
if (((GEN)v[i])[2])
s = gmul(s, gpowgs(aut(pk,z,((GEN)v[i])[1]), ((GEN)v[i])[2]));
return s;
}
static GEN
get_AL(GEN N, int pk)
{
const int r = smodis(N, pk);
GEN AL, E = (GEN)tabE[pkfalse];
int i, LE = lg(E);
AL = cgetg(LE,t_VEC);
for (i=1; i<LE; i++)
{
GEN p1 = cgetg(3,t_VECSMALL); AL[i] = (long)p1;
p1[1] = E[i];
p1[2] = (r*E[i]) / pk;
}
return AL;
}
static int
look_eta(int pk, GEN s)
{
GEN etats = (GEN)tabeta[pkfalse];
long i;
for (i=1; i<=pk; i++)
if (gegal(s, (GEN)etats[i])) return i-1;
return pk;
}
static int
step4a(GEN N, ulong q, int p, int k, GEN jpq)
{
int pk,pk2,ind;
GEN AL,s1,s2,s3;
red_t R;
pk = u_pow(p,k); pk2=pkfalse;
if (dotime) (void)timer2();
if (!jpq)
{
GEN tabf, tabg;
compute_fg(q,1, &tabf,&tabg);
jpq = get_jac(q,p,k,tabf,tabg);
if (dotime) sgtjac+=timer2();
}
R.N = N;
R.C = tabcyc[pk2];
AL = get_AL(N, pk);
s1 = autvec(pk,jpq,(GEN)tabTH[pk2]);
if (dotime) sgtaut+=timer2();
s2 = powpolmod(&R,k,pk,s1);
if (dotime) {sgt[pk2]+=timer2();ctsgt[pk2]++;};
s3 = autvec(pk,jpq,AL);
if (dotime) sgtaut+=timer2();
s3 = _red(gmul(lift(s3),s2), &R);
if (dotime) sgt[pk2]+=timer2();
ind = look_eta(pk, s3);
if (ind == pk) return -1;
return (ind%p) != 0;
}
/* x == -1 mod N ? */
static int
is_m1(GEN x, GEN N)
{
return egalii(addis(x,1), N);
}
/* p=2, k>=3 */
static int
step4b(GEN N, ulong q, int k)
{
const int pk = u_pow(2,k);
int ind,pk2;
GEN AL,s1,s2,s3, j2q,j3q;
red_t R;
if (dotime) (void)timer2();
(void)get_jac2(N,q,k, &j2q,&j3q);
if (dotime) sgtjac+=timer2();
pk2 = pkfalse;
R.N = N;
R.C = tabcyc[pk2];
AL = get_AL(N, pk);
s1 = autvec(pk,j3q,(GEN)tabTH[pk2]);
if (dotime) sgtaut+=timer2();
s2 = powpolmod(&R, k,pk,s1);
if (dotime) {sgt[pk2]+=timer2();ctsgt[pk2]++;}
s3 = autvec(pk,j3q,AL);
if (dotime) sgtaut+=timer2();
s3 = _red(gmul(lift(s3),s2), &R);
if (mod8(N) >= 5) s3 = _red(gmul(j2q, s3), &R);
if (dotime) sgt[pk2]+=timer2();
ind = look_eta(pk, s3);
if (ind == pk) return -1;
if ((ind&1)==0) return 0;
else
{
s3 = powmodulo(stoi(q), shifti(N,-1), N);
if (dotime) {sgt[pk2]+=timer2();ctsgt[pk2]++;}
return is_m1(s3, N);
}
}
/* p=2, k=2 */
static int
step4c(GEN N, ulong q)
{
const int pk = 4;
int ind,pk2;
GEN s0,s1,s3, jpq = get_jac2(N,q,2, NULL,NULL);
red_t R;
pk2 = pkfalse;
R.N = N;
R.C = tabcyc[pk2];
if (dotime) (void)timer2();
s0 = sqrmod4(jpq, &R);
s1 = gmulsg(q,s0);
s3 = powpolmod(&R, 2,pk,s1);
if (mod4(N) == 3) s3 = _red(gmul(s0,s3), &R);
if (dotime) {sgt[pk2]+=timer2();ctsgt[pk2]++;}
ind = look_eta(pk, s3);
if (ind == pk) return -1;
if ((ind&1)==0) return 0;
else
{
s3 = powmodulo(stoi(q), shifti(N,-1), N);
if (dotime) sgt[pk2]+=timer2();
return is_m1(s3,N);
}
}
/* p=2, k=1 */
static int
step4d(GEN N, ulong q)
{
GEN s1;
if (dotime) (void)timer2();
s1 = powmodulo(negi(stoi(q)), shifti(N,-1), N);
if (dotime) {sgt[2]+=timer2();ctsgt[2]++;}
if (gcmp1(s1)) return 0;
if (is_m1(s1,N)) return (mod4(N) == 1);
return -1;
}
/* return 1 [OK so far],0 [APRCL fails] or < 0 [not a prime] */
static int
step5(GEN N, int p, GEN et)
{
ulong ct = 0, q;
gpmem_t av;
int k, fl = -1;
byteptr d = diffptr+2;
const ulong ltab = (NBITSN/kglob)+2;
av = avma;
for (q = 3; *d; )
{
if (q%p != 1 || smodis(et,q) == 0) goto repeat;
if (smodis(N,q) == 0) return -1;
k = u_val(q-1, p);
filltabs(N,p,k,ltab);
if (p>=3) fl = step4a(N,q,p,k, NULL);
else if (k >= 3) fl = step4b(N,q,k);
else if (k == 2) fl = step4c(N,q);
else fl = step4d(N,q);
if (ishack) tabmatvite[1] = gzero;
if (fl == -1) return (int)(-q);
if (fl == 1) {ctglob = max(ctglob,ct); return 1;}
ct++;
avma = av;
repeat:
NEXT_PRIME_VIADIFF(q,d);
}
return 0;
}
static GEN
step6(GEN N, ulong t, GEN et)
{
GEN N1,r,p1;
ulong i;
gpmem_t av;
if (dotime) (void)timer2();
N1 = resii(N, et);
r = gun; av = avma;
for (i=1; i<t; i++)
{
r = resii(mulii(r,N1), et);
if (!signe(resii(N,r)) && !gcmp1(r) && !egalii(r,N))
{
p1 = cgetg(3,t_VEC);
p1[1] = (long)r;
p1[2] = zero; return p1;
}
if ((i & 0x1f) == 0) r = gerepileuptoint(av, r);
}
if (dotime) sgt6 = timer2();
return gun;
}
static GEN
_res(long a, long b)
{
GEN z;
if (b) { z=cgetg(3, t_VEC); z[1]=lstoi(a); z[2]=lstoi(b); }
else { z=cgetg(2, t_VEC); z[1]=lstoi(a); }
return z;
}
void
seetimes()
{
int i;
ulong s,sc;
s = sc = 0; for (i=1; i<lg(sgt); i++) {s+=sgt[i]; sc+=ctsgt[i];}
printf("Timings in ms:\nJacobi sums = %lu, Galois Automorphisms = %lu\n",sgtjac,sgtaut);
printf("Global Fermat powerings = %lu\n",s);
printf("Number of Fermat powerings = %lu\n",sc);
printf("Individual Fermat powerings = "); output(gtovec(sgt));
printf("Number of individual Fermat powerings = "); output(gtovec(ctsgt));
printf("Final trial divisions = %lu\n",sgt6);
printf("Maximal number of nondeterministic steps = %lu\n",ctglob);
}
GEN
aprcl(GEN N)
{
GEN et,fat,flaglp,faq,faqpr,faqex,res;
long fl;
ulong lfat,p,q,lfaq,k,t,i,j,l;
gpmem_t av,av2,avinit = avma;
if (cmpis(N,12) <= 0)
switch(itos(N))
{
case 2: case 3: case 5: case 7: case 11: return gun;
default: return _res(0,0);
}
ctglob = 0;
if (dotime) (void)timer2();
t = compt(N);
if (DEBUGLEVEL) fprintferr("Choosing t = %ld\n",t);
et = e(t);
if (cmpii(sqri(et),N) < 0) err(bugparier,"aprcl: e(t) too small");
if (!gcmp1(mppgcd(N,mulsi(t,et)))) return _res(1,0);
isinstep5 = 0;
errfill = gzero;
fat = calcglobs(N,t); lfat = lg(fat);
if (signe(errfill))
{
GEN z;
avma = avinit;
if (signe(errfill)<0) return _res(1,0);
else {z=cgetg(3,t_VEC); z[1]=(long)errfill; z[2]=zero; return z;}
}
p = itos((GEN)fat[lfat-1]); /* largest p | t */
flaglp = cgetg(p+1, t_VECSMALL);
flaglp[2] = 0;
for (i=2; i<lfat; i++)
{
p = itos((GEN)fat[i]); q = p*p;
flaglp[p] = (powuumod(smodis(N,q),p-1,q) != 1);
}
calcjac(et);
if (dotime)
{
sgtjac = timer2();
fprintferr("Jacobi sums and tables computed\nq-values: ");
}
sgtaut = 0;
av = avma; l = lg(globfa);
for (i=1; i<l; i++)
{
avma = av;
q = itos((GEN)globfa[i]);
if (dotime) fprintferr("%ld ",q);
faq = (GEN)tabfaq[i];
faqpr = (GEN)faq[1];
faqex = (GEN)faq[2]; lfaq = lg(faqpr);
av2 = avma;
for (j=1; j<lfaq; j++)
{
avma = av2;
p = itos((GEN)faqpr[j]);
k = itos((GEN)faqex[j]);
if (p >= 3) fl = step4a(N,q,p,k, gmael(tabj,i,j));
else if (k >= 3) fl = step4b(N,q,k);
else if (k == 2) fl = step4c(N,q);
else fl = step4d(N,q);
if (fl == -1) return _res(q,p);
if (fl == 1) flaglp[p] = 1;
}
}
if (dotime) fprintferr("\nNormal test done; testing conditions lp\n");
isinstep5 = 1;
for (i=1; i<lfat; i++)
{
p = itos((GEN)fat[i]);
if (flaglp[p] == 0)
{
errfill = gzero;
fl = step5(N,p,et);
if (signe(errfill))
{
GEN z;
avma = avinit;
if (signe(errfill)<0) return _res(1,0);
else {z=cgetg(3,t_VEC); z[1]=(long)errfill; z[2]=zero; return z;}
}
if (fl < 0) return _res(fl,p);
if (fl==0) err(talker,"aprcl test fails! this is highly improbable");
}
}
if (dotime) fprintferr("Conditions lp done, doing step6\n");
res = step6(N,t,et);
if (dotime && (typ(res) == t_INT)) seetimes();
return res;
}
long
isprimeAPRCL(GEN N)
{
gpmem_t av = avma;
GEN res = aprcl(N);
avma = av; return (typ(res) == t_INT);
}
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to aprcl.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 / modules