Annotation of OpenXM_contrib/pari-2.2/src/modules/aprcl.c, Revision 1.1
1.1 ! noro 1: /* $Id: aprcl.c,v 1.23 2002/07/15 13:26:20 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: /* THE APRCL PRIMALITY TEST */
! 19: /* */
! 20: /*******************************************************************/
! 21: #include "pari.h"
! 22:
! 23: static ulong kglob,ctglob,NBITSN,sgtaut,sgt6,sgtjac;
! 24: static int isinstep5,ishack;
! 25: static GEN *tabaall,*tabtall,*tabcyc,tabE,tabTH,tabeta,*tabmatvite,*tabmatinvvite,globfa,tabfaq,tabj,sgt,ctsgt,tabavite,tabpkvite,errfill;
! 26: #define pkfalse (ishack ? 1 : pk)
! 27: #define dotime (DEBUGLEVEL>=3)
! 28:
! 29: typedef struct red_t {
! 30: long n;
! 31: GEN C; /* polcyclo(n) */
! 32: GEN N; /* prime we are certifying */
! 33: GEN (*red)(GEN x, struct red_t *);
! 34: } red_t;
! 35:
! 36: /* powering t_POLMOD mod N, flexible window */
! 37:
! 38: static GEN
! 39: makepoldeg1(GEN c, GEN d)
! 40: {
! 41: GEN res;
! 42: if (signe(c)) {
! 43: res = cgetg(4,t_POL); res[1] = evalsigne(1)|evallgef(4);
! 44: res[2] = (long)d; res[3] = (long)c;
! 45: } else if (signe(d)) {
! 46: res = cgetg(3,t_POL); res[1] = evalsigne(1)|evallgef(3);
! 47: res[2] = (long)d;
! 48: } else {
! 49: res = cgetg(2,t_POL); res[1] = evalsigne(0)|evallgef(2);
! 50: }
! 51: return res;
! 52: }
! 53:
! 54: /* T mod polcyclo(p), assume deg(T) < 2p */
! 55: static GEN
! 56: red_mod_cyclo(GEN T, long p)
! 57: {
! 58: long i, d;
! 59: GEN y, *z;
! 60:
! 61: d = degpol(T) - p; /* < p */
! 62: if (d <= -2) return T;
! 63:
! 64: /* reduce mod (x^p - 1) */
! 65: y = dummycopy(T); z = (GEN*)(y+2);
! 66: for (i = 0; i<=d; i++) z[i] = addii(z[i], z[i+p]);
! 67:
! 68: /* reduce mod x^(p-1) + ... + 1 */
! 69: d = p-1;
! 70: if (signe(z[d]))
! 71: for (i=0; i<d; i++) z[i] = subii(z[i], z[d]);
! 72: return normalizepol_i(y, d+2);
! 73: }
! 74:
! 75: /* x t_VECSMALL, n > 0. Return x mod polcyclo(2^n) = (x^(2^(n-1)) + 1). DESTROY x */
! 76: static GEN
! 77: u_red_mod_cyclo2n(GEN x, int n)
! 78: {
! 79: long i, pow2 = 1<<(n-1);
! 80: GEN z;
! 81:
! 82: for (i = lg(x)-1; i>pow2; i--) x[i-pow2] -= x[i];
! 83: for (; i>0; i--)
! 84: if (x[i]) break;
! 85: i += 2;
! 86: z = cgetg(i, t_POL); z[1] = evalsigne(1)|evallgef(i);
! 87: for (i--; i>=2; i--) z[i] = lstoi(x[i-1]);
! 88: return z;
! 89: }
! 90: /* x t_POL, n > 0. Return x mod polcyclo(2^n) = (x^(2^(n-1)) + 1). DESTROY x */
! 91: static GEN
! 92: red_mod_cyclo2n(GEN x, int n)
! 93: {
! 94: long i, pow2 = 1<<(n-1);
! 95:
! 96: for (i = lgef(x)-1; i>pow2+1; i--)
! 97: if (signe(x[i])) x[i-pow2] = lsubii((GEN)x[i-pow2], (GEN)x[i]);
! 98: for (; i>1; i--)
! 99: if (signe(x[i])) break;
! 100: i += 1; setlgef(x,i);
! 101: return x;
! 102: }
! 103:
! 104: /* x a non-zero VECSMALL */
! 105: static GEN
! 106: smallpolrev(GEN x)
! 107: {
! 108: long i,j, lx = lg(x);
! 109: GEN y;
! 110:
! 111: while (lx-- && x[lx]==0) /* empty */;
! 112: i = lx+2; y = cgetg(i,t_POL);
! 113: y[1] = evallgef(i) | evalsigne(1);
! 114: for (j=2; j<i; j++) y[j] = lstoi(x[j-1]);
! 115: return y;
! 116: }
! 117:
! 118: /* x polynomial in t_VECSMALL form, T t_POL return x mod T */
! 119: static GEN
! 120: u_red(GEN x, GEN T)
! 121: {
! 122: return gres(smallpolrev(x), T);
! 123: }
! 124:
! 125: /* special case R->C = polcyclo(p) */
! 126: static GEN
! 127: _red2(GEN x, red_t *R) {
! 128: return FpX_red(red_mod_cyclo(x, R->n), R->N);
! 129: }
! 130: static GEN
! 131: _red(GEN x, red_t *R) {
! 132: return FpX_red(gres(x, R->C), R->N);
! 133: }
! 134:
! 135: static GEN
! 136: _redsimple(GEN x, red_t *R) {
! 137: return modii(x, R->N);
! 138: }
! 139:
! 140: static GEN
! 141: sqrmod(GEN x, red_t *R) {
! 142: return R->red(gsqr(x), R);
! 143: }
! 144:
! 145: static GEN
! 146: sqrconst(GEN pol, red_t *R)
! 147: {
! 148: GEN p1 = cgetg(3,t_POL);
! 149: p1[2] = (long)modii(sqri((GEN)pol[2]), R->N);
! 150: p1[1] = pol[1]; return p1;
! 151: }
! 152:
! 153: /* pol^2 mod (x^2+x+1, N) */
! 154: static GEN
! 155: sqrmod3(GEN pol, red_t *R)
! 156: {
! 157: GEN a,b,bma,A,B;
! 158: long lv=lgef(pol);
! 159:
! 160: if (lv==2) return pol;
! 161: if (lv==3) return sqrconst(pol, R);
! 162: a = (GEN)pol[3];
! 163: b = (GEN)pol[2]; bma = subii(b,a);
! 164: A = modii(mulii(a,addii(b,bma)), R->N);
! 165: B = modii(mulii(bma,addii(a,b)), R->N);
! 166: return makepoldeg1(A,B);
! 167: }
! 168:
! 169: /* pol^2 mod (x^2+1,N) */
! 170: static GEN
! 171: sqrmod4(GEN pol, red_t *R)
! 172: {
! 173: GEN a,b,A,B;
! 174: long lv=lgef(pol);
! 175:
! 176: if (lv==2) return pol;
! 177: if (lv==3) return sqrconst(pol, R);
! 178: a = (GEN)pol[3];
! 179: b = (GEN)pol[2];
! 180: A = modii(mulii(a, shifti(b,1)), R->N);
! 181: B = modii(mulii(subii(b,a),addii(b,a)), R->N);
! 182: return makepoldeg1(A,B);
! 183: }
! 184:
! 185: /* pol^2 mod (polcyclo(5),N) */
! 186: static GEN
! 187: sqrmod5(GEN pol, red_t *R)
! 188: {
! 189: GEN c2,b,c,d,A,B,C,D;
! 190: long lv=lgef(pol);
! 191:
! 192: if (lv==2) return pol;
! 193: if (lv==3) return sqrconst(pol, R);
! 194: b = (GEN)pol[4];
! 195: c = (GEN)pol[3]; c2 = shifti(c,1);
! 196: d = (GEN)pol[2];
! 197: if (lv==4)
! 198: {
! 199: A = mulii(b, subii(c2,b));
! 200: B = addii(sqri(c), mulii(b, subii(shifti(d,1),b)));
! 201: C = subii(mulii(c2,d), sqri(b));
! 202: D = mulii(subii(d,b), addii(d,b));
! 203: }
! 204: else
! 205: {
! 206: GEN a = (GEN)pol[5], a2 = shifti(a,1);
! 207: /* 2a(d - c) + b(2c - b) */
! 208: A = addii(mulii(a2, subii(d,c)), mulii(b, subii(c2,b)));
! 209: /* c(c - 2a) + b(2d - b) */
! 210: B = addii(mulii(c, subii(c,a2)), mulii(b, subii(shifti(d,1),b)));
! 211: /* (a-b)(a+b) + 2c(d - a) */
! 212: C = addii(mulii(subii(a,b),addii(a,b)), mulii(c2,subii(d,a)));
! 213: /* 2a(b - c) + (d-b)(d+b) */
! 214: D = addii(mulii(a2, subii(b,c)), mulii(subii(d,b), addii(d,b)));
! 215: }
! 216: A = modii(A, R->N);
! 217: B = modii(B, R->N);
! 218: C = modii(C, R->N);
! 219: D = modii(D, R->N);
! 220: return coefs_to_pol(4,A,B,C,D);
! 221: }
! 222:
! 223: static GEN
! 224: _mul(GEN x, GEN y, red_t *R) {
! 225: return R->red(gmul(x,y), R);
! 226: }
! 227:
! 228: static GEN
! 229: _powpolmod(int pk, GEN jac, red_t *R, GEN (*_sqr)(GEN, red_t *))
! 230: {
! 231: const GEN taba = tabaall[pkfalse];
! 232: const GEN tabt = tabtall[pkfalse];
! 233: const int efin = lg(taba)-1;
! 234: GEN res,pol2, *vz;
! 235: int lv,tf,f,i;
! 236: gpmem_t av;
! 237:
! 238: lv = 1 << (kglob-1);
! 239: vz = (GEN*)cgetg(lv+1,t_VEC);
! 240: res = lift(jac); pol2 = _sqr(res, R);
! 241: vz[1] = res;
! 242: for (i=2; i<=lv; i++) vz[i] = _mul(vz[i-1],pol2, R);
! 243: av = avma;
! 244: for (f=efin; f>=1; f--)
! 245: {
! 246: res = f==efin ? vz[taba[f]]
! 247: : _mul(vz[taba[f]],res, R);
! 248: tf = tabt[f];
! 249: for (i=1; i<=tf; i++) res = _sqr(res, R);
! 250: if ((f&7) == 0) res = gerepilecopy(av, res);
! 251: }
! 252: return res;
! 253: }
! 254:
! 255: static GEN
! 256: _powpolmodsimple(red_t *R, int pk, GEN jac)
! 257: {
! 258: GEN vres,w,wpow,wnew,ma;
! 259: long lvres;
! 260: int ph,j;
! 261:
! 262: vres = gtovec(lift(jac)); settyp(vres,t_COL);
! 263: lvres = lg(vres);
! 264: ma = tabmatvite[pkfalse]; ph = lg(ma);
! 265: w = cgetg(lvres,t_MAT);
! 266: for (j=1; j<lvres; j++) w[j] = ma[j+ph-lvres];
! 267: w = gmul(w,vres);
! 268: R->red = &_redsimple;
! 269: wpow = cgetg(ph,t_COL);
! 270: for (j=1; j<ph; j++) wpow[j] = (long)_powpolmod(pk,(GEN)w[j],R,&sqrmod);
! 271: wnew = lift(gmul(tabmatinvvite[pkfalse],wpow));
! 272: return gtopoly(wnew,0);
! 273: }
! 274:
! 275: GEN
! 276: powpolmod(red_t *R, int k, int pk, GEN jac)
! 277: {
! 278: GEN (*_sqr)(GEN, red_t *);
! 279:
! 280: if (!gcmp0(tabmatvite[pkfalse])) return _powpolmodsimple(R, pk, jac);
! 281: if (pk == 3) {
! 282: R->red = &_red; _sqr = &sqrmod3;
! 283: } else if (pk == 4) {
! 284: R->red = &_red; _sqr = &sqrmod4;
! 285: } else if (pk == 5) {
! 286: R->red = &_red; _sqr = &sqrmod5;
! 287: } else if (k == 1 && pk >= 7) {
! 288: R->n = pk;
! 289: R->red = &_red2; _sqr = &sqrmod;
! 290: } else {
! 291: R->red = &_red; _sqr = &sqrmod;
! 292: }
! 293: return _powpolmod(pk, jac, R, _sqr);
! 294: }
! 295:
! 296: static GEN
! 297: e(ulong t)
! 298: {
! 299: GEN fa,pr,ex,s;
! 300: ulong nbd,m,k,d;
! 301: int lfa,i,j;
! 302:
! 303: fa = decomp(stoi(t));
! 304: pr = (GEN)fa[1]; settyp(pr, t_VECSMALL);
! 305: ex = (GEN)fa[2]; settyp(ex, t_VECSMALL); lfa = lg(pr);
! 306: nbd = 1;
! 307: for (i=1; i<lfa; i++)
! 308: {
! 309: ex[i] = itos((GEN)ex[i]) + 1;
! 310: pr[i] = itos((GEN)pr[i]);
! 311: nbd *= ex[i];
! 312: }
! 313: s = gdeux; /* nbd = number of divisors */
! 314: for (k=0; k<nbd; k++)
! 315: {
! 316: m = k; d = 1;
! 317: for (j=1; m; j++)
! 318: {
! 319: d *= u_pow(pr[j], m % ex[j]);
! 320: m /= ex[j];
! 321: }
! 322: /* d runs through the divisors of t */
! 323: if (isprime(stoi(++d))) s = muliu(s, u_pow(d, 1 + u_val(t,d)));
! 324: }
! 325: return s;
! 326: }
! 327:
! 328: static ulong
! 329: compt(GEN N)
! 330: {
! 331: gpmem_t av0 = avma;
! 332: ulong Bint,t;
! 333: GEN B;
! 334:
! 335: B = mulsr(100, divrr(glog(N,DEFAULTPREC), dbltor(log(10.))));
! 336: Bint = itos(gceil(B));
! 337: avma = av0;
! 338: /* Bint < [200*log_10 e(t)] ==> return t. For e(t) < 10^529, N < 10^1058 */
! 339: if (Bint < 540) return 6;
! 340: if (Bint < 963) return 12;
! 341: if (Bint < 1023) return 24;
! 342: if (Bint < 1330) return 48;
! 343: if (Bint < 1628) return 36;
! 344: if (Bint < 1967) return 60;
! 345: if (Bint < 2349) return 120;
! 346: if (Bint < 3083) return 180;
! 347: if (Bint < 3132) return 240;
! 348: if (Bint < 3270) return 504;
! 349: if (Bint < 3838) return 360;
! 350: if (Bint < 4115) return 420;
! 351: if (Bint < 4621) return 720;
! 352: if (Bint < 4987) return 840;
! 353: if (Bint < 5079) return 1440;
! 354: if (Bint < 6212) return 1260;
! 355: if (Bint < 6686) return 1680;
! 356: if (Bint < 8137) return 2520;
! 357: if (Bint < 8415) return 3360;
! 358: if (Bint < 10437) return 5040;
! 359: if (Bint < 11643) return 13860;
! 360: if (Bint < 12826) return 10080;
! 361: if (Bint < 11643) return 13860;
! 362: if (Bint < 12826) return 10080;
! 363: if (Bint < 13369) return 16380;
! 364: if (Bint < 13540) return 21840;
! 365: if (Bint < 15060) return 18480;
! 366: if (Bint < 15934) return 27720;
! 367: if (Bint < 17695) return 32760;
! 368: if (Bint < 18816) return 36960;
! 369: if (Bint < 21338) return 55440;
! 370: if (Bint < 23179) return 65520;
! 371: if (Bint < 23484) return 98280;
! 372: if (Bint < 27465) return 110880;
! 373: if (Bint < 30380) return 131040;
! 374: if (Bint < 31369) return 166320;
! 375: if (Bint < 33866) return 196560;
! 376: if (Bint < 34530) return 262080;
! 377: if (Bint < 36195) return 277200;
! 378: if (Bint < 37095) return 360360;
! 379: if (Bint < 38179) return 480480;
! 380: if (Bint < 41396) return 332640;
! 381: if (Bint < 43301) return 554400;
! 382: if (Bint < 47483) return 720720;
! 383: if (Bint < 47742) return 665280;
! 384: if (Bint < 50202) return 831600;
! 385: if (Bint < 52502) return 1113840;
! 386: if (Bint < 60245) return 1441440;
! 387: if (Bint < 63112) return 1663200;
! 388: if (Bint < 65395) return 2227680;
! 389: if (Bint < 69895) return 2162160;
! 390: if (Bint < 71567) return 2827440;
! 391: if (Bint < 75708) return 3326400;
! 392: if (Bint < 79377) return 3603600;
! 393: if (Bint < 82703) return 6126120;
! 394: if (Bint < 91180) return 4324320;
! 395: if (Bint < 93978) return 6683040;
! 396: if (Bint < 98840) return 7207200;
! 397: if (Bint < 99282) return 11138400;
! 398: if (Bint < 105811) return 8648640;
! 399:
! 400: B = racine(N);
! 401: for (t = 8648640+840;; t+=840)
! 402: {
! 403: gpmem_t av = avma;
! 404: if (cmpii(e(t),B) > 0) break;
! 405: avma = av;
! 406: }
! 407: avma = av0; return t;
! 408: }
! 409:
! 410: extern ulong u_gener(ulong p);
! 411:
! 412: /* tabdl[i] = discrete log of i+1 in (Z/q)^*, q odd prime */
! 413: static GEN
! 414: computetabdl(ulong q)
! 415: {
! 416: GEN v = cgetg(q-1,t_VECSMALL), w = v-1; /* w[i] = dl(i) */
! 417: ulong g,qm3s2,qm1s2,a,i;
! 418:
! 419: g = u_gener(q);
! 420: qm3s2 = (q-3)>>1;
! 421: qm1s2 = qm3s2+1;
! 422: w[q-1] = qm1s2; a = 1;
! 423: for (i=1; i<=qm3s2; i++)
! 424: {
! 425: a = mulssmod(g,a,q);
! 426: w[a] = i;
! 427: w[q-a] = i+qm1s2;
! 428: }
! 429: return v;
! 430: }
! 431:
! 432: static void
! 433: compute_fg(ulong q, int half, GEN *tabf, GEN *tabg)
! 434: {
! 435: const long qm3s2 = (q-3)>>1;
! 436: const long l = half? qm3s2: q-2;
! 437: ulong x;
! 438: *tabf = computetabdl(q);
! 439: *tabg = cgetg(l+1,t_VECSMALL);
! 440: for (x=1; x<=l; x++) (*tabg)[x] = (*tabf)[x] + (*tabf)[q-x-1];
! 441: }
! 442:
! 443: /* p odd prime */
! 444: static GEN
! 445: get_jac(ulong q, ulong p, int k, GEN tabf, GEN tabg)
! 446: {
! 447: GEN ze, vpk;
! 448: long i, qm3s2;
! 449: ulong x, pk;
! 450:
! 451: pk = u_pow(p,k);
! 452: vpk = cgetg(pk+1,t_VECSMALL);
! 453: for (i=1; i<=pk; i++) vpk[i] = 0;
! 454: ze = tabcyc[pkfalse];
! 455:
! 456: qm3s2 = (q-3)>>1;
! 457: for (x=1; x<=qm3s2; x++) vpk[ tabg[x]%pk + 1 ] += 2;
! 458: vpk[ (2*tabf[qm3s2+1])%pk + 1 ]++;
! 459: return u_red(vpk,ze);
! 460: }
! 461:
! 462: /* p = 2 */
! 463: static GEN
! 464: get_jac2(GEN N, ulong q, int k, GEN *j2q, GEN *j3q)
! 465: {
! 466: GEN jpq, vpk, tabf, tabg;
! 467: long i, qm3s2;
! 468: ulong x, pk;
! 469:
! 470: if (k == 1) return NULL;
! 471:
! 472: compute_fg(q,0, &tabf,&tabg);
! 473:
! 474: pk = u_pow(2,k);
! 475: vpk = cgetg(pk+1,t_VECSMALL);
! 476: for (i=1; i<=pk; i++) vpk[i] = 0;
! 477:
! 478: qm3s2 = (q-3)>>1;
! 479: for (x=1; x<=qm3s2; x++) vpk[ tabg[x]%pk + 1 ] += 2;
! 480: vpk[ (2*tabf[qm3s2+1])%pk + 1 ]++;
! 481: jpq = u_red_mod_cyclo2n(vpk, k);
! 482:
! 483: if (k == 2) return jpq;
! 484:
! 485: if (mod8(N) >= 5)
! 486: {
! 487: GEN v8 = cgetg(9,t_VECSMALL);
! 488: for (x=1; x<=8; x++) v8[x] = 0;
! 489: for (x=1; x<=q-2; x++) v8[ ((2*tabf[x]+tabg[x])&7) + 1 ]++;
! 490: *j2q = polinflate(gsqr(u_red_mod_cyclo2n(v8,3)), pk>>3);
! 491: *j2q = red_mod_cyclo2n(*j2q, k);
! 492: }
! 493:
! 494: for (i=1; i<=pk; i++) vpk[i] = 0;
! 495: for (x=1; x<=q-2; x++) vpk[ (tabf[x]+tabg[x])%pk + 1 ]++;
! 496: *j3q = gmul(jpq, u_red_mod_cyclo2n(vpk,k));
! 497: *j3q = red_mod_cyclo2n(*j3q, k);
! 498: return jpq;
! 499: }
! 500:
! 501: static void
! 502: calcjac(GEN et)
! 503: {
! 504: gpmem_t av;
! 505: ulong q, l;
! 506: int lfaq,p,k,i,j;
! 507: GEN J,tabf,tabg,faq,faqpr,faqex;
! 508:
! 509: globfa = (GEN)decomp(shifti(et,-vali(et)))[1];
! 510: l = lg(globfa);
! 511: tabfaq= cgetg(l,t_VEC);
! 512: tabj = cgetg(l,t_VEC);
! 513: for (i=1; i<l; i++)
! 514: {
! 515: q = itos((GEN)globfa[i]); /* odd prime */
! 516: faq = decomp(stoi(q-1)); tabfaq[i] = (long)faq;
! 517: faqpr = (GEN)faq[1];
! 518: faqex = (GEN)faq[2]; lfaq = lg(faqpr);
! 519:
! 520: av = avma;
! 521: compute_fg(q, 1, &tabf, &tabg);
! 522:
! 523: J = cgetg(lfaq,t_VEC);
! 524: J[1] = lgetg(1,t_STR); /* dummy */
! 525: for (j=2; j<lfaq; j++) /* skip p = faqpr[1] = 2 */
! 526: {
! 527: p = itos((GEN)faqpr[j]);
! 528: k = itos((GEN)faqex[j]);
! 529: J[j] = (long)get_jac(q,p,k,tabf,tabg);
! 530: }
! 531: tabj[i] = (long)gerepilecopy(av, J);
! 532: }
! 533: }
! 534:
! 535: static void
! 536: inittabs(int lv)
! 537: {
! 538: int i;
! 539: tabaall = (GEN*)cgetg(lv,t_VEC);
! 540: tabtall = (GEN*)cgetg(lv,t_VEC);
! 541: tabcyc = (GEN*)cgetg(lv,t_VEC);
! 542: for (i=1; i<lv; i++) tabcyc[i] = gzero;
! 543: tabE = cgetg(lv,t_VEC);
! 544: tabTH= cgetg(lv,t_VEC);
! 545: tabeta=cgetg(lv,t_VEC);
! 546: tabmatvite = (GEN*)cgetg(lv,t_VEC);
! 547: tabmatinvvite = (GEN*)cgetg(lv,t_VEC);
! 548: tabavite = cgetg(lv,t_VEC);
! 549: tabpkvite = cgetg(lv,t_VECSMALL);
! 550: for (i=1; i<lv; i++)
! 551: {
! 552: tabE[i] = tabTH[i] = tabeta[i] = tabavite[i] = zero;
! 553: tabmatvite[i] = tabmatinvvite[i] = gzero;
! 554: }
! 555: sgt = cgetg(lv,t_VECSMALL);
! 556: ctsgt= cgetg(lv,t_VECSMALL);
! 557: for (i=1; i<lv; i++) tabpkvite[i] = sgt[i] = ctsgt[i] = 0;
! 558: }
! 559:
! 560: static GEN
! 561: finda(GEN N, int pk, int p)
! 562: {
! 563: int ph,k;
! 564: ulong u;
! 565: GEN a,b,N1;
! 566:
! 567: if (!isinstep5 && tabpkvite[p])
! 568: return gpowgs((GEN)tabavite[p],tabpkvite[p]/pk);
! 569: else
! 570: {
! 571: k = pvaluation(addis(N,-1),stoi(p),&N1);
! 572: ph = u_pow(p,k-1); tabpkvite[p] = p*ph;
! 573: u = 2;
! 574: if (p>2)
! 575: {
! 576: for(;;u++)
! 577: {
! 578: a = powgi(gmodulcp(stoi(u),N),N1);
! 579: b = gpowgs(a,ph);
! 580: if (!gcmp1(b)) break;
! 581: }
! 582: }
! 583: else
! 584: {
! 585: while(krosg(u,N)!=-1) u++;
! 586: a = powgi(gmodulcp(stoi(u),N),N1);
! 587: b = gpowgs(a,ph);
! 588: }
! 589: if (!isinstep5) tabavite[p] = (long)a;
! 590: /* Checking b^p = 1 mod N is done economically in filltabs */
! 591: b = mppgcd(addis(lift(b),-1),N);
! 592: if (!gcmp1(b)) {errfill = b; return NULL;}
! 593: return gpowgs(a,(p*ph)/pk);
! 594: }
! 595: }
! 596:
! 597: static void
! 598: filltabs(GEN N, int p, int k, ulong ltab)
! 599: {
! 600: const ulong mask = (1<<kglob)-1;
! 601: gpmem_t av;
! 602: int pk,pk2,i,j,LE=0;
! 603: ulong s,e;
! 604: GEN tabt,taba,m,E=gzero,p1,a=gzero,a2=gzero;
! 605:
! 606: pk = u_pow(p,k);
! 607: ishack = isinstep5 && (pk >= lg((GEN)tabcyc) || !signe(tabcyc[pk]));
! 608: pk2 = pkfalse;
! 609: tabcyc[pk2] = cyclo(pk,0);
! 610:
! 611: p1 = cgetg(pk+1,t_VEC);
! 612: for (i=1; i<=pk; i++)
! 613: p1[i] = (long)FpX_res(gpowgs(polx[0],i-1), tabcyc[pk2], N);
! 614: tabeta[pk2] = (long)p1;
! 615:
! 616: if (p > 2)
! 617: {
! 618: LE = pk - pk/p + 1; E = cgetg(LE,t_VECSMALL);
! 619: for (i=1,j=0; i<pk; i++)
! 620: if (i%p) E[++j] = i;
! 621: }
! 622: else if (k >= 3)
! 623: {
! 624: LE = (pk>>2) + 1; E = cgetg(LE,t_VECSMALL);
! 625: for (i=1,j=0; i<pk; i++)
! 626: if ((i%8)==1 || (i%8)==3) E[++j] = i;
! 627: }
! 628: if (p>2 || k>=3)
! 629: {
! 630: tabE[pk2] = (long)E;
! 631: p1 = cgetg(LE,t_VEC);
! 632: for (i=1; i<LE; i++)
! 633: {
! 634: GEN p2 = cgetg(3,t_VECSMALL);
! 635: p2[1] = p2[2] = E[i];
! 636: p1[i] = (long)p2;
! 637: }
! 638: tabTH[pk2] = (long)p1;
! 639: }
! 640:
! 641: if (pk > 2 && smodis(N,pk) == 1)
! 642: {
! 643: int ph,jj;
! 644: GEN vpa,p1,p2,p3;
! 645:
! 646: if (DEBUGLEVEL&& !isinstep5) fprintferr("%ld ",pk);
! 647: a = finda(N,pk,p); if (!a) return;
! 648: ph = pk - pk/p;
! 649: vpa = cgetg(ph+1,t_COL); vpa[1] = (long)a;
! 650: if (k>1) a2 = gsqr(a);
! 651: jj = 1;
! 652: for (i=2; i<pk; i++) /* vpa[i] = a^i */
! 653: if (i%p) { jj++; vpa[jj] = lmul((i%p==1) ? a2 : a, (GEN)vpa[jj-1]); }
! 654: if (jj!=ph) err(bugparier,"filltabs1");
! 655: if (!gcmp1(gmul(a,(GEN)vpa[ph])))
! 656: {
! 657: if (signe(errfill)<=0) errfill = stoi(-1);
! 658: return;
! 659: }
! 660: p1 = cgetg(ph+1,t_MAT);
! 661: p2 = cgetg(ph+1,t_COL); p1[ph] = (long)p2;
! 662: for (i=1; i<=ph; i++) p2[i] = un;
! 663: p1[ph-1] = (long)vpa; p3 = vpa;
! 664: for (j=3; j<=ph; j++)
! 665: {
! 666: p2 = cgetg(ph+1,t_COL); p1[ph+1-j] = (long)p2;
! 667: for (i=1; i<=ph; i++) p2[i] = lmul((GEN)vpa[i],(GEN)p3[i]);
! 668: p3 = p2;
! 669: }
! 670: tabmatvite[pk2] = p1;
! 671: tabmatinvvite[pk2] = invmat(p1);
! 672: }
! 673:
! 674: tabt = cgetg(ltab+1,t_VECSMALL);
! 675: taba = cgetg(ltab+1,t_VECSMALL);
! 676: av = avma;
! 677: m = divis(N,pk);
! 678: for (e=1; e<=ltab && signe(m); e++)
! 679: {
! 680: s = vali(m); m = shifti(m,-s);
! 681: tabt[e] = e==1? s: s+kglob;
! 682: taba[e] = signe(m)? ((modBIL(m) & mask)+1)>>1: 0;
! 683: m = shifti(m, -kglob);
! 684: }
! 685: avma = av;
! 686: if (e > ltab) err(bugparier,"filltabs");
! 687: setlg(taba, e); tabaall[pk2] = taba;
! 688: setlg(tabt, e); tabtall[pk2] = tabt;
! 689: }
! 690:
! 691:
! 692: static GEN
! 693: calcglobs(GEN N, ulong t)
! 694: {
! 695: GEN fat,fapr,faex;
! 696: int lfa,pk,p,ex,i,k;
! 697: ulong ltab, ltaball;
! 698:
! 699: NBITSN = bit_accuracy(lgefint(N)) - 1;
! 700: while (bittest(N,NBITSN)==0) NBITSN--;
! 701: NBITSN++;
! 702: kglob=3; while (((kglob+1)*(kglob+2) << (kglob-1)) < NBITSN) kglob++;
! 703: ltab = (NBITSN/kglob) + 2;
! 704:
! 705: fat = decomp(stoi(t));
! 706: fapr = (GEN)fat[1];
! 707: faex = (GEN)fat[2];
! 708: lfa = lg(fapr); ltaball = 1;
! 709: for (i=1; i<lfa; i++)
! 710: {
! 711: pk = u_pow(itos((GEN)fapr[i]), itos((GEN)faex[i]));
! 712: if (pk > ltaball) ltaball = pk;
! 713: }
! 714: inittabs(ltaball+1);
! 715: if (DEBUGLEVEL) fprintferr("Fast pk-values: ");
! 716: for (i=1; i<lfa; i++)
! 717: {
! 718: p = itos((GEN)fapr[i]);
! 719: ex= itos((GEN)faex[i]);
! 720: for (k=1; k<=ex; k++)
! 721: {
! 722: filltabs(N,p,k,ltab);
! 723: if (signe(errfill)) break;
! 724: }
! 725: }
! 726: if (DEBUGLEVEL) fprintferr("\n");
! 727: return fapr;
! 728: }
! 729:
! 730: /* Calculer sig_a^{-1}(z) pour z dans Q(ze) et sig_a: ze -> ze^a */
! 731: static GEN
! 732: aut(int pk, GEN z,int a)
! 733: {
! 734: GEN v = cgetg(pk+1,t_VEC);
! 735: int i;
! 736: for (i=1; i<=pk; i++)
! 737: v[i] = (long)polcoeff0(z, (a*(i-1))%pk, 0);
! 738: return gmodulcp(gtopolyrev(v,0), tabcyc[pkfalse]);
! 739: }
! 740:
! 741: /* calculer z^v pour v dans Z[G] represente par des couples [sig_x^{-1},y] */
! 742: static GEN
! 743: autvec(int pk, GEN z, GEN v)
! 744: {
! 745: int i,lv = lg(v);
! 746: GEN s = gun;
! 747: for (i=1; i<lv; i++)
! 748: if (((GEN)v[i])[2])
! 749: s = gmul(s, gpowgs(aut(pk,z,((GEN)v[i])[1]), ((GEN)v[i])[2]));
! 750: return s;
! 751: }
! 752:
! 753: static GEN
! 754: get_AL(GEN N, int pk)
! 755: {
! 756: const int r = smodis(N, pk);
! 757: GEN AL, E = (GEN)tabE[pkfalse];
! 758: int i, LE = lg(E);
! 759:
! 760: AL = cgetg(LE,t_VEC);
! 761: for (i=1; i<LE; i++)
! 762: {
! 763: GEN p1 = cgetg(3,t_VECSMALL); AL[i] = (long)p1;
! 764: p1[1] = E[i];
! 765: p1[2] = (r*E[i]) / pk;
! 766: }
! 767: return AL;
! 768: }
! 769:
! 770: static int
! 771: look_eta(int pk, GEN s)
! 772: {
! 773: GEN etats = (GEN)tabeta[pkfalse];
! 774: long i;
! 775: for (i=1; i<=pk; i++)
! 776: if (gegal(s, (GEN)etats[i])) return i-1;
! 777: return pk;
! 778: }
! 779:
! 780: static int
! 781: step4a(GEN N, ulong q, int p, int k, GEN jpq)
! 782: {
! 783: int pk,pk2,ind;
! 784: GEN AL,s1,s2,s3;
! 785: red_t R;
! 786:
! 787: pk = u_pow(p,k); pk2=pkfalse;
! 788: if (dotime) (void)timer2();
! 789: if (!jpq)
! 790: {
! 791: GEN tabf, tabg;
! 792:
! 793: compute_fg(q,1, &tabf,&tabg);
! 794: jpq = get_jac(q,p,k,tabf,tabg);
! 795: if (dotime) sgtjac+=timer2();
! 796: }
! 797: R.N = N;
! 798: R.C = tabcyc[pk2];
! 799: AL = get_AL(N, pk);
! 800:
! 801: s1 = autvec(pk,jpq,(GEN)tabTH[pk2]);
! 802: if (dotime) sgtaut+=timer2();
! 803: s2 = powpolmod(&R,k,pk,s1);
! 804: if (dotime) {sgt[pk2]+=timer2();ctsgt[pk2]++;};
! 805: s3 = autvec(pk,jpq,AL);
! 806: if (dotime) sgtaut+=timer2();
! 807: s3 = _red(gmul(lift(s3),s2), &R);
! 808: if (dotime) sgt[pk2]+=timer2();
! 809:
! 810: ind = look_eta(pk, s3);
! 811: if (ind == pk) return -1;
! 812: return (ind%p) != 0;
! 813: }
! 814:
! 815: /* x == -1 mod N ? */
! 816: static int
! 817: is_m1(GEN x, GEN N)
! 818: {
! 819: return egalii(addis(x,1), N);
! 820: }
! 821:
! 822: /* p=2, k>=3 */
! 823: static int
! 824: step4b(GEN N, ulong q, int k)
! 825: {
! 826: const int pk = u_pow(2,k);
! 827: int ind,pk2;
! 828: GEN AL,s1,s2,s3, j2q,j3q;
! 829: red_t R;
! 830:
! 831: if (dotime) (void)timer2();
! 832: (void)get_jac2(N,q,k, &j2q,&j3q);
! 833: if (dotime) sgtjac+=timer2();
! 834:
! 835: pk2 = pkfalse;
! 836: R.N = N;
! 837: R.C = tabcyc[pk2];
! 838: AL = get_AL(N, pk);
! 839:
! 840: s1 = autvec(pk,j3q,(GEN)tabTH[pk2]);
! 841: if (dotime) sgtaut+=timer2();
! 842: s2 = powpolmod(&R, k,pk,s1);
! 843: if (dotime) {sgt[pk2]+=timer2();ctsgt[pk2]++;}
! 844: s3 = autvec(pk,j3q,AL);
! 845: if (dotime) sgtaut+=timer2();
! 846: s3 = _red(gmul(lift(s3),s2), &R);
! 847: if (mod8(N) >= 5) s3 = _red(gmul(j2q, s3), &R);
! 848: if (dotime) sgt[pk2]+=timer2();
! 849:
! 850: ind = look_eta(pk, s3);
! 851: if (ind == pk) return -1;
! 852: if ((ind&1)==0) return 0;
! 853: else
! 854: {
! 855: s3 = powmodulo(stoi(q), shifti(N,-1), N);
! 856: if (dotime) {sgt[pk2]+=timer2();ctsgt[pk2]++;}
! 857: return is_m1(s3, N);
! 858: }
! 859: }
! 860:
! 861: /* p=2, k=2 */
! 862: static int
! 863: step4c(GEN N, ulong q)
! 864: {
! 865: const int pk = 4;
! 866: int ind,pk2;
! 867: GEN s0,s1,s3, jpq = get_jac2(N,q,2, NULL,NULL);
! 868: red_t R;
! 869:
! 870: pk2 = pkfalse;
! 871: R.N = N;
! 872: R.C = tabcyc[pk2];
! 873:
! 874: if (dotime) (void)timer2();
! 875: s0 = sqrmod4(jpq, &R);
! 876: s1 = gmulsg(q,s0);
! 877: s3 = powpolmod(&R, 2,pk,s1);
! 878: if (mod4(N) == 3) s3 = _red(gmul(s0,s3), &R);
! 879: if (dotime) {sgt[pk2]+=timer2();ctsgt[pk2]++;}
! 880:
! 881: ind = look_eta(pk, s3);
! 882: if (ind == pk) return -1;
! 883: if ((ind&1)==0) return 0;
! 884: else
! 885: {
! 886: s3 = powmodulo(stoi(q), shifti(N,-1), N);
! 887: if (dotime) sgt[pk2]+=timer2();
! 888: return is_m1(s3,N);
! 889: }
! 890: }
! 891:
! 892: /* p=2, k=1 */
! 893: static int
! 894: step4d(GEN N, ulong q)
! 895: {
! 896: GEN s1;
! 897:
! 898: if (dotime) (void)timer2();
! 899: s1 = powmodulo(negi(stoi(q)), shifti(N,-1), N);
! 900: if (dotime) {sgt[2]+=timer2();ctsgt[2]++;}
! 901: if (gcmp1(s1)) return 0;
! 902: if (is_m1(s1,N)) return (mod4(N) == 1);
! 903: return -1;
! 904: }
! 905:
! 906: /* return 1 [OK so far],0 [APRCL fails] or < 0 [not a prime] */
! 907: static int
! 908: step5(GEN N, int p, GEN et)
! 909: {
! 910: ulong ct = 0, q;
! 911: gpmem_t av;
! 912: int k, fl = -1;
! 913: byteptr d = diffptr+2;
! 914: const ulong ltab = (NBITSN/kglob)+2;
! 915:
! 916: av = avma;
! 917: for (q = 3; *d; )
! 918: {
! 919: if (q%p != 1 || smodis(et,q) == 0) goto repeat;
! 920:
! 921: if (smodis(N,q) == 0) return -1;
! 922: k = u_val(q-1, p);
! 923: filltabs(N,p,k,ltab);
! 924:
! 925: if (p>=3) fl = step4a(N,q,p,k, NULL);
! 926: else if (k >= 3) fl = step4b(N,q,k);
! 927: else if (k == 2) fl = step4c(N,q);
! 928: else fl = step4d(N,q);
! 929: if (ishack) tabmatvite[1] = gzero;
! 930: if (fl == -1) return (int)(-q);
! 931: if (fl == 1) {ctglob = max(ctglob,ct); return 1;}
! 932: ct++;
! 933: avma = av;
! 934: repeat:
! 935: NEXT_PRIME_VIADIFF(q,d);
! 936: }
! 937: return 0;
! 938: }
! 939:
! 940: static GEN
! 941: step6(GEN N, ulong t, GEN et)
! 942: {
! 943: GEN N1,r,p1;
! 944: ulong i;
! 945: gpmem_t av;
! 946:
! 947: if (dotime) (void)timer2();
! 948: N1 = resii(N, et);
! 949: r = gun; av = avma;
! 950: for (i=1; i<t; i++)
! 951: {
! 952: r = resii(mulii(r,N1), et);
! 953: if (!signe(resii(N,r)) && !gcmp1(r) && !egalii(r,N))
! 954: {
! 955: p1 = cgetg(3,t_VEC);
! 956: p1[1] = (long)r;
! 957: p1[2] = zero; return p1;
! 958: }
! 959: if ((i & 0x1f) == 0) r = gerepileuptoint(av, r);
! 960: }
! 961: if (dotime) sgt6 = timer2();
! 962: return gun;
! 963: }
! 964:
! 965: static GEN
! 966: _res(long a, long b)
! 967: {
! 968: GEN z;
! 969: if (b) { z=cgetg(3, t_VEC); z[1]=lstoi(a); z[2]=lstoi(b); }
! 970: else { z=cgetg(2, t_VEC); z[1]=lstoi(a); }
! 971: return z;
! 972: }
! 973:
! 974: void
! 975: seetimes()
! 976: {
! 977: int i;
! 978: ulong s,sc;
! 979:
! 980: s = sc = 0; for (i=1; i<lg(sgt); i++) {s+=sgt[i]; sc+=ctsgt[i];}
! 981: printf("Timings in ms:\nJacobi sums = %lu, Galois Automorphisms = %lu\n",sgtjac,sgtaut);
! 982: printf("Global Fermat powerings = %lu\n",s);
! 983: printf("Number of Fermat powerings = %lu\n",sc);
! 984: printf("Individual Fermat powerings = "); output(gtovec(sgt));
! 985: printf("Number of individual Fermat powerings = "); output(gtovec(ctsgt));
! 986: printf("Final trial divisions = %lu\n",sgt6);
! 987: printf("Maximal number of nondeterministic steps = %lu\n",ctglob);
! 988: }
! 989:
! 990: GEN
! 991: aprcl(GEN N)
! 992: {
! 993: GEN et,fat,flaglp,faq,faqpr,faqex,res;
! 994: long fl;
! 995: ulong lfat,p,q,lfaq,k,t,i,j,l;
! 996: gpmem_t av,av2,avinit = avma;
! 997:
! 998: if (cmpis(N,12) <= 0)
! 999: switch(itos(N))
! 1000: {
! 1001: case 2: case 3: case 5: case 7: case 11: return gun;
! 1002: default: return _res(0,0);
! 1003: }
! 1004: ctglob = 0;
! 1005: if (dotime) (void)timer2();
! 1006: t = compt(N);
! 1007: if (DEBUGLEVEL) fprintferr("Choosing t = %ld\n",t);
! 1008: et = e(t);
! 1009: if (cmpii(sqri(et),N) < 0) err(bugparier,"aprcl: e(t) too small");
! 1010: if (!gcmp1(mppgcd(N,mulsi(t,et)))) return _res(1,0);
! 1011:
! 1012: isinstep5 = 0;
! 1013: errfill = gzero;
! 1014: fat = calcglobs(N,t); lfat = lg(fat);
! 1015: if (signe(errfill))
! 1016: {
! 1017: GEN z;
! 1018:
! 1019: avma = avinit;
! 1020: if (signe(errfill)<0) return _res(1,0);
! 1021: else {z=cgetg(3,t_VEC); z[1]=(long)errfill; z[2]=zero; return z;}
! 1022: }
! 1023: p = itos((GEN)fat[lfat-1]); /* largest p | t */
! 1024: flaglp = cgetg(p+1, t_VECSMALL);
! 1025: flaglp[2] = 0;
! 1026: for (i=2; i<lfat; i++)
! 1027: {
! 1028: p = itos((GEN)fat[i]); q = p*p;
! 1029: flaglp[p] = (powuumod(smodis(N,q),p-1,q) != 1);
! 1030: }
! 1031: calcjac(et);
! 1032: if (dotime)
! 1033: {
! 1034: sgtjac = timer2();
! 1035: fprintferr("Jacobi sums and tables computed\nq-values: ");
! 1036: }
! 1037: sgtaut = 0;
! 1038: av = avma; l = lg(globfa);
! 1039: for (i=1; i<l; i++)
! 1040: {
! 1041: avma = av;
! 1042: q = itos((GEN)globfa[i]);
! 1043: if (dotime) fprintferr("%ld ",q);
! 1044: faq = (GEN)tabfaq[i];
! 1045: faqpr = (GEN)faq[1];
! 1046: faqex = (GEN)faq[2]; lfaq = lg(faqpr);
! 1047: av2 = avma;
! 1048: for (j=1; j<lfaq; j++)
! 1049: {
! 1050: avma = av2;
! 1051: p = itos((GEN)faqpr[j]);
! 1052: k = itos((GEN)faqex[j]);
! 1053: if (p >= 3) fl = step4a(N,q,p,k, gmael(tabj,i,j));
! 1054: else if (k >= 3) fl = step4b(N,q,k);
! 1055: else if (k == 2) fl = step4c(N,q);
! 1056: else fl = step4d(N,q);
! 1057: if (fl == -1) return _res(q,p);
! 1058: if (fl == 1) flaglp[p] = 1;
! 1059: }
! 1060: }
! 1061: if (dotime) fprintferr("\nNormal test done; testing conditions lp\n");
! 1062: isinstep5 = 1;
! 1063: for (i=1; i<lfat; i++)
! 1064: {
! 1065: p = itos((GEN)fat[i]);
! 1066: if (flaglp[p] == 0)
! 1067: {
! 1068: errfill = gzero;
! 1069: fl = step5(N,p,et);
! 1070: if (signe(errfill))
! 1071: {
! 1072: GEN z;
! 1073:
! 1074: avma = avinit;
! 1075: if (signe(errfill)<0) return _res(1,0);
! 1076: else {z=cgetg(3,t_VEC); z[1]=(long)errfill; z[2]=zero; return z;}
! 1077: }
! 1078: if (fl < 0) return _res(fl,p);
! 1079: if (fl==0) err(talker,"aprcl test fails! this is highly improbable");
! 1080: }
! 1081: }
! 1082: if (dotime) fprintferr("Conditions lp done, doing step6\n");
! 1083: res = step6(N,t,et);
! 1084: if (dotime && (typ(res) == t_INT)) seetimes();
! 1085: return res;
! 1086: }
! 1087:
! 1088: long
! 1089: isprimeAPRCL(GEN N)
! 1090: {
! 1091: gpmem_t av = avma;
! 1092: GEN res = aprcl(N);
! 1093: avma = av; return (typ(res) == t_INT);
! 1094: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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