Annotation of OpenXM_contrib2/asir2018/engine/Hgfs.c, Revision 1.2
1.2 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2018/engine/Hgfs.c,v 1.1 2018/09/19 05:45:07 noro Exp $ */
1.1 noro 2:
3: #include "ca.h"
4: #include "inline.h"
5:
6: int debug_sfbfctr;
7:
8: void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1);
9: void extractcoefbm(BM f,int dx,UM r);
10:
11: int comp_dum(DUM a,DUM b)
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:
21: void ufctrsf(P p,DCP *dcp)
22: {
23: int n,i,j,k;
24: DCP dc,dc0;
25: P lc;
26: UM mp;
27: UM *tl;
28: Obj obj;
29: struct oDUM *udc,*udc1;
30:
31: simp_ff((Obj)p,&obj); p = (P)obj;
32: if ( !p ) {
33: NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE;
34: NEXT(dc) = 0; *dcp = dc;
35: return;
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;
68: for ( i = 0; udc[i].f; i++ );
69: qsort(udc,i,sizeof(struct oDUM),
70: (int (*)(const void *,const void *))comp_dum);
71:
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);
1.2 ! noro 75: STOZ(udc[n].n,DEG(dc)); sfumtop(VR(p),udc[n].f,&COEF(dc));
1.1 noro 76: }
77: NEXT(dc) = 0; *dcp = dc0;
78: }
79:
80: void gensqfrsfum(UM p,struct oDUM *dc)
81: {
82: int n,i,j,d,mod;
83: UM t,s,g,f,f1,b;
84: GFS u,v;
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;
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: }
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:
137: void randsfum(int d,UM p)
138: {
139: int i;
140:
141: for ( i = 0; i < d; i++ )
142: COEF(p)[i] = _randomsf();
143: for ( i = d-1; i >= 0 && !COEF(p)[i]; i-- );
144: p->d = i;
145: }
146:
147: void pwrmodsfum(UM p,int e,UM f,UM pr)
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:
175: void spwrsfum(UM m,UM f,Z e,UM r)
176: {
177: UM t,s,q;
178: Z e1,rem,two;
179:
180: if ( !e ) {
181: DEG(r) = 0; COEF(r)[0] = _onesf();
182: } else if ( UNIZ(e) )
183: cpyum(f,r);
184: else {
1.2 ! noro 185: STOZ(2,two);
1.1 noro 186: divqrz(e,two,&e1,&rem);
187: t = W_UMALLOC(2*DEG(m)); spwrsfum(m,f,e1,t);
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 ( rem ) {
191: mulsfum(s,f,t); DEG(t) = divsfum(t,m,q); cpyum(t,r);
192: } else
193: cpyum(s,r);
194: }
195: }
196:
197: void tracemodsfum(UM m,UM f,int e,UM r)
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:
217: void make_qmatsf(UM p,UM *tab,int ***mp)
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:
238: void nullsf(int **mat,int n,int *ind)
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:
270: void null_to_solsf(int **mat,int *ind,int n,UM *r)
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:
299: void czsfum(UM f,UM *r)
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:
353: int berlekampsf(UM p,int df,UM *tab,UM *r)
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: }
399: /* NOT REACHED */
400: error("berlekampsf : cannot happen");
401: return 0;
402: }
403:
404: void minipolysf(UM f,UM p,UM mp)
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:
439: void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1)
440: {
441: int h,d1;
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:
458: int find_rootsf(UM p,int *root)
459: {
460: UM *r;
461: int i,n;
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:
471: void canzassf(UM f,int d,UM *r)
472: {
473: UM t,s,u,w,g,o;
474: Z n1,n2,n3,n4,n5,z;
475: UM *b;
476: int n,q,ed;
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();
490: if ( q % 2 ) {
1.2 ! noro 491: STOZ(q,n1); STOZ(d,z); pwrz(n1,z,&n2); subz(n2,ONE,&n3);
! 492: STOZ(2,n4); divsz(n3,n4,&n5);
1.1 noro 493: } else
494: ed = d*extdeg_sf();
495: while ( 1 ) {
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);
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:
513: /* Hensel related functions */
514:
515: int sfberle(V,V,P,int,GFS *,DCP *);
516: void sfgcdgen(P,ML,ML *);
517: void sfhenmain2(BM,UM,UM,int,BM *);
518: void ptosfbm(int,P,BM);
519: void sfhensel(int count,P f,V x,V y,int degbound,GFS *evp,P *sfp,ML *listp);
520:
521: /* f = f(x,y) */
522:
523: void sfhensel(int count,P f,V x,V y,int degbound,GFS *evp,P *sfp,ML *listp)
524: {
525: int i;
526: int fn;
527: ML rlist;
528: BM fl;
529: VL vl,nvl;
530: int dx,dy,bound;
531: GFS ev;
532: P f1,t,c,sf;
533: DCP dc,dct,dc0;
534: UM q,fm,hm;
535: UM *gm;
536: struct oEGT tmp0,tmp1,eg_hensel,eg_hensel_t;
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: }
543: if ( vl->next )
544: y = vl->next->v;
545: dx = getdeg(x,f);
546: dy = getdeg(y,f);
547: if ( dx == 1 ) {
548: *listp = rlist = MLALLOC(1); rlist->n = 1; rlist->c[0] = 0;
549: return;
550: }
551: fn = sfberle(x,y,f,count,&ev,&dc);
552: if ( fn <= 1 ) {
553: /* fn == 0 => short of evaluation points */
554: *listp = rlist = MLALLOC(1); rlist->n = fn; rlist->c[0] = 0;
555: return;
556: }
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: }
586:
587: /* convert mod y-a factors into UM */
588: gm = (UM *)ALLOCA(fn*sizeof(UM));
589: for ( i = 0; i < fn; i++, dc = NEXT(dc) ) {
590: gm[i] = W_UMALLOC(UDEG(dc->c));
591: ptosfum(dc->c,gm[i]);
592: }
593:
594: /* set bound */
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 */
597: bound = dy+2;
598:
599: /* f(x,y) -> f(x,y+ev) */
600: fl = BMALLOC(dx,bound);
601: ptosfbm(bound,f,fl);
602: if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev)));
603:
604: /* sf = f(x+ev) */
605: sfbmtop(fl,x,y,&sf);
606:
607: /* fm = fl mod y */
608: fm = W_UMALLOC(dx);
609: cpyum(COEF(fl)[0],fm);
610: hm = W_UMALLOC(dx);
611:
612: q = W_UMALLOC(dx);
613: rlist = MLALLOC(fn); rlist->n = fn; rlist->bound = bound;
614: if ( debug_sfbfctr )
615: fprintf(asir_out,"%d candidates\n",fn);
616: init_eg(&eg_hensel);
617: for ( i = 0; i < fn-1; i++ ) {
618: if ( debug_sfbfctr )
619: fprintf(asir_out,"deg(fm) = %d, deg(gm[%d]) = %d\n",
620: DEG(fm),i,DEG(gm[i]));
621: init_eg(&eg_hensel_t);
622: get_eg(&tmp0);
623: /* fl = gm[i]*hm mod y */
624: divsfum(fm,gm[i],hm);
625: /* fl is replaced by the cofactor of gk mod y^bound */
626: /* rlist->c[i] = gk */
627: sfhenmain2(fl,gm[i],hm,bound,(BM *)&rlist->c[i]);
628: cpyum(hm,fm);
629: get_eg(&tmp1); add_eg(&eg_hensel_t,&tmp0,&tmp1);
630: add_eg(&eg_hensel,&tmp0,&tmp1);
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);
638: fprintf(asir_out,"\n");
639: }
640: /* finally, fl must be the lift of gm[fn-1] */
641: rlist->c[i] = fl;
642:
643: #if 0
644: /* y -> y-a */
645: mev = _chsgnsf(FTOIF(CONT(ev)));
646: for ( i = 0; i < fn; i++ )
647: shiftsfbm((BM)(rlist->c[i]),mev);
648: #endif
649: *evp = ev;
650: *sfp = sf;
651: *listp = rlist;
652: }
653:
654: /* main variable of f = x */
655:
656: int sfberle(V x,V y,P f,int count,GFS *ev,DCP *dcp)
657: {
658: UM wf,wf1,wf2,wfs,gcd;
659: int fn,n;
660: GFS m,fm;
661: DCP dc,dct,dc0;
662: VL vl;
663: P lc,lc0,f0;
664: Obj obj;
665: int j,q,index,i;
666:
667: NEWVL(vl); vl->v = x;
668: NEWVL(NEXT(vl)); NEXT(vl)->v = y;
669: NEXT(NEXT(vl)) =0;
670: simp_ff((Obj)f,&obj); f = (P)obj;
1.2 ! noro 671: n = ZTOS(DEG(DC(f)));
1.1 noro 672: wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); wf2 = W_UMALLOC(n);
673: wfs = W_UMALLOC(n); gcd = W_UMALLOC(n);
674: q = field_order_sf();
675: lc = DC(f)->c;
676: for ( j = 0, fn = n + 1, index = 0;
677: index < q && j < count && fn > 1; index++ ) {
678: indextogfs(index,&m);
679: substp(vl,lc,y,(P)m,&lc0);
680: if ( lc0 ) {
681: substp(vl,f,y,(P)m,&f0);
682: ptosfum(f0,wf); cpyum(wf,wf1);
683: diffsfum(wf1,wf2); gcdsfum(wf1,wf2,gcd);
684: if ( DEG(gcd) == 0 ) {
685: ufctrsf(f0,&dc);
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: }
694: if ( index == q )
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:
705: void sfgcdgen(P f,ML blist,ML *clistp)
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:
732: /* f = g0*h0 mod y -> f = gk*hk mod y^(dy+1), f is replaced by hk */
733:
734: void sfhenmain2(BM f,UM g0,UM h0,int dy,BM *gp)
735: {
736: int i,k;
737: int dx;
738: UM wt,wa,wb,q,w1,w2,wh1,wg1,ws;
739: UM wc,wd,we,wz;
740: BM wb0,wb1;
741: int dg,dh;
742: BM fk,gk,hk;
743:
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);
754:
755: /* fk = gk*hk mod y^k */
756: W_BMALLOC(dx,dy,fk);
757: cpyum(COEF(f)[0],COEF(fk)[0]);
758: gk = BMALLOC(dg,dy);
759: cpyum(g0,COEF(gk)[0]);
760: W_BMALLOC(dh,dy,hk);
761: cpyum(h0,COEF(hk)[0]);
762:
763: wc = W_UMALLOC(2*dx); wd = W_UMALLOC(2*dx);
764: we = W_UMALLOC(2*dx); wz = W_UMALLOC(2*dx);
765:
766: /* compute wa,wb s.t. wa*g0+wb*h0 = 1 mod y */
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 */
770: eucsfum(w1,w2,wa,wb);
771:
772: if ( debug_sfbfctr)
773: fprintf(stderr,"dy=%d\n",dy);
774: for ( k = 1; k <= dy; k++ ) {
775: if ( debug_sfbfctr)
776: fprintf(stderr,".");
777:
778: /* at this point, f = gk*hk mod y^k */
779:
780: /* clear wt */
781: clearum(wt,dx);
782:
783: /* wt = (f-gk*hk)/y^k */
784: subsfum(COEF(f)[k],COEF(fk)[k],wt);
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 )
793: error("henmain2 : cannot happen(adj)");
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:
801: /* fk += ((wg1*hk+wh1*gk)*y^k+wg1*wh1*y^(2*k) mod y^(dy+1) */
802: /* wb0 = wh1*y^k */
803: clearbm(dx,wb0);
804: cpyum(wh1,COEF(wb0)[k]);
805:
806: /* wb1 = gk*wb0 mod y^(dy+1) */
807: clearbm(dx,wb1);
808: mulsfbm(gk,wb0,wb1);
809: /* fk += wb1 */
810: addtosfbm(wb1,fk);
811:
812: /* wb0 = wg1*y^k */
813: clearbm(dx,wb0);
814: cpyum(wg1,COEF(wb0)[k]);
815:
816: /* wb1 = hk*wb0 mod y^(dy+1) */
817: clearbm(dx,wb1);
818: mulsfbm(hk,wb0,wb1);
819: /* fk += wb1 */
820: addtosfbm(wb1,fk);
821:
822: /* fk += wg1*wh1*y^(2*k) mod y^(dy+1) */
823: if ( 2*k <= dy ) {
824: mulsfum(wg1,wh1,wt); addsfum(COEF(fk)[2*k],wt,ws);
825: cpyum(ws,COEF(fk)[2*k]);
826: }
827:
828: /* gk += wg1*y^k, hk += wh1*y^k */
829: cpyum(wg1,COEF(gk)[k]);
830: cpyum(wh1,COEF(hk)[k]);
831: }
832: if ( debug_sfbfctr)
833: fprintf(stderr,"\n");
834: *gp = gk;
835: DEG(f) = dy;
836: for ( i = 0; i <= dy; i++ )
837: cpyum(COEF(hk)[i],COEF(f)[i]);
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:
879: if ( debug_sfbfctr)
880: fprintf(stderr,"dy=%d\n",dy);
881: for ( k = 1; k <= dy; k++ ) {
882: if ( debug_sfbfctr)
883: fprintf(stderr,".");
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: }
922: if ( debug_sfbfctr)
923: fprintf(stderr,"\n");
924: DEG(a) = dy;
925: DEG(b) = dy;
926: *ap = a;
927: *bp = b;
928: }
929:
930: /* fl->c[i] = coef_y(f,i) */
931:
932: void ptosfbm(int dy,P f,BM fl)
933: {
934: DCP dc;
935: int d,i,dx;
936: UM t;
937:
1.2 ! noro 938: dx = ZTOS(DEG(DC(f)));
1.1 noro 939: if ( DEG(fl) < dy )
940: error("ptosfbm : invalid input");
941: DEG(fl) = dy;
942: clearbm(dx,fl);
943: t = UMALLOC(dy);
944: for ( dc = DC(f); dc; dc = NEXT(dc) ) {
1.2 ! noro 945: d = ZTOS(DEG(dc));
1.1 noro 946: ptosfum(COEF(dc),t);
947: for ( i = 0; i <= DEG(t); i++ )
948: COEF(COEF(fl)[i])[d] = COEF(t)[i];
949: }
950: for ( i = 0; i <= dy; i++ )
951: degum(COEF(fl)[i],dx);
952: }
953:
954: /* x : main variable */
955:
956: void sfbmtop(BM f,V x,V y,P *fp)
957: {
958: UM *c;
959: int i,j,d,a,dy;
960: GFS b;
961: DCP dc0,dc,dct;
962:
963: dy = DEG(f);
964: c = COEF(f);
965: d = degbm(f);
966:
967: dc0 = 0;
968: for ( i = 0; i <= d; i++ ) {
969: dc = 0;
970: for ( j = 0; j <= dy; j++ ) {
971: if ( DEG(c[j]) >= i && (a = COEF(c[j])[i]) ) {
972: NEWDC(dct);
1.2 ! noro 973: STOZ(j,DEG(dct));
1.1 noro 974: iftogfs(a,&b);
975: COEF(dct) = (P)b;
976: NEXT(dct) = dc;
977: dc = dct;
978: }
979: }
980: if ( dc ) {
981: NEWDC(dct);
1.2 ! noro 982: STOZ(i,DEG(dct));
1.1 noro 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:
994: void sfsqfr(P f,DCP *dcp)
995: {
996: Obj obj;
997: DCP dc;
998: VL vl;
999:
1000: simp_ff((Obj)f,&obj); f = (P)obj;
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
1008: sqfrsf(vl,f,dcp);
1009: }
1010:
1011: void sfusqfr(P f,DCP *dcp)
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++ ) {
1.2 ! noro 1031: NEWDC(dct); STOZ(udc[i].n,DEG(dct));
1.1 noro 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:
1039: #if 0
1040: void sfbsqfrmain(P f,V x,V y,DCP *dcp)
1041: {
1042: /* XXX*/
1043: }
1044:
1045: /* f is bivariate */
1046:
1047: void sfbsqfr(P f,V x,V y,DCP *dcp)
1048: {
1049: P t,rf,cx,cy;
1050: VL vl,rvl;
1051: DCP dcx,dcy,dct,dc;
1052: struct oVL vl0,vl1;
1053:
1054: /* cy(y) = cont(f,x), f /= cy */
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 */
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: }
1082: }
1083: #endif
1084:
1085: void sfdtest(P,ML,V,V,DCP *);
1086:
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)
1090: {
1091: ML list;
1092: P sf;
1093: GFS ev;
1094: DCP dc,dct;
1095: BM fl;
1096: int dx,dy;
1097:
1098: /* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */
1099: sfhensel(5,f,x,y,degbound,&ev,&sf,&list);
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: }
1108: sfdtest(sf,list,x,y,&dc);
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: }
1118: }
1119: *dcp = dc;
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]*... */
1133: sfhensel(5,f,x,y,degbound,&ev,&sf,&list);
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: }
1147: }
1148:
1149: /* f = f(x,y) = list->c[0]*list->c[1]*... mod y^(list->bound+1) */
1150:
1151: void sfdtest(P f,ML list,V x,V y,DCP *dcp)
1152: {
1153: int np,dx,dy;
1154: int i,j,k,bound;
1155: int *win;
1156: P g,lcg,factor,cofactor,lcyx;
1157: P csum;
1158: DCP dcf,dcf0,dc;
1159: BM *c;
1160: BM lcy;
1161: UM lcg0,lcy0,w;
1162: UM *d1c;
1163: ML wlist;
1164: struct oVL vl1,vl0;
1165: VL vl;
1166: int z,dt,dtok;
1167:
1168: /* vl = [x,y] */
1169: vl0.v = x; vl0.next = &vl1; vl1.v = y; vl1.next = 0; vl = &vl0;
1170:
1171: /* setup various structures and arrays */
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;
1178: bound = wlist->bound = list->bound;
1179: c = (BM *)COEF(wlist);
1180: bcopy((char *)COEF(list),(char *)c,(int)(sizeof(BM)*np));
1181:
1182: lcg0 = W_UMALLOC(2*dy);
1183:
1184: /* initialize g by f */
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);
1195:
1196: /* initialize lcy by LC(f) */
1197: W_BMALLOC(0,dy,lcy);
1198: NEWDC(dc); COEF(dc) = COEF(DC(g)); DEG(dc) = 0;
1199: NEWP(lcyx); VR(lcyx) = x; DC(lcyx) = dc;
1200: ptosfbm(dy,lcyx,lcy);
1201:
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:
1221: if ( debug_sfbfctr)
1222: fprintf(stderr,"np = %d\n",np);
1223: dtok = 0;
1224: for ( g = f, k = 1, dcf = dcf0 = 0, win[0] = 1, --np, z = 0; ; z++ ) {
1225: if ( debug_sfbfctr && !(z % 1000) ) fprintf(stderr,".");
1226: dt = sfdegtest(dy,bound,d1c,k,win);
1227: if ( dt )
1228: dtok++;
1229: if ( dt && sfdtestmain(vl,lcg,lcg0,lcy,csum,wlist,
1230: k,win,&factor,&cofactor) ) {
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:
1237: /* update lcg0 */
1238: const_term(lcg,lcg0);
1239:
1240: /* update csum */
1241: sfcsump(vl,lcg,&csum);
1242:
1243: /* update dy */
1244: dy = getdeg(y,g);
1245:
1246: /* update lcy */
1247: clearbm(0,lcy);
1248: COEF(dc) = COEF(DC(g));
1249: ptosfbm(dy,lcyx,lcy);
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;
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: }
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: }
1290: }
1291: if ( debug_sfbfctr )
1292: fprintf(stderr,"total %d, omitted by degtest %d\n",z,z-dtok);
1293: NEXTDC(dcf0,dcf); COEF(dcf) = g;
1294: DEG(dcf) = ONE; NEXT(dcf) = 0; *dcp = dcf0;
1295: }
1296:
1297: void extractcoefbm(BM f,int dx,UM r)
1298: {
1299: int j;
1300: UM fj;
1301:
1302: for ( j = DEG(f); j >= 0; j-- ) {
1303: fj = COEF(f)[j];
1304: if ( fj && DEG(fj) >= dx ) {
1305: COEF(r)[j] = COEF(fj)[dx];
1306: } else
1307: COEF(r)[j] = 0;
1308: }
1309: degum(r,DEG(f));
1310: }
1311:
1312: /* deg_y(prod mod y^(bound+1)) <= dy ? */
1313:
1314: int sfdegtest(int dy,int bound,UM *d1c,int k,int *in)
1315: {
1316: int i,j;
1317: UM w,w1,wt;
1318: BM t;
1319:
1320: w = W_UMALLOC(bound);
1321: w1 = W_UMALLOC(bound);
1322: clearum(w,bound);
1323: for ( i = 0; i < k; i++ ) {
1324: addsfum(w,d1c[in[i]],w1); wt = w; w = w1; w1 = wt;
1325: }
1326: return DEG(w) <= dy ? 1 : 0;
1327: }
1328:
1329: /* lcy = LC(g), lcg = lcy*g, lcg0 = const part of lcg */
1330:
1331: int sfdtestmain(VL vl,P lcg,UM lcg0,BM lcy,P csum,ML list,
1332: int k,int *in,P *fp,P *cofp)
1333: {
1334: P fmul,csumg,q,cont;
1335: V x,y;
1336:
1337: x = vl->v;
1338: y = vl->next->v;
1339: if (!sfctest(lcg0,lcy,list,k,in))
1340: return 0;
1341: mulsfbmarray(UDEG(lcg),lcy,list,k,in,x,y,&fmul);
1342: if ( csum ) {
1343: sfcsump(vl,fmul,&csumg);
1344: if ( csumg ) {
1345: if ( !divtp(vl,csum,csumg,&q) )
1346: return 0;
1347: }
1348: }
1349: if ( divtp_by_sfbm(vl,lcg,fmul,&q) ) {
1350: cont_pp_sfp(vl,fmul,&cont,fp);
1351: cont_pp_sfp(vl,q,&cont,cofp);
1352: return 1;
1353: } else
1354: return 0;
1355: }
1356:
1357: void const_term(P f,UM c)
1358: {
1359: DCP dc;
1360:
1361: for ( dc = DC(f); dc && DEG(dc); dc = NEXT(dc) );
1362: if ( dc )
1363: ptosfum(COEF(dc),c);
1364: else
1365: DEG(c) = -1;
1366: }
1367:
1368: void const_term_sfbm(BM f,UM c)
1369: {
1370: int i,dy;
1371:
1372: dy = DEG(f);
1373: for ( i = 0; i <= dy; i++ )
1374: if ( DEG(COEF(f)[i]) >= 0 )
1375: COEF(c)[i] = COEF(COEF(f)[i])[0];
1376: else
1377: COEF(c)[i] = 0;
1378: degum(c,dy);
1379: }
1380:
1381: /* lcy*(product of const part) | lcg0 ? */
1382:
1383: int sfctest(UM lcg0,BM lcy,ML list,int k,int *in)
1384: {
1385: int dy,i,dr;
1386: UM t,s,u,w;
1387: BM *l;
1388:
1389: dy = list->bound;
1390: t = W_UMALLOC(2*dy);
1391: s = W_UMALLOC(2*dy);
1392: u = W_UMALLOC(2*dy);
1393: const_term_sfbm(lcy,t);
1394: if ( DEG(t) < 0 )
1395: return 1;
1396:
1397: l = (BM *)list->c;
1398: for ( i = 0; i < k; i++ ) {
1399: const_term_sfbm(l[in[i]],s);
1400: mulsfum(t,s,u);
1401: if ( DEG(u) > dy )
1402: degum(u,dy);
1403: w = t; t = u; u = w;
1404: }
1405: cpyum(lcg0,s);
1406: dr = divsfum(s,t,u);
1407: if ( dr >= 0 )
1408: return 0;
1409: else
1410: return 1;
1411: }
1412:
1413: /* main var of f is x */
1414:
1415: void mulsfbmarray(int dx,BM lcy,ML list,int k,int *in,V x,V y,P *g)
1416: {
1417: int dy,i;
1418: BM wb0,wb1,t;
1419: BM *l;
1420:
1421: dy = list->bound;
1422: W_BMALLOC(dx,dy,wb0); W_BMALLOC(dx,dy,wb1);
1423: l = (BM *)list->c;
1424: clearbm(dx,wb0);
1425: mulsfbm(lcy,l[in[0]],wb0);
1426: for ( i = 1; i < k; i++ ) {
1427: clearbm(dx,wb1);
1428: mulsfbm(l[in[i]],wb0,wb1);
1429: t = wb0; wb0 = wb1; wb1 = t;
1430: }
1431: sfbmtop(wb0,x,y,g);
1432: }
1433:
1434: void sfcsump(VL vl,P f,P *s)
1435: {
1436: P t,u;
1437: DCP dc;
1438:
1439: for ( dc = DC(f), t = 0; dc; dc = NEXT(dc) ) {
1440: addp(vl,COEF(dc),t,&u); t = u;
1441: }
1442: *s = t;
1443: }
1444:
1445: /* *fp = primitive part of f w.r.t. x */
1446:
1447: void cont_pp_sfp(VL vl,P f,P *cp,P *fp)
1448: {
1449: V x,y;
1450: int d;
1451: UM t,s,gcd;
1452: DCP dc;
1453: GFS g;
1454:
1455: x = vl->v;
1456: y = vl->next->v;
1457: d = getdeg(y,f);
1458: if ( d == 0 ) {
1459: itogfs(1,&g);
1460: *cp = (P)g;
1461: *fp = f; /* XXX */
1462: } else {
1463: t = W_UMALLOC(2*d);
1464: s = W_UMALLOC(2*d);
1465: gcd = W_UMALLOC(2*d);
1466: dc = DC(f);
1467: ptosfum(COEF(dc),gcd);
1468: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
1469: ptosfum(COEF(dc),t);
1470: gcdsfum(gcd,t,s);
1471: cpyum(s,gcd);
1472: }
1473: sfumtop(y,gcd,cp);
1474: divsp(vl,f,*cp,fp);
1475: }
1476: }
1477:
1478: int divtp_by_sfbm(VL vl,P f,P g,P *qp)
1479: {
1480: V x,y;
1481: int fx,fy,gx,gy;
1482: BM fl,gl,ql;
1483: UM *cf,*cg,*cq;
1484: UM hg,q,t,s;
1485: int i,j,dr;
1486:
1487: x = vl->v; y = vl->next->v;
1488: fx = getdeg(x,f); fy = getdeg(y,f);
1489: gx = getdeg(x,g); gy = getdeg(y,g);
1490:
1491: if ( fx < gx || fy < gy )
1492: return 0;
1493: W_BMALLOC(fx,fy,fl); ptosfbm(fy,f,fl); cf = COEF(fl);
1494: W_BMALLOC(gx,gy,gl); ptosfbm(gy,g,gl); cg = COEF(gl);
1495: W_BMALLOC(fx-gx,fy-gy,ql); cq = COEF(ql);
1496:
1497: hg = cg[gy];
1498: q = W_UMALLOC(fx); t = W_UMALLOC(fx); s = W_UMALLOC(fx);
1499:
1500: for ( i = fy; i >= gy; i-- ) {
1501: if ( DEG(cf[i]) < 0 )
1502: continue;
1503: dr = divsfum(cf[i],hg,q);
1504: if ( dr >= 0 )
1505: return 0;
1506: if ( DEG(q) > fx-gx )
1507: return 0;
1508: cpyum(q,cq[i-gy]);
1509: for ( j = 0; j <= gy; j++ ) {
1510: mulsfum(cg[j],q,t);
1511: subsfum(cf[j+i-gy],t,s);
1512: cpyum(s,cf[j+i-gy]);
1513: }
1514: }
1515: for ( j = gy-1; j >= 0 && DEG(cf[j]) < 0; j-- );
1516: if ( j >= 0 )
1517: return 0;
1518: sfbmtop(ql,x,y,qp);
1519: return 1;
1520: }
1521:
1522: /* XXX generate an irreducible poly of degree n */
1523:
1524: extern int current_gfs_q1;
1525: extern int *current_gfs_ntoi;
1526:
1527: void generate_defpoly_sfum(int n,UM *dp)
1528: {
1529: UM r,dr,t,g;
1530: UM *f;
1531: int *c,*w;
1532: int max,i,j;
1533:
1534: *dp = r = UMALLOC(n);
1535: DEG(r) = n;
1536: c = COEF(r);
1537: c[n] = _onesf();
1538: max = current_gfs_q1;
1539: w = (int *)ALLOCA(n*sizeof(int));
1540: bzero(w,n*sizeof(int));
1541:
1542: dr = W_UMALLOC(n); t = W_UMALLOC(n); g = W_UMALLOC(n);
1543: f = (UM *)ALLOCA((n+1)*sizeof(UM));
1544: while ( 1 ) {
1545: for ( i = 0; i < n && w[i] == max; i++ );
1546: if ( i == n ) {
1547: /* XXX cannot happen */
1548: error("generate_defpoly_sfum : cannot happen");
1549: }
1550: for ( j = 0; j < i; j++ )
1551: w[j] = 0;
1552: w[i]++;
1553: if ( !current_gfs_ntoi )
1554: for ( i = 0; i < n; i++ )
1555: c[i] = w[i]?FTOIF(w[i]):0;
1556: else
1557: for ( i = 0; i < n; i++ )
1558: c[i] = w[i]?FTOIF(w[i]-1):0;
1559: if ( !c[0] )
1560: continue;
1561: diffsfum(r,dr); cpyum(r,t); gcdsfum(t,dr,g);
1562: if ( DEG(g) > 0 )
1563: continue;
1564:
1565: czsfum(r,f);
1566: for ( i = 0; f[i]; i++ );
1567: if ( i == 1 )
1568: return;
1569: }
1570: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Hgfs.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine