![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to H.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to H.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/H.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3: #include "inline.h"
4: #include "base.h"
5: #include <math.h>
6:
7: void modfctrp(P,int,int,DCP *);
8: void gensqfrum(int,UM,struct oDUM *);
9: void srchum(int,UM,UM,UM);
10: UM *resberle(int,UM,UM *);
11: int substum(int,UM,int);
12: void ddd(int,UM,UM *);
13: void canzas(int,UM,int,UM *,UM *);
14: int divum(int,UM,UM,UM);
15: void randum(int,int,UM);
16: void pwrmodum(int,UM,int,UM,UM);
17: void spwrum(int,UM,UM *,UM,N,UM);
18: void spwrum0(int,UM,UM,N,UM);
19: void mulum(int,UM,UM,UM);
20: void mulsum(int,UM,int,UM);
21: void gcdum(int,UM,UM,UM);
22: void mult_mod_tab(UM,int,UM *,UM,int);
23: void make_qmat(UM,int,UM *,int ***);
24: void null_mod(int **,int,int,int *);
25: void null_to_sol(int **,int *,int,int,UM *);
26: void newddd(int,UM,UM *);
27: int irred_check(UM,int);
28: int berlekamp(UM,int,int,UM *,UM *);
29: int nfctr_mod(UM,int);
30: void minipoly_mod(int,UM,UM,UM);
31: int find_root(int,UM,int *);
32: void showum(UM);
33: void showumat(int **,int);
34: #if 1
35: #define Mulum mulum
36: #define Divum divum
37: #define Mulsum mulsum
38: #define Gcdum gcdum
39: #endif
40:
41: #define FCTR 0
42: #define SQFR 1
43: #define DDD 2
44: #define NEWDDD 3
45:
46: LUM LUMALLOC();
47:
48: struct p_pair {
49: UM p0;
50: UM p1;
51: struct p_pair *next;
52: };
53:
54: void lnf_mod(int,int,UM,UM,struct p_pair *,UM,UM);
55:
56: void berle(index,count,f,listp)
57: int index,count;
58: P f;
59: ML *listp;
60: {
61: UM wf,wf1,wf2,wfs,gcd;
62: ML flist;
63: int fn,fn1,fm,m,n,fhd;
64: register int i,j,inv,hd,*ptr,*ptr1;
65:
66: n = UDEG(f);
67: wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); wf2 = W_UMALLOC(n);
68: wfs = W_UMALLOC(n); gcd = W_UMALLOC(n);
69: for ( j = 0, fn = n + 1; (j < count) && (fn > 1); ) {
70: m = sprime[index++];
71: if ( !rem(NM((Q)UCOEF(f)),m) )
72: continue;
73: ptoum(m,f,wf); cpyum(wf,wf1);
74: diffum(m,wf1,wf2); gcdum(m,wf1,wf2,gcd);
75: if ( DEG(gcd) > 0 )
76: continue;
77: hd = COEF(wf)[n]; inv = invm(hd,m);
78: for ( i = n, ptr = COEF(wf); i >= 0; i-- )
79: ptr[i] = ( ptr[i] * inv ) % m;
80: fn1 = berlecnt(m,wf);
81: if ( fn1 < fn ) {
82: fn = fn1; fm = m; fhd = hd;
83: for ( i = n, ptr = COEF(wf), ptr1 = COEF(wfs); i >= 0; i-- )
84: ptr1[i] = ptr[i];
85: }
86: j++;
87: }
88: DEG(wfs) = n;
89: *listp = flist = MLALLOC(fn); flist->n = fn; flist->mod = fm;
90: /* berlemain(fm,wfs,(UM *)flist->c); */
91: if ( fm == 2 )
92: berlemain(fm,wfs,(UM *)flist->c);
93: else
94: newddd(fm,wfs,(UM *)flist->c);
95: for ( i = DEG((UM)(flist->c[0])),
96: ptr = COEF((UM)(flist->c[0])),
97: hd = fhd, m = fm; i >= 0; i-- )
98: ptr[i] = ( ptr[i] * hd ) % m;
99: }
100:
101: int berlecnt(mod,f)
102: register int mod;
103: UM f;
104: {
105: register int i,j,**c;
106: int d,dr,n;
107: UM w,q;
108: int **almat();
109:
110: n = DEG(f); c = almat(n,n);
111: w = W_UMALLOC(mod + n); q = W_UMALLOC(mod + n);
112: for ( i = 1; ( d = ( mod * i ) ) < n; i++ )
113: c[d][i - 1] = 1;
114: DEG(w) = d; COEF(w)[d] = 1;
115: for ( j = d - 1; j >= 0; j-- )
116: COEF(w)[j] = 0;
117: for ( ; ( i < n ) && ( ( dr = divum(mod,w,f,q) ) != -1 ); i++ ) {
118: for ( j = dr; j >= 0; j-- )
119: COEF(w)[j + mod] = c[j][i - 1] = COEF(w)[j];
120: for ( j = mod - 1; j >= 0; j-- )
121: COEF(w)[j] = 0;
122: DEG(w) = dr + mod;
123: }
124: for ( i = 1; i < n; i++ )
125: c[i][i - 1] = ( c[i][i - 1] + mod - 1 ) % mod;
126: return berlecntmain(mod,n,n-1,c);
127: }
128:
129: /* XXX berlecntmain should not be used for large mod */
130:
131: int berlecntmain(mod,n,m,c)
132: register int mod;
133: int n,m;
134: register int **c;
135: {
136: register int *p1,*p2,i,j,k,l,a;
137: int *tmp,inv;
138: int cfs;
139:
140: for ( cfs = 1, j = k = 0; j < m; j++ ) {
141: for ( i = k; ( n > i ) && ( c[i][j] == 0 ); i++ );
142: if ( i == n ) {
143: cfs++; continue;
144: }
145: if ( i != k ) {
146: tmp = c[i]; c[i] = c[k]; c[k] = tmp;
147: }
148: p1 = c[k]; inv = invm((p1[j] + mod) % mod,mod);
149: for ( l = j; l < m; l++ )
150: p1[l] = ( p1[l] * inv ) % mod;
151: for ( i = k + 1; i < n; c[i][j] = 0, i++ )
152: if ( i != k && ( a = -c[i][j] ) )
153: for ( l = j + 1, p2 = c[i]; l < m; l++ )
154: p2[l] = (a*p1[l] + p2[l]) % mod;
155: k++;
156: }
157: return ( cfs );
158: }
159:
160: UM *berlemain(mod,f,fp)
161: register int mod;
162: UM f;
163: UM *fp;
164: {
165: UM wg,ws,wf,f0,gcd,q;
166: int n;
167: register int i;
168:
169: n = DEG(f); wg = W_UMALLOC(n); mini(mod,f,wg);
170: if ( DEG(wg) <= 0 ) {
171: f0 = UMALLOC(n); cpyum(f,f0); *fp++ = f0;
172: return ( fp );
173: }
174: f0 = W_UMALLOC(n); cpyum(f,f0);
175: ws = W_UMALLOC(n); wf = W_UMALLOC(n);
176: q = W_UMALLOC(n); gcd = W_UMALLOC(n);
177: for ( i = 0; i < mod; i++ ) {
178: cpyum(f0,wf); cpyum(wg,ws);
179: COEF(ws)[0] = ( COEF(ws)[0] + mod - i ) % mod;
180: gcdum(mod,wf,ws,gcd);
181: if ( DEG(gcd) > 0 ) {
182: if ( DEG(gcd) < n ) {
183: divum(mod,f0,gcd,q); f0 = q; fp = berlemain(mod,gcd,fp);
184: }
185: break;
186: }
187: }
188: fp = berlemain(mod,f0,fp);
189: return ( fp );
190: }
191:
192: void hensel(index,count,f,listp)
193: int index,count;
194: P f;
195: ML *listp;
196: {
197: register int i,j;
198: int q,n,bound;
199: int *p;
200: int **pp;
201: ML blist,clist,bqlist,cqlist,rlist;
202: UM *b;
203: LUM fl,tl;
204: LUM *l;
205:
206: if ( UDEG(f) == 1 ) {
207: *listp = blist = MLALLOC(1); blist->n = 1; blist->c[0] = 0;
208: return;
209: }
210: berle(index,count,f,&blist);
211: if ( blist->n == 1 ) {
212: *listp = blist = MLALLOC(1); blist->n = 1; blist->c[0] = 0;
213: return;
214: }
215: gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
216: n = bqlist->n; q = bqlist->mod;
217: bqlist->bound = cqlist->bound = bound = mignotte(q,f);
218: if ( bound == 1 ) {
219: *listp = rlist = MLALLOC(n);
220: rlist->n = n; rlist->mod = q; rlist->bound = bound;
221: for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
222: tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
223: for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
224: pp[j][0] = p[j];
225: }
226: } else {
227: W_LUMALLOC((int)UDEG(f),bound,fl);
228: ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
229: }
230: }
231:
232: void hsq(index,count,f,nindex,dcp)
233: int index,count,*nindex;
234: P f;
235: DCP *dcp;
236: {
237: register int i,j,k;
238: register int **pp,**fpp;
239: register int *px,*py;
240: int **wpp;
241: int n,dr,tmp,m,b,e,np,dt;
242: LUM fpa,wb0,wb1,lcpa,tpa,tlum;
243: struct oDUM *dct;
244: UM wt,wq0,wq,wr,wm,wm0,wa,ws,wb;
245: LUM *llist,*ll;
246: UM *dlist,*l,*c;
247: ML list,fp,cfp;
248: DCP dc;
249:
250: sqfrum(index,count,f,nindex,&dct,&fp);
251: np = fp->n; m = fp->mod;
252: if ( ( np == 1 ) && ( dct[0].n == 1 ) ) {
253: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0; *dcp = dc;
254: return;
255: }
256: for ( i = 0, dt = 0; i < np; i++ )
257: dt = MAX(DEG(dct[i].f),dt);
258: b = mig(m,dt,f); fp->bound = b;
259: if ( np == 1 ) {
260: nthrootchk(f,dct,fp,dcp);
261: return;
262: }
263: list = W_MLALLOC(np); list->n = np; list->mod = m; list->bound = 1;
264: for ( i = 0, ll = (LUM *)list->c; i < np; i++ ) {
265: W_LUMALLOC(DEG(dct[i].f),b,ll[i]);
266: for ( j = 0, px = COEF(dct[i].f), pp = COEF(ll[i]);
267: j <= DEG(ll[i]); j++ )
268: pp[j][0] = px[j];
269: }
270: dtestsql(f,list,dct,&dc);
271: if ( dc ) {
272: *dcp = dc;
273: return;
274: }
275: n = UDEG(f);
276: W_LUMALLOC(n,b,fpa); W_LUMALLOC(0,b,lcpa);
277: W_LUMALLOC(n,b,wb0); W_LUMALLOC(n,b,wb1);
278: W_LUMALLOC(n,b,tpa);
279: wt = W_UMALLOC(n); ws = W_UMALLOC(n);
280: wr = W_UMALLOC(n);
281: wq = W_UMALLOC(2*n); wq0 = W_UMALLOC(n);
282: wm = W_UMALLOC(2*n); wm0 = W_UMALLOC(2*n);
283: wa = W_UMALLOC(2*n);
284: ptolum(m,b,f,fpa); DEG(lcpa) = 0;
285: for ( i = 0, pp = COEF(lcpa), fpp = COEF(fpa); i < b; i++ )
286: pp[0][i] = fpp[n][i];
287: gcdgen(f,fp,&cfp);
288: llist = (LUM *) ALLOCA(np*sizeof(LUM));
289: dlist = (UM *) ALLOCA(np*sizeof(UM));
290: l = (UM *)fp->c; c = (UM *)cfp->c;
291: for ( i = 0; i < np; i++ ) {
292: W_LUMALLOC(DEG(l[i]),b,llist[i]);
293: for ( j = DEG(l[i]), pp = COEF(llist[i]), px = COEF(l[i]); j >= 0; j-- )
294: pp[j][0] = px[j];
295: if ( ( e = dct[i].n ) != 1 ) {
296: wb = dct[i].f;
297: dlist[i] = W_UMALLOC(DEG(wb)*e); cpyum(l[i],dlist[i]);
298: divum(m,dlist[i],wb,wq); DEG(dlist[i])= DEG(wq);
299: for ( k = 0; k <= DEG(wq); k++ )
300: COEF(dlist[i])[k] = dmb(m,COEF(wq)[k],e,&tmp);
301: }
302: }
303: for ( i = 1; i < b; i++ ) {
304: mullum(m,i+1,lcpa,llist[0],wb0);
305: for ( j = 1; j < np; j++ ) {
306: mullum(m,i+1,llist[j],wb0,wb1);
307: tlum = wb0; wb0 = wb1; wb1 = tlum;
308: }
309: for ( j = n, px = COEF(wt), pp = COEF(fpa), wpp = COEF(wb0);
310: j >= 0; j-- )
311: px[j] = ( pp[j][i] - wpp[j][i] + m ) % m;
312: degum(wt,n);
313: for ( j = n, px = COEF(wq0); j >= 0; j-- )
314: px[j] = 0;
315: for ( j = 1; j < np; j++ ) {
316: mulum(m,wt,c[j],wm); dr = divum(m,wm,l[j],wq);
317: for ( k = DEG(wq), px = COEF(wq0), py = COEF(wq); k >= 0; k-- )
318: px[k] = ( px[k] + py[k] ) % m;
319: for ( k = dr, pp = COEF(llist[j]), px = COEF(wm); k >= 0; k-- )
320: pp[k][i] = px[k];
321: }
322: degum(wq0,n); mulum(m,wq0,l[0],wm);
323: mulum(m,wt,c[0],wm0); addum(m,wm,wm0,wa);
324: for ( j = DEG(wa), pp = COEF(llist[0]), px = COEF(wa); j >= 0; j-- )
325: pp[j][i] = px[j];
326: for ( j = n, px = COEF(wq0); j >= 0; j-- )
327: px[j] = 0;
328: for ( j = 0; j < np; j++ )
329: if ( dct[j].n == 1 )
330: for ( k = 0,
331: pp = COEF(llist[j]),
332: wpp = COEF(((LUM *)list->c)[j]);
333: k <= DEG(llist[j]); k++ )
334: wpp[k][i] = pp[k][i];
335: else {
336: pwrlum(m,i+1,((LUM *)list->c)[j],dct[j].n,tpa);
337: for ( k = 0,
338: pp = COEF(llist[j]),
339: wpp = COEF(tpa);
340: k <= DEG(l[j]); k++ )
341: COEF(wt)[k] = (pp[k][i]-wpp[k][i]+m)%m;
342: degum(wt,DEG(l[j])); dr = divum(m,wt,dlist[j],ws);
343: if ( dr >= 0 ) {
344: *dcp = 0;
345: return;
346: } else
347: for ( k = 0,
348: pp = COEF(((LUM *)list->c)[j]);
349: k <= DEG(ws); k++ )
350: pp[k][i] = COEF(ws)[k];
351: }
352: list->bound = i+1; dtestsql(f,list,dct,&dc);
353: if ( dc ) {
354: *dcp = dc;
355: return;
356: }
357: }
358: *dcp = 0;
359: }
360:
361: void gcdgen(f,blist,clistp)
362: P f;
363: ML blist,*clistp;
364: {
365: register int i;
366: int n,d,mod,np;
367: UM wf,wm,wx,wy,wu,wv,wa,wb,wg,q,tum;
368: UM *in,*out;
369: ML clist;
370:
371: n = UDEG(f); mod = blist->mod; np = blist->n;
372: d = 2*n;
373: q = W_UMALLOC(d); wf = W_UMALLOC(d);
374: wm = W_UMALLOC(d); wx = W_UMALLOC(d);
375: wy = W_UMALLOC(d); wu = W_UMALLOC(d);
376: wv = W_UMALLOC(d); wg = W_UMALLOC(d);
377: wa = W_UMALLOC(d); wb = W_UMALLOC(d);
378: ptoum(mod,f,wf); DEG(wg) = 0; COEF(wg)[0] = 1;
379: *clistp = clist = MLALLOC(np); clist->mod = mod; clist->n = np;
380: for ( i = 0, in = (UM *)blist->c, out = (UM *)clist->c; i < np; i++ ) {
381: divum(mod,wf,in[i],q); tum = wf; wf = q; q = tum;
382: cpyum(wf,wx); cpyum(in[i],wy);
383: eucum(mod,wx,wy,wa,wb); mulum(mod,wa,wg,wm);
384: DEG(wm) = divum(mod,wm,in[i],q); out[i] = UMALLOC(DEG(wm));
385: cpyum(wm,out[i]); mulum(mod,q,wf,wu);
386: mulum(mod,wg,wb,wv); addum(mod,wu,wv,wg);
387: }
388: }
389:
390: /*
391: henprep(f,&blist,&clist,&bqlist,&cqlist);
392: */
393:
394: void henprep(f,blist,clist,bqlistp,cqlistp)
395: P f;
396: ML blist,clist,*bqlistp,*cqlistp;
397: {
398: register int i,j,k,*px,*py,*pz;
399: int n,pmax,dr,tmp,p,p1,mod,np,b,q;
400: UM w,wm,wn,wa,wt,wq,wf,quot,tum,*in,*inc,*out,*outc;
401: ML bqlist,cqlist;
402:
403: n = UDEG(f); p = mod = blist->mod; np = blist->n;
404: /* for ( b = 1, q = mod; q <= (unsigned int)(LBASE / (L)mod); q *= mod, b++ ); */
405: for ( b = 1, q = mod; q <= ((1<<27) / mod); q *= mod, b++ );
406: w = W_UMALLOC(n); ptoum(q,f,w);
407: wm = W_UMALLOC(2*n); wn = W_UMALLOC(2*n);
408: wa = W_UMALLOC(2*n); wt = W_UMALLOC(2*n);
409: wq = W_UMALLOC(2*n); wf = W_UMALLOC(2*n);
410: quot = W_UMALLOC(2*n);
411: *bqlistp = bqlist = MLALLOC(np); *cqlistp = cqlist = MLALLOC(np);
412: for ( i = 0; i < n+2; i++ )
413: COEF(wq)[i] = 0;
414: for ( i = 0,
415: in = (UM *)blist->c, inc = (UM *)clist->c,
416: out = (UM *)bqlist->c, outc = (UM *)cqlist->c;
417: i < np; i++ ) {
418: out[i] = C_UMALLOC(n+1); cpyum(in[i],out[i]);
419: outc[i] = C_UMALLOC(n+1); cpyum(inc[i],outc[i]);
420: }
421: for ( pmax = 1, i = b; i > 0; i-- )
422: pmax *= mod;
423: for ( i = 1; i < b; i++, p = p1 ) {
424: cpyum(out[0],wm);
425: for ( j = 1; j < np; j++ ) {
426: mulum(pmax,wm,out[j],wn);
427: tum = wm; wm = wn; wn = tum;
428: }
429: for ( j = n, px = COEF(w), py = COEF(wm), pz = COEF(wt); j >= 0; j-- ) {
430: tmp = ( ( px[j] - py[j] ) / p ) % mod;
431: pz[j] = ( tmp >= 0? tmp : tmp + mod );
432: }
433: degum(wt,n);
434: for ( j = 1; j < np; j++ ) {
435: mulum(mod,wt,inc[j],wm); dr = divum(mod,wm,in[j],quot);
436: for ( k = DEG(quot); k >= 0; k-- )
437: COEF(wq)[k] = ( COEF(wq)[k] + COEF(quot)[k] ) % mod;
438: for ( k = dr, px = COEF(out[j]), py = COEF(wm); k >= 0; k-- )
439: px[k] += p * py[k];
440: }
441: degum(wq,n); mulum(mod,wq,in[0],wm);
442: mulum(mod,wt,inc[0],wn); addum(mod,wm,wn,wa);
443: for ( j = DEG(wa), px = COEF(out[0]), py = COEF(wa); j >= 0; j-- )
444: px[j] += p * py[j];
445: for ( j = n, px = COEF(wq); j >= 0; j-- )
446: px[j] = 0;
447: p1 = p * mod;
448: for ( j = n, px = COEF(wt); j >= 1; j-- )
449: px[j] = 0;
450: px[0] = 1;
451: for ( j = 0; j < np; j++ ) {
452: cpyum(w,wf);
453: for ( k = DEG(wf), px = COEF(wf); k >= 0; k-- )
454: px[k] %= p1;
455: divum(p1,wf,out[j],quot); mulum(p1,outc[j],quot,wm);
456: for ( k = DEG(wm), px = COEF(wt), py = COEF(wm); k >= 0; k-- )
457: px[k] = ( px[k] - py[k] ) % p1;
458: }
459: degum(wt,n);
460: for ( j = DEG(wt), px = COEF(wt); j >= 0; j-- )
461: px[j] = ((tmp=(px[j]/p)%mod)>= 0?tmp:tmp + mod);
462: for ( j = 0; j < np; j++ ) {
463: mulum(mod,wt,outc[j],wm); dr = divum(mod,wm,in[j],quot);
464: for ( k = dr, px = COEF(outc[j]), py = COEF(wm); k >= 0; k-- )
465: px[k] += p * py[k];
466: degum(outc[j],MAX(DEG(outc[j]),dr));
467: }
468: }
469: bqlist->n = cqlist->n = np;
470: bqlist->mod = cqlist->mod = q;
471: }
472:
473: /*
474: henmain(fl,bqlist,cqlist,listp)
475: */
476:
477: void henmain(f,bqlist,cqlist,listp)
478: LUM f;
479: ML bqlist,cqlist,*listp;
480: {
481: register int i,j,k;
482: int *px,*py;
483: int **pp,**pp1;
484: int n,np,mod,bound,dr,tmp;
485: UM wt,wq0,wq,wr,wm,wm0,wa,q;
486: LUM wb0,wb1,tlum;
487: UM *b,*c;
488: LUM *l;
489: ML list;
490:
491: n = DEG(f); np = bqlist->n; mod = bqlist->mod; bound = bqlist->bound;
492: *listp = list = MLALLOC(n);
493: list->n = np; list->mod = mod; list->bound = bound;
494: W_LUMALLOC(n,bound,wb0); W_LUMALLOC(n,bound,wb1);
495: wt = W_UMALLOC(n); wq0 = W_UMALLOC(n); wq = W_UMALLOC(n);
496: wr = W_UMALLOC(n); wm = W_UMALLOC(2*n); wm0 = W_UMALLOC(2*n);
497: wa = W_UMALLOC(2*n); q = W_UMALLOC(2*n);
498: b = (UM *)bqlist->c; c = (UM *)cqlist->c; l = (LUM *)list->c;
499: for ( i = 0; i < np; i++ ) {
500: l[i] = LUMALLOC(DEG(b[i]),bound);
501: for ( j = DEG(b[i]), pp = COEF(l[i]), px = COEF(b[i]); j >= 0; j-- )
502: pp[j][0] = px[j];
503: }
504: for ( i = 1; i < bound; i++ ) {
505: mullum(mod,i+1,l[0],l[1],wb0);
506: for ( j = 2; j < np; j++ ) {
507: mullum(mod,i+1,l[j],wb0,wb1);
508: tlum = wb0; wb0 = wb1; wb1 = tlum;
509: }
510: for ( j = n, px = COEF(wt); j >= 0; j-- )
511: px[j] = 0;
512: for ( j = n, pp = COEF(f), pp1 = COEF(wb0); j >= 0; j-- ) {
513: tmp = ( pp[j][i] - pp1[j][i] ) % mod;
514: COEF(wt)[j] = ( tmp < 0 ? tmp + mod : tmp );
515: }
516: degum(wt,n);
517: for ( j = n, px = COEF(wq0); j >= 0; j-- )
518: px[j] = 0;
519: for ( j = 1; j < np; j++ ) {
520: mulum(mod,wt,c[j],wm); dr = divum(mod,wm,b[j],q);
521: for ( k = DEG(q), px = COEF(wq0), py = COEF(q); k >= 0; k-- )
522: px[k] = ( px[k] + py[k] ) % mod;
523: for ( k = dr, pp = COEF(l[j]), px = COEF(wm); k >= 0; k-- )
524: pp[k][i] = px[k];
525: }
526: degum(wq0,n); mulum(mod,wq0,b[0],wm);
527: mulum(mod,wt,c[0],wm0); addum(mod,wm,wm0,wa);
528: for ( j = DEG(wa), pp = COEF(l[0]), px = COEF(wa); j >= 0; j-- )
529: pp[j][i] = px[j];
530: for ( j = n, px = COEF(wq0); j >= 0; j-- )
531: px[j] = 0;
532: }
533: }
534:
535: static double M;
536: static int E;
537:
538: int mignotte(q,f)
539: int q;
540: P f;
541: {
542: int p;
543: unsigned int *b;
544: N c;
545: DCP dc;
546:
547: for ( dc = DC(f), M = 0, E = 0; dc; dc = NEXT(dc) ) {
548: c = NM((Q)COEF(dc)); p = PL(c); b = BD(c);
549: sqad(b[p-1],(p-1)*BSH);
550: }
551: if (E % 2) M *= 2; M = ceil(sqrt(M)); E /= 2;
552: c = NM((Q)COEF(DC(f))); p = PL(c); M *= ((double)BD(c)[p-1]+1.0); E += (p-1) * BSH;
553: return (int)ceil( (0.31*(E+UDEG(f)+1)+log10((double)M)) / log10((double)q) );
554: }
555:
556: int mig(q,d,f)
557: int q,d;
558: P f;
559: {
560: int p;
561: unsigned int *b;
562: N c;
563: DCP dc;
564:
565: for ( dc = DC(f), M = 0, E = 0; dc; dc = NEXT(dc) ) {
566: c = NM((Q)COEF(dc)); p = PL(c); b = BD(c);
567: sqad(b[p-1],(p-1)*BSH);
568: }
569: if (E % 2) M *= 2; M = ceil(sqrt(M)); E /= 2;
570: c = NM((Q)COEF(DC(f))); p = PL(c);
571: M *= (BD(c)[p-1]+1); E += (p-1) * BSH;
572: return (int)ceil( (0.31*(E+d+1)+log10((double)M)) / log10((double)q) );
573: }
574:
575: void sqad(man,exp)
576: unsigned int man;
577: int exp;
578: {
579: int e,sqe;
580: unsigned int t;
581: double man1,d,sqm;
582: int diff;
583:
584: if ( man == BMASK ) {
585: e = BSH; man1 = 1.0;
586: } else {
587: man += 1;
588: for ( e = 0, t = man; t; e++, t >>= 1 );
589: e--; d = (double)(1<<e);
590: man1 = ((double)man)/d;
591: }
592: exp += e; sqm = man1 * man1; sqe = 2 * exp;
593: if ( sqm >= 2.0 ) {
594: sqm /= 2.0; sqe++;
595: }
596: diff = E - sqe;
597: if ( diff > 18 )
598: return;
599: if ( diff < -18 ) {
600: M = sqm; E = sqe;
601: return;
602: }
603: if ( diff >= 0 )
604: M += (sqm / (double)(1<<diff));
605: else {
606: M = ( ( M / (double)(1<<-diff)) + sqm ); E = sqe;
607: }
608: if ( M >= 2.0 ) {
609: M /= 2.0; E++;
610: }
611: }
612:
613: void ptolum(q,bound,f,fl)
614: int q,bound;
615: P f;
616: LUM fl;
617: {
618: DCP dc;
619: int i,j;
620: int **pp;
621: int d,br,s;
622: unsigned int r;
623: int *c;
624: unsigned int *m,*w;
625:
626: for ( dc = DC(f), pp = COEF(fl); dc; dc = NEXT(dc) ) {
627: d = PL(NM((Q)COEF(dc))); m = BD(NM((Q)COEF(dc)));
628: c = pp[QTOS(DEG(dc))]; w = (unsigned int *)W_ALLOC(d);
629: for ( i = 0; i < d; i++ )
630: w[i] = m[i];
631: for ( i = 0; d >= 1; ) {
632: for ( j = d - 1, r = 0; j >= 0; j-- ) {
633: DSAB(q,r,w[j],w[j],r)
634: }
635: c[i++] = (int)r;
636: if ( !w[d-1] )
637: d--;
638: }
639: if ( SGN((Q)COEF(dc)) < 0 )
640: for (i = 0, br = 0; i < bound; i++ )
641: if ( ( s = -(c[i] + br) ) < 0 ) {
642: c[i] = s + q; br = 1;
643: } else {
644: c[i] = 0; br = 0;
645: }
646: }
647: }
648:
649: void modfctrp(p,mod,flag,dcp)
650: P p;
651: int mod,flag;
652: DCP *dcp;
653: {
654: int cm,n,i,j,k;
655: DCP dc,dc0;
656: P zp;
657: Q c,q;
658: UM mp;
659: UM *tl;
660: struct oDUM *udc,*udc1;
661:
662: if ( !p ) {
663: *dcp = 0; return;
664: }
665: ptozp(p,1,&c,&zp);
666: if ( DN(c) || !(cm = rem(NM(c),mod)) ) {
667: *dcp = 0; return;
668: }
669: mp = W_UMALLOC(UDEG(p));
670: ptoum(mod,zp,mp);
671: if ( (n = DEG(mp)) < 0 ) {
672: *dcp = 0; return;
673: } else if ( n == 0 ) {
674: cm = dmar(cm,COEF(mp)[0],0,mod); STOQ(cm,q);
675: NEWDC(dc); COEF(dc) = (P)q; DEG(dc) = ONE;
676: NEXT(dc) = 0; *dcp = dc;
677: return;
678: }
679: if ( COEF(mp)[n] != 1 ) {
680: cm = dmar(cm,COEF(mp)[n],0,mod);
681: i = invm(COEF(mp)[n],mod);
682: for ( j = 0; j <= n; j++ )
683: COEF(mp)[j] = dmar(COEF(mp)[j],i,0,mod);
684: }
685: W_CALLOC(n+1,struct oDUM,udc);
686: gensqfrum(mod,mp,udc);
687: switch ( flag ) {
688: case FCTR:
689: tl = (UM *)ALLOCA((n+1)*sizeof(UM));
690: W_CALLOC(DEG(mp)+1,struct oDUM,udc1);
691: for ( i = 0,j = 0; udc[i].f; i++ )
692: if ( DEG(udc[i].f) == 1 ) {
693: udc1[j].f = udc[i].f; udc1[j].n = udc[i].n; j++;
694: } else {
695: bzero((char *)tl,(n+1)*sizeof(UM));
696: berlemain(mod,udc[i].f,tl);
697: for ( k = 0; tl[k]; k++, j++ ) {
698: udc1[j].f = tl[k]; udc1[j].n = udc[i].n;
699: }
700: }
701: udc = udc1; break;
702: case SQFR:
703: break;
704: case DDD:
705: tl = (UM *)ALLOCA((n+1)*sizeof(UM));
706: W_CALLOC(DEG(mp)+1,struct oDUM,udc1);
707: for ( i = 0,j = 0; udc[i].f; i++ )
708: if ( DEG(udc[i].f) == 1 ) {
709: udc1[j].f = udc[i].f; udc1[j].n = udc[i].n; j++;
710: } else {
711: bzero((char *)tl,(n+1)*sizeof(UM));
712: ddd(mod,udc[i].f,tl);
713: for ( k = 0; tl[k]; k++, j++ ) {
714: udc1[j].f = tl[k]; udc1[j].n = udc[i].n;
715: }
716: }
717: udc = udc1; break;
718: case NEWDDD:
719: tl = (UM *)ALLOCA((n+1)*sizeof(UM));
720: W_CALLOC(DEG(mp)+1,struct oDUM,udc1);
721: for ( i = 0,j = 0; udc[i].f; i++ )
722: if ( DEG(udc[i].f) == 1 ) {
723: udc1[j].f = udc[i].f; udc1[j].n = udc[i].n; j++;
724: } else {
725: bzero((char *)tl,(n+1)*sizeof(UM));
726: if ( mod == 2 )
727: berlemain(mod,udc[i].f,tl);
728: else
729: newddd(mod,udc[i].f,tl);
730: for ( k = 0; tl[k]; k++, j++ ) {
731: udc1[j].f = tl[k]; udc1[j].n = udc[i].n;
732: }
733: }
734: udc = udc1; break;
735: }
736: NEWDC(dc0); STOQ(cm,q); COEF(dc0) = (P)q; DEG(dc0) = ONE; dc = dc0;
737: for ( n = 0; udc[n].f; n++ ) {
738: NEWDC(NEXT(dc)); dc = NEXT(dc);
739: STOQ(udc[n].n,DEG(dc)); umtop(VR(p),udc[n].f,&COEF(dc));
740: }
741: NEXT(dc) = 0; *dcp = dc0;
742: }
743:
744: void gensqfrum(mod,p,dc)
745: int mod;
746: UM p;
747: struct oDUM *dc;
748: {
749: int n,i,j,d;
750: UM t,s,g,f,f1,b;
751:
752: if ( (n = DEG(p)) == 1 ) {
753: dc[0].f = UMALLOC(DEG(p)); cpyum(p,dc[0].f); dc[0].n = 1;
754: return;
755: }
756: t = W_UMALLOC(n); s = W_UMALLOC(n); g = W_UMALLOC(n);
757: f = W_UMALLOC(n); f1 = W_UMALLOC(n); b = W_UMALLOC(n);
758: diffum(mod,p,t); cpyum(p,s); Gcdum(mod,t,s,g);
759: if ( !DEG(g) ) {
760: dc[0].f = UMALLOC(DEG(p)); cpyum(p,dc[0].f); dc[0].n = 1;
761: return;
762: }
763: cpyum(p,b); cpyum(p,t); Divum(mod,t,g,f);
764: for ( i = 0, d = 0; DEG(f); i++ ) {
765: while ( 1 ) {
766: cpyum(b,t);
767: if ( Divum(mod,t,f,s) >= 0 )
768: break;
769: else {
770: cpyum(s,b); d++;
771: }
772: }
773: cpyum(b,t); cpyum(f,s); Gcdum(mod,t,s,f1);
774: Divum(mod,f,f1,s); cpyum(f1,f);
775: dc[i].f = UMALLOC(DEG(s)); cpyum(s,dc[i].f); dc[i].n = d;
776: }
777: if ( DEG(b) > 0 ) {
778: d = 1;
779: while ( 1 ) {
780: cpyum(b,t);
781: for ( j = DEG(t); j >= 0; j-- )
782: if ( COEF(t)[j] && (j % mod) )
783: break;
784: if ( j >= 0 )
785: break;
786: else {
787: DEG(s) = DEG(t)/mod;
788: for ( j = 0; j <= DEG(t); j++ )
789: COEF(s)[j] = COEF(t)[j*mod];
790: cpyum(s,b); d *= mod;
791: }
792: }
793: gensqfrum(mod,b,dc+i);
794: for ( j = i; dc[j].f; j++ )
795: dc[j].n *= d;
796: }
797: }
798:
799: #if 0
800: void srchum(mod,p1,p2,gr)
801: int mod;
802: UM p1,p2,gr;
803: {
804: UM m,m1,m2,q,r,t,g1,g2;
805: int lc,d,d1,d2,i,j,k,l,l1,l2,l3,tmp,adj;
806: V v;
807:
808: d = MAX(DEG(p1),DEG(p2));
809: g1 = W_UMALLOC(d); g2 = W_UMALLOC(d);
810: bzero((char *)g1,(d+2)*sizeof(int)); bzero((char *)g2,(d+2)*sizeof(int));
811: if ( d == DEG(p1) ) {
812: cpyum(p1,g1); cpyum(p2,g2);
813: } else {
814: cpyum(p1,g2); cpyum(p2,g1);
815: }
816: if ( ( d1 = DEG(g1) ) > ( d2 = DEG(g2) ) ) {
817: j = d1 - 1; adj = 1;
818: } else
819: j = d2;
820: lc = 1;
821: r = W_UMALLOC(d1+d2); q = W_UMALLOC(d1+d2);
822: m1 = W_UMALLOC(d1+d2); t = W_UMALLOC(d1+d2);
823: bzero((char *)r,(d1+d2+2)*sizeof(int)); bzero((char *)q,(d1+d2+2)*sizeof(int));
824: bzero((char *)m1,(d1+d2+2)*sizeof(int)); bzero((char *)t,(d1+d2+2)*sizeof(int));
825: m = W_UMALLOC(0); bzero((char *)m,2*sizeof(int));
826: adj = pwrm(mod,COEF(g2)[DEG(g2)],DEG(g1));
827: DEG(m) = 0; COEF(m)[0] = invm(COEF(g2)[DEG(g2)],mod);
828: Mulum(mod,g2,m,r); cpyum(r,g2);
829: while ( 1 ) {
830: if ( ( k = DEG(g2) ) < 0 ) {
831: DEG(gr) = -1;
832: return;
833: }
834: if ( k == j ) {
835: if ( k == 0 ) {
836: DEG(m) = 0; COEF(m)[0] = adj;
837: Mulum(mod,g2,m,gr);
838: return;
839: } else {
840: DEG(m) = 0;
841: COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);
842: Mulum(mod,g1,m,r); DEG(r) = Divum(mod,r,g2,t);
843: DEG(m) = 0; COEF(m)[0] = dmb(mod,lc,lc,&tmp);
844: Divum(mod,r,m,q); cpyum(g2,g1); cpyum(q,g2);
845: lc = COEF(g1)[DEG(g1)]; j = k - 1;
846: }
847: } else {
848: d = j - k;
849: DEG(m) = 0; COEF(m)[0] = pwrm(mod,COEF(g2)[DEG(g2)],d);
850: Mulum(mod,g2,m,m1); l = pwrm(mod,lc,d);
851: DEG(m) = 0; COEF(m)[0] = l; Divum(mod,m1,m,t);
852: if ( k == 0 ) {
853: DEG(m) = 0; COEF(m)[0] = adj;
854: Mulum(mod,t,m,gr);
855: return;
856: } else {
857: DEG(m) = 0;
858: COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);
859: Mulum(mod,g1,m,r); DEG(r) = Divum(mod,r,g2,q);
860: l1 = dmb(mod,lc,lc,&tmp); l2 = dmb(mod,l,l1,&tmp);
861: DEG(m) = 0; COEF(m)[0] = l2;
862: Divum(mod,r,m,q); cpyum(t,g1); cpyum(q,g2);
863: if ( d % 2 )
864: for ( i = DEG(g2); i >= 0; i-- )
865: COEF(g2)[i] = ( mod - COEF(g2)[i] ) % mod;
866: lc = COEF(g1)[DEG(g1)]; j = k - 1;
867: }
868: }
869: }
870: }
871:
872: UM *resberle(mod,f,fp)
873: register int mod;
874: UM f;
875: UM *fp;
876: {
877: UM w,wg,ws,wf,f0,gcd,q,res;
878: int n;
879: register int i;
880:
881: n = DEG(f); wg = W_UMALLOC(n); mini(mod,f,wg);
882: if ( DEG(wg) <= 0 ) {
883: f0 = UMALLOC(n); cpyum(f,f0); *fp++ = f0;
884: return ( fp );
885: }
886: f0 = W_UMALLOC(n); cpyum(f,f0);
887: ws = W_UMALLOC(n); wf = W_UMALLOC(n);
888: q = W_UMALLOC(n); gcd = W_UMALLOC(n);
889: res = W_UMALLOC(2*n);
890: srchum(mod,f,wg,res);
891: for ( i = 0; i < mod; i++ ) {
892: if ( substum(mod,res,i) )
893: continue;
894: cpyum(f0,wf); cpyum(wg,ws);
895: COEF(ws)[0] = ( COEF(ws)[0] + mod - i ) % mod;
896: Gcdum(mod,wf,ws,gcd);
897: if ( DEG(gcd) > 0 ) {
898: if ( DEG(gcd) < n ) {
899: Divum(mod,f0,gcd,q); f0 = q; fp = resberle(mod,gcd,fp);
900: }
901: break;
902: }
903: }
904: fp = resberle(mod,f0,fp);
905: return ( fp );
906: }
907:
908: int substum(mod,p,a)
909: int mod;
910: UM p;
911: int a;
912: {
913: int i,j,s;
914: int *c;
915:
916: if ( DEG(p) < 0 )
917: return 0;
918: if ( DEG(p) == 0 )
919: return COEF(p)[0];
920: for ( i = DEG(p), c = COEF(p), s = c[i]; i >= 0; ) {
921: for ( j = i--; (i>=0) && !c[i]; i-- );
922: if ( i >= 0 )
923: s = (s*pwrm(mod,a,j-i)%mod+c[i])%mod;
924: else
925: s = s*pwrm(mod,a,j)%mod;
926: }
927: return s;
928: }
929: #endif
930:
931: void ddd(mod,f,r)
932: int mod;
933: UM f,*r;
934: {
935: register int i,j;
936: int d,n;
937: UM q,s,t,u,v,w,g,x,m;
938: UM *base;
939:
940: n = DEG(f);
941: if ( n == 1 ) {
942: r[0] = UMALLOC(1); cpyum(f,r[0]); r[1] = 0; return;
943: }
944: base = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)base,n*sizeof(UM));
945: w = W_UMALLOC(2*n); q = W_UMALLOC(2*n); m = W_UMALLOC(2*n);
946: base[0] = W_UMALLOC(0); DEG(base[0]) = 0; COEF(base[0])[0] = 1;
947: t = W_UMALLOC(1); DEG(t) = 1; COEF(t)[0] = 0; COEF(t)[1] = 1;
948: pwrmodum(mod,t,mod,f,w); base[1] = W_UMALLOC(DEG(w)); cpyum(w,base[1]);
949: for ( i = 2; i < n; i++ ) {
950: mulum(mod,base[i-1],base[1],m); DEG(m) = divum(mod,m,f,q);
951: base[i] = W_UMALLOC(DEG(m)); cpyum(m,base[i]);
952: }
953: v = W_UMALLOC(n); cpyum(f,v);
954: DEG(w) = 1; COEF(w)[0] = 0; COEF(w)[1] = 1;
955: x = W_UMALLOC(1); DEG(x) = 1; COEF(x)[0] = 0; COEF(x)[1] = 1;
956: t = W_UMALLOC(n); s = W_UMALLOC(n); u = W_UMALLOC(n); g = W_UMALLOC(n);
957: for ( j = 0, d = 1; 2*d <= DEG(v); d++ ) {
958: for ( DEG(t) = -1, i = 0; i <= DEG(w); i++ )
959: if ( COEF(w)[i] ) {
960: Mulsum(mod,base[i],COEF(w)[i],s);
961: addum(mod,s,t,u); cpyum(u,t);
962: }
963: cpyum(t,w); cpyum(v,s); subum(mod,w,x,t); Gcdum(mod,s,t,g);
964: if ( DEG(g) >= 1 ) {
965: canzas(mod,g,d,base,r+j); j += DEG(g)/d;
966: Divum(mod,v,g,q); cpyum(q,v);
967: DEG(w) = Divum(mod,w,v,q);
968: for ( i = 0; i < DEG(v); i++ )
969: DEG(base[i]) = Divum(mod,base[i],v,q);
970: }
971: }
972: if ( DEG(v) ) {
973: r[j] = UMALLOC(DEG(v)); cpyum(v,r[j]); j++;
974: }
975: r[j] = 0;
976: }
977:
978: #if 0
979: void canzas(mod,f,d,base,r)
980: int mod;
981: UM f,*base,*r;
982: {
983: UM t,s,u,w,g,o,q;
984: N n1,n2,n3,n4,n5;
985: UM *b;
986: int n,m,i;
987:
988: if ( DEG(f) == d ) {
989: r[0] = UMALLOC(d); cpyum(f,r[0]);
990: return;
991: } else {
992: n = DEG(f); b = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)b,n*sizeof(UM));
993: for ( i = 0, m = 0; i < n; i++ )
994: m = MAX(m,DEG(base[i]));
995: q = W_UMALLOC(m);
996: for ( i = 0; i < n; i++ ) {
997: b[i] = W_UMALLOC(DEG(base[i])); cpyum(base[i],b[i]);
998: DEG(b[i]) = Divum(mod,b[i],f,q);
999: }
1000: t = W_UMALLOC(2*d);
1001: s = W_UMALLOC(DEG(f)); u = W_UMALLOC(DEG(f));
1002: w = W_UMALLOC(DEG(f)); g = W_UMALLOC(DEG(f));
1003: o = W_UMALLOC(0); DEG(o) = 0; COEF(o)[0] = 1;
1004: STON(mod,n1); pwrn(n1,d,&n2); subn(n2,ONEN,&n3);
1005: STON(2,n4); divsn(n3,n4,&n5);
1006: while ( 1 ) {
1007: randum(mod,2*d,t); spwrum(mod,f,b,t,n5,s);
1008: subum(mod,s,o,u); cpyum(f,w); Gcdum(mod,w,u,g);
1009: if ( (DEG(g) >= 1) && (DEG(g) < DEG(f)) ) {
1010: canzas(mod,g,d,b,r);
1011: cpyum(f,w); Divum(mod,w,g,s);
1012: canzas(mod,s,d,b,r+DEG(g)/d);
1013: return;
1014: }
1015: }
1016: }
1017: }
1018: #else
1019: void canzas(mod,f,d,base,r)
1020: int mod;
1021: UM f,*base,*r;
1022: {
1023: UM t,s,u,w,g,o,q;
1024: N n1,n2,n3,n4,n5;
1025: UM *b;
1026: int n,m,i;
1027:
1028: if ( DEG(f) == d ) {
1029: r[0] = UMALLOC(d); cpyum(f,r[0]);
1030: return;
1031: } else {
1032: n = DEG(f); b = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)b,n*sizeof(UM));
1033: for ( i = 0, m = 0; i < n; i++ )
1034: m = MAX(m,DEG(base[i]));
1035: q = W_UMALLOC(m);
1036: for ( i = 0; i < n; i++ ) {
1037: b[i] = W_UMALLOC(DEG(base[i])); cpyum(base[i],b[i]);
1038: DEG(b[i]) = Divum(mod,b[i],f,q);
1039: }
1040: t = W_UMALLOC(2*d);
1041: s = W_UMALLOC(DEG(f)); u = W_UMALLOC(DEG(f));
1042: w = W_UMALLOC(DEG(f)); g = W_UMALLOC(DEG(f));
1043: o = W_UMALLOC(0); DEG(o) = 0; COEF(o)[0] = 1;
1044: STON(mod,n1); pwrn(n1,d,&n2); subn(n2,ONEN,&n3);
1045: STON(2,n4); divsn(n3,n4,&n5);
1046: while ( 1 ) {
1047: randum(mod,2*d,t); spwrum0(mod,f,t,n5,s);
1048: subum(mod,s,o,u); cpyum(f,w); Gcdum(mod,w,u,g);
1049: if ( (DEG(g) >= 1) && (DEG(g) < DEG(f)) ) {
1050: canzas(mod,g,d,b,r);
1051: cpyum(f,w); Divum(mod,w,g,s);
1052: canzas(mod,s,d,b,r+DEG(g)/d);
1053: return;
1054: }
1055: }
1056: }
1057: }
1058: #endif
1059:
1060: void randum(mod,d,p)
1061: int mod,d;
1062: UM p;
1063: {
1064: unsigned int n;
1065: int i;
1066:
1067: n = ((unsigned int)random()) % d; DEG(p) = n; COEF(p)[n] = 1;
1068: for ( i = 0; i < (int)n; i++ )
1069: COEF(p)[i] = ((unsigned int)random()) % mod;
1070: }
1071:
1072: void pwrmodum(mod,p,e,f,pr)
1073: int mod,e;
1074: UM p,f,pr;
1075: {
1076: UM wt,ws,q;
1077:
1078: if ( e == 0 ) {
1079: DEG(pr) = 0; COEF(pr)[0] = 1;
1080: } else if ( DEG(p) < 0 )
1081: DEG(pr) = -1;
1082: else if ( e == 1 ) {
1083: q = W_UMALLOC(DEG(p)); cpyum(p,pr);
1084: DEG(pr) = divum(mod,pr,f,q);
1085: } else if ( DEG(p) == 0 ) {
1086: DEG(pr) = 0; COEF(pr)[0] = pwrm(mod,COEF(p)[0],e);
1087: } else {
1088: wt = W_UMALLOC(2*DEG(f)); ws = W_UMALLOC(2*DEG(f));
1089: q = W_UMALLOC(2*DEG(f));
1090: pwrmodum(mod,p,e/2,f,wt);
1091: if ( !(e%2) ) {
1092: mulum(mod,wt,wt,pr); DEG(pr) = divum(mod,pr,f,q);
1093: } else {
1094: mulum(mod,wt,wt,ws); DEG(ws) = divum(mod,ws,f,q);
1095: mulum(mod,ws,p,pr); DEG(pr) = divum(mod,pr,f,q);
1096: }
1097: }
1098: }
1099:
1100: void spwrum(mod,m,base,f,e,r)
1101: int mod;
1102: UM f,m,r;
1103: UM *base;
1104: N e;
1105: {
1106: int a,n,i;
1107: N e1,an;
1108: UM t,s,u,q,r1,r2;
1109:
1110: if ( !e ) {
1111: DEG(r) = 0; COEF(r)[0] = 1;
1112: } else if ( UNIN(e) )
1113: cpyum(f,r);
1114: else if ( (PL(e) == 1) && (BD(e)[0] < (unsigned int)mod) )
1115: spwrum0(mod,m,f,e,r);
1116: else {
1117: a = divin(e,mod,&e1); STON(a,an);
1118: n = DEG(m); t = W_UMALLOC(n); s = W_UMALLOC(n);
1119: u = W_UMALLOC(2*n); q = W_UMALLOC(2*n);
1120: for ( DEG(t) = -1, i = 0; i <= DEG(f); i++ )
1121: if ( COEF(f)[i] ) {
1122: Mulsum(mod,base[i],COEF(f)[i],s);
1123: addum(mod,s,t,u); cpyum(u,t);
1124: }
1125: r1 = W_UMALLOC(n); spwrum0(mod,m,f,an,r1);
1126: r2 = W_UMALLOC(n); spwrum(mod,m,base,t,e1,r2);
1127: Mulum(mod,r1,r2,u); DEG(u) = Divum(mod,u,m,q);
1128: cpyum(u,r);
1129: }
1130: }
1131:
1132: void spwrum0(mod,m,f,e,r)
1133: int mod;
1134: UM f,m,r;
1135: N e;
1136: {
1137: UM t,s,q;
1138: N e1;
1139: int a;
1140:
1141: if ( !e ) {
1142: DEG(r) = 0; COEF(r)[0] = 1;
1143: } else if ( UNIN(e) )
1144: cpyum(f,r);
1145: else {
1146: a = divin(e,2,&e1);
1147: t = W_UMALLOC(2*DEG(m)); spwrum0(mod,m,f,e1,t);
1148: s = W_UMALLOC(2*DEG(m)); q = W_UMALLOC(2*DEG(m));
1149: Mulum(mod,t,t,s); DEG(s) = Divum(mod,s,m,q);
1150: if ( a ) {
1151: Mulum(mod,s,f,t); DEG(t) = Divum(mod,t,m,q); cpyum(t,r);
1152: } else
1153: cpyum(s,r);
1154: }
1155: }
1156:
1157: #if 0
1158: void Mulum(mod,p1,p2,pr)
1159: register int mod;
1160: UM p1,p2,pr;
1161: {
1162: register int *pc1,*pcr;
1163: register int mul,i,j,d1,d2;
1164: int *c1,*c2,*cr;
1165:
1166: if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) {
1167: DEG(pr) = -1;
1168: return;
1169: }
1170: c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
1171: bzero((char *)cr,(d1+d2+1)*sizeof(int));
1172: for ( i = 0; i <= d2; i++, cr++ )
1173: if ( mul = *c2++ )
1174: for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ )
1175: *pcr = (*pc1 * mul + *pcr) % mod;
1176: DEG(pr) = d1 + d2;
1177: }
1178:
1179: void Mulsum(mod,p,n,pr)
1180: register int mod,n;
1181: UM p,pr;
1182: {
1183: register int *sp,*dp;
1184: register int i;
1185:
1186: for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i;
1187: i >= 0; i--, dp--, sp-- )
1188: *dp = (*sp * n) % mod;
1189: }
1190:
1191: int Divum(mod,p1,p2,pq)
1192: register int mod;
1193: UM p1,p2,pq;
1194: {
1195: register int *pc1,*pct;
1196: register int tmp,i,j,inv;
1197: int *c1,*c2,*ct;
1198: int d1,d2,dd,hd;
1199:
1200: if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) {
1201: DEG(pq) = -1;
1202: return( d1 );
1203: }
1204: c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2;
1205: if ( ( hd = c2[d2] ) != 1 ) {
1206: inv = invm(hd,mod);
1207: for ( pc1 = c2 + d2; pc1 >= c2; pc1-- )
1208: *pc1 = (*pc1 * inv) % mod;
1209: } else
1210: inv = 1;
1211: for ( i = dd, ct = c1+d1; i >= 0; i-- )
1212: if ( tmp = *ct-- ) {
1213: tmp = mod - tmp;
1214: for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- )
1215: *pct = (*pc1 * tmp + *pct) % mod;
1216: }
1217: if ( inv != 1 ) {
1218: for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ )
1219: *pc1 = (*pc1 * inv) % mod;
1220: for ( pc1 = c2, pct = c2+d2, inv = hd; pc1 <= pct; pc1++ )
1221: *pc1 = (*pc1 * inv) % mod;
1222: }
1223: for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- );
1224: for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- )
1225: *pct-- = *pc1--;
1226: return( i );
1227: }
1228:
1229: void Gcdum(mod,p1,p2,pr)
1230: register int mod;
1231: UM p1,p2,pr;
1232: {
1233: register int *sp,*dp;
1234: register int i,inv;
1235: UM t1,t2,q,tum;
1236: int drem;
1237:
1238: if ( DEG(p1) < 0 )
1239: cpyum(p2,pr);
1240: else if ( DEG(p2) < 0 )
1241: cpyum(p1,pr);
1242: else {
1243: if ( DEG(p1) >= DEG(p2) ) {
1244: t1 = p1; t2 = p2;
1245: } else {
1246: t1 = p2; t2 = p1;
1247: }
1248: q = W_UMALLOC(DEG(t1));
1249: while ( ( drem = Divum(mod,t1,t2,q) ) >= 0 ) {
1250: tum = t1; t1 = t2; t2 = tum; DEG(t2) = drem;
1251: }
1252: inv = invm(COEF(t2)[DEG(t2)],mod);
1253: Mulsum(mod,t2,inv,pr);
1254: }
1255: }
1256: #endif
1257:
1258: void mult_mod_tab(p,mod,tab,r,d)
1259: UM p,r;
1260: UM *tab;
1261: int mod,d;
1262: {
1263: UM w,w1,c;
1264: int n,i;
1265: int *pc;
1266:
1267: w = W_UMALLOC(d); w1 = W_UMALLOC(d);
1268: c = W_UMALLOC(1); DEG(c) = 0;
1269: n = DEG(p); DEG(r) = -1;
1270: for ( i = 0, pc = COEF(p); i <= n; i++ )
1271: if ( pc[i] ) {
1272: COEF(c)[0] = pc[i];
1273: mulum(mod,tab[i],c,w);
1274: addum(mod,r,w,w1);
1275: cpyum(w1,r);
1276: }
1277: }
1278:
1279: void make_qmat(p,mod,tab,mp)
1280: UM p;
1281: int mod;
1282: UM *tab;
1283: int ***mp;
1284: {
1285: int n,i,j;
1286: int *c;
1287: UM q,r;
1288: int **mat;
1289:
1290: n = DEG(p);
1291: *mp = mat = almat(n,n);
1292: for ( j = 0; j < n; j++ ) {
1293: r = W_UMALLOC(DEG(tab[j])); q = W_UMALLOC(DEG(tab[j]));
1294: cpyum(tab[j],r); DEG(r) = divum(mod,r,p,q);
1295: for ( i = 0, c = COEF(r); i <= DEG(r); i++ )
1296: mat[i][j] = c[i];
1297: }
1298: for ( i = 0; i < n; i++ )
1299: mat[i][i] = (mat[i][i]+mod-1) % mod;
1300: }
1301:
1302: void null_mod(mat,mod,n,ind)
1303: int **mat;
1304: int *ind;
1305: int mod,n;
1306: {
1307: int i,j,l,s,h,inv;
1308: int *t,*u;
1309:
1310: bzero((char *)ind,n*sizeof(int));
1311: ind[0] = 0;
1312: for ( i = j = 0; j < n; i++, j++ ) {
1313: for ( ; j < n; j++ ) {
1314: for ( l = i; l < n; l++ )
1315: if ( mat[l][j] )
1316: break;
1317: if ( l < n ) {
1318: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
1319: } else
1320: ind[j] = 1;
1321: }
1322: if ( j == n )
1323: break;
1324: inv = invm(mat[i][j],mod);
1325: for ( s = j, t = mat[i]; s < n; s++ )
1326: t[s] = dmar(t[s],inv,0,mod);
1327: for ( l = 0; l < n; l++ ) {
1328: if ( l == i )
1329: continue;
1330: for ( s = j, u = mat[l], h = (mod-u[j])%mod; s < n; s++ )
1331: u[s] = dmar(h,t[s],u[s],mod);
1332: }
1333: }
1334: }
1335:
1336: void null_to_sol(mat,ind,mod,n,r)
1337: int **mat;
1338: int *ind;
1339: int mod,n;
1340: UM *r;
1341: {
1342: int i,j,k,l;
1343: int *c;
1344: UM w;
1345:
1346: for ( i = 0, l = 0; i < n; i++ ) {
1347: if ( !ind[i] )
1348: continue;
1349: w = UMALLOC(n);
1350: for ( j = k = 0, c = COEF(w); j < n; j++ )
1351: if ( ind[j] )
1352: c[j] = 0;
1353: else
1354: c[j] = mat[k++][i];
1355: c[i] = mod-1;
1356: for ( j = n; j >= 0; j-- )
1357: if ( c[j] )
1358: break;
1359: DEG(w) = j;
1360: r[l++] = w;
1361: }
1362: }
1363: /*
1364: make_qmat(p,mod,tab,mp)
1365: null_mod(mat,mod,n,ind)
1366: null_to_sol(mat,ind,mod,n,r)
1367: */
1368:
1369: void newddd(mod,f,r)
1370: int mod;
1371: UM f,*r;
1372: {
1373: register int i,j;
1374: int d,n;
1375: UM q,s,t,u,v,w,g,x,m;
1376: UM *base;
1377:
1378: n = DEG(f); base = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)base,n*sizeof(UM));
1379: w = W_UMALLOC(2*n); q = W_UMALLOC(2*n); m = W_UMALLOC(2*n);
1380: base[0] = W_UMALLOC(0); DEG(base[0]) = 0; COEF(base[0])[0] = 1;
1381: t = W_UMALLOC(1); DEG(t) = 1; COEF(t)[0] = 0; COEF(t)[1] = 1;
1382: pwrmodum(mod,t,mod,f,w); base[1] = W_UMALLOC(DEG(w)); cpyum(w,base[1]);
1383: for ( i = 2; i < n; i++ ) {
1384: /* fprintf(stderr,"i=%d\n",i); */
1385: mulum(mod,base[i-1],base[1],m); DEG(m) = divum(mod,m,f,q);
1386: base[i] = W_UMALLOC(DEG(m)); cpyum(m,base[i]);
1387: }
1388: v = W_UMALLOC(n); cpyum(f,v);
1389: DEG(w) = 1; COEF(w)[0] = 0; COEF(w)[1] = 1;
1390: x = W_UMALLOC(1); DEG(x) = 1; COEF(x)[0] = 0; COEF(x)[1] = 1;
1391: t = W_UMALLOC(n); s = W_UMALLOC(n); u = W_UMALLOC(n); g = W_UMALLOC(n);
1392: for ( j = 0, d = 1; 2*d <= DEG(v); d++ ) {
1393: /* fprintf(stderr,"d=%d\n",d); */
1394: for ( DEG(t) = -1, i = 0; i <= DEG(w); i++ )
1395: if ( COEF(w)[i] ) {
1396: Mulsum(mod,base[i],COEF(w)[i],s);
1397: addum(mod,s,t,u); cpyum(u,t);
1398: }
1399: cpyum(t,w); cpyum(v,s); subum(mod,w,x,t); Gcdum(mod,s,t,g);
1400: if ( DEG(g) >= 1 ) {
1401: berlekamp(g,mod,d,base,r+j); j += DEG(g)/d;
1402: Divum(mod,v,g,q); cpyum(q,v);
1403: DEG(w) = Divum(mod,w,v,q);
1404: for ( i = 0; i < DEG(v); i++ )
1405: DEG(base[i]) = Divum(mod,base[i],v,q);
1406: }
1407: }
1408: if ( DEG(v) ) {
1409: r[j] = UMALLOC(DEG(v)); cpyum(v,r[j]); j++;
1410: }
1411: r[j] = 0;
1412: }
1413:
1414: int nfctr_mod(f,mod)
1415: int mod;
1416: UM f;
1417: {
1418: register int i,j;
1419: int d,n;
1420: UM q,s,t,u,v,w,g,x,m;
1421: UM *base;
1422:
1423: n = DEG(f); base = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)base,n*sizeof(UM));
1424: w = W_UMALLOC(2*n); q = W_UMALLOC(2*n); m = W_UMALLOC(2*n);
1425: base[0] = W_UMALLOC(0); DEG(base[0]) = 0; COEF(base[0])[0] = 1;
1426: t = W_UMALLOC(1); DEG(t) = 1; COEF(t)[0] = 0; COEF(t)[1] = 1;
1427: pwrmodum(mod,t,mod,f,w); base[1] = W_UMALLOC(DEG(w)); cpyum(w,base[1]);
1428: for ( i = 2; i < n; i++ ) {
1429: /* fprintf(stderr,"i=%d\n",i); */
1430: mulum(mod,base[i-1],base[1],m); DEG(m) = divum(mod,m,f,q);
1431: base[i] = W_UMALLOC(DEG(m)); cpyum(m,base[i]);
1432: }
1433: v = W_UMALLOC(n); cpyum(f,v);
1434: DEG(w) = 1; COEF(w)[0] = 0; COEF(w)[1] = 1;
1435: x = W_UMALLOC(1); DEG(x) = 1; COEF(x)[0] = 0; COEF(x)[1] = 1;
1436: t = W_UMALLOC(n); s = W_UMALLOC(n); u = W_UMALLOC(n); g = W_UMALLOC(n);
1437: for ( j = 0, d = 1; 2*d <= DEG(v); d++ ) {
1438: /* fprintf(stderr,"d=%d\n",d); */
1439: for ( DEG(t) = -1, i = 0; i <= DEG(w); i++ )
1440: if ( COEF(w)[i] ) {
1441: Mulsum(mod,base[i],COEF(w)[i],s);
1442: addum(mod,s,t,u); cpyum(u,t);
1443: }
1444: cpyum(t,w); cpyum(v,s); subum(mod,w,x,t); Gcdum(mod,s,t,g);
1445: if ( DEG(g) >= 1 ) {
1446: j += DEG(g)/d;
1447: Divum(mod,v,g,q); cpyum(q,v);
1448: DEG(w) = Divum(mod,w,v,q);
1449: for ( i = 0; i < DEG(v); i++ )
1450: DEG(base[i]) = Divum(mod,base[i],v,q);
1451: }
1452: }
1453: if ( DEG(v) ) j++;
1454: return j;
1455: }
1456:
1457: int irred_check(f,mod)
1458: UM f;
1459: int mod;
1460: {
1461: register int i,j;
1462: int d,n;
1463: UM q,s,t,u,v,w,g,x,m,f1,b;
1464: UM *base;
1465:
1466: if ( (n = DEG(f)) == 1 )
1467: return 1;
1468: t = W_UMALLOC(n); s = W_UMALLOC(n); g = W_UMALLOC(n);
1469: f1 = W_UMALLOC(n); b = W_UMALLOC(n);
1470: diffum(mod,f,t); cpyum(f,s); Gcdum(mod,t,s,g);
1471: if ( DEG(g) )
1472: return 0;
1473:
1474: base = (UM *)ALLOCA(n*sizeof(UM)); bzero((char *)base,n*sizeof(UM));
1475: w = W_UMALLOC(2*n); q = W_UMALLOC(2*n); m = W_UMALLOC(2*n);
1476: base[0] = W_UMALLOC(0); DEG(base[0]) = 0; COEF(base[0])[0] = 1;
1477: t = W_UMALLOC(1); DEG(t) = 1; COEF(t)[0] = 0; COEF(t)[1] = 1;
1478: pwrmodum(mod,t,mod,f,w); base[1] = W_UMALLOC(DEG(w)); cpyum(w,base[1]);
1479: for ( i = 2; i < n; i++ ) {
1480: /* fprintf(stderr,"i=%d\n",i); */
1481: mulum(mod,base[i-1],base[1],m); DEG(m) = divum(mod,m,f,q);
1482: base[i] = W_UMALLOC(DEG(m)); cpyum(m,base[i]);
1483: }
1484: v = W_UMALLOC(n); cpyum(f,v);
1485: DEG(w) = 1; COEF(w)[0] = 0; COEF(w)[1] = 1;
1486: x = W_UMALLOC(1); DEG(x) = 1; COEF(x)[0] = 0; COEF(x)[1] = 1;
1487: t = W_UMALLOC(n); s = W_UMALLOC(n); u = W_UMALLOC(n); g = W_UMALLOC(n);
1488: for ( j = 0, d = 1; 2*d <= n; d++ ) {
1489: /* fprintf(stderr,"d=%d\n",d); */
1490: for ( DEG(t) = -1, i = 0; i <= DEG(w); i++ )
1491: if ( COEF(w)[i] ) {
1492: Mulsum(mod,base[i],COEF(w)[i],s);
1493: addum(mod,s,t,u); cpyum(u,t);
1494: }
1495: cpyum(t,w); cpyum(v,s); subum(mod,w,x,t); Gcdum(mod,s,t,g);
1496: if ( DEG(g) >= 1 )
1497: return 0;
1498: }
1499: return 1;
1500: }
1501:
1502: int berlekamp(p,mod,df,tab,r)
1503: UM p;
1504: int mod,df;
1505: UM *tab,*r;
1506: {
1507: int n,i,j,k,nf,d,nr;
1508: int **mat;
1509: int *ind;
1510: UM mp,w,q,gcd,w1,w2;
1511: UM *u;
1512: int *root;
1513:
1514: n = DEG(p);
1515: ind = ALLOCA(n*sizeof(int));
1516: make_qmat(p,mod,tab,&mat);
1517: null_mod(mat,mod,n,ind);
1518: for ( i = 0, d = 0; i < n; i++ )
1519: if ( ind[i] )
1520: d++;
1521: if ( d == 1 ) {
1522: r[0] = UMALLOC(n); cpyum(p,r[0]); return 1;
1523: }
1524: u = ALLOCA(d*sizeof(UM *));
1525: r[0] = UMALLOC(n); cpyum(p,r[0]);
1526: null_to_sol(mat,ind,mod,n,u);
1527: root = ALLOCA(d*sizeof(int));
1528: w = W_UMALLOC(n); mp = W_UMALLOC(d);
1529: w1 = W_UMALLOC(n); w2 = W_UMALLOC(n);
1530: for ( i = 1, nf = 1; i < d; i++ ) {
1531: minipoly_mod(mod,u[i],p,mp);
1532: nr = find_root(mod,mp,root);
1533: for ( j = 0; j < nf; j++ ) {
1534: if ( DEG(r[j]) == df )
1535: continue;
1536: for ( k = 0; k < nr; k++ ) {
1537: cpyum(u[i],w1); cpyum(r[j],w2);
1538: COEF(w1)[0] = (mod-root[k]) % mod;
1539: gcdum(mod,w1,w2,w);
1540: if ( DEG(w) > 0 && DEG(w) < DEG(r[j]) ) {
1541: gcd = UMALLOC(DEG(w));
1542: q = UMALLOC(DEG(r[j])-DEG(w));
1543: cpyum(w,gcd); divum(mod,r[j],w,q);
1544: r[j] = q; r[nf++] = gcd;
1545: }
1546: if ( nf == d )
1547: return d;
1548: }
1549: }
1550: }
1551: }
1552:
1553: void minipoly_mod(mod,f,p,mp)
1554: int mod;
1555: UM f,p,mp;
1556: {
1557: struct p_pair *list,*l,*l1,*lprev;
1558: int n,d;
1559: UM u,p0,p1,np0,np1,q,w;
1560:
1561: list = (struct p_pair *)MALLOC(sizeof(struct p_pair));
1562: list->p0 = u = W_UMALLOC(0); DEG(u) = 0; COEF(u)[0] = 1;
1563: list->p1 = W_UMALLOC(0); cpyum(list->p0,list->p1);
1564: list->next = 0;
1565: n = DEG(p); w = UMALLOC(2*n);
1566: p0 = UMALLOC(2*n); cpyum(list->p0,p0);
1567: p1 = UMALLOC(2*n); cpyum(list->p1,p1);
1568: q = W_UMALLOC(2*n);
1569: while ( 1 ) {
1570: COEF(p0)[DEG(p0)] = 0; DEG(p0)++; COEF(p0)[DEG(p0)] = 1;
1571: mulum(mod,f,p1,w); DEG(w) = divum(mod,w,p,q); cpyum(w,p1);
1572: np0 = UMALLOC(n); np1 = UMALLOC(n);
1573: lnf_mod(mod,n,p0,p1,list,np0,np1);
1574: if ( DEG(np1) < 0 ) {
1575: cpyum(np0,mp); return;
1576: } else {
1577: l1 = (struct p_pair *)MALLOC(sizeof(struct p_pair));
1578: l1->p0 = np0; l1->p1 = np1;
1579: for ( l = list, lprev = 0, d = DEG(np1);
1580: l && (DEG(l->p1) > d); lprev = l, l = l->next );
1581: if ( lprev ) {
1582: lprev->next = l1; l1->next = l;
1583: } else {
1584: l1->next = list; list = l1;
1585: }
1586: }
1587: }
1588: }
1589:
1590: void lnf_mod(mod,n,p0,p1,list,np0,np1)
1591: int mod,n;
1592: UM p0,p1;
1593: struct p_pair *list;
1594: UM np0,np1;
1595: {
1596: int inv,h,d1;
1597: UM t0,t1,s0,s1;
1598: struct p_pair *l;
1599:
1600: cpyum(p0,np0); cpyum(p1,np1);
1601: t0 = W_UMALLOC(n); t1 = W_UMALLOC(n);
1602: s0 = W_UMALLOC(n); s1 = W_UMALLOC(n);
1603: for ( l = list; l; l = l->next ) {
1604: d1 = DEG(np1);
1605: if ( d1 == DEG(l->p1) ) {
1606: inv = invm((mod-COEF(l->p1)[d1])%mod,mod);
1607: h = dmar(COEF(np1)[d1],inv,0,mod);
1608: mulsum(mod,l->p0,h,t0); addum(mod,np0,t0,s0); cpyum(s0,np0);
1609: mulsum(mod,l->p1,h,t1); addum(mod,np1,t1,s1); cpyum(s1,np1);
1610: }
1611: }
1612: }
1613:
1614: int find_root(mod,p,root)
1615: int mod;
1616: UM p;
1617: int *root;
1618: {
1619: UM *r;
1620: int i,j;
1621:
1622: r = ALLOCA((DEG(p)+1)*sizeof(UM));
1623: ddd(mod,p,r);
1624: for ( i = 0, j = 0; r[i]; i++ )
1625: if ( DEG(r[i]) == 1 )
1626: root[j++] = (mod - COEF(r[i])[0]) % mod;
1627: return j;
1628: }
1629:
1630: void showum(p)
1631: UM p;
1632: {
1633: int i;
1634: int *c;
1635:
1636: for ( i = DEG(p), c = COEF(p); i >= 0; i-- )
1637: if ( c[i] )
1638: printf("+%dx^%d",c[i],i);
1639: printf("\n");
1640: }
1641:
1642: void showumat(mat,n)
1643: int **mat;
1644: int n;
1645: {
1646: int i,j;
1647:
1648: for ( i = 0; i < n; i++ ) {
1649: for ( j = 0; j < n; j++ )
1650: printf("%d ",mat[i][j]);
1651: printf("\n");
1652: }
1653: }