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