=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/Hgfs.c,v retrieving revision 1.15 retrieving revision 1.28 diff -u -p -r1.15 -r1.28 --- OpenXM_contrib2/asir2000/engine/Hgfs.c 2001/08/02 03:59:15 1.15 +++ OpenXM_contrib2/asir2000/engine/Hgfs.c 2002/11/22 07:32:10 1.28 @@ -1,29 +1,12 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.14 2001/07/03 01:41:25 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.27 2002/11/01 05:43:35 noro Exp $ */ #include "ca.h" +#include "inline.h" -struct p_pair { - UM p0; - UM p1; - struct p_pair *next; -}; +void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1); +void extractcoefbm(BM f,int dx,UM r); -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(BM,V,V,P *); -void cont_pp_sfp(VL,P,P *,P *); -void sfcsump(VL,P,P *); -void mulsfbmarray(int,BM,ML,int,int *,V,V,P *); -void const_term(P,UM); - -void sfbsqfr(P,V,V,DCP *); -void sfusqfr(P,DCP *); - -int comp_dum(a,b) -DUM a,b; +int comp_dum(DUM a,DUM b) { if ( DEG(a->f) > DEG(b->f) ) return -1; @@ -33,21 +16,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); @@ -92,9 +75,7 @@ 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; @@ -147,9 +128,7 @@ struct oDUM *dc; } } -void randsfum(d,p) -int d; -UM p; +void randsfum(int d,UM p) { int i; @@ -159,9 +138,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; @@ -189,9 +166,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; @@ -213,9 +188,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; @@ -235,10 +208,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; @@ -259,10 +229,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; @@ -294,11 +261,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; @@ -327,8 +290,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; @@ -382,10 +344,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; @@ -431,10 +390,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; @@ -469,13 +430,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; @@ -492,12 +449,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)); @@ -507,14 +462,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]); @@ -553,32 +506,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 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; + 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; @@ -587,22 +534,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)); @@ -611,9 +585,14 @@ ML *listp; ptosfum(dc->c,gm[i]); } + /* 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,dy); - ptosfbm(dy,f,fl); + fl = BMALLOC(dx,bound); + ptosfbm(bound,f,fl); if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev))); /* sf = f(x+ev) */ @@ -625,7 +604,7 @@ ML *listp; hm = W_UMALLOC(dx); q = W_UMALLOC(dx); - rlist = MLALLOC(fn); rlist->n = fn; rlist->bound = dy; + rlist = MLALLOC(fn); rlist->n = fn; rlist->bound = bound; fprintf(asir_out,"%d candidates\n",fn); init_eg(&eg_hensel); for ( i = 0; i < fn-1; i++ ) { @@ -635,9 +614,9 @@ ML *listp; get_eg(&tmp0); /* fl = gm[i]*hm mod y */ divsfum(fm,gm[i],hm); - /* fl is replaced by the cofactor of gk mod y^dy */ + /* fl is replaced by the cofactor of gk mod y^bound */ /* rlist->c[i] = gk */ - sfhenmain2(fl,gm[i],hm,dy,(BM *)&rlist->c[i]); + sfhenmain2(fl,gm[i],hm,bound,(BM *)&rlist->c[i]); cpyum(hm,fm); get_eg(&tmp1); add_eg(&eg_hensel_t,&tmp0,&tmp1); add_eg(&eg_hensel,&tmp0,&tmp1); @@ -662,41 +641,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 +679,7 @@ DCP *dcp; } } } - if ( index == q1 ) + if ( index == q ) return 0; else if ( fn == 1 ) return 1; @@ -716,9 +690,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; @@ -747,17 +719,11 @@ ML blist,*clistp; /* f = g0*h0 mod y -> f = gk*hk mod y^(dy+1), f is replaced by hk */ -void sfhenmain2(f,g0,h0,dy,gp) -BM f; -UM g0,h0; -int dy; -BM *gp; +void sfhenmain2(BM f,UM g0,UM h0,int dy,BM *gp) { - int i,j,k,l; - int *px,*py; - int **pp,**pp1; - int dx,np,dr,tmp; - UM wt,wa,wb,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 dg,dh; @@ -856,12 +822,96 @@ BM *gp; cpyum(COEF(hk)[i],COEF(f)[i]); } +/* 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) +{ + 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; + + 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; + + fprintf(stderr,"dy=%d\n",dy); + for ( k = 1; k <= dy; k++ ) { + 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]); + } + fprintf(stderr,"\n"); + DEG(a) = dy; + DEG(b) = dy; + *ap = a; + *bp = b; +} + /* fl->c[i] = coef_y(f,i) */ -void ptosfbm(dy,f,fl) -int dy; -P f; -BM fl; +void ptosfbm(int dy,P f,BM fl) { DCP dc; int d,i,dx; @@ -885,10 +935,7 @@ BM fl; /* x : main variable */ -void sfbmtop(f,x,y,fp) -BM f; -V x,y; -P *fp; +void sfbmtop(BM f,V x,V y,P *fp) { UM *c; int i,j,d,a,dy; @@ -906,7 +953,7 @@ P *fp; 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; @@ -926,31 +973,24 @@ P *fp; *fp = 0; } -void sfsqfr(f,dcp) -P f; -DCP *dcp; +void sfsqfr(P f,DCP *dcp) { - P g; - V x; + Obj obj; DCP dc; VL vl; - simp_ff(f,&g); f = g; + 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 if ( !NEXT(NEXT(vl)) ) - sfbsqfr(f,x,NEXT(vl)->v,dcp); else - error("sfsqfr : not implemented yet"); + sqfrsf(f,dcp); } -void sfusqfr(f,dcp) -P f; -DCP *dcp; +void sfusqfr(P f,DCP *dcp) { DCP dc,dct; struct oDUM *udc; @@ -978,19 +1018,22 @@ DCP *dcp; *dcp = dct; } -void sfbsqfr(f,x,y,dcp) -P f; -V x,y; -DCP *dcp; +#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; + DCP dcx,dcy,dct,dc; struct oVL vl0,vl1; - /* vl = [x,y] */ - vl0.v = x; vl0.next = &vl1; vl1.v = y; vl1.next = 0; vl = &vl0; - /* cy(y) = cont(f,x), f /= cx */ + /* 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); @@ -999,16 +1042,33 @@ DCP *dcp; reorderp(vl,rvl,rf,&f); /* f -> cx*cy*f */ - sfusqfr(cx,&dcx); - sfusqfr(cy,&dcy); + 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; @@ -1018,7 +1078,7 @@ DCP *dcp; 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 ) { @@ -1041,27 +1101,51 @@ DCP *dcp; *dcp = dc; } +/* returns shifted f, shifted factors and the eval pt */ + +void sfbfctr_shift(P f,V x,V y,int degbound,GFS *evp,P *sfp,DCP *dcp) +{ + 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(f,list,x,y,dcp) -P f; -ML list; -V x,y; -DCP *dcp; +void sfdtest(P f,ML list,V x,V y,DCP *dcp) { int np,dx,dy; - int i,j,k; + 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; @@ -1073,7 +1157,7 @@ DCP *dcp; win = W_ALLOC(np+1); wlist = W_MLALLOC(np); wlist->n = list->n; - wlist->bound = dy; + bound = wlist->bound = list->bound; c = (BM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(BM)*np)); @@ -1097,10 +1181,34 @@ DCP *dcp; NEWP(lcyx); VR(lcyx) = x; DC(lcyx) = dc; ptosfbm(dy,lcyx,lcy); + /* 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); + } + } + 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) ) { + 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; @@ -1113,6 +1221,9 @@ DCP *dcp; /* update csum */ sfcsump(vl,lcg,&csum); + /* update dy */ + dy = getdeg(y,g); + /* update lcy */ clearbm(0,lcy); COEF(dc) = COEF(DC(g)); @@ -1134,6 +1245,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; @@ -1141,22 +1268,47 @@ DCP *dcp; for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1; } - fprintf(stderr,"\n"); + 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,cont; V x,y; @@ -1181,9 +1333,7 @@ P *fp,*cofp; return 0; } -void const_term(f,c) -P f; -UM c; +void const_term(P f,UM c) { DCP dc; @@ -1194,9 +1344,7 @@ UM c; DEG(c) = -1; } -void const_term_sfbm(f,c) -BM f; -UM c; +void const_term_sfbm(BM f,UM c) { int i,dy; @@ -1211,14 +1359,8 @@ UM c; /* 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 dy,i,dr; UM t,s,u,w; BM *l; @@ -1249,17 +1391,10 @@ int *in; /* main var of f is x */ -void mulsfbmarray(dx,lcy,list,k,in,x,y,g) -int dx; -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 dy,i; - BM wb0,wb1,t,lcbm; + BM wb0,wb1,t; BM *l; dy = list->bound; @@ -1275,10 +1410,7 @@ P *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; @@ -1291,10 +1423,7 @@ P *s; /* *fp = primitive part of f w.r.t. x */ -void cont_pp_sfp(vl,f,cp,fp) -VL vl; -P f; -P *cp,*fp; +void cont_pp_sfp(VL vl,P f,P *cp,P *fp) { V x,y; int d; @@ -1306,7 +1435,7 @@ P *cp,*fp; y = vl->next->v; d = getdeg(y,f); if ( d == 0 ) { - MKGFS(0,g); + itogfs(1,&g); *cp = (P)g; *fp = f; /* XXX */ } else { @@ -1325,10 +1454,7 @@ P *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; @@ -1369,4 +1495,55 @@ P *qp; if ( j >= 0 ) return 0; 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; + } }