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