![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gf.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / builtin|
Return to gf.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/builtin/gf.c,v 1.1.1.1 1999/11/10 08:12:25 noro Exp $ */
! 2: #include "ca.h"
! 3: #include "parse.h"
! 4:
! 5: struct resf_dlist {
! 6: int ib,id;
! 7: };
! 8:
! 9: int resf_degtest(int,int *,int,struct resf_dlist *);
! 10: void uhensel(P,NODE,int,int,NODE *);
! 11: void resf_hensel(int,P,int,P *,ML *);
! 12: void resf_dtest(P,ML,int,int *,int *,DCP *);
! 13: void resf_dtest_special(P,ML,int,int *,int *,DCP *);
! 14: void dtest_special(P,ML,P *);
! 15: void hensel_special(int,int,P,P *,ML *);
! 16:
! 17: void nullspace(UM **,UM,int,int,int *);
! 18: void nullspace_lm(LM **,int,int *);
! 19: void nullspace_gf2n(GF2N **,int,int *);
! 20: void nullspace_gfpn(GFPN **,int,int *);
! 21: void null_to_sol(int **,int *,int,int,UM *);
! 22:
! 23: void showgfmat(UM **,int);
! 24: void pwr_mod(P,P,V,P,int,N,P *);
! 25: void rem_mod(P,P,V,P,int,P *);
! 26:
! 27: void Pnullspace(),Pgcda_mod(),Pftest(),Presfmain(),Ppwr_mod(),Puhensel();
! 28:
! 29: void Pnullspace_ff();
! 30:
! 31: void Psolve_linear_equation_gf2n();
! 32: void Plinear_form_to_vect(),Pvect_to_linear_form();
! 33:
! 34: void solve_linear_equation_gf2n(GF2N **,int,int,int *);
! 35: void linear_form_to_array(P,VL,int,Num *);
! 36: void array_to_linear_form(Num *,VL,int,P *);
! 37:
! 38: extern int current_ff;
! 39:
! 40: struct ftab gf_tab[] = {
! 41: {"solve_linear_equation_gf2n",Psolve_linear_equation_gf2n,1},
! 42: {"nullspace",Pnullspace,3},
! 43: {"nullspace_ff",Pnullspace_ff,1},
! 44: /* {"gcda_mod",Pgcda_mod,5}, */
! 45: {"ftest",Pftest,2},
! 46: {"resfmain",Presfmain,4},
! 47: {"pwr_mod",Ppwr_mod,6},
! 48: {"uhensel",Puhensel,4},
! 49: {0,0,0},
! 50: };
! 51:
! 52: int resf_degtest(k,in,nb,dlist)
! 53: int k;
! 54: int *in;
! 55: int nb;
! 56: struct resf_dlist *dlist;
! 57: {
! 58: int i,d0,d;
! 59: int dl_i;
! 60: struct resf_dlist *t;
! 61:
! 62: if ( k < nb )
! 63: return 0;
! 64: if ( nb == 1 )
! 65: return 1;
! 66: d0 = 0; d = 0; dl_i = 0;
! 67: for ( i = 0; i < k; i++ ) {
! 68: t = &dlist[in[i]];
! 69: if ( t->ib > dl_i + 1 )
! 70: return 0;
! 71: else if ( t->ib == dl_i )
! 72: d += t->id;
! 73: else if ( !d || (dl_i && d0 != d) )
! 74: return 0;
! 75: else {
! 76: d0 = d; dl_i++; d = t->id;
! 77: }
! 78: }
! 79: if ( dl_i != nb - 1 || d0 != d )
! 80: return 0;
! 81: else
! 82: return 1;
! 83: }
! 84:
! 85: void Puhensel(arg,rp)
! 86: NODE arg;
! 87: LIST *rp;
! 88: {
! 89: P f;
! 90: NODE mfl,r;
! 91: int mod,bound;
! 92:
! 93: f = (P)ARG0(arg);
! 94: mfl = BDY((LIST)ARG1(arg));
! 95: mod = QTOS((Q)ARG2(arg));
! 96: bound = QTOS((Q)ARG3(arg));
! 97: uhensel(f,mfl,mod,bound,&r);
! 98: MKLIST(*rp,r);
! 99: }
! 100:
! 101: void uhensel(f,mfl,mod,bound,rp)
! 102: P f;
! 103: NODE mfl;
! 104: int mod,bound;
! 105: NODE *rp;
! 106: {
! 107: ML blist,clist,rlist;
! 108: LUM fl;
! 109: int nf,i;
! 110: P s;
! 111: V v;
! 112: NODE t,top;
! 113:
! 114: nf = length(mfl);
! 115: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
! 116: for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {
! 117: blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));
! 118: ptoum(mod,(P)BDY(t),blist->c[i]);
! 119: }
! 120: gcdgen(f,blist,&clist);
! 121: blist->bound = clist->bound = bound;
! 122: W_LUMALLOC((int)UDEG(f),bound,fl);
! 123: ptolum(mod,bound,f,fl);
! 124: henmain(fl,blist,clist,&rlist);
! 125: v = VR(f);
! 126: for ( i = nf-1, top = 0; i >= 0; i-- ) {
! 127: lumtop(v,mod,bound,rlist->c[i],&s);
! 128: MKNODE(t,s,top); top = t;
! 129: }
! 130: *rp = top;
! 131: }
! 132:
! 133: void Presfmain(arg,rp)
! 134: NODE arg;
! 135: LIST *rp;
! 136: {
! 137: P f;
! 138: int mod,n,nb,i,j,k;
! 139: int *nf,*md;
! 140: P *mf;
! 141: NODE mfl,mdl,t,s,u;
! 142: ML list;
! 143: DCP dc;
! 144: int sflag;
! 145:
! 146: f = (P)ARG0(arg);
! 147: mod = QTOS((Q)ARG1(arg));
! 148: mfl = BDY((LIST)ARG2(arg));
! 149: mdl = BDY((LIST)ARG3(arg));
! 150: for ( n = nb = 0, t = mfl; t; nb++, t = NEXT(t) )
! 151: for ( s = BDY((LIST)BDY(t)); s; n++, s = NEXT(s) );
! 152: if ( n == nb ) {
! 153: /* f must be irreducible! */
! 154: NEWDC(dc);
! 155: dc->c = f; dc->d = ONE;
! 156: } else {
! 157: nf = W_ALLOC(nb); md = W_ALLOC(nb); mf = (P *)ALLOCA(n*sizeof(P));
! 158: for ( k = i = 0, t = mfl, u = mdl, sflag = 1; t;
! 159: i++, t = NEXT(t), u = NEXT(u) ) {
! 160: if ( sflag && length(BDY((LIST)BDY(t))) != 2 )
! 161: sflag = 0;
! 162: for ( j = 0, s = BDY((LIST)BDY(t)); s; j++, s = NEXT(s) )
! 163: mf[k++] = (P)BDY(s);
! 164: nf[i] = j; md[i] = QTOS((Q)BDY(u));
! 165: }
! 166: resf_hensel(mod,f,n,mf,&list);
! 167: if ( sflag )
! 168: resf_dtest_special(f,list,nb,nf,md,&dc);
! 169: else
! 170: resf_dtest(f,list,nb,nf,md,&dc);
! 171: }
! 172: dcptolist(dc,rp);
! 173: }
! 174:
! 175: void resf_hensel(mod,f,nf,mfl,listp)
! 176: int mod;
! 177: P f;
! 178: int nf;
! 179: P *mfl;
! 180: ML *listp;
! 181: {
! 182: register int i,j;
! 183: int q,n,bound;
! 184: int *p;
! 185: int **pp;
! 186: ML blist,clist,bqlist,cqlist,rlist;
! 187: UM *b;
! 188: LUM fl,tl;
! 189: LUM *l;
! 190:
! 191: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
! 192: for ( i = 0; i < nf; i++ ) {
! 193: blist->c[i] = (pointer)UMALLOC(UDEG(mfl[i]));
! 194: ptoum(mod,mfl[i],blist->c[i]);
! 195: }
! 196: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 197: n = bqlist->n; q = bqlist->mod;
! 198: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 199: if ( bound == 1 ) {
! 200: *listp = rlist = MLALLOC(n);
! 201: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 202: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 203: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 204: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 205: pp[j][0] = p[j];
! 206: }
! 207: } else {
! 208: W_LUMALLOC((int)UDEG(f),bound,fl);
! 209: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 210: }
! 211: }
! 212:
! 213: void resf_dtest(f,list,nb,nfl,mdl,dcp)
! 214: P f;
! 215: ML list;
! 216: int nb;
! 217: int *nfl,*mdl;
! 218: DCP *dcp;
! 219: {
! 220: int n,np,bound,q;
! 221: int i,j,k;
! 222: int *win;
! 223: P g,factor,cofactor;
! 224: Q csum,csumt;
! 225: DCP dcf,dcf0;
! 226: LUM *c;
! 227: ML wlist;
! 228: struct resf_dlist *dlist;
! 229: int z;
! 230:
! 231: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 232: win = W_ALLOC(np+1);
! 233: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 234: wlist = W_MLALLOC(np); wlist->n = list->n;
! 235: wlist->mod = list->mod; wlist->bound = list->bound;
! 236: c = (LUM *)COEF(wlist);
! 237: bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 238: dlist = (struct resf_dlist *)ALLOCA(np*sizeof(struct resf_dlist));
! 239: for ( i = 0, j = 0; j < nb; j++ )
! 240: for ( k = 0; k < nfl[j]; k++, i++ ) {
! 241: dlist[i].ib = j; dlist[i].id = DEG(c[i])/mdl[j];
! 242: }
! 243: k = nb;
! 244: for ( i = 0; i < nb; i++ )
! 245: win[i] = i+1;
! 246: for ( g = f, dcf = dcf0 = 0, --np, z = 0; ; ) {
! 247: #if 0
! 248: if ( !(z++ % 10000) )
! 249: fputc('.',stderr);
! 250: #endif
! 251: if ( resf_degtest(k,win,nb,dlist) &&
! 252: dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
! 253: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
! 254: g = cofactor;
! 255: ucsump(g,&csum); mulq(csum,(Q)COEF(DC(g)),&csumt); csum = csumt;
! 256: for ( i = 0; i < k - 1; i++ )
! 257: for ( j = win[i] + 1; j < win[i + 1]; j++ ) {
! 258: c[j-i-1] = c[j];
! 259: dlist[j-i-1] = dlist[j];
! 260: }
! 261: for ( j = win[k-1] + 1; j <= np; j++ ) {
! 262: c[j-k] = c[j];
! 263: dlist[j-k] = dlist[j];
! 264: }
! 265: if ( ( np -= k ) < k )
! 266: break;
! 267: if ( np - win[0] + 1 < k )
! 268: if ( ++k > np )
! 269: break;
! 270: else
! 271: for ( i = 0; i < k; i++ )
! 272: win[i] = i + 1;
! 273: else
! 274: for ( i = 1; i < k; i++ )
! 275: win[i] = win[0] + i;
! 276: } else if ( !ncombi(1,np,k,win) )
! 277: if ( k == np )
! 278: break;
! 279: else
! 280: for ( i = 0, ++k; i < k; i++ )
! 281: win[i] = i + 1;
! 282: }
! 283: NEXTDC(dcf0,dcf); COEF(dcf) = g;
! 284: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
! 285: }
! 286:
! 287: void resf_dtest_special(f,list,nb,nfl,mdl,dcp)
! 288: P f;
! 289: ML list;
! 290: int nb;
! 291: int *nfl,*mdl;
! 292: DCP *dcp;
! 293: {
! 294: int n,np,bound,q;
! 295: int i,j;
! 296: int *win;
! 297: P g,factor,cofactor;
! 298: Q csum,csumt;
! 299: DCP dcf,dcf0;
! 300: LUM *c;
! 301: ML wlist;
! 302: int max;
! 303:
! 304: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 305: win = W_ALLOC(np+1);
! 306: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 307: wlist = W_MLALLOC(np); wlist->n = list->n;
! 308: wlist->mod = list->mod; wlist->bound = list->bound;
! 309: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 310: max = 1<<nb;
! 311: for ( g = f, dcf = dcf0 = 0, i = 0; i < max; i++ ) {
! 312: #if 0
! 313: if ( !(i % 1000) )
! 314: fprintf(stderr,"i=%d\n",i);
! 315: #endif
! 316: for ( j = 0; j < nb; j++ )
! 317: win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
! 318: if ( dtestmain(g,csum,wlist,nb,win,&factor,&cofactor) ) {
! 319: #if 0
! 320: fprintf(stderr,"success : i=%d\n",i);
! 321: #endif
! 322: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
! 323: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = cofactor;
! 324: NEXT(dcf) = 0;*dcp = dcf0;
! 325: return;
! 326: }
! 327: }
! 328: NEXTDC(dcf0,dcf); COEF(dcf) = g;
! 329: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
! 330: }
! 331:
! 332: #if 0
! 333: void Pftest(arg,rp)
! 334: NODE arg;
! 335: P *rp;
! 336: {
! 337: ML list;
! 338: DCP dc;
! 339: P p;
! 340: P *mfl;
! 341:
! 342: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
! 343: hensel_special(4,1,p,mfl,&list);
! 344: dtest_special(p,list,rp);
! 345: }
! 346:
! 347: void dtest_special(f,list,pr)
! 348: P f;
! 349: ML list;
! 350: P *pr;
! 351: {
! 352: int n,np,bound,q;
! 353: int i,j,k;
! 354: int *win;
! 355: P g,factor,cofactor;
! 356: Q csum,csumt;
! 357: DCP dc,dcf,dcf0;
! 358: LUM *c;
! 359: ML wlist;
! 360:
! 361: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 362: win = W_ALLOC(np+1);
! 363: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 364: wlist = W_MLALLOC(np); wlist->n = list->n;
! 365: wlist->mod = list->mod; wlist->bound = list->bound;
! 366: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 367: for ( g = f, i = 130000; i < (1<<20); i++ ) {
! 368: #if 0
! 369: if ( !(i % 1000) )
! 370: fprintf(stderr,"i=%d\n",i);
! 371: #endif
! 372: for ( j = 0; j < 20; j++ )
! 373: win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
! 374: if ( dtestmain(g,csum,wlist,20,win,&factor,&cofactor) ) {
! 375: #if 0
! 376: fprintf(stderr,"success : i=%d\n",i);
! 377: #endif
! 378: *pr = factor; return;
! 379: }
! 380: }
! 381: }
! 382:
! 383: void hensel_special(index,count,f,mfl,listp)
! 384: int index,count;
! 385: P f;
! 386: P *mfl;
! 387: ML *listp;
! 388: {
! 389: register int i,j;
! 390: int q,n,t,d,r,u,br,tmp,bound;
! 391: int *c,*p,*m,*w;
! 392: int **pp;
! 393: DCP dc;
! 394: ML blist,clist,bqlist,cqlist,rlist;
! 395: UM *b;
! 396: LUM fl,tl;
! 397: LUM *l;
! 398:
! 399: blist = MLALLOC(40); blist->n = 40; blist->mod = 11;
! 400: for ( i = 0; i < 40; i++ ) {
! 401: blist->c[i] = (pointer)UMALLOC(6);
! 402: ptoum(11,mfl[i],blist->c[i]);
! 403: }
! 404: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 405: n = bqlist->n; q = bqlist->mod;
! 406: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 407: if ( bound == 1 ) {
! 408: *listp = rlist = MLALLOC(n);
! 409: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 410: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 411: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 412: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 413: pp[j][0] = p[j];
! 414: }
! 415: } else {
! 416: W_LUMALLOC(UDEG(f),bound,fl);
! 417: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 418: }
! 419: }
! 420: #endif
! 421:
! 422: #if 0
! 423: void Pftest(arg,rp)
! 424: NODE arg;
! 425: P *rp;
! 426: {
! 427: ML list;
! 428: DCP dc;
! 429: P p;
! 430: P *mfl;
! 431:
! 432: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
! 433: hensel_special(2,1,p,mfl,&list);
! 434: dtest_special(p,list,rp);
! 435: }
! 436:
! 437: void dtest_special(f,list,pr)
! 438: P f;
! 439: ML list;
! 440: P *pr;
! 441: {
! 442: int n,np,bound,q;
! 443: int i,j,k,t,b0;
! 444: int *win;
! 445: P g,factor,cofactor;
! 446: Q csum,csumt;
! 447: DCP dc,dcf,dcf0;
! 448: LUM *c;
! 449: ML wlist;
! 450: static int nbits[16] = {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
! 451:
! 452: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 453: win = W_ALLOC(np+1);
! 454: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 455: wlist = W_MLALLOC(np); wlist->n = list->n;
! 456: wlist->mod = list->mod; wlist->bound = list->bound;
! 457: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 458: for ( g = f, i = 0; i < (1<<23); i++ ) {
! 459: #if 0
! 460: if ( !(i % 1000) )
! 461: fprintf(stderr,"i=%d\n",i);
! 462: #endif
! 463: t = i>>20; b0 = nbits[t];
! 464: if ( !b0 )
! 465: continue;
! 466: for ( j = 1; j < 6; j++ ) {
! 467: t = (i>>(20-4*j))&0xf;
! 468: if ( nbits[t] != b0 )
! 469: break;
! 470: }
! 471: if ( j != 6 )
! 472: continue;
! 473: for ( j = k = 0; j < 24; j++ )
! 474: if ( i & (1<<(23-j)) )
! 475: win[k++] = j;
! 476: if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
! 477: #if 0
! 478: fprintf(stderr,"success : i=%d\n",i);
! 479: #endif
! 480: *pr = factor; return;
! 481: }
! 482: }
! 483: *pr = f;
! 484: }
! 485:
! 486: void hensel_special(index,count,f,mfl,listp)
! 487: int index,count;
! 488: P f;
! 489: P *mfl;
! 490: ML *listp;
! 491: {
! 492: register int i,j;
! 493: int q,n,t,d,r,u,br,tmp,bound;
! 494: int *c,*p,*m,*w;
! 495: int **pp;
! 496: DCP dc;
! 497: ML blist,clist,bqlist,cqlist,rlist;
! 498: UM *b;
! 499: LUM fl,tl;
! 500: LUM *l;
! 501:
! 502: blist = MLALLOC(24); blist->n = 24; blist->mod = 5;
! 503: for ( i = 0; i < 24; i++ ) {
! 504: blist->c[i] = (pointer)UMALLOC(7);
! 505: ptoum(5,mfl[i],blist->c[i]);
! 506: }
! 507: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 508: n = bqlist->n; q = bqlist->mod;
! 509: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 510: if ( bound == 1 ) {
! 511: *listp = rlist = MLALLOC(n);
! 512: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 513: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 514: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 515: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 516: pp[j][0] = p[j];
! 517: }
! 518: } else {
! 519: W_LUMALLOC(UDEG(f),bound,fl);
! 520: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 521: }
! 522: }
! 523: #endif
! 524:
! 525: void Pftest(arg,rp)
! 526: NODE arg;
! 527: P *rp;
! 528: {
! 529: ML list;
! 530: P p;
! 531: P *mfl;
! 532:
! 533: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
! 534: hensel_special(5,1,p,mfl,&list);
! 535: dtest_special(p,list,rp);
! 536: }
! 537:
! 538: int nbits(a)
! 539: int a;
! 540: {
! 541: int i,s;
! 542:
! 543: for ( i = 0, s = 0; a && (i < 20); i++, a >>= 1 )
! 544: if ( a & 1 ) s++;
! 545: return s;
! 546: }
! 547:
! 548: void dtest_special(f,list,pr)
! 549: P f;
! 550: ML list;
! 551: P *pr;
! 552: {
! 553: int n,np,bound,q;
! 554: int i,j,k,b0;
! 555: int *win;
! 556: P g,factor,cofactor;
! 557: Q csum,csumt;
! 558: LUM *c;
! 559: ML wlist;
! 560:
! 561: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
! 562: win = W_ALLOC(np+1);
! 563: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
! 564: wlist = W_MLALLOC(np); wlist->n = list->n;
! 565: wlist->mod = list->mod; wlist->bound = list->bound;
! 566: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
! 567: for ( g = f, i = 0; i < (1<<19); i++ ) {
! 568: #if 0
! 569: if ( !(i % 10000) )
! 570: fprintf(stderr,"i=%d\n",i);
! 571: #endif
! 572: b0 = nbits(i>>10);
! 573: if ( !b0 || (nbits(i&0x3ff) != b0) )
! 574: continue;
! 575: for ( j = k = 0; j < 20; j++ )
! 576: if ( i & (1<<(19-j)) )
! 577: win[k++] = j;
! 578: if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
! 579: #if 0
! 580: fprintf(stderr,"success : i=%d\n",i);
! 581: #endif
! 582: *pr = factor; return;
! 583: }
! 584: }
! 585: *pr = f;
! 586: }
! 587:
! 588: void hensel_special(index,count,f,mfl,listp)
! 589: int index,count;
! 590: P f;
! 591: P *mfl;
! 592: ML *listp;
! 593: {
! 594: register int i,j;
! 595: int q,n,bound;
! 596: int *p;
! 597: int **pp;
! 598: ML blist,clist,bqlist,cqlist,rlist;
! 599: UM *b;
! 600: LUM fl,tl;
! 601: LUM *l;
! 602:
! 603: blist = MLALLOC(20); blist->n = 20; blist->mod = 11;
! 604: for ( i = 0; i < 20; i++ ) {
! 605: blist->c[i] = (pointer)UMALLOC(10);
! 606: ptoum(11,mfl[i],blist->c[i]);
! 607: }
! 608: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
! 609: n = bqlist->n; q = bqlist->mod;
! 610: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
! 611: if ( bound == 1 ) {
! 612: *listp = rlist = MLALLOC(n);
! 613: rlist->n = n; rlist->mod = q; rlist->bound = bound;
! 614: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
! 615: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
! 616: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
! 617: pp[j][0] = p[j];
! 618: }
! 619: } else {
! 620: W_LUMALLOC((int)UDEG(f),bound,fl);
! 621: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
! 622: }
! 623: }
! 624:
! 625: void Pnullspace(arg,rp)
! 626: NODE arg;
! 627: LIST *rp;
! 628: {
! 629: int i,j,n,mod;
! 630: MAT mat,r;
! 631: VECT u;
! 632: V v;
! 633: P p,z;
! 634: Q q;
! 635: UM **w;
! 636: UM mp;
! 637: P *t;
! 638: UM *s;
! 639: int *ind;
! 640: NODE n0,n1;
! 641:
! 642: mat = (MAT)ARG0(arg);
! 643: p = (P)ARG1(arg);
! 644: v = VR(p);
! 645: mod = QTOS((Q)ARG2(arg));
! 646: n = mat->row;
! 647: w = (UM **)almat_pointer(n,n);
! 648: for ( i = 0; i < n; i++ )
! 649: for ( j = 0, t = (P *)mat->body[i], s = w[i]; j < n; j++ ) {
! 650: ptomp(mod,t[j],&z);
! 651: s[j] = W_UMALLOC((z&&!NUM(z))?UDEG(z):0);
! 652: mptoum(z,s[j]);
! 653: }
! 654: mp = W_UMALLOC(UDEG(p)); ptoum(mod,p,mp);
! 655: ind = (int *)ALLOCA(n*sizeof(int));
! 656: nullspace(w,mp,mod,n,ind);
! 657: MKMAT(r,n,n);
! 658: for ( i = 0; i < n; i++ )
! 659: for ( j = 0, t = (P *)r->body[i], s = w[i]; j < n; j++ )
! 660: umtop(v,s[j],&t[j]);
! 661: MKVECT(u,n);
! 662: for ( i = 0; i < n; i++ ) {
! 663: STOQ(ind[i],q); u->body[i] = (pointer)q;
! 664: }
! 665: MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
! 666: }
! 667:
! 668: void nullspace(mat,p,mod,n,ind)
! 669: UM **mat;
! 670: UM p;
! 671: int mod,n;
! 672: int *ind;
! 673: {
! 674: int i,j,l,s,d;
! 675: UM q,w,w1,h,inv;
! 676: UM *t,*u;
! 677:
! 678: d = DEG(p); inv = W_UMALLOC(d); q = W_UMALLOC(2*d);
! 679: w = W_UMALLOC(2*d); w1 = W_UMALLOC(2*d); h = W_UMALLOC(d);
! 680: bzero(ind,n*sizeof(int));
! 681: ind[0] = 0;
! 682: for ( i = j = 0; j < n; i++, j++ ) {
! 683: for ( ; j < n; j++ ) {
! 684: for ( l = i; l < n; l++ )
! 685: if ( DEG(mat[l][j])>=0 )
! 686: break;
! 687: if ( l < n ) {
! 688: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 689: } else
! 690: ind[j] = 1;
! 691: }
! 692: if ( j == n )
! 693: break;
! 694: invum(mod,p,mat[i][j],inv);
! 695: for ( s = j, t = mat[i]; s < n; s++ ) {
! 696: mulum(mod,t[s],inv,w);
! 697: DEG(w) = divum(mod,w,p,q);
! 698: cpyum(w,t[s]);
! 699: }
! 700: for ( l = 0; l < n; l++ ) {
! 701: if ( l == i )
! 702: continue;
! 703: u = mat[l]; DEG(w) = -1; subum(mod,w,u[j],h);
! 704: for ( s = j; s < n; s++ ) {
! 705: mulum(mod,h,t[s],w); addum(mod,w,u[s],w1);
! 706: DEG(w1) = divum(mod,w1,p,q); cpyum(w1,u[s]);
! 707: }
! 708: }
! 709: }
! 710: }
! 711:
! 712: void Pnullspace_ff(arg,rp)
! 713: NODE arg;
! 714: LIST *rp;
! 715: {
! 716: int i,j,n;
! 717: Q mod;
! 718: MAT mat,r;
! 719: VECT u;
! 720: Q q;
! 721: Obj **w;
! 722: Obj *t;
! 723: Obj *s;
! 724: int *ind;
! 725: NODE n0,n1;
! 726:
! 727: mat = (MAT)ARG0(arg);
! 728: n = mat->row;
! 729: w = (Obj **)almat_pointer(n,n);
! 730: for ( i = 0; i < n; i++ )
! 731: for ( j = 0, t = (Obj *)mat->body[i], s = w[i]; j < n; j++ )
! 732: s[j] = t[j];
! 733: ind = (int *)ALLOCA(n*sizeof(int));
! 734: switch ( current_ff ) {
! 735: case FF_GFP:
! 736: nullspace_lm((LM **)w,n,ind); break;
! 737: case FF_GF2N:
! 738: nullspace_gf2n((GF2N **)w,n,ind); break;
! 739: case FF_GFPN:
! 740: nullspace_gfpn((GFPN **)w,n,ind); break;
! 741: default:
! 742: error("nullspace_ff : current_ff is not set");
! 743: }
! 744: MKMAT(r,n,n);
! 745: for ( i = 0; i < n; i++ )
! 746: for ( j = 0, t = (Obj *)r->body[i], s = w[i]; j < n; j++ )
! 747: t[j] = s[j];
! 748: MKVECT(u,n);
! 749: for ( i = 0; i < n; i++ ) {
! 750: STOQ(ind[i],q); u->body[i] = (pointer)q;
! 751: }
! 752: MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
! 753: }
! 754:
! 755: void nullspace_lm(mat,n,ind)
! 756: LM **mat;
! 757: int n;
! 758: int *ind;
! 759: {
! 760: int i,j,l,s;
! 761: Q q,mod;
! 762: N lm;
! 763: LM w,w1,h,inv;
! 764: LM *t,*u;
! 765:
! 766: getmod_lm(&lm); NTOQ(lm,1,mod);
! 767:
! 768: bzero(ind,n*sizeof(int));
! 769: ind[0] = 0;
! 770: for ( i = j = 0; j < n; i++, j++ ) {
! 771: for ( ; j < n; j++ ) {
! 772: for ( l = i; l < n; l++ ) {
! 773: simplm(mat[l][j],&w); mat[l][j] = w;
! 774: if ( mat[l][j] )
! 775: break;
! 776: }
! 777: if ( l < n ) {
! 778: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 779: } else
! 780: ind[j] = 1;
! 781: }
! 782: if ( j == n )
! 783: break;
! 784: NTOQ(mat[i][j]->body,1,q); invl(q,mod,(Q *)&inv);
! 785: for ( s = j, t = mat[i]; s < n; s++ ) {
! 786: mullm(t[s],inv,&w); t[s] = w;
! 787: }
! 788: for ( l = 0; l < n; l++ ) {
! 789: if ( l == i )
! 790: continue;
! 791: u = mat[l]; chsgnlm(u[j],&h);
! 792: for ( s = j; s < n; s++ ) {
! 793: mullm(h,t[s],&w); addlm(w,u[s],&w1); u[s] = w1;
! 794: }
! 795: }
! 796: }
! 797: }
! 798:
! 799: void nullspace_gf2n(mat,n,ind)
! 800: GF2N **mat;
! 801: int n;
! 802: int *ind;
! 803: {
! 804: int i,j,l,s;
! 805: GF2N w,w1,h,inv;
! 806: GF2N *t,*u;
! 807: extern gf2n_lazy;
! 808:
! 809: bzero(ind,n*sizeof(int));
! 810: ind[0] = 0;
! 811: for ( i = j = 0; j < n; i++, j++ ) {
! 812: for ( ; j < n; j++ ) {
! 813: for ( l = i; l < n; l++ ) {
! 814: simpgf2n(mat[l][j],&w); mat[l][j] = w;
! 815: if ( mat[l][j] )
! 816: break;
! 817: }
! 818: if ( l < n ) {
! 819: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 820: } else
! 821: ind[j] = 1;
! 822: }
! 823: if ( j == n )
! 824: break;
! 825: invgf2n(mat[i][j],&inv);
! 826: for ( s = j, t = mat[i]; s < n; s++ ) {
! 827: mulgf2n(t[s],inv,&w); t[s] = w;
! 828: }
! 829: for ( l = 0; l < n; l++ ) {
! 830: if ( l == i )
! 831: continue;
! 832: u = mat[l]; h = u[j];
! 833: for ( s = j; s < n; s++ ) {
! 834: mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
! 835: }
! 836: }
! 837: }
! 838: }
! 839:
! 840: void nullspace_gfpn(mat,n,ind)
! 841: GFPN **mat;
! 842: int n;
! 843: int *ind;
! 844: {
! 845: int i,j,l,s;
! 846: GFPN w,w1,h,inv;
! 847: GFPN *t,*u;
! 848:
! 849: bzero(ind,n*sizeof(int));
! 850: ind[0] = 0;
! 851: for ( i = j = 0; j < n; i++, j++ ) {
! 852: for ( ; j < n; j++ ) {
! 853: for ( l = i; l < n; l++ ) {
! 854: simpgfpn(mat[l][j],&w); mat[l][j] = w;
! 855: if ( mat[l][j] )
! 856: break;
! 857: }
! 858: if ( l < n ) {
! 859: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 860: } else
! 861: ind[j] = 1;
! 862: }
! 863: if ( j == n )
! 864: break;
! 865: divgfpn((GFPN)ONE,(GFPN)mat[i][j],&inv);
! 866: for ( s = j, t = mat[i]; s < n; s++ ) {
! 867: mulgfpn(t[s],inv,&w); t[s] = w;
! 868: }
! 869: for ( l = 0; l < n; l++ ) {
! 870: if ( l == i )
! 871: continue;
! 872: u = mat[l]; chsgngfpn(u[j],&h);
! 873: for ( s = j; s < n; s++ ) {
! 874: mulgfpn(h,t[s],&w); addgfpn(w,u[s],&w1); u[s] = w1;
! 875: }
! 876: }
! 877: }
! 878: }
! 879: /* p = a(0)vl[0]+a(1)vl[1]+...+a(m-1)vl[m-1]+a(m) -> array = [a(0) a(1) ... a(m)] */
! 880:
! 881: void linear_form_to_array(p,vl,m,array)
! 882: P p;
! 883: VL vl;
! 884: int m;
! 885: Num *array;
! 886: {
! 887: int i;
! 888: DCP dc;
! 889:
! 890: bzero((char *)array,(m+1)*sizeof(Num *));
! 891: for ( i = 0; p && vl; vl = NEXT(vl), i++ ) {
! 892: if ( ID(p) == O_N )
! 893: break;
! 894: else if ( VR(p) == vl->v ) {
! 895: dc = DC(p);
! 896: array[i] = (Num)COEF(dc);
! 897: dc = NEXT(dc);
! 898: p = dc ? COEF(dc) : 0;
! 899: }
! 900: }
! 901: array[m] = (Num)p;
! 902: }
! 903:
! 904: void array_to_linear_form(array,vl,m,r)
! 905: Num *array;
! 906: VL vl;
! 907: int m;
! 908: P *r;
! 909: {
! 910: P t;
! 911: DCP dc0,dc1;
! 912:
! 913: if ( !m )
! 914: *r = (P)array[0];
! 915: else {
! 916: array_to_linear_form(array+1,NEXT(vl),m-1,&t);
! 917: if ( !array[0] )
! 918: *r = t;
! 919: else {
! 920: NEWDC(dc0); DEG(dc0) = ONE; COEF(dc0) = (P)array[0];
! 921: if ( !t )
! 922: NEXT(dc0) = 0;
! 923: else {
! 924: NEWDC(dc1); DEG(dc1) = 0; COEF(dc1) = t;
! 925: NEXT(dc1) = 0;
! 926: NEXT(dc0) = dc1;
! 927: }
! 928: MKP(vl->v,dc0,*r);
! 929: }
! 930: }
! 931: }
! 932:
! 933: void Psolve_linear_equation_gf2n(arg,rp)
! 934: NODE arg;
! 935: LIST *rp;
! 936: {
! 937: NODE eqs,tn;
! 938: VL vars,tvl;
! 939: int i,j,n,m,dim,codim;
! 940: GF2N **w;
! 941: int *ind;
! 942: NODE n0,n1;
! 943:
! 944: get_vars(ARG0(arg),&vars);
! 945: eqs = BDY((LIST)ARG0(arg));
! 946: for ( n = 0, tn = eqs; tn; tn = NEXT(tn), n++);
! 947: for ( m = 0, tvl = vars; tvl; tvl = NEXT(tvl), m++);
! 948: w = (GF2N **)almat_pointer(n,m+1);
! 949: for ( i = 0, tn = eqs; i < n; i++, tn = NEXT(tn) )
! 950: linear_form_to_array(BDY(tn),vars,m,(Num *)w[i]);
! 951: ind = (int *)ALLOCA(m*sizeof(int));
! 952: solve_linear_equation_gf2n(w,n,m,ind);
! 953: for ( j = 0, dim = 0; j < m; j++ )
! 954: if ( ind[j] )
! 955: dim++;
! 956: codim = m-dim;
! 957: for ( i = codim; i < n; i++ )
! 958: if ( w[i][m] ) {
! 959: MKLIST(*rp,0); return;
! 960: }
! 961: for ( i = 0, n0 = 0; i < codim; i++ ) {
! 962: NEXTNODE(n0,n1);
! 963: array_to_linear_form((Num *)w[i],vars,m,(P *)&BDY(n1));
! 964: }
! 965: if ( n0 )
! 966: NEXT(n1) = 0;
! 967: MKLIST(*rp,n0);
! 968: }
! 969:
! 970: void solve_linear_equation_gf2n(mat,n,m,ind)
! 971: GF2N **mat;
! 972: int n;
! 973: int *ind;
! 974: {
! 975: int i,j,l,s;
! 976: GF2N w,w1,h,inv;
! 977: GF2N *t,*u;
! 978: extern gf2n_lazy;
! 979:
! 980: bzero(ind,m*sizeof(int));
! 981: ind[0] = 0;
! 982: for ( i = j = 0; j < m; i++, j++ ) {
! 983: for ( ; j < m; j++ ) {
! 984: for ( l = i; l < n; l++ ) {
! 985: simpgf2n(mat[l][j],&w); mat[l][j] = w;
! 986: if ( mat[l][j] )
! 987: break;
! 988: }
! 989: if ( l < n ) {
! 990: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
! 991: } else
! 992: ind[j] = 1;
! 993: }
! 994: if ( j == m )
! 995: break;
! 996: invgf2n(mat[i][j],&inv);
! 997: for ( s = j, t = mat[i]; s <= m; s++ ) {
! 998: mulgf2n(t[s],inv,&w); t[s] = w;
! 999: }
! 1000: for ( l = 0; l < n; l++ ) {
! 1001: if ( l == i )
! 1002: continue;
! 1003: u = mat[l]; h = u[j];
! 1004: for ( s = j; s <= m; s++ ) {
! 1005: mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
! 1006: }
! 1007: }
! 1008: }
! 1009: }
! 1010:
! 1011: /*
! 1012: void null_to_sol(mat,ind,mod,n,r)
! 1013: int **mat;
! 1014: int *ind;
! 1015: int mod,n;
! 1016: UM *r;
! 1017: {
! 1018: int i,j,k,l;
! 1019: int *c;
! 1020: UM w;
! 1021:
! 1022: for ( i = 0, l = 0; i < n; i++ ) {
! 1023: if ( !ind[i] )
! 1024: continue;
! 1025: w = UMALLOC(n);
! 1026: for ( j = k = 0, c = COEF(w); j < n; j++ )
! 1027: if ( ind[j] )
! 1028: c[j] = 0;
! 1029: else
! 1030: c[j] = mat[k++][i];
! 1031: c[i] = mod-1;
! 1032: for ( j = n; j >= 0; j-- )
! 1033: if ( c[j] )
! 1034: break;
! 1035: DEG(w) = j;
! 1036: r[l++] = w;
! 1037: }
! 1038: }
! 1039: */
! 1040:
! 1041: void showgfmat(mat,n)
! 1042: UM **mat;
! 1043: int n;
! 1044: {
! 1045: int i,j,k;
! 1046: int *c;
! 1047: UM p;
! 1048:
! 1049: for ( i = 0; i < n; i++ ) {
! 1050: for ( j = 0; j < n; j++ ) {
! 1051: p = mat[i][j];
! 1052: if ( DEG(p) < 0 )
! 1053: fprintf(asir_out,"0");
! 1054: else
! 1055: for ( p = mat[i][j], k = DEG(p), c = COEF(p); k >= 0; k-- ) {
! 1056: if ( c[k] )
! 1057: fprintf(asir_out,"+%d",c[k]);
! 1058: if ( k > 1 )
! 1059: fprintf(asir_out,"a^%d",k);
! 1060: else if ( k == 1 )
! 1061: fprintf(asir_out,"a",k);
! 1062: }
! 1063: fprintf(asir_out," ");
! 1064: }
! 1065: fprintf(asir_out,"\n");
! 1066: }
! 1067: }
! 1068:
! 1069: #if 0
! 1070: void Pgcda_mod(arg,rp)
! 1071: NODE arg;
! 1072: P *rp;
! 1073: {
! 1074: p1 = (P)ARG0(arg);
! 1075: p2 = (P)ARG1(arg);
! 1076: v = VR((P)ARG2(arg));
! 1077: d = (P)ARG3(arg);
! 1078: m = QTOS((Q)ARG4(arg));
! 1079: reordvar(CO,v,&vl);
! 1080: reorderp(vl,CO,p1,&t); ptomp(m,t,&m1);
! 1081: reorderp(vl,CO,p2,&t); ptomp(m,t,&m2);
! 1082: if ( NUM(m1) || NUM(m2) || VR(m1) != v || VR(m2) != v ) {
! 1083: *rp = ONE; return;
! 1084: }
! 1085: if ( deg(v,m1) >= deg(v,m2) ) {
! 1086: t = m1; m1 = m2; m2 = t;
! 1087: }
! 1088: while ( 1 ) {
! 1089: inva_mod(COEF(DC(m2)),d,m,&inv);
! 1090: NEWDC(dc); NEXT(dc) = 0; MKP(v,dc,h);
! 1091: d0 = deg(v,m1)-deg(v,m2); STOQ(d0,DEG(dc));
! 1092: mulgq(m,d,inv,COEF(DC(m1)),&COEF(dc));
! 1093: mulgp(vl,m,d,m2,h,&t); subgp(vl,m,d,m1,t,&s);
! 1094: }
! 1095: }
! 1096: #endif
! 1097:
! 1098: void Ppwr_mod(arg,rp)
! 1099: NODE arg;
! 1100: P *rp;
! 1101: {
! 1102: P p,a,d,r;
! 1103: int m;
! 1104: Q q;
! 1105: N n;
! 1106:
! 1107: m = QTOS((Q)ARG4(arg)); q = (Q)ARG5(arg); n = q ? NM(q) : 0;
! 1108: ptomp(m,(P)ARG0(arg),&p); ptomp(m,(P)ARG1(arg),&a);
! 1109: ptomp(m,(P)ARG3(arg),&d);
! 1110: pwr_mod(p,a,VR((P)ARG2(arg)),d,m,n,&r);
! 1111: mptop(r,rp);
! 1112: }
! 1113:
! 1114: void pwr_mod(p,a,v,d,m,n,rp)
! 1115: P p,a,d;
! 1116: V v;
! 1117: int m;
! 1118: N n;
! 1119: P *rp;
! 1120: {
! 1121: int b;
! 1122: P t,s,r;
! 1123: N n1;
! 1124:
! 1125: if ( !n )
! 1126: *rp = (P)ONEM;
! 1127: else if ( UNIN(n) )
! 1128: *rp = p;
! 1129: else {
! 1130: b = divin(n,2,&n1);
! 1131: pwr_mod(p,a,v,d,m,n1,&t);
! 1132: mulmp(CO,m,t,t,&s); rem_mod(s,a,v,d,m,&r);
! 1133: if ( b ) {
! 1134: mulmp(CO,m,r,p,&t); rem_mod(t,a,v,d,m,rp);
! 1135: } else
! 1136: *rp = r;
! 1137: }
! 1138: }
! 1139:
! 1140: void rem_mod(p,a,v,d,m,rp)
! 1141: P p,a,d;
! 1142: V v;
! 1143: int m;
! 1144: P *rp;
! 1145: {
! 1146: P tmp,r1,r2;
! 1147:
! 1148: divsrmp(CO,m,p,d,&tmp,&r1);
! 1149: divsrmp(CO,m,r1,a,&tmp,&r2);
! 1150: divsrmp(CO,m,r2,d,&tmp,rp);
! 1151: }