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