Annotation of OpenXM_contrib2/asir2000/engine/Hgfs.c, Revision 1.35
1.35 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.34 2018/03/29 01:32:51 noro Exp $ */
1.1 noro 2:
3: #include "ca.h"
1.22 noro 4: #include "inline.h"
1.1 noro 5:
1.30 noro 6: int debug_sfbfctr;
7:
1.18 noro 8: void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1);
1.20 noro 9: void extractcoefbm(BM f,int dx,UM r);
1.35 ! noro 10: void sqfrsf(VL vl, P f, DCP *dcp);
1.1 noro 11:
1.18 noro 12: int comp_dum(DUM a,DUM b)
1.12 noro 13: {
1.34 noro 14: if ( DEG(a->f) > DEG(b->f) )
15: return -1;
16: else if ( DEG(a->f) < DEG(b->f) )
17: return 1;
18: else
19: return 0;
1.12 noro 20: }
21:
1.25 noro 22: void ufctrsf(P p,DCP *dcp)
1.1 noro 23: {
1.34 noro 24: int n,i,j,k;
25: DCP dc,dc0;
26: P lc;
27: UM mp;
28: UM *tl;
29: Obj obj;
30: struct oDUM *udc,*udc1;
31:
32: simp_ff((Obj)p,&obj); p = (P)obj;
33: if ( !p ) {
34: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE;
35: NEXT(dc) = 0; *dcp = dc;
36: return;
37: }
38: mp = W_UMALLOC(UDEG(p));
39: ptosfum(p,mp);
40: if ( (n = DEG(mp)) < 0 ) {
41: *dcp = 0; return;
42: } else if ( n == 0 ) {
43: NEWDC(dc); COEF(dc) = p; DEG(dc) = ONE;
44: NEXT(dc) = 0; *dcp = dc;
45: return;
46: }
47: lc = COEF(DC(p));
48: if ( !_isonesf(COEF(mp)[n]) ) {
49: monicsfum(mp);
50: }
51:
52: W_CALLOC(n+1,struct oDUM,udc);
53: gensqfrsfum(mp,udc);
54:
55: tl = (UM *)ALLOCA((n+1)*sizeof(UM));
56: W_CALLOC(DEG(mp)+1,struct oDUM,udc1);
57:
58: for ( i = 0,j = 0; udc[i].f; i++ )
59: if ( DEG(udc[i].f) == 1 ) {
60: udc1[j].f = udc[i].f; udc1[j].n = udc[i].n; j++;
61: } else {
62: bzero((char *)tl,(n+1)*sizeof(UM));
63: czsfum(udc[i].f,tl);
64: for ( k = 0; tl[k]; k++, j++ ) {
65: udc1[j].f = tl[k]; udc1[j].n = udc[i].n;
66: }
67: }
68: udc = udc1;
69: for ( i = 0; udc[i].f; i++ );
70: qsort(udc,i,sizeof(struct oDUM),
71: (int (*)(const void *,const void *))comp_dum);
72:
73: NEWDC(dc0); COEF(dc0) = lc; DEG(dc0) = ONE; dc = dc0;
74: for ( n = 0; udc[n].f; n++ ) {
75: NEWDC(NEXT(dc)); dc = NEXT(dc);
76: STOQ(udc[n].n,DEG(dc)); sfumtop(VR(p),udc[n].f,&COEF(dc));
77: }
78: NEXT(dc) = 0; *dcp = dc0;
1.1 noro 79: }
80:
1.18 noro 81: void gensqfrsfum(UM p,struct oDUM *dc)
1.1 noro 82: {
1.34 noro 83: int n,i,j,d,mod;
84: UM t,s,g,f,f1,b;
85: GFS u,v;
86:
87: if ( (n = DEG(p)) == 1 ) {
88: dc[0].f = UMALLOC(DEG(p)); cpyum(p,dc[0].f); dc[0].n = 1;
89: return;
90: }
91: t = W_UMALLOC(n); s = W_UMALLOC(n); g = W_UMALLOC(n);
92: f = W_UMALLOC(n); f1 = W_UMALLOC(n); b = W_UMALLOC(n);
93: diffsfum(p,t); cpyum(p,s); gcdsfum(t,s,g);
94: if ( !DEG(g) ) {
95: dc[0].f = UMALLOC(DEG(p)); cpyum(p,dc[0].f); dc[0].n = 1;
96: return;
97: }
98: cpyum(p,b); cpyum(p,t); divsfum(t,g,f);
99: for ( i = 0, d = 0; DEG(f); i++ ) {
100: while ( 1 ) {
101: cpyum(b,t);
102: if ( divsfum(t,f,s) >= 0 )
103: break;
104: else {
105: cpyum(s,b); d++;
106: }
107: }
108: cpyum(b,t); cpyum(f,s); gcdsfum(t,s,f1);
109: divsfum(f,f1,s); cpyum(f1,f);
110: dc[i].f = UMALLOC(DEG(s)); cpyum(s,dc[i].f); dc[i].n = d;
111: }
112: mod = characteristic_sf();
113: if ( DEG(b) > 0 ) {
114: d = 1;
115: while ( 1 ) {
116: cpyum(b,t);
117: for ( j = DEG(t); j >= 0; j-- )
118: if ( COEF(t)[j] && (j % mod) )
119: break;
120: if ( j >= 0 )
121: break;
122: else {
123: DEG(s) = DEG(t)/mod;
124: for ( j = 0; j <= DEG(t); j++ ) {
125: iftogfs(COEF(t)[j*mod],&u);
126: pthrootgfs(u,&v);
127: COEF(s)[j] = v?FTOIF(CONT(v)):0;
128: }
129: cpyum(s,b); d *= mod;
130: }
131: }
132: gensqfrsfum(b,dc+i);
133: for ( j = i; dc[j].f; j++ )
134: dc[j].n *= d;
135: }
1.1 noro 136: }
137:
1.18 noro 138: void randsfum(int d,UM p)
1.1 noro 139: {
1.34 noro 140: int i;
1.1 noro 141:
1.34 noro 142: for ( i = 0; i < d; i++ )
143: COEF(p)[i] = _randomsf();
144: for ( i = d-1; i >= 0 && !COEF(p)[i]; i-- );
145: p->d = i;
1.1 noro 146: }
147:
1.18 noro 148: void pwrmodsfum(UM p,int e,UM f,UM pr)
1.1 noro 149: {
1.34 noro 150: UM wt,ws,q;
1.1 noro 151:
1.34 noro 152: if ( e == 0 ) {
153: DEG(pr) = 0; COEF(pr)[0] = _onesf();
154: } else if ( DEG(p) < 0 )
155: DEG(pr) = -1;
156: else if ( e == 1 ) {
157: q = W_UMALLOC(DEG(p)); cpyum(p,pr);
158: DEG(pr) = divsfum(pr,f,q);
159: } else if ( DEG(p) == 0 ) {
160: DEG(pr) = 0; COEF(pr)[0] = _pwrsf(COEF(p)[0],e);
161: } else {
162: wt = W_UMALLOC(2*DEG(f)); ws = W_UMALLOC(2*DEG(f));
163: q = W_UMALLOC(2*DEG(f));
164: pwrmodsfum(p,e/2,f,wt);
165: if ( !(e%2) ) {
166: mulsfum(wt,wt,pr); DEG(pr) = divsfum(pr,f,q);
167: } else {
168: mulsfum(wt,wt,ws);
169: DEG(ws) = divsfum(ws,f,q);
170: mulsfum(ws,p,pr);
171: DEG(pr) = divsfum(pr,f,q);
172: }
173: }
1.1 noro 174: }
175:
1.18 noro 176: void spwrsfum(UM m,UM f,N e,UM r)
1.1 noro 177: {
1.34 noro 178: UM t,s,q;
179: N e1;
180: int a;
181:
182: if ( !e ) {
183: DEG(r) = 0; COEF(r)[0] = _onesf();
184: } else if ( UNIN(e) )
185: cpyum(f,r);
186: else {
187: a = divin(e,2,&e1);
188: t = W_UMALLOC(2*DEG(m)); spwrsfum(m,f,e1,t);
189: s = W_UMALLOC(2*DEG(m)); q = W_UMALLOC(2*DEG(m));
190: mulsfum(t,t,s); DEG(s) = divsfum(s,m,q);
191: if ( a ) {
192: mulsfum(s,f,t); DEG(t) = divsfum(t,m,q); cpyum(t,r);
1.1 noro 193: } else
1.34 noro 194: cpyum(s,r);
195: }
1.1 noro 196: }
197:
1.18 noro 198: void tracemodsfum(UM m,UM f,int e,UM r)
1.2 noro 199: {
1.34 noro 200: UM t,s,q,u;
201: int i;
1.2 noro 202:
1.34 noro 203: q = W_UMALLOC(2*DEG(m)+DEG(f)); /* XXX */
204: t = W_UMALLOC(2*DEG(m));
205: s = W_UMALLOC(2*DEG(m));
206: u = W_UMALLOC(2*DEG(m));
207: DEG(f) = divsfum(f,m,q);
208: cpyum(f,s);
209: cpyum(f,t);
210: for ( i = 1; i < e; i++ ) {
211: mulsfum(t,t,u);
212: DEG(u) = divsfum(u,m,q); cpyum(u,t);
213: addsfum(t,s,u); cpyum(u,s);
214: }
215: cpyum(s,r);
1.2 noro 216: }
217:
1.18 noro 218: void make_qmatsf(UM p,UM *tab,int ***mp)
1.1 noro 219: {
1.34 noro 220: int n,i,j;
221: int *c;
222: UM q,r;
223: int **mat;
224: int one;
225:
226: n = DEG(p);
227: *mp = mat = almat(n,n);
228: for ( j = 0; j < n; j++ ) {
229: r = W_UMALLOC(DEG(tab[j])); q = W_UMALLOC(DEG(tab[j]));
230: cpyum(tab[j],r); DEG(r) = divsfum(r,p,q);
231: for ( i = 0, c = COEF(r); i <= DEG(r); i++ )
232: mat[i][j] = c[i];
233: }
234: one = _onesf();
235: for ( i = 0; i < n; i++ )
236: mat[i][i] = _subsf(mat[i][i],one);
1.1 noro 237: }
238:
1.18 noro 239: void nullsf(int **mat,int n,int *ind)
1.1 noro 240: {
1.34 noro 241: int i,j,l,s,h,inv;
242: int *t,*u;
1.1 noro 243:
1.34 noro 244: bzero((char *)ind,n*sizeof(int));
245: ind[0] = 0;
246: for ( i = j = 0; j < n; i++, j++ ) {
247: for ( ; j < n; j++ ) {
248: for ( l = i; l < n; l++ )
249: if ( mat[l][j] )
250: break;
251: if ( l < n ) {
252: t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
253: } else
254: ind[j] = 1;
255: }
256: if ( j == n )
257: break;
258: inv = _invsf(mat[i][j]);
259: for ( s = j, t = mat[i]; s < n; s++ )
260: t[s] = _mulsf(t[s],inv);
261: for ( l = 0; l < n; l++ ) {
262: if ( l == i )
263: continue;
264: u = mat[l]; h = _chsgnsf(u[j]);
265: for ( s = j; s < n; s++ )
266: u[s] = _addsf(_mulsf(h,t[s]),u[s]);
267: }
268: }
1.1 noro 269: }
270:
1.18 noro 271: void null_to_solsf(int **mat,int *ind,int n,UM *r)
1.1 noro 272: {
1.34 noro 273: int i,j,k,l;
274: int *c;
275: UM w;
276:
277: for ( i = 0, l = 0; i < n; i++ ) {
278: if ( !ind[i] )
279: continue;
280: w = UMALLOC(n);
281: for ( j = k = 0, c = COEF(w); j < n; j++ )
282: if ( ind[j] )
283: c[j] = 0;
284: else
285: c[j] = mat[k++][i];
286: c[i] = _chsgnsf(_onesf());
287: for ( j = n; j >= 0; j-- )
288: if ( c[j] )
289: break;
290: DEG(w) = j;
291: r[l++] = w;
292: }
1.1 noro 293: }
294: /*
295: make_qmatsf(p,tab,mp)
296: nullsf(mat,n,ind)
297: null_to_solsf(ind,n,r)
298: */
299:
1.18 noro 300: void czsfum(UM f,UM *r)
1.1 noro 301: {
1.34 noro 302: int i,j;
303: int d,n,ord;
304: UM s,t,u,v,w,g,x,m,q;
305: UM *base;
306:
307: n = DEG(f); base = (UM *)ALLOCA(n*sizeof(UM));
308: bzero((char *)base,n*sizeof(UM));
309:
310: w = W_UMALLOC(2*n); q = W_UMALLOC(2*n); m = W_UMALLOC(2*n);
311:
312: base[0] = W_UMALLOC(0); DEG(base[0]) = 0; COEF(base[0])[0] = _onesf();
313:
314: t = W_UMALLOC(1); DEG(t) = 1; COEF(t)[0] = 0; COEF(t)[1] = _onesf();
315:
316: ord = field_order_sf();
317: pwrmodsfum(t,ord,f,w);
318: base[1] = W_UMALLOC(DEG(w));
319: cpyum(w,base[1]);
320:
321: for ( i = 2; i < n; i++ ) {
322: mulsfum(base[i-1],base[1],m);
323: DEG(m) = divsfum(m,f,q);
324: base[i] = W_UMALLOC(DEG(m)); cpyum(m,base[i]);
325: }
326:
327: v = W_UMALLOC(n); cpyum(f,v);
328: DEG(w) = 1; COEF(w)[0] = 0; COEF(w)[1] = _onesf();
329: x = W_UMALLOC(1); DEG(x) = 1; COEF(x)[0] = 0; COEF(x)[1] = _onesf();
330: t = W_UMALLOC(n); s = W_UMALLOC(n); u = W_UMALLOC(n); g = W_UMALLOC(n);
331:
332: for ( j = 0, d = 1; 2*d <= DEG(v); d++ ) {
333: for ( DEG(t) = -1, i = 0; i <= DEG(w); i++ )
334: if ( COEF(w)[i] ) {
335: mulssfum(base[i],COEF(w)[i],s);
336: addsfum(s,t,u); cpyum(u,t);
337: }
338: cpyum(t,w); cpyum(v,s); subsfum(w,x,t);
339: gcdsfum(s,t,g);
340: if ( DEG(g) >= 1 ) {
341: berlekampsf(g,d,base,r+j); j += DEG(g)/d;
342: divsfum(v,g,q); cpyum(q,v);
343: DEG(w) = divsfum(w,v,q);
344: for ( i = 0; i < DEG(v); i++ )
345: DEG(base[i]) = divsfum(base[i],v,q);
346: }
347: }
348: if ( DEG(v) ) {
349: r[j] = UMALLOC(DEG(v)); cpyum(v,r[j]); j++;
350: }
351: r[j] = 0;
1.1 noro 352: }
353:
1.18 noro 354: int berlekampsf(UM p,int df,UM *tab,UM *r)
1.1 noro 355: {
1.34 noro 356: int n,i,j,k,nf,d,nr;
357: int **mat;
358: int *ind;
359: UM mp,w,q,gcd,w1,w2;
360: UM *u;
361: int *root;
362:
363: n = DEG(p);
364: ind = ALLOCA(n*sizeof(int));
365: make_qmatsf(p,tab,&mat);
366: nullsf(mat,n,ind);
367: for ( i = 0, d = 0; i < n; i++ )
368: if ( ind[i] )
369: d++;
370: if ( d == 1 ) {
371: r[0] = UMALLOC(n); cpyum(p,r[0]); return 1;
372: }
373: u = ALLOCA(d*sizeof(UM *));
374: r[0] = UMALLOC(n); cpyum(p,r[0]);
375: null_to_solsf(mat,ind,n,u);
376: root = ALLOCA(d*sizeof(int));
377: w = W_UMALLOC(n); mp = W_UMALLOC(d);
378: w1 = W_UMALLOC(n); w2 = W_UMALLOC(n);
379: for ( i = 1, nf = 1; i < d; i++ ) {
380: minipolysf(u[i],p,mp);
381: nr = find_rootsf(mp,root);
382: for ( j = 0; j < nf; j++ ) {
383: if ( DEG(r[j]) == df )
384: continue;
385: for ( k = 0; k < nr; k++ ) {
386: cpyum(u[i],w1); cpyum(r[j],w2);
387: COEF(w1)[0] = _chsgnsf(root[k]);
388: gcdsfum(w1,w2,w);
389: if ( DEG(w) > 0 && DEG(w) < DEG(r[j]) ) {
390: gcd = UMALLOC(DEG(w));
391: q = UMALLOC(DEG(r[j])-DEG(w));
392: cpyum(w,gcd); divsfum(r[j],w,q);
393: r[j] = q; r[nf++] = gcd;
394: }
395: if ( nf == d )
396: return d;
397: }
398: }
399: }
400: /* NOT REACHED */
401: error("berlekampsf : cannot happen");
402: return 0;
1.1 noro 403: }
404:
1.18 noro 405: void minipolysf(UM f,UM p,UM mp)
1.1 noro 406: {
1.34 noro 407: struct p_pair *list,*l,*l1,*lprev;
408: int n,d;
409: UM u,p0,p1,np0,np1,q,w;
410:
411: list = (struct p_pair *)MALLOC(sizeof(struct p_pair));
412: list->p0 = u = W_UMALLOC(0); DEG(u) = 0; COEF(u)[0] = _onesf();
413: list->p1 = W_UMALLOC(0); cpyum(list->p0,list->p1);
414: list->next = 0;
415: n = DEG(p); w = UMALLOC(2*n);
416: p0 = UMALLOC(2*n); cpyum(list->p0,p0);
417: p1 = UMALLOC(2*n); cpyum(list->p1,p1);
418: q = W_UMALLOC(2*n);
419: while ( 1 ) {
420: COEF(p0)[DEG(p0)] = 0; DEG(p0)++; COEF(p0)[DEG(p0)] = _onesf();
421: mulsfum(f,p1,w); DEG(w) = divsfum(w,p,q); cpyum(w,p1);
422: np0 = UMALLOC(n); np1 = UMALLOC(n);
423: lnfsf(n,p0,p1,list,np0,np1);
424: if ( DEG(np1) < 0 ) {
425: cpyum(np0,mp); return;
426: } else {
427: l1 = (struct p_pair *)MALLOC(sizeof(struct p_pair));
428: l1->p0 = np0; l1->p1 = np1;
429: for ( l = list, lprev = 0, d = DEG(np1);
430: l && (DEG(l->p1) > d); lprev = l, l = l->next );
431: if ( lprev ) {
432: lprev->next = l1; l1->next = l;
433: } else {
434: l1->next = list; list = l1;
435: }
436: }
437: }
1.1 noro 438: }
439:
1.18 noro 440: void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1)
1.1 noro 441: {
1.34 noro 442: int h,d1;
443: UM t0,t1,s0,s1;
444: struct p_pair *l;
445:
446: cpyum(p0,np0); cpyum(p1,np1);
447: t0 = W_UMALLOC(n); t1 = W_UMALLOC(n);
448: s0 = W_UMALLOC(n); s1 = W_UMALLOC(n);
449: for ( l = list; l; l = l->next ) {
450: d1 = DEG(np1);
451: if ( d1 == DEG(l->p1) ) {
452: h = _divsf(COEF(np1)[d1],_chsgnsf(COEF(l->p1)[d1]));
453: mulssfum(l->p0,h,t0); addsfum(np0,t0,s0); cpyum(s0,np0);
454: mulssfum(l->p1,h,t1); addsfum(np1,t1,s1); cpyum(s1,np1);
455: }
456: }
1.1 noro 457: }
458:
1.18 noro 459: int find_rootsf(UM p,int *root)
1.1 noro 460: {
1.34 noro 461: UM *r;
462: int i,n;
1.1 noro 463:
1.34 noro 464: n = DEG(p);
465: r = ALLOCA((DEG(p))*sizeof(UM));
466: canzassf(p,1,r);
467: for ( i = 0; i < n; i++ )
468: root[i] = _chsgnsf(COEF(r[i])[0]);
469: return n;
1.1 noro 470: }
471:
1.18 noro 472: void canzassf(UM f,int d,UM *r)
1.1 noro 473: {
1.34 noro 474: UM t,s,u,w,g,o;
475: N n1,n2,n3,n4,n5;
476: UM *b;
477: int n,q,ed;
478:
479: if ( DEG(f) == d ) {
480: r[0] = UMALLOC(d); cpyum(f,r[0]);
481: return;
482: } else {
483: n = DEG(f); b = (UM *)ALLOCA(n*sizeof(UM));
484: bzero((char *)b,n*sizeof(UM));
485:
486: t = W_UMALLOC(2*d);
487: s = W_UMALLOC(DEG(f)); u = W_UMALLOC(DEG(f));
488: w = W_UMALLOC(DEG(f)); g = W_UMALLOC(DEG(f));
489: o = W_UMALLOC(0); DEG(o) = 0; COEF(o)[0] = _onesf();
490: q = field_order_sf();
491: if ( q % 2 ) {
492: STON(q,n1); pwrn(n1,d,&n2); subn(n2,ONEN,&n3);
493: STON(2,n4); divsn(n3,n4,&n5);
494: } else
495: ed = d*extdeg_sf();
496: while ( 1 ) {
497: randsfum(2*d,t);
498: if ( q % 2 ) {
499: spwrsfum(f,t,n5,s); subsfum(s,o,u);
500: } else
501: tracemodsfum(f,t,ed,u);
502: cpyum(f,w);
503: gcdsfum(w,u,g);
504: if ( (DEG(g) >= 1) && (DEG(g) < DEG(f)) ) {
505: canzassf(g,d,r);
506: cpyum(f,w); divsfum(w,g,s);
507: canzassf(s,d,r+DEG(g)/d);
508: return;
509: }
510: }
511: }
1.1 noro 512: }
513:
1.3 noro 514: /* Hensel related functions */
515:
1.28 noro 516: int sfberle(V,V,P,int,GFS *,DCP *);
1.3 noro 517: void sfgcdgen(P,ML,ML *);
1.5 noro 518: void sfhenmain2(BM,UM,UM,int,BM *);
519: void ptosfbm(int,P,BM);
1.28 noro 520: void sfhensel(int count,P f,V x,V y,int degbound,GFS *evp,P *sfp,ML *listp);
1.3 noro 521:
522: /* f = f(x,y) */
523:
1.28 noro 524: void sfhensel(int count,P f,V x,V y,int degbound,GFS *evp,P *sfp,ML *listp)
1.3 noro 525: {
1.34 noro 526: int i;
527: int fn;
528: ML rlist;
529: BM fl;
530: VL vl,nvl;
531: int dx,dy,bound;
532: GFS ev;
533: P f1,t,c,sf;
534: DCP dc,dct,dc0;
535: UM q,fm,hm;
536: UM *gm;
537: struct oEGT tmp0,tmp1,eg_hensel,eg_hensel_t;
538:
539: clctv(CO,f,&vl);
540: if ( vl->v != x ) {
541: reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&f1);
542: vl = nvl; f = f1;
543: }
544: if ( vl->next )
545: y = vl->next->v;
546: dx = getdeg(x,f);
547: dy = getdeg(y,f);
548: if ( dx == 1 ) {
549: *listp = rlist = MLALLOC(1); rlist->n = 1; rlist->c[0] = 0;
550: return;
551: }
552: fn = sfberle(x,y,f,count,&ev,&dc);
553: if ( fn <= 1 ) {
554: /* fn == 0 => short of evaluation points */
555: *listp = rlist = MLALLOC(1); rlist->n = fn; rlist->c[0] = 0;
556: return;
557: }
558: if ( degbound >= 0 ) {
559: /*
560: * reconstruct dc so that
561: * dc[1],... : factors satisfy degree bound
562: * dc[0] : product of others
563: */
564: c = dc->c; dc = NEXT(dc);
565: dc0 = 0;
566: fn = 0;
567: while ( dc ) {
568: if ( getdeg(x,COEF(dc)) <= degbound ) {
569: dct = NEXT(dc); NEXT(dc) = dc0; dc0 = dc; dc = dct;
570: fn++;
571: } else {
572: mulp(vl,COEF(dc),c,&t); c = t;
573: dc = NEXT(dc);
574: }
575: }
576: if ( OID(c) == O_P ) {
577: NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = dc0;
578: fn++;
579: } else {
580: mulp(vl,dc0->c,c,&t); dc0->c = t; dc = dc0;
581: }
582: } else {
583: /* pass the the leading coeff. to the first element */
584: c = dc->c; dc = NEXT(dc);
585: mulp(vl,dc->c,c,&t); dc->c = t;
586: }
587:
588: /* convert mod y-a factors into UM */
589: gm = (UM *)ALLOCA(fn*sizeof(UM));
590: for ( i = 0; i < fn; i++, dc = NEXT(dc) ) {
591: gm[i] = W_UMALLOC(UDEG(dc->c));
592: ptosfum(dc->c,gm[i]);
593: }
594:
595: /* set bound */
596: /* g | f, lc_y(g) = lc_y(f) => deg_y(g) <= deg_y(f) */
597: /* so, bound = dy is sufficient, but we use slightly large value */
598: bound = dy+2;
599:
600: /* f(x,y) -> f(x,y+ev) */
601: fl = BMALLOC(dx,bound);
602: ptosfbm(bound,f,fl);
603: if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev)));
604:
605: /* sf = f(x+ev) */
606: sfbmtop(fl,x,y,&sf);
607:
608: /* fm = fl mod y */
609: fm = W_UMALLOC(dx);
610: cpyum(COEF(fl)[0],fm);
611: hm = W_UMALLOC(dx);
612:
613: q = W_UMALLOC(dx);
614: rlist = MLALLOC(fn); rlist->n = fn; rlist->bound = bound;
615: if ( debug_sfbfctr )
616: fprintf(asir_out,"%d candidates\n",fn);
617: init_eg(&eg_hensel);
618: for ( i = 0; i < fn-1; i++ ) {
619: if ( debug_sfbfctr )
620: fprintf(asir_out,"deg(fm) = %d, deg(gm[%d]) = %d\n",
621: DEG(fm),i,DEG(gm[i]));
622: init_eg(&eg_hensel_t);
623: get_eg(&tmp0);
624: /* fl = gm[i]*hm mod y */
625: divsfum(fm,gm[i],hm);
626: /* fl is replaced by the cofactor of gk mod y^bound */
627: /* rlist->c[i] = gk */
628: sfhenmain2(fl,gm[i],hm,bound,(BM *)&rlist->c[i]);
629: cpyum(hm,fm);
630: get_eg(&tmp1); add_eg(&eg_hensel_t,&tmp0,&tmp1);
631: add_eg(&eg_hensel,&tmp0,&tmp1);
632: if ( debug_sfbfctr) {
633: print_eg("Hensel",&eg_hensel_t);
634: fprintf(asir_out,"\n");
635: }
636: }
637: if ( debug_sfbfctr) {
638: print_eg("Hensel total",&eg_hensel);
639: fprintf(asir_out,"\n");
640: }
641: /* finally, fl must be the lift of gm[fn-1] */
642: rlist->c[i] = fl;
1.4 noro 643:
1.8 noro 644: #if 0
1.34 noro 645: /* y -> y-a */
646: mev = _chsgnsf(FTOIF(CONT(ev)));
647: for ( i = 0; i < fn; i++ )
648: shiftsfbm((BM)(rlist->c[i]),mev);
1.8 noro 649: #endif
1.34 noro 650: *evp = ev;
651: *sfp = sf;
652: *listp = rlist;
1.3 noro 653: }
654:
655: /* main variable of f = x */
656:
1.28 noro 657: int sfberle(V x,V y,P f,int count,GFS *ev,DCP *dcp)
1.3 noro 658: {
1.34 noro 659: UM wf,wf1,wf2,wfs,gcd;
660: int fn,n;
661: GFS m,fm;
662: DCP dc,dct,dc0;
663: VL vl;
664: P lc,lc0,f0;
665: Obj obj;
666: int j,q,index,i;
667:
668: NEWVL(vl); vl->v = x;
669: NEWVL(NEXT(vl)); NEXT(vl)->v = y;
670: NEXT(NEXT(vl)) =0;
671: simp_ff((Obj)f,&obj); f = (P)obj;
672: n = QTOS(DEG(DC(f)));
673: wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); wf2 = W_UMALLOC(n);
674: wfs = W_UMALLOC(n); gcd = W_UMALLOC(n);
675: q = field_order_sf();
676: lc = DC(f)->c;
677: for ( j = 0, fn = n + 1, index = 0;
678: index < q && j < count && fn > 1; index++ ) {
679: indextogfs(index,&m);
680: substp(vl,lc,y,(P)m,&lc0);
681: if ( lc0 ) {
682: substp(vl,f,y,(P)m,&f0);
683: ptosfum(f0,wf); cpyum(wf,wf1);
684: diffsfum(wf1,wf2); gcdsfum(wf1,wf2,gcd);
685: if ( DEG(gcd) == 0 ) {
686: ufctrsf(f0,&dc);
687: for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ );
688: if ( i < fn ) {
689: dc0 = dc; fn = i; fm = m;
690: }
691: j++;
692: }
693: }
694: }
695: if ( index == q )
696: return 0;
697: else if ( fn == 1 )
698: return 1;
699: else {
700: *dcp = dc0;
701: *ev = fm;
702: return fn;
703: }
1.3 noro 704: }
705:
1.18 noro 706: void sfgcdgen(P f,ML blist,ML *clistp)
1.3 noro 707: {
1.34 noro 708: int i;
709: int n,d,np;
710: UM wf,wm,wx,wy,wu,wv,wa,wb,wg,q,tum;
711: UM *in,*out;
712: ML clist;
713:
714: n = UDEG(f); np = blist->n;
715: d = 2*n;
716: q = W_UMALLOC(d); wf = W_UMALLOC(d);
717: wm = W_UMALLOC(d); wx = W_UMALLOC(d);
718: wy = W_UMALLOC(d); wu = W_UMALLOC(d);
719: wv = W_UMALLOC(d); wg = W_UMALLOC(d);
720: wa = W_UMALLOC(d); wb = W_UMALLOC(d);
721: ptosfum(f,wf); DEG(wg) = 0; COEF(wg)[0] = _onesf();
722: *clistp = clist = MLALLOC(np); clist->n = np;
723: for ( i = 0, in = (UM *)blist->c, out = (UM *)clist->c; i < np; i++ ) {
724: divsfum(wf,in[i],q); tum = wf; wf = q; q = tum;
725: cpyum(wf,wx); cpyum(in[i],wy);
726: eucsfum(wx,wy,wa,wb); mulsfum(wa,wg,wm);
727: DEG(wm) = divsfum(wm,in[i],q); out[i] = UMALLOC(DEG(wm));
728: cpyum(wm,out[i]); mulsfum(q,wf,wu);
729: mulsfum(wg,wb,wv); addsfum(wu,wv,wg);
730: }
1.3 noro 731: }
732:
1.14 noro 733: /* f = g0*h0 mod y -> f = gk*hk mod y^(dy+1), f is replaced by hk */
1.3 noro 734:
1.18 noro 735: void sfhenmain2(BM f,UM g0,UM h0,int dy,BM *gp)
1.4 noro 736: {
1.34 noro 737: int i,k;
738: int dx;
739: UM wt,wa,wb,q,w1,w2,wh1,wg1,ws;
740: UM wc,wd,we,wz;
741: BM wb0,wb1;
742: int dg,dh;
743: BM fk,gk,hk;
744:
745: if ( DEG(f) < dy )
746: error("sfhenmain2 : invalid input");
747:
748: dx = degbm(f);
749: dg = DEG(g0);
750: dh = DEG(h0);
751:
752: W_BMALLOC(dx,dy,wb0); W_BMALLOC(dx,dy,wb1);
753: wt = W_UMALLOC(dx); ws = W_UMALLOC(dx); q = W_UMALLOC(2*dx);
754: wg1 = W_UMALLOC(2*dx); wh1 = W_UMALLOC(2*dx);
755:
756: /* fk = gk*hk mod y^k */
757: W_BMALLOC(dx,dy,fk);
758: cpyum(COEF(f)[0],COEF(fk)[0]);
759: gk = BMALLOC(dg,dy);
760: cpyum(g0,COEF(gk)[0]);
761: W_BMALLOC(dh,dy,hk);
762: cpyum(h0,COEF(hk)[0]);
763:
764: wc = W_UMALLOC(2*dx); wd = W_UMALLOC(2*dx);
765: we = W_UMALLOC(2*dx); wz = W_UMALLOC(2*dx);
766:
767: /* compute wa,wb s.t. wa*g0+wb*h0 = 1 mod y */
768: w1 = W_UMALLOC(dg); cpyum(g0,w1);
769: w2 = W_UMALLOC(dh); cpyum(h0,w2);
770: wa = W_UMALLOC(2*dx); wb = W_UMALLOC(2*dx); /* XXX */
771: eucsfum(w1,w2,wa,wb);
772:
773: if ( debug_sfbfctr)
774: fprintf(stderr,"dy=%d\n",dy);
775: for ( k = 1; k <= dy; k++ ) {
776: if ( debug_sfbfctr)
777: fprintf(stderr,".");
778:
779: /* at this point, f = gk*hk mod y^k */
780:
781: /* clear wt */
782: clearum(wt,dx);
783:
784: /* wt = (f-gk*hk)/y^k */
785: subsfum(COEF(f)[k],COEF(fk)[k],wt);
786:
787: /* compute wf1,wg1 s.t. wh1*g0+wg1*h0 = wt */
788: mulsfum(wa,wt,wh1); DEG(wh1) = divsfum(wh1,h0,q);
789: mulsfum(wh1,g0,wc); subsfum(wt,wc,wd); DEG(wd) = divsfum(wd,h0,wg1);
1.4 noro 790:
1.34 noro 791: /* check */
1.4 noro 792: #if 0
1.34 noro 793: if ( DEG(wd) >= 0 || DEG(wg1) > ng )
794: error("henmain2 : cannot happen(adj)");
1.4 noro 795:
1.34 noro 796: mulsfum(wg1,h0,wc); mulsfum(wh1,g0,wd); addsfum(wc,wd,we);
797: subsfum(we,wt,wz);
798: if ( DEG(wz) >= 0 )
799: error("henmain2 : cannot happen");
1.4 noro 800: #endif
801:
1.34 noro 802: /* fk += ((wg1*hk+wh1*gk)*y^k+wg1*wh1*y^(2*k) mod y^(dy+1) */
803: /* wb0 = wh1*y^k */
804: clearbm(dx,wb0);
805: cpyum(wh1,COEF(wb0)[k]);
806:
807: /* wb1 = gk*wb0 mod y^(dy+1) */
808: clearbm(dx,wb1);
809: mulsfbm(gk,wb0,wb1);
810: /* fk += wb1 */
811: addtosfbm(wb1,fk);
812:
813: /* wb0 = wg1*y^k */
814: clearbm(dx,wb0);
815: cpyum(wg1,COEF(wb0)[k]);
816:
817: /* wb1 = hk*wb0 mod y^(dy+1) */
818: clearbm(dx,wb1);
819: mulsfbm(hk,wb0,wb1);
820: /* fk += wb1 */
821: addtosfbm(wb1,fk);
822:
823: /* fk += wg1*wh1*y^(2*k) mod y^(dy+1) */
824: if ( 2*k <= dy ) {
825: mulsfum(wg1,wh1,wt); addsfum(COEF(fk)[2*k],wt,ws);
826: cpyum(ws,COEF(fk)[2*k]);
827: }
828:
829: /* gk += wg1*y^k, hk += wh1*y^k */
830: cpyum(wg1,COEF(gk)[k]);
831: cpyum(wh1,COEF(hk)[k]);
832: }
833: if ( debug_sfbfctr)
834: fprintf(stderr,"\n");
835: *gp = gk;
836: DEG(f) = dy;
837: for ( i = 0; i <= dy; i++ )
838: cpyum(COEF(hk)[i],COEF(f)[i]);
1.27 noro 839: }
840:
841: /* a0*g+b0*h = 1 mod y -> a*g+b*h = 1 mod y^(dy+1) */
842:
843: void sfexgcd_by_hensel(BM g,BM h,int dy,BM *ap,BM *bp)
844: {
1.34 noro 845: int i,k;
846: int dx;
847: UM wt,wa,wb,q,w1,w2,ws;
848: UM wc,wd,we,wz,wa1,wb1;
849: BM wz0,wz1;
850: int dg,dh;
851: BM a,b,c;
852:
853: dg = degbm(g);
854: dh = degbm(h);
855: dx = dg+dh;
856:
857: a = BMALLOC(dh,dy);
858: b = BMALLOC(dg,dy);
859: /* c holds a*g+b*h-1 */
860: c = BMALLOC(dg+dh,dy);
861:
862: W_BMALLOC(dx,dy,wz0); W_BMALLOC(dx,dy,wz1);
863:
864: wt = W_UMALLOC(dx); ws = W_UMALLOC(dx); q = W_UMALLOC(2*dx);
865: wa1 = W_UMALLOC(2*dx); wb1 = W_UMALLOC(2*dx);
866: wc = W_UMALLOC(2*dx); wd = W_UMALLOC(2*dx);
867: we = W_UMALLOC(2*dx); wz = W_UMALLOC(2*dx);
868:
869: /* compute wa,wb s.t. wa*g0+wb*h0 = 1 mod y */
870: w1 = W_UMALLOC(dg); cpyum(COEF(g)[0],w1);
871: w2 = W_UMALLOC(dh); cpyum(COEF(h)[0],w2);
872: wa = W_UMALLOC(2*dx); wb = W_UMALLOC(2*dx); /* XXX */
873: eucsfum(w1,w2,wa,wb);
874: cpyum(wa,COEF(a)[0]); cpyum(wb,COEF(b)[0]);
875:
876: /* initialize c to a*g+b*h-1 */
877: mulsfbm(a,g,c); mulsfbm(b,h,wz0); addtosfbm(wz0,c);
878: COEF(COEF(c)[0])[0] = 0;
879:
880: if ( debug_sfbfctr)
881: fprintf(stderr,"dy=%d\n",dy);
882: for ( k = 1; k <= dy; k++ ) {
883: if ( debug_sfbfctr)
884: fprintf(stderr,".");
885:
886: /* at this point, a*g+b*h = 1 mod y^k, c = a*g+b*h-1 */
887:
888: /* wt = -((a*g+b*h-1)/y^k) */
889: cpyum(COEF(c)[k],wt);
890: for ( i = DEG(wt); i >= 0; i-- )
891: COEF(wt)[i] = _chsgnsf(COEF(wt)[i]);
892:
893: /* compute wa1,wb1 s.t. wa1*g0+wb1*h0 = wt */
894: mulsfum(wa,wt,wa1); DEG(wa1) = divsfum(wa1,COEF(h)[0],q);
895: mulsfum(wa1,COEF(g)[0],wc); subsfum(wt,wc,wd);
896: DEG(wd) = divsfum(wd,COEF(h)[0],wb1);
897:
898: /* c += ((wa1*g+wb1*h)*y^k mod y^(dy+1) */
899: /* wz0 = wa1*y^k */
900: clearbm(dx,wz0);
901: cpyum(wa1,COEF(wz0)[k]);
902:
903: /* wz1 = wz0*g mod y^(dy+1) */
904: clearbm(dx,wz1);
905: mulsfbm(g,wz0,wz1);
906: /* c += wz1 */
907: addtosfbm(wz1,c);
908:
909: /* wz0 = wb1*y^k */
910: clearbm(dx,wz0);
911: cpyum(wb1,COEF(wz0)[k]);
912:
913: /* wz1 = wz0*h mod y^(dy+1) */
914: clearbm(dx,wz1);
915: mulsfbm(h,wz0,wz1);
916: /* c += wz1 */
917: addtosfbm(wz1,c);
918:
919: /* a += wa1*y^k, b += wb1*y^k */
920: cpyum(wa1,COEF(a)[k]);
921: cpyum(wb1,COEF(b)[k]);
922: }
923: if ( debug_sfbfctr)
924: fprintf(stderr,"\n");
925: DEG(a) = dy;
926: DEG(b) = dy;
927: *ap = a;
928: *bp = b;
1.3 noro 929: }
930:
1.5 noro 931: /* fl->c[i] = coef_y(f,i) */
932:
1.18 noro 933: void ptosfbm(int dy,P f,BM fl)
1.5 noro 934: {
1.34 noro 935: DCP dc;
936: int d,i,dx;
937: UM t;
938:
939: dx = QTOS(DEG(DC(f)));
940: if ( DEG(fl) < dy )
941: error("ptosfbm : invalid input");
942: DEG(fl) = dy;
943: clearbm(dx,fl);
944: t = UMALLOC(dy);
945: for ( dc = DC(f); dc; dc = NEXT(dc) ) {
946: d = QTOS(DEG(dc));
947: ptosfum(COEF(dc),t);
948: for ( i = 0; i <= DEG(t); i++ )
949: COEF(COEF(fl)[i])[d] = COEF(t)[i];
950: }
951: for ( i = 0; i <= dy; i++ )
952: degum(COEF(fl)[i],dx);
1.5 noro 953: }
954:
955: /* x : main variable */
956:
1.18 noro 957: void sfbmtop(BM f,V x,V y,P *fp)
1.5 noro 958: {
1.34 noro 959: UM *c;
960: int i,j,d,a,dy;
961: GFS b;
962: DCP dc0,dc,dct;
963:
964: dy = DEG(f);
965: c = COEF(f);
966: d = degbm(f);
967:
968: dc0 = 0;
969: for ( i = 0; i <= d; i++ ) {
970: dc = 0;
971: for ( j = 0; j <= dy; j++ ) {
972: if ( DEG(c[j]) >= i && (a = COEF(c[j])[i]) ) {
973: NEWDC(dct);
974: STOQ(j,DEG(dct));
975: iftogfs(a,&b);
976: COEF(dct) = (P)b;
977: NEXT(dct) = dc;
978: dc = dct;
979: }
980: }
981: if ( dc ) {
982: NEWDC(dct);
983: STOQ(i,DEG(dct));
984: MKP(y,dc,COEF(dct));
985: NEXT(dct) = dc0;
986: dc0 = dct;
987: }
988: }
989: if ( dc0 )
990: MKP(x,dc0,*fp);
991: else
992: *fp = 0;
1.8 noro 993: }
994:
1.18 noro 995: void sfsqfr(P f,DCP *dcp)
1.14 noro 996: {
1.34 noro 997: Obj obj;
998: DCP dc;
999: VL vl;
1000:
1001: simp_ff((Obj)f,&obj); f = (P)obj;
1002: clctv(CO,f,&vl);
1003: if ( !vl ) {
1004: /* f is a const */
1005: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0; *dcp = dc;
1006: } else if ( !NEXT(vl) )
1007: sfusqfr(f,dcp);
1008: else
1009: sqfrsf(vl,f,dcp);
1.14 noro 1010: }
1011:
1.18 noro 1012: void sfusqfr(P f,DCP *dcp)
1.14 noro 1013: {
1.34 noro 1014: DCP dc,dct;
1015: struct oDUM *udc;
1016: V x;
1017: P lc;
1018: int n,i;
1019: UM mf;
1020:
1021: x = VR(f);
1022: n = getdeg(x,f);
1023: mf = W_UMALLOC(n);
1024: ptosfum(f,mf);
1025: lc = COEF(DC(f));
1026: if ( !_isonesf(COEF(mf)[n]) ) {
1027: monicsfum(mf);
1028: }
1029: W_CALLOC(n+1,struct oDUM,udc);
1030: gensqfrsfum(mf,udc);
1031: for ( i = 0, dc = 0; udc[i].f; i++ ) {
1032: NEWDC(dct); STOQ(udc[i].n,DEG(dct));
1033: sfumtop(x,udc[i].f,&COEF(dct));
1034: NEXT(dct) = dc; dc = dct;
1035: }
1036: NEWDC(dct); DEG(dct) = ONE; COEF(dct) = (P)lc; NEXT(dct) = dc;
1037: *dcp = dct;
1.14 noro 1038: }
1039:
1.24 noro 1040: #if 0
1.21 noro 1041: void sfbsqfrmain(P f,V x,V y,DCP *dcp)
1042: {
1.34 noro 1043: /* XXX*/
1.21 noro 1044: }
1045:
1046: /* f is bivariate */
1047:
1.18 noro 1048: void sfbsqfr(P f,V x,V y,DCP *dcp)
1.14 noro 1049: {
1.34 noro 1050: P t,rf,cx,cy;
1051: VL vl,rvl;
1052: DCP dcx,dcy,dct,dc;
1053: struct oVL vl0,vl1;
1054:
1055: /* cy(y) = cont(f,x), f /= cy */
1056: cont_pp_sfp(vl,f,&cy,&t); f = t;
1057: /* rvl = [y,x] */
1058: reordvar(vl,y,&rvl); reorderp(rvl,vl,f,&rf);
1059: /* cx(x) = cont(rf,y), Rf /= cy */
1060: cont_pp_sfp(rvl,rf,&cx,&t); rf = t;
1061: reorderp(vl,rvl,rf,&f);
1062:
1063: /* f -> cx*cy*f */
1064: sfsqfr(cx,&dcx); dcx = NEXT(dcx);
1065: sfsqfr(cy,&dcy); dcy = NEXT(dcy);
1066: if ( dcx ) {
1067: for ( dct = dcx; NEXT(dct); dct = NEXT(dct) );
1068: NEXT(dct) = dcy;
1069: } else
1070: dcx = dcy;
1071: if ( OID(f) == O_N )
1072: *dcp = dcx;
1073: else {
1074: /* f must be bivariate */
1075: sfbsqfrmain(f,x,y,&dc);
1076: if ( dcx ) {
1077: for ( dct = dcx; NEXT(dct); dct = NEXT(dct) );
1078: NEXT(dct) = dc;
1079: } else
1080: dcx = dc;
1081: *dcp = dcx;
1082: }
1.14 noro 1083: }
1.24 noro 1084: #endif
1.14 noro 1085:
1.8 noro 1086: void sfdtest(P,ML,V,V,DCP *);
1087:
1.21 noro 1088: /* if degbound >= 0 find factor s.t. deg_x(factor) <= degbound */
1089:
1090: void sfbfctr(P f,V x,V y,int degbound,DCP *dcp)
1.8 noro 1091: {
1.34 noro 1092: ML list;
1093: P sf;
1094: GFS ev;
1095: DCP dc,dct;
1096: BM fl;
1097: int dx,dy;
1098:
1099: /* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */
1100: sfhensel(5,f,x,y,degbound,&ev,&sf,&list);
1101: if ( list->n == 0 )
1102: error("sfbfctr : short of evaluation points");
1103: else if ( list->n == 1 ) {
1104: /* f is irreducible */
1105: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0;
1106: *dcp = dc;
1107: return;
1108: }
1109: sfdtest(sf,list,x,y,&dc);
1110: if ( ev ) {
1111: dx = getdeg(x,sf);
1112: dy = getdeg(y,sf);
1113: W_BMALLOC(dx,dy,fl);
1114: for ( dct = dc; dct; dct = NEXT(dct) ) {
1115: ptosfbm(dy,COEF(dct),fl);
1116: shiftsfbm(fl,_chsgnsf(FTOIF(CONT(ev))));
1117: sfbmtop(fl,x,y,&COEF(dct));
1118: }
1119: }
1120: *dcp = dc;
1.26 noro 1121: }
1122:
1123: /* returns shifted f, shifted factors and the eval pt */
1124:
1125: void sfbfctr_shift(P f,V x,V y,int degbound,GFS *evp,P *sfp,DCP *dcp)
1126: {
1.34 noro 1127: ML list;
1128: P sf;
1129: GFS ev;
1130: DCP dc,dct;
1131: int dx,dy;
1132:
1133: /* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */
1134: sfhensel(5,f,x,y,degbound,&ev,&sf,&list);
1135: if ( list->n == 0 )
1136: error("sfbfctr_shift : short of evaluation points");
1137: else if ( list->n == 1 ) {
1138: /* f is irreducible */
1139: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0;
1140: *evp = 0;
1141: *sfp = f;
1142: *dcp = dc;
1143: } else {
1144: sfdtest(sf,list,x,y,dcp);
1145: *evp = ev;
1146: *sfp = sf;
1147: }
1.8 noro 1148: }
1149:
1.14 noro 1150: /* f = f(x,y) = list->c[0]*list->c[1]*... mod y^(list->bound+1) */
1.8 noro 1151:
1.18 noro 1152: void sfdtest(P f,ML list,V x,V y,DCP *dcp)
1.8 noro 1153: {
1.34 noro 1154: int np,dx,dy;
1155: int i,j,k,bound;
1156: int *win;
1157: P g,lcg,factor,cofactor,lcyx;
1158: P csum;
1159: DCP dcf,dcf0,dc;
1160: BM *c;
1161: BM lcy;
1162: UM lcg0,lcy0,w;
1163: UM *d1c;
1164: ML wlist;
1165: struct oVL vl1,vl0;
1166: VL vl;
1167: int z,dt,dtok;
1168:
1169: /* vl = [x,y] */
1170: vl0.v = x; vl0.next = &vl1; vl1.v = y; vl1.next = 0; vl = &vl0;
1171:
1172: /* setup various structures and arrays */
1173: dx = getdeg(x,f);
1174: dy = getdeg(y,f);
1175: np = list->n;
1176: win = W_ALLOC(np+1);
1177: wlist = W_MLALLOC(np);
1178: wlist->n = list->n;
1179: bound = wlist->bound = list->bound;
1180: c = (BM *)COEF(wlist);
1181: bcopy((char *)COEF(list),(char *)c,(int)(sizeof(BM)*np));
1182:
1183: lcg0 = W_UMALLOC(2*dy);
1184:
1185: /* initialize g by f */
1186: g = f;
1187:
1188: /* initialize lcg */
1189: mulp(vl,g,COEF(DC(g)),&lcg);
1190:
1191: /* initialize lcg0 */
1192: const_term(lcg,lcg0);
1193:
1194: /* initialize csum = lcg(1) */
1195: sfcsump(vl,lcg,&csum);
1196:
1197: /* initialize lcy by LC(f) */
1198: W_BMALLOC(0,dy,lcy);
1199: NEWDC(dc); COEF(dc) = COEF(DC(g)); DEG(dc) = 0;
1200: NEWP(lcyx); VR(lcyx) = x; DC(lcyx) = dc;
1201: ptosfbm(dy,lcyx,lcy);
1202:
1203: /* initialize lcy0 by LC(f) */
1204: lcy0 = W_UMALLOC(bound);
1205: ptosfum(COEF(DC(g)),lcy0);
1206:
1207: /* ((d-1 coefs)*lcy0 */
1208: d1c = (UM *)W_ALLOC(np*sizeof(UM));
1209: w = W_UMALLOC(2*bound);
1210: for ( i = 1; i < np; i++ ) {
1211: extractcoefbm(c[i],degbm(c[i])-1,w);
1212: d1c[i] = W_UMALLOC(2*bound);
1213: mulsfum(w,lcy0,d1c[i]);
1214: /* d1c[i] = d1c[i] mod y^(bound+1) */
1215: if ( DEG(d1c[i]) > bound ) {
1216: for ( j = DEG(d1c[i]); j > bound; j-- )
1217: COEF(d1c[i])[j] = 0;
1218: degum(d1c[i],bound);
1219: }
1220: }
1221:
1222: if ( debug_sfbfctr)
1223: fprintf(stderr,"np = %d\n",np);
1224: dtok = 0;
1225: for ( g = f, k = 1, dcf = dcf0 = 0, win[0] = 1, --np, z = 0; ; z++ ) {
1226: if ( debug_sfbfctr && !(z % 1000) ) fprintf(stderr,".");
1227: dt = sfdegtest(dy,bound,d1c,k,win);
1228: if ( dt )
1229: dtok++;
1230: if ( dt && sfdtestmain(vl,lcg,lcg0,lcy,csum,wlist,
1231: k,win,&factor,&cofactor) ) {
1232: NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
1233: g = cofactor;
1234:
1235: /* update lcg */
1236: mulp(vl,g,COEF(DC(g)),&lcg);
1237:
1238: /* update lcg0 */
1239: const_term(lcg,lcg0);
1240:
1241: /* update csum */
1242: sfcsump(vl,lcg,&csum);
1243:
1244: /* update dy */
1245: dy = getdeg(y,g);
1246:
1247: /* update lcy */
1248: clearbm(0,lcy);
1249: COEF(dc) = COEF(DC(g));
1250: ptosfbm(dy,lcyx,lcy);
1251:
1252: for ( i = 0; i < k - 1; i++ )
1253: for ( j = win[i] + 1; j < win[i + 1]; j++ )
1254: c[j-i-1] = c[j];
1255: for ( j = win[k-1] + 1; j <= np; j++ )
1256: c[j-k] = c[j];
1257: if ( ( np -= k ) < k )
1258: break;
1259: if ( np - win[0] + 1 < k )
1260: if ( ++k > np )
1261: break;
1262: else
1263: for ( i = 0; i < k; i++ )
1264: win[i] = i + 1;
1265: else
1266: for ( i = 1; i < k; i++ )
1267: win[i] = win[0] + i;
1268:
1269:
1270: /* update lcy0 */
1271: ptosfum(COEF(DC(g)),lcy0);
1272:
1273: /* update d-1 coeffs */
1274: for ( i = 1; i <= np; i++ ) {
1275: extractcoefbm(c[i],degbm(c[i])-1,w);
1276: mulsfum(w,lcy0,d1c[i]);
1277: /* d1c[i] = d1c[1] mod y^(bound+1) */
1278: if ( DEG(d1c[i]) > bound ) {
1279: for ( j = DEG(d1c[i]); j > bound; j-- )
1280: COEF(d1c[i])[j] = 0;
1281: degum(d1c[i],bound);
1282: }
1283: }
1284: } else if ( !ncombi(1,np,k,win) )
1285: if ( k == np )
1286: break;
1287: else
1288: for ( i = 0, ++k; i < k; i++ )
1289: win[i] = i + 1;
1290: }
1291: if ( debug_sfbfctr )
1292: fprintf(stderr,"total %d, omitted by degtest %d\n",z,z-dtok);
1293: NEXTDC(dcf0,dcf); COEF(dcf) = g;
1294: DEG(dcf) = ONE; NEXT(dcf) = 0; *dcp = dcf0;
1.8 noro 1295: }
1296:
1.20 noro 1297: void extractcoefbm(BM f,int dx,UM r)
1.19 noro 1298: {
1.34 noro 1299: int j;
1300: UM fj;
1.19 noro 1301:
1.34 noro 1302: for ( j = DEG(f); j >= 0; j-- ) {
1303: fj = COEF(f)[j];
1304: if ( fj && DEG(fj) >= dx ) {
1305: COEF(r)[j] = COEF(fj)[dx];
1306: } else
1307: COEF(r)[j] = 0;
1308: }
1309: degum(r,DEG(f));
1.19 noro 1310: }
1311:
1312: /* deg_y(prod mod y^(bound+1)) <= dy ? */
1313:
1.20 noro 1314: int sfdegtest(int dy,int bound,UM *d1c,int k,int *in)
1.19 noro 1315: {
1.34 noro 1316: int i,j;
1317: UM w,w1,wt;
1318: BM t;
1319:
1320: w = W_UMALLOC(bound);
1321: w1 = W_UMALLOC(bound);
1322: clearum(w,bound);
1323: for ( i = 0; i < k; i++ ) {
1324: addsfum(w,d1c[in[i]],w1); wt = w; w = w1; w1 = wt;
1325: }
1326: return DEG(w) <= dy ? 1 : 0;
1.19 noro 1327: }
1328:
1.9 noro 1329: /* lcy = LC(g), lcg = lcy*g, lcg0 = const part of lcg */
1.19 noro 1330:
1.18 noro 1331: int sfdtestmain(VL vl,P lcg,UM lcg0,BM lcy,P csum,ML list,
1.34 noro 1332: int k,int *in,P *fp,P *cofp)
1.8 noro 1333: {
1.34 noro 1334: P fmul,csumg,q,cont;
1335: V x,y;
1.8 noro 1336:
1.34 noro 1337: x = vl->v;
1338: y = vl->next->v;
1339: if (!sfctest(lcg0,lcy,list,k,in))
1340: return 0;
1341: mulsfbmarray(UDEG(lcg),lcy,list,k,in,x,y,&fmul);
1342: if ( csum ) {
1343: sfcsump(vl,fmul,&csumg);
1344: if ( csumg ) {
1345: if ( !divtp(vl,csum,csumg,&q) )
1346: return 0;
1347: }
1348: }
1349: if ( divtp_by_sfbm(vl,lcg,fmul,&q) ) {
1350: cont_pp_sfp(vl,fmul,&cont,fp);
1351: cont_pp_sfp(vl,q,&cont,cofp);
1352: return 1;
1353: } else
1354: return 0;
1.8 noro 1355: }
1356:
1.18 noro 1357: void const_term(P f,UM c)
1.9 noro 1358: {
1.34 noro 1359: DCP dc;
1.9 noro 1360:
1.34 noro 1361: for ( dc = DC(f); dc && DEG(dc); dc = NEXT(dc) );
1362: if ( dc )
1363: ptosfum(COEF(dc),c);
1364: else
1365: DEG(c) = -1;
1.9 noro 1366: }
1367:
1.18 noro 1368: void const_term_sfbm(BM f,UM c)
1.9 noro 1369: {
1.34 noro 1370: int i,dy;
1.9 noro 1371:
1.34 noro 1372: dy = DEG(f);
1373: for ( i = 0; i <= dy; i++ )
1374: if ( DEG(COEF(f)[i]) >= 0 )
1375: COEF(c)[i] = COEF(COEF(f)[i])[0];
1376: else
1377: COEF(c)[i] = 0;
1378: degum(c,dy);
1.9 noro 1379: }
1380:
1381: /* lcy*(product of const part) | lcg0 ? */
1382:
1.18 noro 1383: int sfctest(UM lcg0,BM lcy,ML list,int k,int *in)
1.8 noro 1384: {
1.34 noro 1385: int dy,i,dr;
1386: UM t,s,u,w;
1387: BM *l;
1388:
1389: dy = list->bound;
1390: t = W_UMALLOC(2*dy);
1391: s = W_UMALLOC(2*dy);
1392: u = W_UMALLOC(2*dy);
1393: const_term_sfbm(lcy,t);
1394: if ( DEG(t) < 0 )
1395: return 1;
1396:
1397: l = (BM *)list->c;
1398: for ( i = 0; i < k; i++ ) {
1399: const_term_sfbm(l[in[i]],s);
1400: mulsfum(t,s,u);
1401: if ( DEG(u) > dy )
1402: degum(u,dy);
1403: w = t; t = u; u = w;
1404: }
1405: cpyum(lcg0,s);
1406: dr = divsfum(s,t,u);
1407: if ( dr >= 0 )
1408: return 0;
1409: else
1410: return 1;
1.8 noro 1411: }
1412:
1413: /* main var of f is x */
1414:
1.18 noro 1415: void mulsfbmarray(int dx,BM lcy,ML list,int k,int *in,V x,V y,P *g)
1.8 noro 1416: {
1.34 noro 1417: int dy,i;
1418: BM wb0,wb1,t;
1419: BM *l;
1420:
1421: dy = list->bound;
1422: W_BMALLOC(dx,dy,wb0); W_BMALLOC(dx,dy,wb1);
1423: l = (BM *)list->c;
1424: clearbm(dx,wb0);
1425: mulsfbm(lcy,l[in[0]],wb0);
1426: for ( i = 1; i < k; i++ ) {
1427: clearbm(dx,wb1);
1428: mulsfbm(l[in[i]],wb0,wb1);
1429: t = wb0; wb0 = wb1; wb1 = t;
1430: }
1431: sfbmtop(wb0,x,y,g);
1.8 noro 1432: }
1433:
1.18 noro 1434: void sfcsump(VL vl,P f,P *s)
1.8 noro 1435: {
1.34 noro 1436: P t,u;
1437: DCP dc;
1.8 noro 1438:
1.34 noro 1439: for ( dc = DC(f), t = 0; dc; dc = NEXT(dc) ) {
1440: addp(vl,COEF(dc),t,&u); t = u;
1441: }
1442: *s = t;
1.8 noro 1443: }
1444:
1445: /* *fp = primitive part of f w.r.t. x */
1446:
1.18 noro 1447: void cont_pp_sfp(VL vl,P f,P *cp,P *fp)
1.8 noro 1448: {
1.34 noro 1449: V x,y;
1450: int d;
1451: UM t,s,gcd;
1452: DCP dc;
1453: GFS g;
1454:
1455: x = vl->v;
1456: y = vl->next->v;
1457: d = getdeg(y,f);
1458: if ( d == 0 ) {
1459: itogfs(1,&g);
1460: *cp = (P)g;
1461: *fp = f; /* XXX */
1462: } else {
1463: t = W_UMALLOC(2*d);
1464: s = W_UMALLOC(2*d);
1465: gcd = W_UMALLOC(2*d);
1466: dc = DC(f);
1467: ptosfum(COEF(dc),gcd);
1468: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
1469: ptosfum(COEF(dc),t);
1470: gcdsfum(gcd,t,s);
1471: cpyum(s,gcd);
1472: }
1473: sfumtop(y,gcd,cp);
1474: divsp(vl,f,*cp,fp);
1475: }
1.11 noro 1476: }
1477:
1.18 noro 1478: int divtp_by_sfbm(VL vl,P f,P g,P *qp)
1.11 noro 1479: {
1.34 noro 1480: V x,y;
1481: int fx,fy,gx,gy;
1482: BM fl,gl,ql;
1483: UM *cf,*cg,*cq;
1484: UM hg,q,t,s;
1485: int i,j,dr;
1486:
1487: x = vl->v; y = vl->next->v;
1488: fx = getdeg(x,f); fy = getdeg(y,f);
1489: gx = getdeg(x,g); gy = getdeg(y,g);
1490:
1491: if ( fx < gx || fy < gy )
1492: return 0;
1493: W_BMALLOC(fx,fy,fl); ptosfbm(fy,f,fl); cf = COEF(fl);
1494: W_BMALLOC(gx,gy,gl); ptosfbm(gy,g,gl); cg = COEF(gl);
1495: W_BMALLOC(fx-gx,fy-gy,ql); cq = COEF(ql);
1496:
1497: hg = cg[gy];
1498: q = W_UMALLOC(fx); t = W_UMALLOC(fx); s = W_UMALLOC(fx);
1499:
1500: for ( i = fy; i >= gy; i-- ) {
1501: if ( DEG(cf[i]) < 0 )
1502: continue;
1503: dr = divsfum(cf[i],hg,q);
1504: if ( dr >= 0 )
1505: return 0;
1506: if ( DEG(q) > fx-gx )
1507: return 0;
1508: cpyum(q,cq[i-gy]);
1509: for ( j = 0; j <= gy; j++ ) {
1510: mulsfum(cg[j],q,t);
1511: subsfum(cf[j+i-gy],t,s);
1512: cpyum(s,cf[j+i-gy]);
1513: }
1514: }
1515: for ( j = gy-1; j >= 0 && DEG(cf[j]) < 0; j-- );
1516: if ( j >= 0 )
1517: return 0;
1518: sfbmtop(ql,x,y,qp);
1519: return 1;
1.16 noro 1520: }
1521:
1522: /* XXX generate an irreducible poly of degree n */
1523:
1.17 noro 1524: extern int current_gfs_q1;
1.22 noro 1525: extern int *current_gfs_ntoi;
1.17 noro 1526:
1.18 noro 1527: void generate_defpoly_sfum(int n,UM *dp)
1.16 noro 1528: {
1.34 noro 1529: UM r,dr,t,g;
1530: UM *f;
1531: int *c,*w;
1532: int max,i,j;
1533:
1534: *dp = r = UMALLOC(n);
1535: DEG(r) = n;
1536: c = COEF(r);
1537: c[n] = _onesf();
1538: max = current_gfs_q1;
1539: w = (int *)ALLOCA(n*sizeof(int));
1540: bzero(w,n*sizeof(int));
1541:
1542: dr = W_UMALLOC(n); t = W_UMALLOC(n); g = W_UMALLOC(n);
1543: f = (UM *)ALLOCA((n+1)*sizeof(UM));
1544: while ( 1 ) {
1545: for ( i = 0; i < n && w[i] == max; i++ );
1546: if ( i == n ) {
1547: /* XXX cannot happen */
1548: error("generate_defpoly_sfum : cannot happen");
1549: }
1550: for ( j = 0; j < i; j++ )
1551: w[j] = 0;
1552: w[i]++;
1553: if ( !current_gfs_ntoi )
1554: for ( i = 0; i < n; i++ )
1555: c[i] = w[i]?FTOIF(w[i]):0;
1556: else
1557: for ( i = 0; i < n; i++ )
1558: c[i] = w[i]?FTOIF(w[i]-1):0;
1559: if ( !c[0] )
1560: continue;
1561: diffsfum(r,dr); cpyum(r,t); gcdsfum(t,dr,g);
1562: if ( DEG(g) > 0 )
1563: continue;
1564:
1565: czsfum(r,f);
1566: for ( i = 0; f[i]; i++ );
1567: if ( i == 1 )
1568: return;
1569: }
1.3 noro 1570: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Hgfs.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine