![[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.3 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.4 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.9 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/builtin/gf.c,v 1.8 2001/06/25 10:01:27 noro Exp $
1.3 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "parse.h"
52:
53: struct resf_dlist {
54: int ib,id;
55: };
56:
57: int resf_degtest(int,int *,int,struct resf_dlist *);
58: void uhensel(P,NODE,int,int,NODE *);
59: void resf_hensel(int,P,int,P *,ML *);
60: void resf_dtest(P,ML,int,int *,int *,DCP *);
61: void resf_dtest_special(P,ML,int,int *,int *,DCP *);
62: void dtest_special(P,ML,P *);
63: void hensel_special(int,int,P,P *,ML *);
64:
65: void nullspace(UM **,UM,int,int,int *);
66: void nullspace_lm(LM **,int,int *);
67: void nullspace_gf2n(GF2N **,int,int *);
68: void nullspace_gfpn(GFPN **,int,int *);
1.5 noro 69: void nullspace_gfs(GFS **,int,int *);
1.1 noro 70: void null_to_sol(int **,int *,int,int,UM *);
71:
72: void showgfmat(UM **,int);
73: void pwr_mod(P,P,V,P,int,N,P *);
74: void rem_mod(P,P,V,P,int,P *);
75:
76: void Pnullspace(),Pgcda_mod(),Pftest(),Presfmain(),Ppwr_mod(),Puhensel();
1.6 noro 77: void Psfuhensel();
1.1 noro 78:
79: void Pnullspace_ff();
80:
81: void Psolve_linear_equation_gf2n();
82: void Plinear_form_to_vect(),Pvect_to_linear_form();
83:
84: void solve_linear_equation_gf2n(GF2N **,int,int,int *);
85: void linear_form_to_array(P,VL,int,Num *);
86: void array_to_linear_form(Num *,VL,int,P *);
1.9 ! noro 87: void sfuhensel(P,NODE,GFS,int,NODE *);
1.1 noro 88:
89: extern int current_ff;
90:
91: struct ftab gf_tab[] = {
92: {"solve_linear_equation_gf2n",Psolve_linear_equation_gf2n,1},
93: {"nullspace",Pnullspace,3},
94: {"nullspace_ff",Pnullspace_ff,1},
95: /* {"gcda_mod",Pgcda_mod,5}, */
96: {"ftest",Pftest,2},
97: {"resfmain",Presfmain,4},
98: {"pwr_mod",Ppwr_mod,6},
99: {"uhensel",Puhensel,4},
1.9 ! noro 100: {"sfuhensel",Psfuhensel,4},
1.1 noro 101: {0,0,0},
102: };
103:
104: int resf_degtest(k,in,nb,dlist)
105: int k;
106: int *in;
107: int nb;
108: struct resf_dlist *dlist;
109: {
110: int i,d0,d;
111: int dl_i;
112: struct resf_dlist *t;
113:
114: if ( k < nb )
115: return 0;
116: if ( nb == 1 )
117: return 1;
118: d0 = 0; d = 0; dl_i = 0;
119: for ( i = 0; i < k; i++ ) {
120: t = &dlist[in[i]];
121: if ( t->ib > dl_i + 1 )
122: return 0;
123: else if ( t->ib == dl_i )
124: d += t->id;
125: else if ( !d || (dl_i && d0 != d) )
126: return 0;
127: else {
128: d0 = d; dl_i++; d = t->id;
129: }
130: }
131: if ( dl_i != nb - 1 || d0 != d )
132: return 0;
133: else
134: return 1;
135: }
136:
137: void Puhensel(arg,rp)
138: NODE arg;
139: LIST *rp;
140: {
141: P f;
142: NODE mfl,r;
143: int mod,bound;
144:
145: f = (P)ARG0(arg);
146: mfl = BDY((LIST)ARG1(arg));
147: mod = QTOS((Q)ARG2(arg));
148: bound = QTOS((Q)ARG3(arg));
149: uhensel(f,mfl,mod,bound,&r);
150: MKLIST(*rp,r);
151: }
152:
153: void uhensel(f,mfl,mod,bound,rp)
154: P f;
155: NODE mfl;
156: int mod,bound;
157: NODE *rp;
158: {
159: ML blist,clist,rlist;
160: LUM fl;
161: int nf,i;
162: P s;
163: V v;
164: NODE t,top;
165:
166: nf = length(mfl);
167: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
168: for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {
169: blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));
170: ptoum(mod,(P)BDY(t),blist->c[i]);
171: }
172: gcdgen(f,blist,&clist);
173: blist->bound = clist->bound = bound;
174: W_LUMALLOC((int)UDEG(f),bound,fl);
175: ptolum(mod,bound,f,fl);
176: henmain(fl,blist,clist,&rlist);
177: v = VR(f);
178: for ( i = nf-1, top = 0; i >= 0; i-- ) {
179: lumtop(v,mod,bound,rlist->c[i],&s);
180: MKNODE(t,s,top); top = t;
1.6 noro 181: }
182: *rp = top;
183: }
184:
185: void Psfuhensel(arg,rp)
186: NODE arg;
187: LIST *rp;
188: {
189: P f;
1.9 ! noro 190: int bound;
! 191: NODE r,mfl;
1.8 noro 192: GFS ev;
1.6 noro 193:
194: f = (P)ARG0(arg);
1.9 ! noro 195: mfl = BDY((LIST)ARG1(arg));
! 196: ev = (GFS)ARG2(arg);
! 197: bound = QTOS((Q)ARG3(arg));
! 198: sfuhensel(f,mfl,ev,bound,&r);
! 199: MKLIST(*rp,r);
1.6 noro 200: }
201:
1.9 ! noro 202: void sfuhensel(f,mfl,ev,bound,rp)
1.6 noro 203: P f;
1.9 ! noro 204: NODE mfl;
! 205: GFS ev;
! 206: int bound;
1.6 noro 207: NODE *rp;
208: {
1.9 ! noro 209: BM fl;
! 210: BM *r;
! 211: VL vl,nvl;
! 212: int i,fn,dx,dy;
1.6 noro 213: NODE t,top;
1.9 ! noro 214: UM fm,hm,q;
! 215: UM *gm;
! 216: V x,y;
! 217: P g,s,u;
! 218:
! 219: clctv(CO,f,&vl);
! 220: if ( !vl || !vl->next || vl->next->next )
! 221: error("sfuhensel : f must be a bivariate poly");
1.8 noro 222:
1.9 ! noro 223: for ( i = 0, t = mfl; t; i++, t = NEXT(t) );
! 224: fn = i;
1.8 noro 225:
1.9 ! noro 226: gm = (UM *)MALLOC(fn*sizeof(UM));
! 227:
! 228: /* XXX : more severe check is necessary */
! 229: x = VR((P)BDY(mfl));
! 230: y = vl->v == x ? vl->next->v : vl->v;
! 231:
! 232: for ( i = 0, t = mfl; i < fn; i++, t = NEXT(t) ) {
! 233: gm[i] = (pointer)UMALLOC(getdeg(x,(P)BDY(t)));
! 234: ptosfum((P)BDY(t),gm[i]);
! 235: }
! 236:
! 237: /* reorder f if necessary */
! 238: if ( vl->v != x ) {
! 239: reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&g);
! 240: vl = nvl; f = g;
! 241: }
! 242: dx = getdeg(x,f);
! 243: dy = getdeg(y,f);
! 244: if ( bound < dy+1 ) bound = dy+1;
! 245: fl = BMALLOC(dx,bound);
! 246: ptosfbm(bound,f,fl);
! 247: shiftsfbm(bound,fl,FTOIF(CONT(ev)));
! 248:
! 249: /* fm = fl mod y */
! 250: fm = W_UMALLOC(dx);
! 251: cpyum(COEF(fl)[0],fm);
! 252: hm = W_UMALLOC(dx);
! 253:
! 254: q = W_UMALLOC(dx);
! 255: r = (BM *)MLALLOC(fn*sizeof(BM));
! 256: for ( i = 0; i < fn-1; i++ ) {
! 257: /* fl = gm[i]*hm mod y */
! 258: divsfum(fm,gm[i],hm);
! 259: /* fl is replaced by the cofactor of gk mod y^bound */
! 260: /* r[i] = gk */
! 261: sfhenmain2(fl,gm[i],hm,bound,r+i);
! 262: cpyum(hm,fm);
! 263: }
! 264: /* finally, fl must be the lift of gm[fn-1] */
! 265: r[i] = fl;
1.6 noro 266:
1.9 ! noro 267: for ( i = fn-1, top = 0; i >= 0; i-- ) {
! 268: sfbmtop(bound,r[i],x,y,&s);
! 269: reorderp(CO,vl,s,&u);
1.8 noro 270: MKNODE(t,u,top); top = t;
1.1 noro 271: }
272: *rp = top;
273: }
274:
275: void Presfmain(arg,rp)
276: NODE arg;
277: LIST *rp;
278: {
279: P f;
280: int mod,n,nb,i,j,k;
281: int *nf,*md;
282: P *mf;
283: NODE mfl,mdl,t,s,u;
284: ML list;
285: DCP dc;
286: int sflag;
287:
288: f = (P)ARG0(arg);
289: mod = QTOS((Q)ARG1(arg));
290: mfl = BDY((LIST)ARG2(arg));
291: mdl = BDY((LIST)ARG3(arg));
292: for ( n = nb = 0, t = mfl; t; nb++, t = NEXT(t) )
293: for ( s = BDY((LIST)BDY(t)); s; n++, s = NEXT(s) );
294: if ( n == nb ) {
295: /* f must be irreducible! */
296: NEWDC(dc);
297: dc->c = f; dc->d = ONE;
298: } else {
299: nf = W_ALLOC(nb); md = W_ALLOC(nb); mf = (P *)ALLOCA(n*sizeof(P));
300: for ( k = i = 0, t = mfl, u = mdl, sflag = 1; t;
301: i++, t = NEXT(t), u = NEXT(u) ) {
302: if ( sflag && length(BDY((LIST)BDY(t))) != 2 )
303: sflag = 0;
304: for ( j = 0, s = BDY((LIST)BDY(t)); s; j++, s = NEXT(s) )
305: mf[k++] = (P)BDY(s);
306: nf[i] = j; md[i] = QTOS((Q)BDY(u));
307: }
308: resf_hensel(mod,f,n,mf,&list);
309: if ( sflag )
310: resf_dtest_special(f,list,nb,nf,md,&dc);
311: else
312: resf_dtest(f,list,nb,nf,md,&dc);
313: }
314: dcptolist(dc,rp);
315: }
316:
317: void resf_hensel(mod,f,nf,mfl,listp)
318: int mod;
319: P f;
320: int nf;
321: P *mfl;
322: ML *listp;
323: {
324: register int i,j;
1.2 noro 325: int q,n,bound,inv,lc;
1.1 noro 326: int *p;
327: int **pp;
328: ML blist,clist,bqlist,cqlist,rlist;
329: UM *b;
330: LUM fl,tl;
331: LUM *l;
1.2 noro 332: UM w;
1.1 noro 333:
1.2 noro 334: w = W_UMALLOC(UDEG(f));
1.1 noro 335: blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
1.2 noro 336:
337: /* c[0] must have lc(f) */
338: blist->c[0] = (pointer)UMALLOC(UDEG(mfl[0]));
339: ptoum(mod,mfl[0],w);
340: inv = invm(w->c[UDEG(mfl[0])],mod);
341: lc = rem(NM((Q)LC(f)),mod);
342: if ( SGN((Q)LC(f)) < 0 )
343: lc = (mod-lc)%mod;
344: lc = dmar(inv,lc,0,mod);
345: if ( lc == 1 )
346: copyum(w,blist->c[0]);
347: else
348: mulsum(mod,w,lc,blist->c[0]);
349:
350: /* c[i] (i=1,...,nf-1) must be monic */
351: for ( i = 1; i < nf; i++ ) {
1.1 noro 352: blist->c[i] = (pointer)UMALLOC(UDEG(mfl[i]));
1.2 noro 353: ptoum(mod,mfl[i],w);
354: inv = invm(w->c[UDEG(mfl[i])],mod);
355: if ( inv == 1 )
356: copyum(w,blist->c[i]);
357: else
358: mulsum(mod,w,inv,blist->c[i]);
1.1 noro 359: }
1.2 noro 360:
1.1 noro 361: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
362: n = bqlist->n; q = bqlist->mod;
363: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
364: if ( bound == 1 ) {
365: *listp = rlist = MLALLOC(n);
366: rlist->n = n; rlist->mod = q; rlist->bound = bound;
367: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
368: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
369: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
370: pp[j][0] = p[j];
371: }
372: } else {
373: W_LUMALLOC((int)UDEG(f),bound,fl);
374: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
375: }
376: }
377:
378: void resf_dtest(f,list,nb,nfl,mdl,dcp)
379: P f;
380: ML list;
381: int nb;
382: int *nfl,*mdl;
383: DCP *dcp;
384: {
385: int n,np,bound,q;
386: int i,j,k;
387: int *win;
388: P g,factor,cofactor;
389: Q csum,csumt;
390: DCP dcf,dcf0;
391: LUM *c;
392: ML wlist;
393: struct resf_dlist *dlist;
394: int z;
395:
396: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
397: win = W_ALLOC(np+1);
398: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
399: wlist = W_MLALLOC(np); wlist->n = list->n;
400: wlist->mod = list->mod; wlist->bound = list->bound;
401: c = (LUM *)COEF(wlist);
402: bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
403: dlist = (struct resf_dlist *)ALLOCA(np*sizeof(struct resf_dlist));
404: for ( i = 0, j = 0; j < nb; j++ )
405: for ( k = 0; k < nfl[j]; k++, i++ ) {
406: dlist[i].ib = j; dlist[i].id = DEG(c[i])/mdl[j];
407: }
408: k = nb;
409: for ( i = 0; i < nb; i++ )
410: win[i] = i+1;
411: for ( g = f, dcf = dcf0 = 0, --np, z = 0; ; ) {
412: #if 0
413: if ( !(z++ % 10000) )
414: fputc('.',stderr);
415: #endif
416: if ( resf_degtest(k,win,nb,dlist) &&
417: dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
418: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
419: g = cofactor;
420: ucsump(g,&csum); mulq(csum,(Q)COEF(DC(g)),&csumt); csum = csumt;
421: for ( i = 0; i < k - 1; i++ )
422: for ( j = win[i] + 1; j < win[i + 1]; j++ ) {
423: c[j-i-1] = c[j];
424: dlist[j-i-1] = dlist[j];
425: }
426: for ( j = win[k-1] + 1; j <= np; j++ ) {
427: c[j-k] = c[j];
428: dlist[j-k] = dlist[j];
429: }
430: if ( ( np -= k ) < k )
431: break;
432: if ( np - win[0] + 1 < k )
433: if ( ++k > np )
434: break;
435: else
436: for ( i = 0; i < k; i++ )
437: win[i] = i + 1;
438: else
439: for ( i = 1; i < k; i++ )
440: win[i] = win[0] + i;
441: } else if ( !ncombi(1,np,k,win) )
442: if ( k == np )
443: break;
444: else
445: for ( i = 0, ++k; i < k; i++ )
446: win[i] = i + 1;
447: }
448: NEXTDC(dcf0,dcf); COEF(dcf) = g;
449: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
450: }
451:
452: void resf_dtest_special(f,list,nb,nfl,mdl,dcp)
453: P f;
454: ML list;
455: int nb;
456: int *nfl,*mdl;
457: DCP *dcp;
458: {
459: int n,np,bound,q;
460: int i,j;
461: int *win;
462: P g,factor,cofactor;
463: Q csum,csumt;
464: DCP dcf,dcf0;
465: LUM *c;
466: ML wlist;
467: int max;
468:
469: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
470: win = W_ALLOC(np+1);
471: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
472: wlist = W_MLALLOC(np); wlist->n = list->n;
473: wlist->mod = list->mod; wlist->bound = list->bound;
474: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
475: max = 1<<nb;
476: for ( g = f, dcf = dcf0 = 0, i = 0; i < max; i++ ) {
477: #if 0
478: if ( !(i % 1000) )
479: fprintf(stderr,"i=%d\n",i);
480: #endif
481: for ( j = 0; j < nb; j++ )
482: win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
483: if ( dtestmain(g,csum,wlist,nb,win,&factor,&cofactor) ) {
484: #if 0
485: fprintf(stderr,"success : i=%d\n",i);
486: #endif
487: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
488: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = cofactor;
489: NEXT(dcf) = 0;*dcp = dcf0;
490: return;
491: }
492: }
493: NEXTDC(dcf0,dcf); COEF(dcf) = g;
494: DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
495: }
496:
497: #if 0
498: void Pftest(arg,rp)
499: NODE arg;
500: P *rp;
501: {
502: ML list;
503: DCP dc;
504: P p;
505: P *mfl;
506:
507: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
508: hensel_special(4,1,p,mfl,&list);
509: dtest_special(p,list,rp);
510: }
511:
512: void dtest_special(f,list,pr)
513: P f;
514: ML list;
515: P *pr;
516: {
517: int n,np,bound,q;
518: int i,j,k;
519: int *win;
520: P g,factor,cofactor;
521: Q csum,csumt;
522: DCP dc,dcf,dcf0;
523: LUM *c;
524: ML wlist;
525:
526: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
527: win = W_ALLOC(np+1);
528: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
529: wlist = W_MLALLOC(np); wlist->n = list->n;
530: wlist->mod = list->mod; wlist->bound = list->bound;
531: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
532: for ( g = f, i = 130000; i < (1<<20); i++ ) {
533: #if 0
534: if ( !(i % 1000) )
535: fprintf(stderr,"i=%d\n",i);
536: #endif
537: for ( j = 0; j < 20; j++ )
538: win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
539: if ( dtestmain(g,csum,wlist,20,win,&factor,&cofactor) ) {
540: #if 0
541: fprintf(stderr,"success : i=%d\n",i);
542: #endif
543: *pr = factor; return;
544: }
545: }
546: }
547:
548: void hensel_special(index,count,f,mfl,listp)
549: int index,count;
550: P f;
551: P *mfl;
552: ML *listp;
553: {
554: register int i,j;
555: int q,n,t,d,r,u,br,tmp,bound;
556: int *c,*p,*m,*w;
557: int **pp;
558: DCP dc;
559: ML blist,clist,bqlist,cqlist,rlist;
560: UM *b;
561: LUM fl,tl;
562: LUM *l;
563:
564: blist = MLALLOC(40); blist->n = 40; blist->mod = 11;
565: for ( i = 0; i < 40; i++ ) {
566: blist->c[i] = (pointer)UMALLOC(6);
567: ptoum(11,mfl[i],blist->c[i]);
568: }
569: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
570: n = bqlist->n; q = bqlist->mod;
571: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
572: if ( bound == 1 ) {
573: *listp = rlist = MLALLOC(n);
574: rlist->n = n; rlist->mod = q; rlist->bound = bound;
575: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
576: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
577: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
578: pp[j][0] = p[j];
579: }
580: } else {
581: W_LUMALLOC(UDEG(f),bound,fl);
582: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
583: }
584: }
585: #endif
586:
587: #if 0
588: void Pftest(arg,rp)
589: NODE arg;
590: P *rp;
591: {
592: ML list;
593: DCP dc;
594: P p;
595: P *mfl;
596:
597: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
598: hensel_special(2,1,p,mfl,&list);
599: dtest_special(p,list,rp);
600: }
601:
602: void dtest_special(f,list,pr)
603: P f;
604: ML list;
605: P *pr;
606: {
607: int n,np,bound,q;
608: int i,j,k,t,b0;
609: int *win;
610: P g,factor,cofactor;
611: Q csum,csumt;
612: DCP dc,dcf,dcf0;
613: LUM *c;
614: ML wlist;
615: static int nbits[16] = {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
616:
617: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
618: win = W_ALLOC(np+1);
619: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
620: wlist = W_MLALLOC(np); wlist->n = list->n;
621: wlist->mod = list->mod; wlist->bound = list->bound;
622: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
623: for ( g = f, i = 0; i < (1<<23); i++ ) {
624: #if 0
625: if ( !(i % 1000) )
626: fprintf(stderr,"i=%d\n",i);
627: #endif
628: t = i>>20; b0 = nbits[t];
629: if ( !b0 )
630: continue;
631: for ( j = 1; j < 6; j++ ) {
632: t = (i>>(20-4*j))&0xf;
633: if ( nbits[t] != b0 )
634: break;
635: }
636: if ( j != 6 )
637: continue;
638: for ( j = k = 0; j < 24; j++ )
639: if ( i & (1<<(23-j)) )
640: win[k++] = j;
641: if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
642: #if 0
643: fprintf(stderr,"success : i=%d\n",i);
644: #endif
645: *pr = factor; return;
646: }
647: }
648: *pr = f;
649: }
650:
651: void hensel_special(index,count,f,mfl,listp)
652: int index,count;
653: P f;
654: P *mfl;
655: ML *listp;
656: {
657: register int i,j;
658: int q,n,t,d,r,u,br,tmp,bound;
659: int *c,*p,*m,*w;
660: int **pp;
661: DCP dc;
662: ML blist,clist,bqlist,cqlist,rlist;
663: UM *b;
664: LUM fl,tl;
665: LUM *l;
666:
667: blist = MLALLOC(24); blist->n = 24; blist->mod = 5;
668: for ( i = 0; i < 24; i++ ) {
669: blist->c[i] = (pointer)UMALLOC(7);
670: ptoum(5,mfl[i],blist->c[i]);
671: }
672: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
673: n = bqlist->n; q = bqlist->mod;
674: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
675: if ( bound == 1 ) {
676: *listp = rlist = MLALLOC(n);
677: rlist->n = n; rlist->mod = q; rlist->bound = bound;
678: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
679: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
680: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
681: pp[j][0] = p[j];
682: }
683: } else {
684: W_LUMALLOC(UDEG(f),bound,fl);
685: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
686: }
687: }
688: #endif
689:
690: void Pftest(arg,rp)
691: NODE arg;
692: P *rp;
693: {
694: ML list;
695: P p;
696: P *mfl;
697:
698: p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
699: hensel_special(5,1,p,mfl,&list);
700: dtest_special(p,list,rp);
701: }
702:
703: int nbits(a)
704: int a;
705: {
706: int i,s;
707:
708: for ( i = 0, s = 0; a && (i < 20); i++, a >>= 1 )
709: if ( a & 1 ) s++;
710: return s;
711: }
712:
713: void dtest_special(f,list,pr)
714: P f;
715: ML list;
716: P *pr;
717: {
718: int n,np,bound,q;
719: int i,j,k,b0;
720: int *win;
721: P g,factor,cofactor;
722: Q csum,csumt;
723: LUM *c;
724: ML wlist;
725:
726: n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
727: win = W_ALLOC(np+1);
728: ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
729: wlist = W_MLALLOC(np); wlist->n = list->n;
730: wlist->mod = list->mod; wlist->bound = list->bound;
731: c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
732: for ( g = f, i = 0; i < (1<<19); i++ ) {
733: #if 0
734: if ( !(i % 10000) )
735: fprintf(stderr,"i=%d\n",i);
736: #endif
737: b0 = nbits(i>>10);
738: if ( !b0 || (nbits(i&0x3ff) != b0) )
739: continue;
740: for ( j = k = 0; j < 20; j++ )
741: if ( i & (1<<(19-j)) )
742: win[k++] = j;
743: if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
744: #if 0
745: fprintf(stderr,"success : i=%d\n",i);
746: #endif
747: *pr = factor; return;
748: }
749: }
750: *pr = f;
751: }
752:
753: void hensel_special(index,count,f,mfl,listp)
754: int index,count;
755: P f;
756: P *mfl;
757: ML *listp;
758: {
759: register int i,j;
760: int q,n,bound;
761: int *p;
762: int **pp;
763: ML blist,clist,bqlist,cqlist,rlist;
764: UM *b;
765: LUM fl,tl;
766: LUM *l;
767:
768: blist = MLALLOC(20); blist->n = 20; blist->mod = 11;
769: for ( i = 0; i < 20; i++ ) {
770: blist->c[i] = (pointer)UMALLOC(10);
771: ptoum(11,mfl[i],blist->c[i]);
772: }
773: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
774: n = bqlist->n; q = bqlist->mod;
775: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
776: if ( bound == 1 ) {
777: *listp = rlist = MLALLOC(n);
778: rlist->n = n; rlist->mod = q; rlist->bound = bound;
779: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
780: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
781: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
782: pp[j][0] = p[j];
783: }
784: } else {
785: W_LUMALLOC((int)UDEG(f),bound,fl);
786: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
787: }
788: }
789:
790: void Pnullspace(arg,rp)
791: NODE arg;
792: LIST *rp;
793: {
794: int i,j,n,mod;
795: MAT mat,r;
796: VECT u;
797: V v;
798: P p,z;
799: Q q;
800: UM **w;
801: UM mp;
802: P *t;
803: UM *s;
804: int *ind;
805: NODE n0,n1;
806:
807: mat = (MAT)ARG0(arg);
808: p = (P)ARG1(arg);
809: v = VR(p);
810: mod = QTOS((Q)ARG2(arg));
811: n = mat->row;
812: w = (UM **)almat_pointer(n,n);
813: for ( i = 0; i < n; i++ )
814: for ( j = 0, t = (P *)mat->body[i], s = w[i]; j < n; j++ ) {
815: ptomp(mod,t[j],&z);
816: s[j] = W_UMALLOC((z&&!NUM(z))?UDEG(z):0);
817: mptoum(z,s[j]);
818: }
819: mp = W_UMALLOC(UDEG(p)); ptoum(mod,p,mp);
820: ind = (int *)ALLOCA(n*sizeof(int));
821: nullspace(w,mp,mod,n,ind);
822: MKMAT(r,n,n);
823: for ( i = 0; i < n; i++ )
824: for ( j = 0, t = (P *)r->body[i], s = w[i]; j < n; j++ )
825: umtop(v,s[j],&t[j]);
826: MKVECT(u,n);
827: for ( i = 0; i < n; i++ ) {
828: STOQ(ind[i],q); u->body[i] = (pointer)q;
829: }
830: MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
831: }
832:
833: void nullspace(mat,p,mod,n,ind)
834: UM **mat;
835: UM p;
836: int mod,n;
837: int *ind;
838: {
839: int i,j,l,s,d;
840: UM q,w,w1,h,inv;
841: UM *t,*u;
842:
843: d = DEG(p); inv = W_UMALLOC(d); q = W_UMALLOC(2*d);
844: w = W_UMALLOC(2*d); w1 = W_UMALLOC(2*d); h = W_UMALLOC(d);
845: bzero(ind,n*sizeof(int));
846: ind[0] = 0;
847: for ( i = j = 0; j < n; i++, j++ ) {
848: for ( ; j < n; j++ ) {
849: for ( l = i; l < n; l++ )
850: if ( DEG(mat[l][j])>=0 )
851: break;
852: if ( l < n ) {
853: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
854: } else
855: ind[j] = 1;
856: }
857: if ( j == n )
858: break;
859: invum(mod,p,mat[i][j],inv);
860: for ( s = j, t = mat[i]; s < n; s++ ) {
861: mulum(mod,t[s],inv,w);
862: DEG(w) = divum(mod,w,p,q);
863: cpyum(w,t[s]);
864: }
865: for ( l = 0; l < n; l++ ) {
866: if ( l == i )
867: continue;
868: u = mat[l]; DEG(w) = -1; subum(mod,w,u[j],h);
869: for ( s = j; s < n; s++ ) {
870: mulum(mod,h,t[s],w); addum(mod,w,u[s],w1);
871: DEG(w1) = divum(mod,w1,p,q); cpyum(w1,u[s]);
872: }
873: }
874: }
875: }
876:
877: void Pnullspace_ff(arg,rp)
878: NODE arg;
879: LIST *rp;
880: {
881: int i,j,n;
882: Q mod;
883: MAT mat,r;
884: VECT u;
885: Q q;
886: Obj **w;
887: Obj *t;
888: Obj *s;
889: int *ind;
890: NODE n0,n1;
891:
892: mat = (MAT)ARG0(arg);
893: n = mat->row;
894: w = (Obj **)almat_pointer(n,n);
895: for ( i = 0; i < n; i++ )
896: for ( j = 0, t = (Obj *)mat->body[i], s = w[i]; j < n; j++ )
897: s[j] = t[j];
898: ind = (int *)ALLOCA(n*sizeof(int));
899: switch ( current_ff ) {
900: case FF_GFP:
901: nullspace_lm((LM **)w,n,ind); break;
902: case FF_GF2N:
903: nullspace_gf2n((GF2N **)w,n,ind); break;
904: case FF_GFPN:
905: nullspace_gfpn((GFPN **)w,n,ind); break;
1.5 noro 906: case FF_GFS:
907: nullspace_gfs((GFS **)w,n,ind); break;
1.1 noro 908: default:
909: error("nullspace_ff : current_ff is not set");
910: }
911: MKMAT(r,n,n);
912: for ( i = 0; i < n; i++ )
913: for ( j = 0, t = (Obj *)r->body[i], s = w[i]; j < n; j++ )
914: t[j] = s[j];
915: MKVECT(u,n);
916: for ( i = 0; i < n; i++ ) {
917: STOQ(ind[i],q); u->body[i] = (pointer)q;
918: }
919: MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
920: }
921:
922: void nullspace_lm(mat,n,ind)
923: LM **mat;
924: int n;
925: int *ind;
926: {
927: int i,j,l,s;
928: Q q,mod;
929: N lm;
930: LM w,w1,h,inv;
931: LM *t,*u;
932:
933: getmod_lm(&lm); NTOQ(lm,1,mod);
934:
935: bzero(ind,n*sizeof(int));
936: ind[0] = 0;
937: for ( i = j = 0; j < n; i++, j++ ) {
938: for ( ; j < n; j++ ) {
939: for ( l = i; l < n; l++ ) {
940: simplm(mat[l][j],&w); mat[l][j] = w;
941: if ( mat[l][j] )
942: break;
943: }
944: if ( l < n ) {
945: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
946: } else
947: ind[j] = 1;
948: }
949: if ( j == n )
950: break;
951: NTOQ(mat[i][j]->body,1,q); invl(q,mod,(Q *)&inv);
952: for ( s = j, t = mat[i]; s < n; s++ ) {
953: mullm(t[s],inv,&w); t[s] = w;
954: }
955: for ( l = 0; l < n; l++ ) {
956: if ( l == i )
957: continue;
958: u = mat[l]; chsgnlm(u[j],&h);
959: for ( s = j; s < n; s++ ) {
960: mullm(h,t[s],&w); addlm(w,u[s],&w1); u[s] = w1;
961: }
962: }
963: }
964: }
965:
966: void nullspace_gf2n(mat,n,ind)
967: GF2N **mat;
968: int n;
969: int *ind;
970: {
971: int i,j,l,s;
972: GF2N w,w1,h,inv;
973: GF2N *t,*u;
974: extern gf2n_lazy;
975:
976: bzero(ind,n*sizeof(int));
977: ind[0] = 0;
978: for ( i = j = 0; j < n; i++, j++ ) {
979: for ( ; j < n; j++ ) {
980: for ( l = i; l < n; l++ ) {
981: simpgf2n(mat[l][j],&w); mat[l][j] = w;
982: if ( mat[l][j] )
983: break;
984: }
985: if ( l < n ) {
986: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
987: } else
988: ind[j] = 1;
989: }
990: if ( j == n )
991: break;
992: invgf2n(mat[i][j],&inv);
993: for ( s = j, t = mat[i]; s < n; s++ ) {
994: mulgf2n(t[s],inv,&w); t[s] = w;
995: }
996: for ( l = 0; l < n; l++ ) {
997: if ( l == i )
998: continue;
999: u = mat[l]; h = u[j];
1000: for ( s = j; s < n; s++ ) {
1001: mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
1002: }
1003: }
1004: }
1005: }
1006:
1007: void nullspace_gfpn(mat,n,ind)
1008: GFPN **mat;
1009: int n;
1010: int *ind;
1011: {
1012: int i,j,l,s;
1013: GFPN w,w1,h,inv;
1014: GFPN *t,*u;
1015:
1016: bzero(ind,n*sizeof(int));
1017: ind[0] = 0;
1018: for ( i = j = 0; j < n; i++, j++ ) {
1019: for ( ; j < n; j++ ) {
1020: for ( l = i; l < n; l++ ) {
1021: simpgfpn(mat[l][j],&w); mat[l][j] = w;
1022: if ( mat[l][j] )
1023: break;
1024: }
1025: if ( l < n ) {
1026: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
1027: } else
1028: ind[j] = 1;
1029: }
1030: if ( j == n )
1031: break;
1032: divgfpn((GFPN)ONE,(GFPN)mat[i][j],&inv);
1033: for ( s = j, t = mat[i]; s < n; s++ ) {
1034: mulgfpn(t[s],inv,&w); t[s] = w;
1035: }
1036: for ( l = 0; l < n; l++ ) {
1037: if ( l == i )
1038: continue;
1039: u = mat[l]; chsgngfpn(u[j],&h);
1040: for ( s = j; s < n; s++ ) {
1041: mulgfpn(h,t[s],&w); addgfpn(w,u[s],&w1); u[s] = w1;
1042: }
1043: }
1044: }
1045: }
1.5 noro 1046:
1047: void nullspace_gfs(mat,n,ind)
1048: GFS **mat;
1049: int n;
1050: int *ind;
1051: {
1052: int i,j,l,s;
1053: GFS w,w1,h,inv;
1054: GFS *t,*u;
1055: GFS one;
1056:
1057: bzero(ind,n*sizeof(int));
1058: ind[0] = 0;
1059: mqtogfs(ONEM,&one);
1060:
1061: for ( i = j = 0; j < n; i++, j++ ) {
1062: for ( ; j < n; j++ ) {
1063: for ( l = i; l < n; l++ )
1064: if ( mat[l][j] )
1065: break;
1066: if ( l < n ) {
1067: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
1068: } else
1069: ind[j] = 1;
1070: }
1071: if ( j == n )
1072: break;
1073: divgfs(one,mat[i][j],&inv);
1074: for ( s = j, t = mat[i]; s < n; s++ ) {
1075: mulgfs(t[s],inv,&w); t[s] = w;
1076: }
1077: for ( l = 0; l < n; l++ ) {
1078: if ( l == i )
1079: continue;
1080: u = mat[l];
1081: chsgngfs(u[j],&h);
1082: for ( s = j; s < n; s++ ) {
1083: mulgfs(h,t[s],&w); addgfs(w,u[s],&w1); u[s] = w1;
1084: }
1085: }
1086: }
1087: }
1088:
1.1 noro 1089: /* p = a(0)vl[0]+a(1)vl[1]+...+a(m-1)vl[m-1]+a(m) -> array = [a(0) a(1) ... a(m)] */
1090:
1091: void linear_form_to_array(p,vl,m,array)
1092: P p;
1093: VL vl;
1094: int m;
1095: Num *array;
1096: {
1097: int i;
1098: DCP dc;
1099:
1100: bzero((char *)array,(m+1)*sizeof(Num *));
1101: for ( i = 0; p && vl; vl = NEXT(vl), i++ ) {
1102: if ( ID(p) == O_N )
1103: break;
1104: else if ( VR(p) == vl->v ) {
1105: dc = DC(p);
1106: array[i] = (Num)COEF(dc);
1107: dc = NEXT(dc);
1108: p = dc ? COEF(dc) : 0;
1109: }
1110: }
1111: array[m] = (Num)p;
1112: }
1113:
1114: void array_to_linear_form(array,vl,m,r)
1115: Num *array;
1116: VL vl;
1117: int m;
1118: P *r;
1119: {
1120: P t;
1121: DCP dc0,dc1;
1122:
1123: if ( !m )
1124: *r = (P)array[0];
1125: else {
1126: array_to_linear_form(array+1,NEXT(vl),m-1,&t);
1127: if ( !array[0] )
1128: *r = t;
1129: else {
1130: NEWDC(dc0); DEG(dc0) = ONE; COEF(dc0) = (P)array[0];
1131: if ( !t )
1132: NEXT(dc0) = 0;
1133: else {
1134: NEWDC(dc1); DEG(dc1) = 0; COEF(dc1) = t;
1135: NEXT(dc1) = 0;
1136: NEXT(dc0) = dc1;
1137: }
1138: MKP(vl->v,dc0,*r);
1139: }
1140: }
1141: }
1142:
1143: void Psolve_linear_equation_gf2n(arg,rp)
1144: NODE arg;
1145: LIST *rp;
1146: {
1147: NODE eqs,tn;
1148: VL vars,tvl;
1149: int i,j,n,m,dim,codim;
1150: GF2N **w;
1151: int *ind;
1152: NODE n0,n1;
1153:
1154: get_vars(ARG0(arg),&vars);
1155: eqs = BDY((LIST)ARG0(arg));
1156: for ( n = 0, tn = eqs; tn; tn = NEXT(tn), n++);
1157: for ( m = 0, tvl = vars; tvl; tvl = NEXT(tvl), m++);
1158: w = (GF2N **)almat_pointer(n,m+1);
1159: for ( i = 0, tn = eqs; i < n; i++, tn = NEXT(tn) )
1160: linear_form_to_array(BDY(tn),vars,m,(Num *)w[i]);
1161: ind = (int *)ALLOCA(m*sizeof(int));
1162: solve_linear_equation_gf2n(w,n,m,ind);
1163: for ( j = 0, dim = 0; j < m; j++ )
1164: if ( ind[j] )
1165: dim++;
1166: codim = m-dim;
1167: for ( i = codim; i < n; i++ )
1168: if ( w[i][m] ) {
1169: MKLIST(*rp,0); return;
1170: }
1171: for ( i = 0, n0 = 0; i < codim; i++ ) {
1172: NEXTNODE(n0,n1);
1173: array_to_linear_form((Num *)w[i],vars,m,(P *)&BDY(n1));
1174: }
1175: if ( n0 )
1176: NEXT(n1) = 0;
1177: MKLIST(*rp,n0);
1178: }
1179:
1180: void solve_linear_equation_gf2n(mat,n,m,ind)
1181: GF2N **mat;
1182: int n;
1183: int *ind;
1184: {
1185: int i,j,l,s;
1186: GF2N w,w1,h,inv;
1187: GF2N *t,*u;
1188: extern gf2n_lazy;
1189:
1190: bzero(ind,m*sizeof(int));
1191: ind[0] = 0;
1192: for ( i = j = 0; j < m; i++, j++ ) {
1193: for ( ; j < m; j++ ) {
1194: for ( l = i; l < n; l++ ) {
1195: simpgf2n(mat[l][j],&w); mat[l][j] = w;
1196: if ( mat[l][j] )
1197: break;
1198: }
1199: if ( l < n ) {
1200: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
1201: } else
1202: ind[j] = 1;
1203: }
1204: if ( j == m )
1205: break;
1206: invgf2n(mat[i][j],&inv);
1207: for ( s = j, t = mat[i]; s <= m; s++ ) {
1208: mulgf2n(t[s],inv,&w); t[s] = w;
1209: }
1210: for ( l = 0; l < n; l++ ) {
1211: if ( l == i )
1212: continue;
1213: u = mat[l]; h = u[j];
1214: for ( s = j; s <= m; s++ ) {
1215: mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
1216: }
1217: }
1218: }
1219: }
1220:
1221: /*
1222: void null_to_sol(mat,ind,mod,n,r)
1223: int **mat;
1224: int *ind;
1225: int mod,n;
1226: UM *r;
1227: {
1228: int i,j,k,l;
1229: int *c;
1230: UM w;
1231:
1232: for ( i = 0, l = 0; i < n; i++ ) {
1233: if ( !ind[i] )
1234: continue;
1235: w = UMALLOC(n);
1236: for ( j = k = 0, c = COEF(w); j < n; j++ )
1237: if ( ind[j] )
1238: c[j] = 0;
1239: else
1240: c[j] = mat[k++][i];
1241: c[i] = mod-1;
1242: for ( j = n; j >= 0; j-- )
1243: if ( c[j] )
1244: break;
1245: DEG(w) = j;
1246: r[l++] = w;
1247: }
1248: }
1249: */
1250:
1251: void showgfmat(mat,n)
1252: UM **mat;
1253: int n;
1254: {
1255: int i,j,k;
1256: int *c;
1257: UM p;
1258:
1259: for ( i = 0; i < n; i++ ) {
1260: for ( j = 0; j < n; j++ ) {
1261: p = mat[i][j];
1262: if ( DEG(p) < 0 )
1263: fprintf(asir_out,"0");
1264: else
1265: for ( p = mat[i][j], k = DEG(p), c = COEF(p); k >= 0; k-- ) {
1266: if ( c[k] )
1267: fprintf(asir_out,"+%d",c[k]);
1268: if ( k > 1 )
1269: fprintf(asir_out,"a^%d",k);
1270: else if ( k == 1 )
1271: fprintf(asir_out,"a",k);
1272: }
1273: fprintf(asir_out," ");
1274: }
1275: fprintf(asir_out,"\n");
1276: }
1277: }
1278:
1279: #if 0
1280: void Pgcda_mod(arg,rp)
1281: NODE arg;
1282: P *rp;
1283: {
1284: p1 = (P)ARG0(arg);
1285: p2 = (P)ARG1(arg);
1286: v = VR((P)ARG2(arg));
1287: d = (P)ARG3(arg);
1288: m = QTOS((Q)ARG4(arg));
1289: reordvar(CO,v,&vl);
1290: reorderp(vl,CO,p1,&t); ptomp(m,t,&m1);
1291: reorderp(vl,CO,p2,&t); ptomp(m,t,&m2);
1292: if ( NUM(m1) || NUM(m2) || VR(m1) != v || VR(m2) != v ) {
1293: *rp = ONE; return;
1294: }
1295: if ( deg(v,m1) >= deg(v,m2) ) {
1296: t = m1; m1 = m2; m2 = t;
1297: }
1298: while ( 1 ) {
1299: inva_mod(COEF(DC(m2)),d,m,&inv);
1300: NEWDC(dc); NEXT(dc) = 0; MKP(v,dc,h);
1301: d0 = deg(v,m1)-deg(v,m2); STOQ(d0,DEG(dc));
1302: mulgq(m,d,inv,COEF(DC(m1)),&COEF(dc));
1303: mulgp(vl,m,d,m2,h,&t); subgp(vl,m,d,m1,t,&s);
1304: }
1305: }
1306: #endif
1307:
1308: void Ppwr_mod(arg,rp)
1309: NODE arg;
1310: P *rp;
1311: {
1312: P p,a,d,r;
1313: int m;
1314: Q q;
1315: N n;
1316:
1317: m = QTOS((Q)ARG4(arg)); q = (Q)ARG5(arg); n = q ? NM(q) : 0;
1318: ptomp(m,(P)ARG0(arg),&p); ptomp(m,(P)ARG1(arg),&a);
1319: ptomp(m,(P)ARG3(arg),&d);
1320: pwr_mod(p,a,VR((P)ARG2(arg)),d,m,n,&r);
1321: mptop(r,rp);
1322: }
1323:
1324: void pwr_mod(p,a,v,d,m,n,rp)
1325: P p,a,d;
1326: V v;
1327: int m;
1328: N n;
1329: P *rp;
1330: {
1331: int b;
1332: P t,s,r;
1333: N n1;
1334:
1335: if ( !n )
1336: *rp = (P)ONEM;
1337: else if ( UNIN(n) )
1338: *rp = p;
1339: else {
1340: b = divin(n,2,&n1);
1341: pwr_mod(p,a,v,d,m,n1,&t);
1342: mulmp(CO,m,t,t,&s); rem_mod(s,a,v,d,m,&r);
1343: if ( b ) {
1344: mulmp(CO,m,r,p,&t); rem_mod(t,a,v,d,m,rp);
1345: } else
1346: *rp = r;
1347: }
1348: }
1349:
1350: void rem_mod(p,a,v,d,m,rp)
1351: P p,a,d;
1352: V v;
1353: int m;
1354: P *rp;
1355: {
1356: P tmp,r1,r2;
1357:
1358: divsrmp(CO,m,p,d,&tmp,&r1);
1359: divsrmp(CO,m,r1,a,&tmp,&r2);
1360: divsrmp(CO,m,r2,d,&tmp,rp);
1361: }