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