=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/Hgfs.c,v retrieving revision 1.12 retrieving revision 1.30 diff -u -p -r1.12 -r1.30 --- OpenXM_contrib2/asir2000/engine/Hgfs.c 2001/06/27 04:07:57 1.12 +++ OpenXM_contrib2/asir2000/engine/Hgfs.c 2003/01/15 10:21:17 1.30 @@ -1,26 +1,14 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.11 2001/06/26 09:47:05 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.29 2002/11/22 08:44:57 noro Exp $ */ #include "ca.h" +#include "inline.h" -struct p_pair { - UM p0; - UM p1; - struct p_pair *next; -}; +int debug_sfbfctr; -void canzassf(UM,int,UM *); -void lnfsf(int,UM,UM,struct p_pair *,UM,UM); -void minipolysf(UM,UM,UM); -void czsfum(UM,UM *); -void gensqfrsfum(UM,DUM); -void sfbmtop(int,BM,V,V,P *); -void pp_sfp(VL,P,P *); -void sfcsump(VL,P,P *); -void mulsfbmarray(int,BM,ML,int,int *,V,V,P *); -void const_term(P,UM); +void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1); +void extractcoefbm(BM f,int dx,UM r); -int comp_dum(a,b) -DUM a,b; +int comp_dum(DUM a,DUM b) { if ( DEG(a->f) > DEG(b->f) ) return -1; @@ -30,21 +18,21 @@ DUM a,b; return 0; } -void fctrsf(p,dcp) -P p; -DCP *dcp; +void ufctrsf(P p,DCP *dcp) { int n,i,j,k; DCP dc,dc0; P lc; - P zp; UM mp; UM *tl; + Obj obj; struct oDUM *udc,*udc1; - simp_ff(p,&zp); p = zp; + simp_ff((Obj)p,&obj); p = (P)obj; if ( !p ) { - *dcp = 0; return; + NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; + NEXT(dc) = 0; *dcp = dc; + return; } mp = W_UMALLOC(UDEG(p)); ptosfum(p,mp); @@ -89,12 +77,11 @@ DCP *dcp; NEXT(dc) = 0; *dcp = dc0; } -void gensqfrsfum(p,dc) -UM p; -struct oDUM *dc; +void gensqfrsfum(UM p,struct oDUM *dc) { int n,i,j,d,mod; UM t,s,g,f,f1,b; + GFS u,v; if ( (n = DEG(p)) == 1 ) { dc[0].f = UMALLOC(DEG(p)); cpyum(p,dc[0].f); dc[0].n = 1; @@ -133,8 +120,11 @@ struct oDUM *dc; break; else { DEG(s) = DEG(t)/mod; - for ( j = 0; j <= DEG(t); j++ ) - COEF(s)[j] = COEF(t)[j*mod]; + for ( j = 0; j <= DEG(t); j++ ) { + iftogfs(COEF(t)[j*mod],&u); + pthrootgfs(u,&v); + COEF(s)[j] = v?FTOIF(CONT(v)):0; + } cpyum(s,b); d *= mod; } } @@ -144,9 +134,7 @@ struct oDUM *dc; } } -void randsfum(d,p) -int d; -UM p; +void randsfum(int d,UM p) { int i; @@ -156,9 +144,7 @@ UM p; p->d = i; } -void pwrmodsfum(p,e,f,pr) -int e; -UM p,f,pr; +void pwrmodsfum(UM p,int e,UM f,UM pr) { UM wt,ws,q; @@ -186,9 +172,7 @@ UM p,f,pr; } } -void spwrsfum(m,f,e,r) -UM f,m,r; -N e; +void spwrsfum(UM m,UM f,N e,UM r) { UM t,s,q; N e1; @@ -210,9 +194,7 @@ N e; } } -void tracemodsfum(m,f,e,r) -UM f,m,r; -int e; +void tracemodsfum(UM m,UM f,int e,UM r) { UM t,s,q,u; int i; @@ -232,10 +214,7 @@ int e; cpyum(s,r); } -void make_qmatsf(p,tab,mp) -UM p; -UM *tab; -int ***mp; +void make_qmatsf(UM p,UM *tab,int ***mp) { int n,i,j; int *c; @@ -256,10 +235,7 @@ int ***mp; mat[i][i] = _subsf(mat[i][i],one); } -void nullsf(mat,n,ind) -int **mat; -int *ind; -int n; +void nullsf(int **mat,int n,int *ind) { int i,j,l,s,h,inv; int *t,*u; @@ -291,11 +267,7 @@ int n; } } -void null_to_solsf(mat,ind,n,r) -int **mat; -int *ind; -int n; -UM *r; +void null_to_solsf(int **mat,int *ind,int n,UM *r) { int i,j,k,l; int *c; @@ -324,8 +296,7 @@ nullsf(mat,n,ind) null_to_solsf(ind,n,r) */ -void czsfum(f,r) -UM f,*r; +void czsfum(UM f,UM *r) { int i,j; int d,n,ord; @@ -379,10 +350,7 @@ UM f,*r; r[j] = 0; } -int berlekampsf(p,df,tab,r) -UM p; -int df; -UM *tab,*r; +int berlekampsf(UM p,int df,UM *tab,UM *r) { int n,i,j,k,nf,d,nr; int **mat; @@ -428,10 +396,12 @@ UM *tab,*r; } } } + /* NOT REACHED */ + error("berlekampsf : cannot happen"); + return 0; } -void minipolysf(f,p,mp) -UM f,p,mp; +void minipolysf(UM f,UM p,UM mp) { struct p_pair *list,*l,*l1,*lprev; int n,d; @@ -466,13 +436,9 @@ UM f,p,mp; } } -void lnfsf(n,p0,p1,list,np0,np1) -int n; -UM p0,p1; -struct p_pair *list; -UM np0,np1; +void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1) { - int inv,h,d1; + int h,d1; UM t0,t1,s0,s1; struct p_pair *l; @@ -489,12 +455,10 @@ UM np0,np1; } } -int find_rootsf(p,root) -UM p; -int *root; +int find_rootsf(UM p,int *root) { UM *r; - int i,j,n; + int i,n; n = DEG(p); r = ALLOCA((DEG(p))*sizeof(UM)); @@ -504,14 +468,12 @@ int *root; return n; } -void canzassf(f,d,r) -UM f,*r; -int d; +void canzassf(UM f,int d,UM *r) { UM t,s,u,w,g,o; N n1,n2,n3,n4,n5; UM *b; - int n,m,i,q,ed; + int n,q,ed; if ( DEG(f) == d ) { r[0] = UMALLOC(d); cpyum(f,r[0]); @@ -550,34 +512,26 @@ int d; /* Hensel related functions */ -int sfberle(VL,P,int,GFS *,DCP *); +int sfberle(V,V,P,int,GFS *,DCP *); void sfgcdgen(P,ML,ML *); -void sfhenmain(LUM,ML,ML,ML *); -void ptosflum(int,P,LUM); void sfhenmain2(BM,UM,UM,int,BM *); void ptosfbm(int,P,BM); +void sfhensel(int count,P f,V x,V y,int degbound,GFS *evp,P *sfp,ML *listp); /* f = f(x,y) */ -void sfhensel(count,f,x,evp,sfp,listp) -int count; -P f; -V x; -GFS *evp; -P *sfp; -ML *listp; +void sfhensel(int count,P f,V x,V y,int degbound,GFS *evp,P *sfp,ML *listp) { - int i,j; - int fn,n,bound; + int i; + int fn; ML rlist; BM fl; VL vl,nvl; - V y; - int dx,dy,mev; + int dx,dy,bound; GFS ev; - P f1,t,yev,c,sf; - DCP dc; - UM w,w1,q,fm,hm; + P f1,t,c,sf; + DCP dc,dct,dc0; + UM q,fm,hm; UM *gm; struct oEGT tmp0,tmp1,eg_hensel,eg_hensel_t; @@ -586,22 +540,49 @@ ML *listp; reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&f1); vl = nvl; f = f1; } - y = vl->next->v; + if ( vl->next ) + y = vl->next->v; dx = getdeg(x,f); dy = getdeg(y,f); if ( dx == 1 ) { *listp = rlist = MLALLOC(1); rlist->n = 1; rlist->c[0] = 0; return; } - fn = sfberle(vl,f,count,&ev,&dc); + fn = sfberle(x,y,f,count,&ev,&dc); if ( fn <= 1 ) { /* fn == 0 => short of evaluation points */ *listp = rlist = MLALLOC(1); rlist->n = fn; rlist->c[0] = 0; return; } - /* pass the the leading coeff. to the first element */ - c = dc->c; dc = NEXT(dc); - mulp(vl,dc->c,c,&t); dc->c = t; + if ( degbound >= 0 ) { + /* + * reconstruct dc so that + * dc[1],... : factors satisfy degree bound + * dc[0] : product of others + */ + c = dc->c; dc = NEXT(dc); + dc0 = 0; + fn = 0; + while ( dc ) { + if ( getdeg(x,COEF(dc)) <= degbound ) { + dct = NEXT(dc); NEXT(dc) = dc0; dc0 = dc; dc = dct; + fn++; + } else { + mulp(vl,COEF(dc),c,&t); c = t; + dc = NEXT(dc); + } + } + if ( OID(c) == O_P ) { + NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = dc0; + fn++; + } else { + mulp(vl,dc0->c,c,&t); dc0->c = t; dc = dc0; + } + } else { + /* pass the the leading coeff. to the first element */ + c = dc->c; dc = NEXT(dc); + mulp(vl,dc->c,c,&t); dc->c = t; + } /* convert mod y-a factors into UM */ gm = (UM *)ALLOCA(fn*sizeof(UM)); @@ -610,14 +591,18 @@ ML *listp; ptosfum(dc->c,gm[i]); } - bound = dy+1; + /* set bound */ + /* g | f, lc_y(g) = lc_y(f) => deg_y(g) <= deg_y(f) */ + /* so, bound = dy is sufficient, but we use slightly large value */ + bound = dy+2; + /* f(x,y) -> f(x,y+ev) */ fl = BMALLOC(dx,bound); ptosfbm(bound,f,fl); - shiftsfbm(bound,fl,FTOIF(CONT(ev))); + if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev))); /* sf = f(x+ev) */ - sfbmtop(bound,fl,x,y,&sf); + sfbmtop(fl,x,y,&sf); /* fm = fl mod y */ fm = W_UMALLOC(dx); @@ -626,10 +611,12 @@ ML *listp; q = W_UMALLOC(dx); rlist = MLALLOC(fn); rlist->n = fn; rlist->bound = bound; - fprintf(asir_out,"%d candidates\n",fn); + if ( debug_sfbfctr ) + fprintf(asir_out,"%d candidates\n",fn); init_eg(&eg_hensel); for ( i = 0; i < fn-1; i++ ) { - fprintf(asir_out,"deg(fm) = %d, deg(gm[%d]) = %d\n", + if ( debug_sfbfctr ) + fprintf(asir_out,"deg(fm) = %d, deg(gm[%d]) = %d\n", DEG(fm),i,DEG(gm[i])); init_eg(&eg_hensel_t); get_eg(&tmp0); @@ -641,11 +628,15 @@ ML *listp; cpyum(hm,fm); get_eg(&tmp1); add_eg(&eg_hensel_t,&tmp0,&tmp1); add_eg(&eg_hensel,&tmp0,&tmp1); - print_eg("Hensel",&eg_hensel_t); + if ( debug_sfbfctr) { + print_eg("Hensel",&eg_hensel_t); + fprintf(asir_out,"\n"); + } + } + if ( debug_sfbfctr) { + print_eg("Hensel total",&eg_hensel); fprintf(asir_out,"\n"); } - print_eg("Hensel total",&eg_hensel); - fprintf(asir_out,"\n"); /* finally, fl must be the lift of gm[fn-1] */ rlist->c[i] = fl; @@ -653,7 +644,7 @@ ML *listp; /* y -> y-a */ mev = _chsgnsf(FTOIF(CONT(ev))); for ( i = 0; i < fn; i++ ) - shiftsfbm(bound,(BM)(rlist->c[i]),mev); + shiftsfbm((BM)(rlist->c[i]),mev); #endif *evp = ev; *sfp = sf; @@ -662,41 +653,36 @@ ML *listp; /* main variable of f = x */ -int sfberle(vl,f,count,ev,dcp) -VL vl; -P f; -int count; -GFS *ev; -DCP *dcp; +int sfberle(V x,V y,P f,int count,GFS *ev,DCP *dcp) { UM wf,wf1,wf2,wfs,gcd; - ML flist; - int fn,fn1,n; + int fn,n; GFS m,fm; DCP dc,dct,dc0; - VL nvl; - V x,y; - P g,lc,lc0,f0; - int j,q1,index,i; + VL vl; + P lc,lc0,f0; + Obj obj; + int j,q,index,i; - clctv(vl,f,&nvl); vl = nvl; - x = vl->v; y = vl->next->v; - simp_ff(f,&g); g = f; + NEWVL(vl); vl->v = x; + NEWVL(NEXT(vl)); NEXT(vl)->v = y; + NEXT(NEXT(vl)) =0; + simp_ff((Obj)f,&obj); f = (P)obj; n = QTOS(DEG(DC(f))); wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); wf2 = W_UMALLOC(n); wfs = W_UMALLOC(n); gcd = W_UMALLOC(n); - q1 = field_order_sf()-1; + q = field_order_sf(); lc = DC(f)->c; for ( j = 0, fn = n + 1, index = 0; - index < q1 && j < count && fn > 1; index++ ) { - MKGFS(index,m); + index < q && j < count && fn > 1; index++ ) { + indextogfs(index,&m); substp(vl,lc,y,(P)m,&lc0); if ( lc0 ) { substp(vl,f,y,(P)m,&f0); ptosfum(f0,wf); cpyum(wf,wf1); diffsfum(wf1,wf2); gcdsfum(wf1,wf2,gcd); if ( DEG(gcd) == 0 ) { - fctrsf(f0,&dc); + ufctrsf(f0,&dc); for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ ); if ( i < fn ) { dc0 = dc; fn = i; fm = m; @@ -705,7 +691,7 @@ DCP *dcp; } } } - if ( index == q1 ) + if ( index == q ) return 0; else if ( fn == 1 ) return 1; @@ -716,9 +702,7 @@ DCP *dcp; } } -void sfgcdgen(f,blist,clistp) -P f; -ML blist,*clistp; +void sfgcdgen(P f,ML blist,ML *clistp) { int i; int n,d,np; @@ -745,146 +729,60 @@ ML blist,*clistp; } } -/* - sfhenmain(fl,bqlist,cqlist,listp) -*/ +/* f = g0*h0 mod y -> f = gk*hk mod y^(dy+1), f is replaced by hk */ -void sfhenmain(f,bqlist,cqlist,listp) -LUM f; -ML bqlist,cqlist,*listp; +void sfhenmain2(BM f,UM g0,UM h0,int dy,BM *gp) { - int i,j,k; - int *px,*py; - int **pp,**pp1; - int n,np,bound,dr,tmp; - UM wt,wq0,wq,wr,wm,wm0,wa,q; - LUM wb0,wb1,tlum; - UM *b,*c; - LUM *l; - ML list; - - n = DEG(f); np = bqlist->n; bound = bqlist->bound; - *listp = list = MLALLOC(n); - list->n = np; list->bound = bound; - W_LUMALLOC(n,bound,wb0); W_LUMALLOC(n,bound,wb1); - wt = W_UMALLOC(n); wq0 = W_UMALLOC(n); wq = W_UMALLOC(n); - wr = W_UMALLOC(n); wm = W_UMALLOC(2*n); wm0 = W_UMALLOC(2*n); - wa = W_UMALLOC(2*n); q = W_UMALLOC(2*n); - b = (UM *)bqlist->c; c = (UM *)cqlist->c; l = (LUM *)list->c; - for ( i = 0; i < np; i++ ) { - l[i] = LUMALLOC(DEG(b[i]),bound); - for ( j = DEG(b[i]), pp = COEF(l[i]), px = COEF(b[i]); j >= 0; j-- ) - pp[j][0] = px[j]; - } - for ( i = 1; i < bound; i++ ) { - fprintf(stderr,"."); - /* at this point, f = l[0]*l[1]*...*l[np-1] mod y^i */ - mulsflum(i+1,l[0],l[1],wb0); - for ( j = 2; j < np; j++ ) { - mulsflum(i+1,l[j],wb0,wb1); - tlum = wb0; wb0 = wb1; wb1 = tlum; - } -#if 0 - /* check */ - for ( j = 0, pp = COEF(f), pp1 = COEF(wb0); j <= n; j++ ) - for ( k = 0; k < i; k++ ) - if ( pp[j][k] != pp1[j][k] ) - error("henmain : cannot happen"); -#endif - for ( j = n, px = COEF(wt); j >= 0; j-- ) - px[j] = 0; - for ( j = n, pp = COEF(f), pp1 = COEF(wb0); j >= 0; j-- ) - COEF(wt)[j] = _subsf(pp[j][i],pp1[j][i]); - degum(wt,n); - for ( j = n, px = COEF(wq0); j >= 0; j-- ) - px[j] = 0; - for ( j = 1; j < np; j++ ) { - mulsfum(wt,c[j],wm); dr = divsfum(wm,b[j],q); - for ( k = DEG(q), px = COEF(wq0), py = COEF(q); k >= 0; k-- ) - px[k] = _addsf(px[k],py[k]); - for ( k = dr, pp = COEF(l[j]), px = COEF(wm); k >= 0; k-- ) - pp[k][i] = px[k]; - } - degum(wq0,n); mulsfum(wq0,b[0],wm); - mulsfum(wt,c[0],wm0); addsfum(wm,wm0,wa); - for ( j = DEG(wa), pp = COEF(l[0]), px = COEF(wa); j >= 0; j-- ) - pp[j][i] = px[j]; - for ( j = n, px = COEF(wq0); j >= 0; j-- ) - px[j] = 0; - } - fprintf(stderr,"\n"); -} - -/* f = g0*h0 mod y -> f = gk*hk mod y^bound, f is replaced by hk */ - -void sfhenmain2(f,g0,h0,bound,gp) -BM f; -UM g0,h0; -int bound; -BM *gp; -{ - int i,j,k,l; - int *px,*py; - int **pp,**pp1; - int n,np,dr,tmp; - UM wt,wa,wb,wq,wm,q,w1,w2,wh1,wg1,ws; + int i,k; + int dx; + UM wt,wa,wb,q,w1,w2,wh1,wg1,ws; UM wc,wd,we,wz; BM wb0,wb1; - int ng,nh; + int dg,dh; BM fk,gk,hk; - n = degsfbm(bound,f); - ng = g0->d; - nh = h0->d; + if ( DEG(f) < dy ) + error("sfhenmain2 : invalid input"); - W_BMALLOC(n,bound,wb0); W_BMALLOC(n,bound,wb1); - wt = W_UMALLOC(n); ws = W_UMALLOC(n); - wq = W_UMALLOC(n); q = W_UMALLOC(2*n); - wg1 = W_UMALLOC(2*n); wh1 = W_UMALLOC(2*n); + dx = degbm(f); + dg = DEG(g0); + dh = DEG(h0); + W_BMALLOC(dx,dy,wb0); W_BMALLOC(dx,dy,wb1); + wt = W_UMALLOC(dx); ws = W_UMALLOC(dx); q = W_UMALLOC(2*dx); + wg1 = W_UMALLOC(2*dx); wh1 = W_UMALLOC(2*dx); + /* fk = gk*hk mod y^k */ - W_BMALLOC(n,bound,fk); + W_BMALLOC(dx,dy,fk); cpyum(COEF(f)[0],COEF(fk)[0]); - gk = BMALLOC(ng,bound); + gk = BMALLOC(dg,dy); cpyum(g0,COEF(gk)[0]); - W_BMALLOC(nh,bound,hk); + W_BMALLOC(dh,dy,hk); cpyum(h0,COEF(hk)[0]); - wc = W_UMALLOC(2*n); wd = W_UMALLOC(2*n); - we = W_UMALLOC(2*n); wz = W_UMALLOC(2*n); + wc = W_UMALLOC(2*dx); wd = W_UMALLOC(2*dx); + we = W_UMALLOC(2*dx); wz = W_UMALLOC(2*dx); /* compute wa,wb s.t. wa*g0+wb*h0 = 1 mod y */ - w1 = W_UMALLOC(ng); cpyum(g0,w1); - w2 = W_UMALLOC(nh); cpyum(h0,w2); - wa = W_UMALLOC(2*n); wb = W_UMALLOC(2*n); /* XXX */ + w1 = W_UMALLOC(dg); cpyum(g0,w1); + w2 = W_UMALLOC(dh); cpyum(h0,w2); + wa = W_UMALLOC(2*dx); wb = W_UMALLOC(2*dx); /* XXX */ eucsfum(w1,w2,wa,wb); -#if 0 - mulsfum(wa,g0,wc); mulsfum(wb,h0,wd); addsfum(wc,wd,we); - if ( DEG(we) != 0 || COEF(we)[0] != _onesf() ) - error("henmain2 : cannot happen(euc)"); -#endif + if ( debug_sfbfctr) + fprintf(stderr,"dy=%d\n",dy); + for ( k = 1; k <= dy; k++ ) { + if ( debug_sfbfctr) + fprintf(stderr,"."); - fprintf(stderr,"bound=%d\n",bound); - for ( k = 1; k < bound; k++ ) { - fprintf(stderr,"."); - /* at this point, f = gk*hk mod y^k */ -#if 0 - for ( j = 0; j < k; j++ ) - if ( !isequalum(COEF(f)[j],COEF(fk)[j]) ) - error("henmain2 : cannot happen(history)"); -#endif /* clear wt */ - bzero(COEF(wt),(n+1)*sizeof(int)); + clearum(wt,dx); /* wt = (f-gk*hk)/y^k */ subsfum(COEF(f)[k],COEF(fk)[k],wt); - /* clear wq */ - bzero(COEF(wq),(n+1)*sizeof(int)); - /* compute wf1,wg1 s.t. wh1*g0+wg1*h0 = wt */ mulsfum(wa,wt,wh1); DEG(wh1) = divsfum(wh1,h0,q); mulsfum(wh1,g0,wc); subsfum(wt,wc,wd); DEG(wd) = divsfum(wd,h0,wg1); @@ -900,31 +798,29 @@ BM *gp; error("henmain2 : cannot happen"); #endif - /* fk += ((wg1*hk+wh1*gk)*y^k+wg1*wh1*y^(2*k) mod y^bound */ + /* fk += ((wg1*hk+wh1*gk)*y^k+wg1*wh1*y^(2*k) mod y^(dy+1) */ /* wb0 = wh1*y^k */ - clearsfbm(bound,n,wb0); - DEG(wb0) = bound; + clearbm(dx,wb0); cpyum(wh1,COEF(wb0)[k]); - /* wb1 = gk*wb0 mod y^bound */ - clearsfbm(bound,n,wb1); - mulsfbm(bound,gk,wb0,wb1); + /* wb1 = gk*wb0 mod y^(dy+1) */ + clearbm(dx,wb1); + mulsfbm(gk,wb0,wb1); /* fk += wb1 */ - addtosfbm(bound,wb1,fk); + addtosfbm(wb1,fk); /* wb0 = wg1*y^k */ - clearsfbm(bound,n,wb0); - DEG(wb0) = bound; + clearbm(dx,wb0); cpyum(wg1,COEF(wb0)[k]); - /* wb1 = hk*wb0 mod y^bound */ - clearsfbm(bound,n,wb1); - mulsfbm(bound,hk,wb0,wb1); + /* wb1 = hk*wb0 mod y^(dy+1) */ + clearbm(dx,wb1); + mulsfbm(hk,wb0,wb1); /* fk += wb1 */ - addtosfbm(bound,wb1,fk); + addtosfbm(wb1,fk); - /* fk += wg1*wh1*y^(2*k) mod y^bound */ - if ( 2*k < bound ) { + /* fk += wg1*wh1*y^(2*k) mod y^(dy+1) */ + if ( 2*k <= dy ) { mulsfum(wg1,wh1,wt); addsfum(COEF(fk)[2*k],wt,ws); cpyum(ws,COEF(fk)[2*k]); } @@ -933,111 +829,149 @@ BM *gp; cpyum(wg1,COEF(gk)[k]); cpyum(wh1,COEF(hk)[k]); } - fprintf(stderr,"\n"); + if ( debug_sfbfctr) + fprintf(stderr,"\n"); *gp = gk; - DEG(f) = bound; - for ( i = 0; i < bound; i++ ) + DEG(f) = dy; + for ( i = 0; i <= dy; i++ ) cpyum(COEF(hk)[i],COEF(f)[i]); } -void ptosflum(bound,f,fl) -int bound; -P f; -LUM fl; +/* a0*g+b0*h = 1 mod y -> a*g+b*h = 1 mod y^(dy+1) */ + +void sfexgcd_by_hensel(BM g,BM h,int dy,BM *ap,BM *bp) { - DCP dc; - int **pp; - int d; - UM t; + int i,k; + int dx; + UM wt,wa,wb,q,w1,w2,ws; + UM wc,wd,we,wz,wa1,wb1; + BM wz0,wz1; + int dg,dh; + BM a,b,c; - t = UMALLOC(bound); - for ( dc = DC(f), pp = COEF(fl); dc; dc = NEXT(dc) ) { - d = QTOS(DEG(dc)); - ptosfum(COEF(dc),t); - bcopy(t->c,pp[d],(t->d+1)*sizeof(int)); + dg = degbm(g); + dh = degbm(h); + dx = dg+dh; + + a = BMALLOC(dh,dy); + b = BMALLOC(dg,dy); + /* c holds a*g+b*h-1 */ + c = BMALLOC(dg+dh,dy); + + W_BMALLOC(dx,dy,wz0); W_BMALLOC(dx,dy,wz1); + + wt = W_UMALLOC(dx); ws = W_UMALLOC(dx); q = W_UMALLOC(2*dx); + wa1 = W_UMALLOC(2*dx); wb1 = W_UMALLOC(2*dx); + wc = W_UMALLOC(2*dx); wd = W_UMALLOC(2*dx); + we = W_UMALLOC(2*dx); wz = W_UMALLOC(2*dx); + + /* compute wa,wb s.t. wa*g0+wb*h0 = 1 mod y */ + w1 = W_UMALLOC(dg); cpyum(COEF(g)[0],w1); + w2 = W_UMALLOC(dh); cpyum(COEF(h)[0],w2); + wa = W_UMALLOC(2*dx); wb = W_UMALLOC(2*dx); /* XXX */ + eucsfum(w1,w2,wa,wb); + cpyum(wa,COEF(a)[0]); cpyum(wb,COEF(b)[0]); + + /* initialize c to a*g+b*h-1 */ + mulsfbm(a,g,c); mulsfbm(b,h,wz0); addtosfbm(wz0,c); + COEF(COEF(c)[0])[0] = 0; + + if ( debug_sfbfctr) + fprintf(stderr,"dy=%d\n",dy); + for ( k = 1; k <= dy; k++ ) { + if ( debug_sfbfctr) + fprintf(stderr,"."); + + /* at this point, a*g+b*h = 1 mod y^k, c = a*g+b*h-1 */ + + /* wt = -((a*g+b*h-1)/y^k) */ + cpyum(COEF(c)[k],wt); + for ( i = DEG(wt); i >= 0; i-- ) + COEF(wt)[i] = _chsgnsf(COEF(wt)[i]); + + /* compute wa1,wb1 s.t. wa1*g0+wb1*h0 = wt */ + mulsfum(wa,wt,wa1); DEG(wa1) = divsfum(wa1,COEF(h)[0],q); + mulsfum(wa1,COEF(g)[0],wc); subsfum(wt,wc,wd); + DEG(wd) = divsfum(wd,COEF(h)[0],wb1); + + /* c += ((wa1*g+wb1*h)*y^k mod y^(dy+1) */ + /* wz0 = wa1*y^k */ + clearbm(dx,wz0); + cpyum(wa1,COEF(wz0)[k]); + + /* wz1 = wz0*g mod y^(dy+1) */ + clearbm(dx,wz1); + mulsfbm(g,wz0,wz1); + /* c += wz1 */ + addtosfbm(wz1,c); + + /* wz0 = wb1*y^k */ + clearbm(dx,wz0); + cpyum(wb1,COEF(wz0)[k]); + + /* wz1 = wz0*h mod y^(dy+1) */ + clearbm(dx,wz1); + mulsfbm(h,wz0,wz1); + /* c += wz1 */ + addtosfbm(wz1,c); + + /* a += wa1*y^k, b += wb1*y^k */ + cpyum(wa1,COEF(a)[k]); + cpyum(wb1,COEF(b)[k]); } + if ( debug_sfbfctr) + fprintf(stderr,"\n"); + DEG(a) = dy; + DEG(b) = dy; + *ap = a; + *bp = b; } /* fl->c[i] = coef_y(f,i) */ -void ptosfbm(bound,f,fl) -int bound; -P f; -BM fl; +void ptosfbm(int dy,P f,BM fl) { DCP dc; - int d,i,n; + int d,i,dx; UM t; - n = QTOS(DEG(DC(f))); - clearsfbm(bound,n,fl); - DEG(fl) = bound; - t = UMALLOC(bound); + dx = QTOS(DEG(DC(f))); + if ( DEG(fl) < dy ) + error("ptosfbm : invalid input"); + DEG(fl) = dy; + clearbm(dx,fl); + t = UMALLOC(dy); for ( dc = DC(f); dc; dc = NEXT(dc) ) { d = QTOS(DEG(dc)); ptosfum(COEF(dc),t); for ( i = 0; i <= DEG(t); i++ ) COEF(COEF(fl)[i])[d] = COEF(t)[i]; } - for ( i = 0; i < bound; i++ ) - degum(COEF(fl)[i],n); + for ( i = 0; i <= dy; i++ ) + degum(COEF(fl)[i],dx); } /* x : main variable */ -void sflumtop(bound,fl,x,y,fp) -int bound; -LUM fl; -V x,y; -P *fp; +void sfbmtop(BM f,V x,V y,P *fp) { - int i,j,n; - int **c; - UM w; - int *coef; - DCP dc,dct; - - n = fl->d; - c = fl->c; - w = W_UMALLOC(bound); - for ( i = 0, dc = 0; i <= n; i++ ) { - coef = c[i]; - for ( j = bound-1; j >= 0 && coef[j] == 0; j-- ); - if ( j < 0 ) - continue; - DEG(w) = j; bcopy(coef,COEF(w),(j+1)*sizeof(int)); - NEWDC(dct); STOQ(i,DEG(dct)); - sfumtop(y,w,&COEF(dct)); NEXT(dct) = dc; dc = dct; - } - MKP(x,dc,*fp); -} - -/* x : main variable */ - -void sfbmtop(bound,f,x,y,fp) -int bound; -BM f; -V x,y; -P *fp; -{ UM *c; - int i,j,d,a; + int i,j,d,a,dy; GFS b; DCP dc0,dc,dct; + dy = DEG(f); c = COEF(f); - d = DEG(c[0]); - for ( i = 1; i < bound; i++ ) - d = MAX(DEG(c[i]),d); + d = degbm(f); dc0 = 0; for ( i = 0; i <= d; i++ ) { dc = 0; - for ( j = 0; j < bound; j++ ) { + for ( j = 0; j <= dy; j++ ) { if ( DEG(c[j]) >= i && (a = COEF(c[j])[i]) ) { NEWDC(dct); STOQ(j,DEG(dct)); - MKGFS(IFTOF(a),b); + iftogfs(a,&b); COEF(dct) = (P)b; NEXT(dct) = dc; dc = dct; @@ -1057,67 +991,195 @@ P *fp; *fp = 0; } +void sfsqfr(P f,DCP *dcp) +{ + Obj obj; + DCP dc; + VL vl; + + simp_ff((Obj)f,&obj); f = (P)obj; + clctv(CO,f,&vl); + if ( !vl ) { + /* f is a const */ + NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0; *dcp = dc; + } else if ( !NEXT(vl) ) + sfusqfr(f,dcp); + else + sqfrsf(f,dcp); +} + +void sfusqfr(P f,DCP *dcp) +{ + DCP dc,dct; + struct oDUM *udc; + V x; + P lc; + int n,i; + UM mf; + + x = VR(f); + n = getdeg(x,f); + mf = W_UMALLOC(n); + ptosfum(f,mf); + lc = COEF(DC(f)); + if ( !_isonesf(COEF(mf)[n]) ) { + monicsfum(mf); + } + W_CALLOC(n+1,struct oDUM,udc); + gensqfrsfum(mf,udc); + for ( i = 0, dc = 0; udc[i].f; i++ ) { + NEWDC(dct); STOQ(udc[i].n,DEG(dct)); + sfumtop(x,udc[i].f,&COEF(dct)); + NEXT(dct) = dc; dc = dct; + } + NEWDC(dct); DEG(dct) = ONE; COEF(dct) = (P)lc; NEXT(dct) = dc; + *dcp = dct; +} + +#if 0 +void sfbsqfrmain(P f,V x,V y,DCP *dcp) +{ + /* XXX*/ +} + +/* f is bivariate */ + +void sfbsqfr(P f,V x,V y,DCP *dcp) +{ + P t,rf,cx,cy; + VL vl,rvl; + DCP dcx,dcy,dct,dc; + struct oVL vl0,vl1; + + /* cy(y) = cont(f,x), f /= cy */ + cont_pp_sfp(vl,f,&cy,&t); f = t; + /* rvl = [y,x] */ + reordvar(vl,y,&rvl); reorderp(rvl,vl,f,&rf); + /* cx(x) = cont(rf,y), Rf /= cy */ + cont_pp_sfp(rvl,rf,&cx,&t); rf = t; + reorderp(vl,rvl,rf,&f); + + /* f -> cx*cy*f */ + sfsqfr(cx,&dcx); dcx = NEXT(dcx); + sfsqfr(cy,&dcy); dcy = NEXT(dcy); + if ( dcx ) { + for ( dct = dcx; NEXT(dct); dct = NEXT(dct) ); + NEXT(dct) = dcy; + } else + dcx = dcy; + if ( OID(f) == O_N ) + *dcp = dcx; + else { + /* f must be bivariate */ + sfbsqfrmain(f,x,y,&dc); + if ( dcx ) { + for ( dct = dcx; NEXT(dct); dct = NEXT(dct) ); + NEXT(dct) = dc; + } else + dcx = dc; + *dcp = dcx; + } +} +#endif + void sfdtest(P,ML,V,V,DCP *); -void sfbfctr(f,x,y,dcp) -P f; -V x,y; -DCP *dcp; +/* if degbound >= 0 find factor s.t. deg_x(factor) <= degbound */ + +void sfbfctr(P f,V x,V y,int degbound,DCP *dcp) { ML list; P sf; GFS ev; DCP dc,dct; BM fl; - int n,bound; + int dx,dy; /* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */ - sfhensel(5,f,x,&ev,&sf,&list); + sfhensel(5,f,x,y,degbound,&ev,&sf,&list); + if ( list->n == 0 ) + error("sfbfctr : short of evaluation points"); + else if ( list->n == 1 ) { + /* f is irreducible */ + NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0; + *dcp = dc; + return; + } sfdtest(sf,list,x,y,&dc); - n = getdeg(x,sf); - bound = list->bound; - W_BMALLOC(n,bound,fl); - for ( dct = dc; dct; dct = NEXT(dct) ) { - ptosfbm(bound,COEF(dct),fl); - shiftsfbm(bound,fl,_chsgnsf(FTOIF(CONT(ev)))); - sfbmtop(bound,fl,x,y,&COEF(dct)); + if ( ev ) { + dx = getdeg(x,sf); + dy = getdeg(y,sf); + W_BMALLOC(dx,dy,fl); + for ( dct = dc; dct; dct = NEXT(dct) ) { + ptosfbm(dy,COEF(dct),fl); + shiftsfbm(fl,_chsgnsf(FTOIF(CONT(ev)))); + sfbmtop(fl,x,y,&COEF(dct)); + } } *dcp = dc; } -/* f = f(x,y) = list->c[0]*list->c[1]*... mod y^list->bound */ +/* returns shifted f, shifted factors and the eval pt */ -void sfdtest(f,list,x,y,dcp) -P f; -ML list; -V x,y; -DCP *dcp; +void sfbfctr_shift(P f,V x,V y,int degbound,GFS *evp,P *sfp,DCP *dcp) { - int n,np,bound; - int i,j,k; + ML list; + P sf; + GFS ev; + DCP dc,dct; + int dx,dy; + + /* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */ + sfhensel(5,f,x,y,degbound,&ev,&sf,&list); + if ( list->n == 0 ) + error("sfbfctr_shift : short of evaluation points"); + else if ( list->n == 1 ) { + /* f is irreducible */ + NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0; + *evp = 0; + *sfp = f; + *dcp = dc; + } else { + sfdtest(sf,list,x,y,dcp); + *evp = ev; + *sfp = sf; + } +} + +/* f = f(x,y) = list->c[0]*list->c[1]*... mod y^(list->bound+1) */ + +void sfdtest(P f,ML list,V x,V y,DCP *dcp) +{ + int np,dx,dy; + int i,j,k,bound; int *win; P g,lcg,factor,cofactor,lcyx; - P t,csum; + P csum; DCP dcf,dcf0,dc; BM *c; BM lcy; - UM lcg0; + UM lcg0,lcy0,w; + UM *d1c; ML wlist; struct oVL vl1,vl0; VL vl; - int z; + int z,dt,dtok; /* vl = [x,y] */ vl0.v = x; vl0.next = &vl1; vl1.v = y; vl1.next = 0; vl = &vl0; /* setup various structures and arrays */ - n = UDEG(f); np = list->n; bound = list->bound; win = W_ALLOC(np+1); - wlist = W_MLALLOC(np); wlist->n = list->n; - wlist->bound = list->bound; + dx = getdeg(x,f); + dy = getdeg(y,f); + np = list->n; + win = W_ALLOC(np+1); + wlist = W_MLALLOC(np); + wlist->n = list->n; + bound = wlist->bound = list->bound; c = (BM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(BM)*np)); - lcg0 = W_UMALLOC(2*bound); + lcg0 = W_UMALLOC(2*dy); /* initialize g by f */ g = f; @@ -1132,15 +1194,40 @@ DCP *dcp; sfcsump(vl,lcg,&csum); /* initialize lcy by LC(f) */ - W_BMALLOC(0,bound,lcy); + W_BMALLOC(0,dy,lcy); NEWDC(dc); COEF(dc) = COEF(DC(g)); DEG(dc) = 0; NEWP(lcyx); VR(lcyx) = x; DC(lcyx) = dc; - ptosfbm(bound,lcyx,lcy); + ptosfbm(dy,lcyx,lcy); - fprintf(stderr,"np = %d\n",np); + /* initialize lcy0 by LC(f) */ + lcy0 = W_UMALLOC(bound); + ptosfum(COEF(DC(g)),lcy0); + + /* ((d-1 coefs)*lcy0 */ + d1c = (UM *)W_ALLOC(np*sizeof(UM)); + w = W_UMALLOC(2*bound); + for ( i = 1; i < np; i++ ) { + extractcoefbm(c[i],degbm(c[i])-1,w); + d1c[i] = W_UMALLOC(2*bound); + mulsfum(w,lcy0,d1c[i]); + /* d1c[i] = d1c[i] mod y^(bound+1) */ + if ( DEG(d1c[i]) > bound ) { + for ( j = DEG(d1c[i]); j > bound; j-- ) + COEF(d1c[i])[j] = 0; + degum(d1c[i],bound); + } + } + + if ( debug_sfbfctr) + fprintf(stderr,"np = %d\n",np); + dtok = 0; for ( g = f, k = 1, dcf = dcf0 = 0, win[0] = 1, --np, z = 0; ; z++ ) { - if ( !(z % 1000) ) fprintf(stderr,"."); - if ( sfdtestmain(vl,lcg,lcg0,lcy,csum,wlist,k,win,&factor,&cofactor) ) { + if ( debug_sfbfctr && !(z % 1000) ) fprintf(stderr,"."); + dt = sfdegtest(dy,bound,d1c,k,win); + if ( dt ) + dtok++; + if ( dt && sfdtestmain(vl,lcg,lcg0,lcy,csum,wlist, + k,win,&factor,&cofactor) ) { NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor; g = cofactor; @@ -1153,9 +1240,13 @@ DCP *dcp; /* update csum */ sfcsump(vl,lcg,&csum); + /* update dy */ + dy = getdeg(y,g); + /* update lcy */ - clearsfbm(bound,0,lcy); - COEF(dc) = COEF(DC(g)); ptosfbm(bound,lcyx,lcy); + clearbm(0,lcy); + COEF(dc) = COEF(DC(g)); + ptosfbm(dy,lcyx,lcy); for ( i = 0; i < k - 1; i++ ) for ( j = win[i] + 1; j < win[i + 1]; j++ ) @@ -1173,6 +1264,22 @@ DCP *dcp; else for ( i = 1; i < k; i++ ) win[i] = win[0] + i; + + + /* update lcy0 */ + ptosfum(COEF(DC(g)),lcy0); + + /* update d-1 coeffs */ + for ( i = 1; i <= np; i++ ) { + extractcoefbm(c[i],degbm(c[i])-1,w); + mulsfum(w,lcy0,d1c[i]); + /* d1c[i] = d1c[1] mod y^(bound+1) */ + if ( DEG(d1c[i]) > bound ) { + for ( j = DEG(d1c[i]); j > bound; j-- ) + COEF(d1c[i])[j] = 0; + degum(d1c[i],bound); + } + } } else if ( !ncombi(1,np,k,win) ) if ( k == np ) break; @@ -1180,24 +1287,50 @@ DCP *dcp; for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1; } - fprintf(stderr,"\n"); + if ( debug_sfbfctr ) + fprintf(stderr,"total %d, omitted by degtest %d\n",z,z-dtok); NEXTDC(dcf0,dcf); COEF(dcf) = g; DEG(dcf) = ONE; NEXT(dcf) = 0; *dcp = dcf0; } +void extractcoefbm(BM f,int dx,UM r) +{ + int j; + UM fj; + + for ( j = DEG(f); j >= 0; j-- ) { + fj = COEF(f)[j]; + if ( fj && DEG(fj) >= dx ) { + COEF(r)[j] = COEF(fj)[dx]; + } else + COEF(r)[j] = 0; + } + degum(r,DEG(f)); +} + +/* deg_y(prod mod y^(bound+1)) <= dy ? */ + +int sfdegtest(int dy,int bound,UM *d1c,int k,int *in) +{ + int i,j; + UM w,w1,wt; + BM t; + + w = W_UMALLOC(bound); + w1 = W_UMALLOC(bound); + clearum(w,bound); + for ( i = 0; i < k; i++ ) { + addsfum(w,d1c[in[i]],w1); wt = w; w = w1; w1 = wt; + } + return DEG(w) <= dy ? 1 : 0; +} + /* lcy = LC(g), lcg = lcy*g, lcg0 = const part of lcg */ -int sfdtestmain(vl,lcg,lcg0,lcy,csum,list,k,in,fp,cofp) -VL vl; -P lcg; -UM lcg0; -BM lcy; -P csum; -ML list; -int k; -int *in; -P *fp,*cofp; + +int sfdtestmain(VL vl,P lcg,UM lcg0,BM lcy,P csum,ML list, + int k,int *in,P *fp,P *cofp) { - P fmul,csumg,q; + P fmul,csumg,q,cont; V x,y; x = vl->v; @@ -1213,16 +1346,14 @@ P *fp,*cofp; } } if ( divtp_by_sfbm(vl,lcg,fmul,&q) ) { - pp_sfp(vl,fmul,fp); - pp_sfp(vl,q,cofp); + cont_pp_sfp(vl,fmul,&cont,fp); + cont_pp_sfp(vl,q,&cont,cofp); return 1; } else return 0; } -void const_term(f,c) -P f; -UM c; +void const_term(P f,UM c) { DCP dc; @@ -1233,49 +1364,41 @@ UM c; DEG(c) = -1; } -void const_term_sfbm(f,bound,c) -BM f; -int bound; -UM c; +void const_term_sfbm(BM f,UM c) { - int i; + int i,dy; - for ( i = 0; i < bound; i++ ) + dy = DEG(f); + for ( i = 0; i <= dy; i++ ) if ( DEG(COEF(f)[i]) >= 0 ) COEF(c)[i] = COEF(COEF(f)[i])[0]; else COEF(c)[i] = 0; - degum(c,bound-1); + degum(c,dy); } /* lcy*(product of const part) | lcg0 ? */ -int sfctest(lcg0,lcy,list,k,in) -UM lcg0; -BM lcy; -ML list; -int k; -int *in; +int sfctest(UM lcg0,BM lcy,ML list,int k,int *in) { - DCP dc; - int bound,i,dr; + int dy,i,dr; UM t,s,u,w; BM *l; - bound = list->bound; - t = W_UMALLOC(2*bound); - s = W_UMALLOC(2*bound); - u = W_UMALLOC(2*bound); - const_term_sfbm(lcy,bound,t); + dy = list->bound; + t = W_UMALLOC(2*dy); + s = W_UMALLOC(2*dy); + u = W_UMALLOC(2*dy); + const_term_sfbm(lcy,t); if ( DEG(t) < 0 ) return 1; l = (BM *)list->c; for ( i = 0; i < k; i++ ) { - const_term_sfbm(l[in[i]],bound,s); + const_term_sfbm(l[in[i]],s); mulsfum(t,s,u); - if ( DEG(u) >= bound ) - degum(u,bound-1); + if ( DEG(u) > dy ) + degum(u,dy); w = t; t = u; u = w; } cpyum(lcg0,s); @@ -1288,36 +1411,26 @@ int *in; /* main var of f is x */ -void mulsfbmarray(n,lcy,list,k,in,x,y,g) -int n; -BM lcy; -ML list; -int k; -int *in; -V x,y; -P *g; +void mulsfbmarray(int dx,BM lcy,ML list,int k,int *in,V x,V y,P *g) { - int bound,i; - BM wb0,wb1,t,lcbm; + int dy,i; + BM wb0,wb1,t; BM *l; - bound = list->bound; - W_BMALLOC(n,bound,wb0); W_BMALLOC(n,bound,wb1); + dy = list->bound; + W_BMALLOC(dx,dy,wb0); W_BMALLOC(dx,dy,wb1); l = (BM *)list->c; - clearsfbm(bound,n,wb0); - mulsfbm(bound,lcy,l[in[0]],wb0); + clearbm(dx,wb0); + mulsfbm(lcy,l[in[0]],wb0); for ( i = 1; i < k; i++ ) { - clearsfbm(bound,n,wb1); - mulsfbm(bound,l[in[i]],wb0,wb1); + clearbm(dx,wb1); + mulsfbm(l[in[i]],wb0,wb1); t = wb0; wb0 = wb1; wb1 = t; } - sfbmtop(bound,wb0,x,y,g); + sfbmtop(wb0,x,y,g); } -void sfcsump(vl,f,s) -VL vl; -P f; -P *s; +void sfcsump(VL vl,P f,P *s) { P t,u; DCP dc; @@ -1330,23 +1443,22 @@ P *s; /* *fp = primitive part of f w.r.t. x */ -void pp_sfp(vl,f,fp) -VL vl; -P f; -P *fp; +void cont_pp_sfp(VL vl,P f,P *cp,P *fp) { V x,y; int d; UM t,s,gcd; DCP dc; - P dvr; + GFS g; x = vl->v; y = vl->next->v; d = getdeg(y,f); - if ( d == 0 ) + if ( d == 0 ) { + itogfs(1,&g); + *cp = (P)g; *fp = f; /* XXX */ - else { + } else { t = W_UMALLOC(2*d); s = W_UMALLOC(2*d); gcd = W_UMALLOC(2*d); @@ -1357,15 +1469,12 @@ P *fp; gcdsfum(gcd,t,s); cpyum(s,gcd); } - sfumtop(y,gcd,&dvr); - divsp(vl,f,dvr,fp); + sfumtop(y,gcd,cp); + divsp(vl,f,*cp,fp); } } -int divtp_by_sfbm(vl,f,g,qp) -VL vl; -P f,g; -P *qp; +int divtp_by_sfbm(VL vl,P f,P g,P *qp) { V x,y; int fx,fy,gx,gy; @@ -1380,9 +1489,9 @@ P *qp; if ( fx < gx || fy < gy ) return 0; - W_BMALLOC(fx,fy+1,fl); ptosfbm(fy+1,f,fl); cf = COEF(fl); - W_BMALLOC(gx,gy+1,gl); ptosfbm(gy+1,g,gl); cg = COEF(gl); - W_BMALLOC(fx-gx,fy-gy+1,ql); cq = COEF(ql); + W_BMALLOC(fx,fy,fl); ptosfbm(fy,f,fl); cf = COEF(fl); + W_BMALLOC(gx,gy,gl); ptosfbm(gy,g,gl); cg = COEF(gl); + W_BMALLOC(fx-gx,fy-gy,ql); cq = COEF(ql); hg = cg[gy]; q = W_UMALLOC(fx); t = W_UMALLOC(fx); s = W_UMALLOC(fx); @@ -1405,5 +1514,56 @@ P *qp; for ( j = gy-1; j >= 0 && DEG(cf[j]) < 0; j-- ); if ( j >= 0 ) return 0; - sfbmtop(DEG(ql),ql,x,y,qp); + sfbmtop(ql,x,y,qp); + return 1; +} + +/* XXX generate an irreducible poly of degree n */ + +extern int current_gfs_q1; +extern int *current_gfs_ntoi; + +void generate_defpoly_sfum(int n,UM *dp) +{ + UM r,dr,t,g; + UM *f; + int *c,*w; + int max,i,j; + + *dp = r = UMALLOC(n); + DEG(r) = n; + c = COEF(r); + c[n] = _onesf(); + max = current_gfs_q1; + w = (int *)ALLOCA(n*sizeof(int)); + bzero(w,n*sizeof(int)); + + dr = W_UMALLOC(n); t = W_UMALLOC(n); g = W_UMALLOC(n); + f = (UM *)ALLOCA((n+1)*sizeof(UM)); + while ( 1 ) { + for ( i = 0; i < n && w[i] == max; i++ ); + if ( i == n ) { + /* XXX cannot happen */ + error("generate_defpoly_sfum : cannot happen"); + } + for ( j = 0; j < i; j++ ) + w[j] = 0; + w[i]++; + if ( !current_gfs_ntoi ) + for ( i = 0; i < n; i++ ) + c[i] = w[i]?FTOIF(w[i]):0; + else + for ( i = 0; i < n; i++ ) + c[i] = w[i]?FTOIF(w[i]-1):0; + if ( !c[0] ) + continue; + diffsfum(r,dr); cpyum(r,t); gcdsfum(t,dr,g); + if ( DEG(g) > 0 ) + continue; + + czsfum(r,f); + for ( i = 0; f[i]; i++ ); + if ( i == 1 ) + return; + } }