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