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