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