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