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