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