![[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: }