=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/dp-supp.c,v retrieving revision 1.18 retrieving revision 1.61 diff -u -p -r1.18 -r1.61 --- OpenXM_contrib2/asir2000/builtin/dp-supp.c 2001/09/17 10:32:40 1.18 +++ OpenXM_contrib2/asir2000/builtin/dp-supp.c 2015/08/14 13:51:54 1.61 @@ -45,7 +45,7 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.17 2001/09/17 02:58:27 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.60 2015/08/08 14:19:41 fujimoto Exp $ */ #include "ca.h" #include "base.h" @@ -62,14 +62,41 @@ extern int NoGCD; extern int GenTrace; extern NODE TraceList; +int show_orderspec; + +void print_composite_order_spec(struct order_spec *spec); + /* * content reduction * */ -void dp_ptozp(p,rp) -DP p,*rp; +static NODE RatDenomList; + +void init_denomlist() { + RatDenomList = 0; +} + +void add_denomlist(P f) +{ + NODE n; + + if ( OID(f)==O_P ) { + MKNODE(n,f,RatDenomList); RatDenomList = n; + } +} + +LIST get_denomlist() +{ + LIST l; + + MKLIST(l,RatDenomList); RatDenomList = 0; + return l; +} + +void dp_ptozp(DP p,DP *rp) +{ MP m,mr,mr0; int i,n; Q *w; @@ -95,9 +122,7 @@ DP p,*rp; } } -void dp_ptozp2(p0,p1,hp,rp) -DP p0,p1; -DP *hp,*rp; +void dp_ptozp2(DP p0,DP p1,DP *hp,DP *rp) { DP t,s,h,r; MP m,mr,mr0,m0; @@ -121,11 +146,34 @@ DP *hp,*rp; *hp = h; *rp = r; } -void dp_idiv(p,c,rp) -DP p; -Q c; -DP *rp; +void dp_ptozp3(DP p,Q *dvr,DP *rp) { + MP m,mr,mr0; + int i,n; + Q *w; + P t; + + if ( !p ) { + *rp = 0; *dvr = 0; + }else { + for ( m =BDY(p), n = 0; m; m = NEXT(m), n++ ); + w = (Q *)ALLOCA(n*sizeof(Q)); + for ( m =BDY(p), i = 0; i < n; m = NEXT(m), i++ ) + if ( NUM(m->c) ) + w[i] = (Q)m->c; + else + ptozp(m->c,1,&w[i],&t); + sortbynm(w,n); + qltozl(w,n,dvr); + for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { + NEXTMP(mr0,mr); divsp(CO,m->c,(P)(*dvr),&mr->c); mr->dl = m->dl; + } + NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar; + } +} + +void dp_idiv(DP p,Q c,DP *rp) +{ Q t; N nm,q; int sgn,s; @@ -154,9 +202,7 @@ DP *rp; } } -void dp_mbase(hlist,mbase) -NODE hlist; -NODE *mbase; +void dp_mbase(NODE hlist,NODE *mbase) { DL *dl; DL d; @@ -170,7 +216,8 @@ NODE *mbase; while ( 1 ) { insert_to_node(d,mbase,nvar); for ( i = nvar-1; i >= 0; ) { - d->d[i]++; d->td++; + d->d[i]++; + d->td += MUL_WEIGHT(1,i); for ( j = 0; j < n; j++ ) { if ( _dl_redble(dl[j],d,nvar) ) break; @@ -179,7 +226,7 @@ NODE *mbase; for ( j = nvar-1; j >= i; j-- ) d->d[j] = 0; for ( j = 0, td = 0; j < i; j++ ) - td += d->d[j]; + td += MUL_WEIGHT(d->d[j],j); d->td = td; i--; } else @@ -190,9 +237,7 @@ NODE *mbase; } } -int _dl_redble(d1,d2,nvar) -DL d1,d2; -int nvar; +int _dl_redble(DL d1,DL d2,int nvar) { int i; @@ -207,10 +252,7 @@ int nvar; return 1; } -void insert_to_node(d,n,nvar) -DL d; -NODE *n; -int nvar; +void insert_to_node(DL d,NODE *n,int nvar) { DL d1; MP m; @@ -239,10 +281,7 @@ int nvar; } } -void dp_vtod(c,p,rp) -Q *c; -DP p; -DP *rp; +void dp_vtod(Q *c,DP p,DP *rp) { MP mr0,m,mr; int i; @@ -261,8 +300,7 @@ DP *rp; extern int mpi_mag; extern int PCoeffs; -void dp_ptozp_d(p,rp) -DP p,*rp; +void dp_ptozp_d(DP p,DP *rp) { int i,j,k,l,n,nsep; MP m; @@ -271,7 +309,6 @@ DP p,*rp; VECT c,cs; VECT qi,ri; LIST *qr; - int s,id; Obj dmy; Q d0,d1,gcd,a,u,u1; Q *q,*r; @@ -286,6 +323,8 @@ DP p,*rp; double t_e,t_d,t_d1,t_c; extern int DP_NFStat; extern LIST Dist; + void Pox_rpc(); + void Pox_pop_local(); if ( !p ) *rp = 0; @@ -362,9 +401,7 @@ DP p,*rp; } } -void dp_ptozp2_d(p0,p1,hp,rp) -DP p0,p1; -DP *hp,*rp; +void dp_ptozp2_d(DP p0,DP p1,DP *hp,DP *rp) { DP t,s,h,r; MP m,mr,mr0,m0; @@ -388,9 +425,46 @@ DP *hp,*rp; *hp = h; *rp = r; } -void dp_prim(p,rp) -DP p,*rp; +int have_sf_coef(P p) { + DCP dc; + + if ( !p ) + return 0; + else if ( NUM(p) ) + return NID((Num)p) == N_GFS ? 1 : 0; + else { + for ( dc = DC(p); dc; dc = NEXT(dc) ) + if ( have_sf_coef(COEF(dc)) ) + return 1; + return 0; + } +} + +void head_coef(P p,Num *c) +{ + if ( !p ) + *c = 0; + else if ( NUM(p) ) + *c = (Num)p; + else + head_coef(COEF(DC(p)),c); +} + +void dp_monic_sf(DP p,DP *rp) +{ + Num c; + + if ( !p ) + *rp = 0; + else { + head_coef(BDY(p)->c,&c); + divsdc(CO,p,(P)c,rp); + } +} + +void dp_prim(DP p,DP *rp) +{ P t,g; DP p1; MP m,mr,mr0; @@ -398,10 +472,34 @@ DP p,*rp; P *w; Q *c; Q dvr; + NODE tn; if ( !p ) *rp = 0; - else if ( dp_fcoeffs ) + else if ( dp_fcoeffs == N_GFS ) { + for ( m = BDY(p); m; m = NEXT(m) ) + if ( OID(m->c) == O_N ) { + /* GCD of coeffs = 1 */ + dp_monic_sf(p,rp); + return; + } else break; + /* compute GCD over the finite fieid */ + for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ ); + w = (P *)ALLOCA(n*sizeof(P)); + for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ ) + w[i] = m->c; + gcdsf(CO,w,n,&g); + if ( NUM(g) ) + dp_monic_sf(p,rp); + else { + for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { + NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl; + } + NEXT(mr) = 0; MKDP(p->nv,mr0,p1); p1->sugar = p->sugar; + dp_monic_sf(p1,rp); + } + return; + } else if ( dp_fcoeffs ) *rp = p; else if ( NoGCD ) dp_ptozp(p,rp); @@ -429,19 +527,19 @@ DP p,*rp; NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl; } NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar; + add_denomlist(g); } } } -void heu_nezgcdnpz(vl,pl,m,pr) -VL vl; -P *pl,*pr; -int m; +void heu_nezgcdnpz(VL vl,P *pl,int m,P *pr) { int i,r; P gcd,t,s1,s2,u; Q rq; - + DCP dc; + extern int DP_Print; + while ( 1 ) { for ( i = 0, s1 = 0; i < m; i++ ) { r = random(); UTOQ(r,rq); @@ -452,6 +550,8 @@ int m; mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u; } ezgcdp(vl,s1,s2,&gcd); + if ( DP_Print > 2 ) + { fprintf(asir_out,"(%d)",nmonop(gcd)); fflush(asir_out); } for ( i = 0; i < m; i++ ) { if ( !divtpz(vl,pl[i],gcd,&t) ) break; @@ -462,9 +562,7 @@ int m; *pr = gcd; } -void dp_prim_mod(p,mod,rp) -int mod; -DP p,*rp; +void dp_prim_mod(DP p,int mod,DP *rp) { P t,g; MP m,mr,mr0; @@ -484,18 +582,14 @@ DP p,*rp; } } -void dp_cont(p,rp) -DP p; -Q *rp; +void dp_cont(DP p,Q *rp) { VECT v; dp_dtov(p,&v); igcdv(v,rp); } -void dp_dtov(dp,rp) -DP dp; -VECT *rp; +void dp_dtov(DP dp,VECT *rp) { MP m,t; int i,n; @@ -515,9 +609,7 @@ VECT *rp; * */ -void dp_sp(p1,p2,rp) -DP p1,p2; -DP *rp; +void dp_sp(DP p1,DP p2,DP *rp) { int i,n,td; int *w; @@ -530,7 +622,7 @@ DP *rp; n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; w = (int *)ALLOCA(n*sizeof(int)); for ( i = 0, td = 0; i < n; i++ ) { - w[i] = MAX(d1->d[i],d2->d[i]); td += w[i]; + w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i); } NEWDL(d,n); d->td = td - d1->td; @@ -559,20 +651,18 @@ DP *rp; LIST hist; NODE node; - node = mknode(4,ONE,0,s1,ONE); + node = mknode(4,ONE,NULLP,s1,ONE); MKLIST(hist,node); MKNODE(TraceList,hist,0); - node = mknode(4,ONE,0,0,ONE); + node = mknode(4,ONE,NULLP,NULLP,ONE); chsgnd(s2,(DP *)&ARG2(node)); MKLIST(hist,node); MKNODE(node,hist,TraceList); TraceList = node; } } -void _dp_sp_dup(p1,p2,rp) -DP p1,p2; -DP *rp; +void _dp_sp_dup(DP p1,DP p2,DP *rp) { int i,n,td; int *w; @@ -585,7 +675,7 @@ DP *rp; n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; w = (int *)ALLOCA(n*sizeof(int)); for ( i = 0, td = 0; i < n; i++ ) { - w[i] = MAX(d1->d[i],d2->d[i]); td += w[i]; + w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i); } _NEWDL(d,n); d->td = td - d1->td; @@ -614,21 +704,18 @@ DP *rp; LIST hist; NODE node; - node = mknode(4,ONE,0,s1,ONE); + node = mknode(4,ONE,NULLP,s1,ONE); MKLIST(hist,node); MKNODE(TraceList,hist,0); - node = mknode(4,ONE,0,0,ONE); + node = mknode(4,ONE,NULLP,NULLP,ONE); chsgnd(s2,(DP *)&ARG2(node)); MKLIST(hist,node); MKNODE(node,hist,TraceList); TraceList = node; } } -void dp_sp_mod(p1,p2,mod,rp) -DP p1,p2; -int mod; -DP *rp; +void dp_sp_mod(DP p1,DP p2,int mod,DP *rp) { int i,n,td; int *w; @@ -639,7 +726,7 @@ DP *rp; n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; w = (int *)ALLOCA(n*sizeof(int)); for ( i = 0, td = 0; i < n; i++ ) { - w[i] = MAX(d1->d[i],d2->d[i]); td += w[i]; + w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i); } NEWDL_NOINIT(d,n); d->td = td - d1->td; for ( i = 0; i < n; i++ ) @@ -654,10 +741,7 @@ DP *rp; submd(CO,mod,t,u,rp); } -void _dp_sp_mod_dup(p1,p2,mod,rp) -DP p1,p2; -int mod; -DP *rp; +void _dp_sp_mod_dup(DP p1,DP p2,int mod,DP *rp) { int i,n,td; int *w; @@ -668,7 +752,7 @@ DP *rp; n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; w = (int *)ALLOCA(n*sizeof(int)); for ( i = 0, td = 0; i < n; i++ ) { - w[i] = MAX(d1->d[i],d2->d[i]); td += w[i]; + w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i); } _NEWDL(d,n); d->td = td - d1->td; for ( i = 0; i < n; i++ ) @@ -683,10 +767,7 @@ DP *rp; _addmd_destructive(mod,t,u,rp); } -void _dp_sp_mod(p1,p2,mod,rp) -DP p1,p2; -int mod; -DP *rp; +void _dp_sp_mod(DP p1,DP p2,int mod,DP *rp) { int i,n,td; int *w; @@ -697,7 +778,7 @@ DP *rp; n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; w = (int *)ALLOCA(n*sizeof(int)); for ( i = 0, td = 0; i < n; i++ ) { - w[i] = MAX(d1->d[i],d2->d[i]); td += w[i]; + w[i] = MAX(d1->d[i],d2->d[i]); td += MUL_WEIGHT(w[i],i); } NEWDL(d,n); d->td = td - d1->td; for ( i = 0; i < n; i++ ) @@ -719,11 +800,7 @@ DP *rp; * */ -void dp_red(p0,p1,p2,head,rest,dnp,multp) -DP p0,p1,p2; -DP *head,*rest; -P *dnp; -DP *multp; +void dp_red(DP p0,DP p1,DP p2,DP *head,DP *rest,P *dnp,DP *multp) { int i,n; DL d1,d2,d; @@ -732,13 +809,18 @@ DP *multp; Q c,c1,c2; N gn,tn; P g,a; + P p[2]; n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; NEWDL(d,n); d->td = d1->td - d2->td; for ( i = 0; i < n; i++ ) d->d[i] = d1->d[i]-d2->d[i]; c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(p2)->c; - if ( dp_fcoeffs ) { + if ( dp_fcoeffs == N_GFS ) { + p[0] = (P)c1; p[1] = (P)c2; + gcdsf(CO,p,2,&g); + divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a; + } else if ( dp_fcoeffs ) { /* do nothing */ } else if ( INT(c1) && INT(c2) ) { gcdn(NM(c1),NM(c2),&gn); @@ -749,6 +831,7 @@ DP *multp; } else { ezgcdpz(CO,(P)c1,(P)c2,&g); divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a; + add_denomlist(g); } NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td; *multp = s; @@ -757,16 +840,91 @@ DP *multp; *head = h; *rest = r; *dnp = (P)c2; } -/* m-reduction over a field */ +/* + * m-reduction by a marked poly + * do content reduction over Z or Q(x,...) + * do nothing over finite fields + * + */ -void dp_red_f(p1,p2,rest) -DP p1,p2; -DP *rest; + +void dp_red_marked(DP p0,DP p1,DP p2,DP hp2,DP *head,DP *rest,P *dnp,DP *multp) { int i,n; DL d1,d2,d; MP m; DP t,s,r,h; + Q c,c1,c2; + N gn,tn; + P g,a; + P p[2]; + + n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(hp2)->dl; + NEWDL(d,n); d->td = d1->td - d2->td; + for ( i = 0; i < n; i++ ) + d->d[i] = d1->d[i]-d2->d[i]; + c1 = (Q)BDY(p1)->c; c2 = (Q)BDY(hp2)->c; + if ( dp_fcoeffs == N_GFS ) { + p[0] = (P)c1; p[1] = (P)c2; + gcdsf(CO,p,2,&g); + divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a; + } else if ( dp_fcoeffs ) { + /* do nothing */ + } else if ( INT(c1) && INT(c2) ) { + gcdn(NM(c1),NM(c2),&gn); + if ( !UNIN(gn) ) { + divsn(NM(c1),gn,&tn); NTOQ(tn,SGN(c1),c); c1 = c; + divsn(NM(c2),gn,&tn); NTOQ(tn,SGN(c2),c); c2 = c; + } + } else { + ezgcdpz(CO,(P)c1,(P)c2,&g); + divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a; + } + NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td; + *multp = s; + muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); subd(CO,s,t,&r); + muldc(CO,p0,(P)c2,&h); + *head = h; *rest = r; *dnp = (P)c2; +} + +void dp_red_marked_mod(DP p0,DP p1,DP p2,DP hp2,int mod,DP *head,DP *rest,P *dnp,DP *multp) +{ + int i,n; + DL d1,d2,d; + MP m; + DP t,s,r,h; + P c1,c2,g,u; + + n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(hp2)->dl; + NEWDL(d,n); d->td = d1->td - d2->td; + for ( i = 0; i < n; i++ ) + d->d[i] = d1->d[i]-d2->d[i]; + c1 = (P)BDY(p1)->c; c2 = (P)BDY(hp2)->c; + gcdprsmp(CO,mod,c1,c2,&g); + divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u; + if ( NUM(c2) ) { + divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM; + } + NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0; + MKDP(n,m,s); s->sugar = d->td; + *multp = s; + mulmd(CO,mod,s,p2,&t); + if ( NUM(c2) ) { + submd(CO,mod,p1,t,&r); h = p0; + } else { + mulmdc(CO,mod,p1,c2,&s); submd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h); + } + *head = h; *rest = r; *dnp = c2; +} + +/* m-reduction over a field */ + +void dp_red_f(DP p1,DP p2,DP *rest) +{ + int i,n; + DL d1,d2,d; + MP m; + DP t,s; Obj a,b; n = p1->nv; @@ -784,11 +942,7 @@ DP *rest; muld(CO,s,p2,&t); addd(CO,p1,t,rest); } -void dp_red_mod(p0,p1,p2,mod,head,rest,dnp) -DP p0,p1,p2; -int mod; -DP *head,*rest; -P *dnp; +void dp_red_mod(DP p0,DP p1,DP p2,int mod,DP *head,DP *rest,P *dnp) { int i,n; DL d1,d2,d; @@ -818,17 +972,13 @@ P *dnp; struct oEGT eg_red_mod; -void _dp_red_mod_destructive(p1,p2,mod,rp) -DP p1,p2; -int mod; -DP *rp; +void _dp_red_mod_destructive(DP p1,DP p2,int mod,DP *rp) { int i,n; DL d1,d2,d; MP m; DP t,s; int c,c1,c2; - struct oEGT t0,t1; extern int do_weyl; n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; @@ -844,7 +994,8 @@ DP *rp; _mulmd_dup(mod,s,p2,&t); _free_dp(s); #else if ( do_weyl ) { - _weyl_mulmdm_dup(mod,p2,m,&t); _FREEMP(m); + _MKDP(n,m,s); s->sugar = d->td; + _mulmd_dup(mod,s,p2,&t); _free_dp(s); } else { _mulmdm_dup(mod,p2,m,&t); _FREEMP(m); } @@ -859,13 +1010,7 @@ DP *rp; * */ -void dp_true_nf(b,g,ps,full,rp,dnp) -NODE b; -DP g; -DP *ps; -int full; -DP *rp; -P *dnp; +void dp_true_nf(NODE b,DP g,DP *ps,int full,DP *rp,P *dnp) { DP u,p,d,s,t,dmy; NODE l; @@ -920,14 +1065,289 @@ P *dnp; *rp = d; *dnp = dn; } +void dp_removecont2(DP p1,DP p2,DP *r1p,DP *r2p,Q *contp) +{ + struct oVECT v; + int i,n1,n2,n; + MP m,m0,t; + Q *w; + Q h; + + if ( p1 ) { + for ( i = 0, m = BDY(p1); m; m = NEXT(m), i++ ); + n1 = i; + } else + n1 = 0; + if ( p2 ) { + for ( i = 0, m = BDY(p2); m; m = NEXT(m), i++ ); + n2 = i; + } else + n2 = 0; + n = n1+n2; + if ( !n ) { + *r1p = 0; *r2p = 0; *contp = ONE; return; + } + w = (Q *)ALLOCA(n*sizeof(Q)); + v.len = n; + v.body = (pointer *)w; + i = 0; + if ( p1 ) + for ( m = BDY(p1); i < n1; m = NEXT(m), i++ ) w[i] = (Q)m->c; + if ( p2 ) + for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) w[i] = (Q)m->c; + h = w[0]; removecont_array((P *)w,n,1); divq(h,w[0],contp); + i = 0; + if ( p1 ) { + for ( m0 = 0, t = BDY(p1); i < n1; i++, t = NEXT(t) ) { + NEXTMP(m0,m); m->c = (P)w[i]; m->dl = t->dl; + } + NEXT(m) = 0; + MKDP(p1->nv,m0,*r1p); (*r1p)->sugar = p1->sugar; + } else + *r1p = 0; + if ( p2 ) { + for ( m0 = 0, t = BDY(p2); i < n; i++, t = NEXT(t) ) { + NEXTMP(m0,m); m->c = (P)w[i]; m->dl = t->dl; + } + NEXT(m) = 0; + MKDP(p2->nv,m0,*r2p); (*r2p)->sugar = p2->sugar; + } else + *r2p = 0; +} + +/* true nf by a marked GB */ + +void dp_true_nf_marked(NODE b,DP g,DP *ps,DP *hps,DP *rp,P *nmp,P *dnp) +{ + DP u,p,d,s,t,dmy,hp; + NODE l; + MP m,mr; + int i,n,hmag; + int *wb; + int sugar,psugar,multiple; + P nm,tnm1,dn,tdn,tdn1; + Q cont; + + multiple = 0; + hmag = multiple*HMAG(g); + nm = (P)ONE; + dn = (P)ONE; + if ( !g ) { + *rp = 0; *dnp = dn; return; + } + for ( n = 0, l = b; l; l = NEXT(l), n++ ); + wb = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0, l = b; i < n; l = NEXT(l), i++ ) + wb[i] = QTOS((Q)BDY(l)); + sugar = g->sugar; + for ( d = 0; g; ) { + for ( u = 0, i = 0; i < n; i++ ) { + if ( dp_redble(g,hp = hps[wb[i]]) ) { + p = ps[wb[i]]; + dp_red_marked(d,g,p,hp,&t,&u,&tdn,&dmy); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + if ( !u ) { + goto last; + } else { + d = t; + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + } + break; + } + } + if ( u ) { + g = u; + if ( multiple && ((d && HMAG(d)>hmag) || (HMAG(g)>hmag)) ) { + dp_removecont2(d,g,&t,&u,&cont); d = t; g = u; + mulp(CO,nm,(P)cont,&tnm1); nm = tnm1; + if ( d ) + hmag = multiple*HMAG(d); + else + hmag = multiple*HMAG(g); + } + } else { + m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c; + NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td; + addd(CO,d,t,&s); d = s; + dp_rest(g,&t); g = t; + } + } +last: + if ( d ) { + dp_removecont2(d,0,&t,&u,&cont); d = t; + mulp(CO,nm,(P)cont,&tnm1); nm = tnm1; + d->sugar = sugar; + } + *rp = d; *nmp = nm; *dnp = dn; +} + +void dp_true_nf_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp) +{ + DP hp,u,p,d,s,t,dmy; + NODE l; + MP m,mr; + int i,n; + int *wb; + int sugar,psugar; + P dn,tdn,tdn1; + + dn = (P)ONEM; + if ( !g ) { + *rp = 0; *dnp = dn; return; + } + for ( n = 0, l = b; l; l = NEXT(l), n++ ); + wb = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0, l = b; i < n; l = NEXT(l), i++ ) + wb[i] = QTOS((Q)BDY(l)); + sugar = g->sugar; + for ( d = 0; g; ) { + for ( u = 0, i = 0; i < n; i++ ) { + if ( dp_redble(g,hp = hps[wb[i]]) ) { + p = ps[wb[i]]; + dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn,&dmy); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + if ( !u ) { + if ( d ) + d->sugar = sugar; + *rp = d; *dnp = dn; return; + } else { + d = t; + mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1; + } + break; + } + } + if ( u ) + g = u; + else { + m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c; + NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td; + addmd(CO,mod,d,t,&s); d = s; + dp_rest(g,&t); g = t; + } + } + if ( d ) + d->sugar = sugar; + *rp = d; *dnp = dn; +} + +/* true nf by a marked GB and collect quotients */ + +DP *dp_true_nf_and_quotient_marked (NODE b,DP g,DP *ps,DP *hps,DP *rp,P *dnp) +{ + DP u,p,d,s,t,dmy,hp,mult; + DP *q; + NODE l; + MP m,mr; + int i,n,j; + int *wb; + int sugar,psugar,multiple; + P nm,tnm1,dn,tdn,tdn1; + Q cont; + + dn = (P)ONE; + if ( !g ) { + *rp = 0; *dnp = dn; return 0; + } + for ( n = 0, l = b; l; l = NEXT(l), n++ ); + wb = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0, l = b; i < n; l = NEXT(l), i++ ) + wb[i] = QTOS((Q)BDY(l)); + q = (DP *)MALLOC(n*sizeof(DP)); + for ( i = 0; i < n; i++ ) q[i] = 0; + sugar = g->sugar; + for ( d = 0; g; ) { + for ( u = 0, i = 0; i < n; i++ ) { + if ( dp_redble(g,hp = hps[wb[i]]) ) { + p = ps[wb[i]]; + dp_red_marked(d,g,p,hp,&t,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + for ( j = 0; j < n; j++ ) { + muldc(CO,q[j],(P)tdn,&dmy); q[j] = dmy; + } + addd(CO,q[wb[i]],mult,&dmy); q[wb[i]] = dmy; + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + d = t; + if ( !u ) goto last; + break; + } + } + if ( u ) { + g = u; + } else { + m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c; + NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td; + addd(CO,d,t,&s); d = s; + dp_rest(g,&t); g = t; + } + } +last: + if ( d ) d->sugar = sugar; + *rp = d; *dnp = dn; + return q; +} + +DP *dp_true_nf_and_quotient_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp) +{ + DP u,p,d,s,t,dmy,hp,mult; + DP *q; + NODE l; + MP m,mr; + int i,n,j; + int *wb; + int sugar,psugar; + P dn,tdn,tdn1; + + for ( n = 0, l = b; l; l = NEXT(l), n++ ); + q = (DP *)MALLOC(n*sizeof(DP)); + for ( i = 0; i < n; i++ ) q[i] = 0; + dn = (P)ONEM; + if ( !g ) { + *rp = 0; *dnp = dn; return 0; + } + wb = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0, l = b; i < n; l = NEXT(l), i++ ) + wb[i] = QTOS((Q)BDY(l)); + sugar = g->sugar; + for ( d = 0; g; ) { + for ( u = 0, i = 0; i < n; i++ ) { + if ( dp_redble(g,hp = hps[wb[i]]) ) { + p = ps[wb[i]]; + dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + for ( j = 0; j < n; j++ ) { + mulmdc(CO,mod,q[j],(P)tdn,&dmy); q[j] = dmy; + } + addmd(CO,mod,q[wb[i]],mult,&dmy); q[wb[i]] = dmy; + mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1; + d = t; + if ( !u ) goto last; + break; + } + } + if ( u ) + g = u; + else { + m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c; + NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td; + addmd(CO,mod,d,t,&s); d = s; + dp_rest(g,&t); g = t; + } + } +last: + if ( d ) + d->sugar = sugar; + *rp = d; *dnp = dn; + return q; +} + /* nf computation over Z */ -void dp_nf_z(b,g,ps,full,multiple,rp) -NODE b; -DP g; -DP *ps; -int full,multiple; -DP *rp; +void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp) { DP u,p,d,s,t,dmy1; P dmy; @@ -998,15 +1418,9 @@ DP *rp; /* nf computation over a field */ -void dp_nf_f(b,g,ps,full,rp) -NODE b; -DP g; -DP *ps; -int full; -DP *rp; +void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp) { DP u,p,d,s,t; - P dmy; NODE l; MP m,mr; int i,n; @@ -1057,12 +1471,7 @@ DP *rp; /* nf computation over GF(mod) (only for internal use) */ -void dp_nf_mod(b,g,ps,mod,full,rp) -NODE b; -DP g; -DP *ps; -int mod,full; -DP *rp; +void dp_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp) { DP u,p,d,s,t; P dmy; @@ -1108,13 +1517,7 @@ DP *rp; *rp = d; } -void dp_true_nf_mod(b,g,ps,mod,full,rp,dnp) -NODE b; -DP g; -DP *ps; -int mod,full; -DP *rp; -P *dnp; +void dp_true_nf_mod(NODE b,DP g,DP *ps,int mod,int full,DP *rp,P *dnp) { DP u,p,d,s,t; NODE l; @@ -1169,17 +1572,12 @@ P *dnp; *rp = d; *dnp = dn; } -void _dp_nf_mod_destructive(b,g,ps,mod,full,rp) -NODE b; -DP g; -DP *ps; -int mod,full; -DP *rp; +void _dp_nf_mod_destructive(NODE b,DP g,DP *ps,int mod,int full,DP *rp) { - DP u,p,d,s,t; + DP u,p,d; NODE l; - MP m,mr,mrd; - int sugar,psugar,n,h_reducible,i; + MP m,mrd; + int sugar,psugar,n,h_reducible; if ( !g ) { *rp = 0; return; @@ -1230,10 +1628,7 @@ DP *rp; /* reduction by linear base over a field */ -void dp_lnf_f(p1,p2,g,r1p,r2p) -DP p1,p2; -NODE g; -DP *r1p,*r2p; +void dp_lnf_f(DP p1,DP p2,NODE g,DP *r1p,DP *r2p) { DP r1,r2,b1,b2,t,s; Obj c,c1,c2; @@ -1262,11 +1657,7 @@ DP *r1p,*r2p; /* reduction by linear base over GF(mod) */ -void dp_lnf_mod(p1,p2,g,mod,r1p,r2p) -DP p1,p2; -NODE g; -int mod; -DP *r1p,*r2p; +void dp_lnf_mod(DP p1,DP p2,NODE g,int mod,DP *r1p,DP *r2p) { DP r1,r2,b1,b2,t,s; P c; @@ -1294,11 +1685,7 @@ DP *r1p,*r2p; *r1p = r1; *r2p = r2; } -void dp_nf_tab_mod(p,tab,mod,rp) -DP p; -LIST *tab; -int mod; -DP *rp; +void dp_nf_tab_mod(DP p,LIST *tab,int mod,DP *rp) { DP s,t,u; MP m; @@ -1319,28 +1706,75 @@ DP *rp; *rp = s; } +void dp_nf_tab_f(DP p,LIST *tab,DP *rp) +{ + DP s,t,u; + MP m; + DL h; + int i,n; + + if ( !p ) { + *rp = p; return; + } + n = p->nv; + for ( s = 0, i = 0, m = BDY(p); m; m = NEXT(m) ) { + h = m->dl; + while ( !dl_equal(n,h,BDY((DP)BDY(BDY(tab[i])))->dl ) ) + i++; + muldc(CO,(DP)BDY(NEXT(BDY(tab[i]))),m->c,&t); + addd(CO,s,t,&u); s = u; + } + *rp = s; +} + /* * setting flags + * call create_order_spec with vl=0 to set old type order. * */ -int create_order_spec(obj,spec) -Obj obj; -struct order_spec *spec; +int create_order_spec(VL vl,Obj obj,struct order_spec **specp) { - int i,j,n,s,row,col; + int i,j,n,s,row,col,ret; + struct order_spec *spec; struct order_pair *l; NODE node,t,tn; MAT m; - pointer **b; + VECT v; + pointer **b,*bv; int **w; + if ( vl && obj && OID(obj) == O_LIST ) { + ret = create_composite_order_spec(vl,(LIST)obj,specp); + if ( show_orderspec ) + print_composite_order_spec(*specp); + return ret; + } + + *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec)); if ( !obj || NUM(obj) ) { spec->id = 0; spec->obj = obj; spec->ord.simple = QTOS((Q)obj); return 1; } else if ( OID(obj) == O_LIST ) { node = BDY((LIST)obj); + if ( !BDY(node) || NUM(BDY(node)) ) { + if ( length(node) < 2 ) + error("create_order_spec : invalid argument"); + create_order_spec(0,(Obj)BDY(NEXT(node)),&spec); + spec->id += 256; spec->obj = obj; + spec->ispot = (BDY(node)!=0); + if ( spec->ispot ) { + n = QTOS((Q)BDY(node)); + if ( n < 0 ) + spec->pot_nelim = -n; + else + spec->pot_nelim = 0; + } + *specp = spec; + return 1; + } + for ( n = 0, t = node; t; t = NEXT(t), n++ ); l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair)); for ( i = 0, t = node, s = 0; i < n; t = NEXT(t), i++ ) { @@ -1366,14 +1800,339 @@ struct order_spec *spec; return 0; } +void print_composite_order_spec(struct order_spec *spec) +{ + int nv,n,len,i,j,k,start; + struct weight_or_block *worb; + + nv = spec->nv; + n = spec->ord.composite.length; + worb = spec->ord.composite.w_or_b; + for ( i = 0; i < n; i++, worb++ ) { + len = worb->length; + printf("[ "); + switch ( worb->type ) { + case IS_DENSE_WEIGHT: + for ( j = 0; j < len; j++ ) + printf("%d ",worb->body.dense_weight[j]); + for ( ; j < nv; j++ ) + printf("0 "); + break; + case IS_SPARSE_WEIGHT: + for ( j = 0, k = 0; j < nv; j++ ) + if ( j == worb->body.sparse_weight[k].pos ) + printf("%d ",worb->body.sparse_weight[k++].value); + else + printf("0 "); + break; + case IS_BLOCK: + start = worb->body.block.start; + for ( j = 0; j < start; j++ ) printf("0 "); + switch ( worb->body.block.order ) { + case 0: + for ( k = 0; k < len; k++, j++ ) printf("R "); + break; + case 1: + for ( k = 0; k < len; k++, j++ ) printf("G "); + break; + case 2: + for ( k = 0; k < len; k++, j++ ) printf("L "); + break; + } + for ( ; j < nv; j++ ) printf("0 "); + break; + } + printf("]\n"); + } +} + +struct order_spec *append_block(struct order_spec *spec, + int nv,int nalg,int ord) +{ + MAT m,mat; + int i,j,row,col,n; + Q **b,**wp; + int **w; + NODE t,s,s0; + struct order_pair *l,*l0; + int n0,nv0; + LIST list0,list1,list; + Q oq,nq; + struct order_spec *r; + + r = (struct order_spec *)MALLOC(sizeof(struct order_spec)); + switch ( spec->id ) { + case 0: + STOQ(spec->ord.simple,oq); STOQ(nv,nq); + t = mknode(2,oq,nq); MKLIST(list0,t); + STOQ(ord,oq); STOQ(nalg,nq); + t = mknode(2,oq,nq); MKLIST(list1,t); + t = mknode(2,list0,list1); MKLIST(list,t); + l = (struct order_pair *)MALLOC_ATOMIC(2*sizeof(struct order_pair)); + l[0].order = spec->ord.simple; l[0].length = nv; + l[1].order = ord; l[1].length = nalg; + r->id = 1; r->obj = (Obj)list; + r->ord.block.order_pair = l; + r->ord.block.length = 2; + r->nv = nv+nalg; + break; + case 1: + if ( spec->nv != nv ) + error("append_block : number of variables mismatch"); + l0 = spec->ord.block.order_pair; + n0 = spec->ord.block.length; + nv0 = spec->nv; + list0 = (LIST)spec->obj; + n = n0+1; + l = (struct order_pair *)MALLOC_ATOMIC(n*sizeof(struct order_pair)); + for ( i = 0; i < n0; i++ ) + l[i] = l0[i]; + l[i].order = ord; l[i].length = nalg; + for ( t = BDY(list0), s0 = 0; t; t = NEXT(t) ) { + NEXTNODE(s0,s); BDY(s) = BDY(t); + } + STOQ(ord,oq); STOQ(nalg,nq); + t = mknode(2,oq,nq); MKLIST(list,t); + NEXTNODE(s0,s); BDY(s) = (pointer)list; NEXT(s) = 0; + MKLIST(list,s0); + r->id = 1; r->obj = (Obj)list; + r->ord.block.order_pair = l; + r->ord.block.length = n; + r->nv = nv+nalg; + break; + case 2: + if ( spec->nv != nv ) + error("append_block : number of variables mismatch"); + m = (MAT)spec->obj; + row = m->row; col = m->col; b = (Q **)BDY(m); + w = almat(row+nalg,col+nalg); + MKMAT(mat,row+nalg,col+nalg); wp = (Q **)BDY(mat); + for ( i = 0; i < row; i++ ) + for ( j = 0; j < col; j++ ) { + w[i][j] = QTOS(b[i][j]); + wp[i][j] = b[i][j]; + } + for ( i = 0; i < nalg; i++ ) { + w[i+row][i+col] = 1; + wp[i+row][i+col] = ONE; + } + r->id = 2; r->obj = (Obj)mat; + r->nv = col+nalg; r->ord.matrix.row = row+nalg; + r->ord.matrix.matrix = w; + break; + case 3: + default: + /* XXX */ + error("append_block : not implemented yet"); + } + return r; +} + +int comp_sw(struct sparse_weight *a, struct sparse_weight *b) +{ + if ( a->pos > b->pos ) return 1; + else if ( a->pos < b->pos ) return -1; + else return 0; +} + +/* order = [w_or_b, w_or_b, ... ] */ +/* w_or_b = w or b */ +/* w = [1,2,...] or [x,1,y,2,...] */ +/* b = [@lex,x,y,...,z] etc */ + +int create_composite_order_spec(VL vl,LIST order,struct order_spec **specp) +{ + NODE wb,t,p; + struct order_spec *spec; + VL tvl; + int n,i,j,k,l,start,end,len,w; + int *dw; + struct sparse_weight *sw; + struct weight_or_block *w_or_b; + Obj a0; + NODE a; + V v,sv,ev; + SYMBOL sym; + int *top; + + /* l = number of vars in vl */ + for ( l = 0, tvl = vl; tvl; tvl = NEXT(tvl), l++ ); + /* n = number of primitives in order */ + wb = BDY(order); + n = length(wb); + *specp = spec = (struct order_spec *)MALLOC(sizeof(struct order_spec)); + spec->id = 3; + spec->obj = (Obj)order; + spec->nv = l; + spec->ord.composite.length = n; + w_or_b = spec->ord.composite.w_or_b = (struct weight_or_block *) + MALLOC(sizeof(struct weight_or_block)*(n+1)); + + /* top : register the top variable in each w_or_b specification */ + top = (int *)ALLOCA(l*sizeof(int)); + for ( i = 0; i < l; i++ ) top[i] = 0; + + for ( t = wb, i = 0; t; t = NEXT(t), i++ ) { + if ( !BDY(t) || OID((Obj)BDY(t)) != O_LIST ) + error("a list of lists must be specified for the key \"order\""); + a = BDY((LIST)BDY(t)); + len = length(a); + a0 = (Obj)BDY(a); + if ( !a0 || OID(a0) == O_N ) { + /* a is a dense weight vector */ + dw = (int *)MALLOC(sizeof(int)*len); + for ( j = 0, p = a; j < len; p = NEXT(p), j++ ) { + if ( !INT((Q)BDY(p)) ) + error("a dense weight vector must be specified as a list of integers"); + dw[j] = QTOS((Q)BDY(p)); + } + w_or_b[i].type = IS_DENSE_WEIGHT; + w_or_b[i].length = len; + w_or_b[i].body.dense_weight = dw; + + /* find the top */ + for ( k = 0; k < len && !dw[k]; k++ ); + if ( k < len ) top[k] = 1; + + } else if ( OID(a0) == O_P ) { + /* a is a sparse weight vector */ + len >>= 1; + sw = (struct sparse_weight *) + MALLOC(sizeof(struct sparse_weight)*len); + for ( j = 0, p = a; j < len; j++ ) { + if ( !BDY(p) || OID((P)BDY(p)) != O_P ) + error("a sparse weight vector must be specified as [var1,weight1,...]"); + v = VR((P)BDY(p)); p = NEXT(p); + for ( tvl = vl, k = 0; tvl && tvl->v != v; + k++, tvl = NEXT(tvl) ); + if ( !tvl ) + error("invalid variable name in a sparse weight vector"); + sw[j].pos = k; + if ( !INT((Q)BDY(p)) ) + error("a sparse weight vector must be specified as [var1,weight1,...]"); + sw[j].value = QTOS((Q)BDY(p)); p = NEXT(p); + } + qsort(sw,len,sizeof(struct sparse_weight), + (int (*)(const void *,const void *))comp_sw); + w_or_b[i].type = IS_SPARSE_WEIGHT; + w_or_b[i].length = len; + w_or_b[i].body.sparse_weight = sw; + + /* find the top */ + for ( k = 0; k < len && !sw[k].value; k++ ); + if ( k < len ) top[sw[k].pos] = 1; + } else if ( OID(a0) == O_RANGE ) { + /* [range(v1,v2),w] */ + sv = VR((P)(((RANGE)a0)->start)); + ev = VR((P)(((RANGE)a0)->end)); + for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) ); + if ( !tvl ) + error("invalid range"); + for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) ); + if ( !tvl ) + error("invalid range"); + len = end-start+1; + sw = (struct sparse_weight *) + MALLOC(sizeof(struct sparse_weight)*len); + w = QTOS((Q)BDY(NEXT(a))); + for ( tvl = vl, k = 0; k < start; k++, tvl = NEXT(tvl) ); + for ( j = 0 ; k <= end; k++, tvl = NEXT(tvl), j++ ) { + sw[j].pos = k; + sw[j].value = w; + } + w_or_b[i].type = IS_SPARSE_WEIGHT; + w_or_b[i].length = len; + w_or_b[i].body.sparse_weight = sw; + + /* register the top */ + if ( w ) top[start] = 1; + } else if ( OID(a0) == O_SYMBOL ) { + /* a is a block */ + sym = (SYMBOL)a0; a = NEXT(a); len--; + if ( OID((Obj)BDY(a)) == O_RANGE ) { + sv = VR((P)(((RANGE)BDY(a))->start)); + ev = VR((P)(((RANGE)BDY(a))->end)); + for ( tvl = vl, start = 0; tvl && tvl->v != sv; start++, tvl = NEXT(tvl) ); + if ( !tvl ) + error("invalid range"); + for ( end = start; tvl && tvl->v != ev; end++, tvl = NEXT(tvl) ); + if ( !tvl ) + error("invalid range"); + len = end-start+1; + } else { + for ( start = 0, tvl = vl; tvl->v != VR((P)BDY(a)); + tvl = NEXT(tvl), start++ ); + for ( p = NEXT(a), tvl = NEXT(tvl); p; + p = NEXT(p), tvl = NEXT(tvl) ) { + if ( !BDY(p) || OID((P)BDY(p)) != O_P ) + error("a block must be specified as [ordsymbol,var1,var2,...]"); + if ( tvl->v != VR((P)BDY(p)) ) break; + } + if ( p ) + error("a block must be contiguous in the variable list"); + } + w_or_b[i].type = IS_BLOCK; + w_or_b[i].length = len; + w_or_b[i].body.block.start = start; + if ( !strcmp(sym->name,"@grlex") ) + w_or_b[i].body.block.order = 0; + else if ( !strcmp(sym->name,"@glex") ) + w_or_b[i].body.block.order = 1; + else if ( !strcmp(sym->name,"@lex") ) + w_or_b[i].body.block.order = 2; + else + error("invalid ordername"); + /* register the tops */ + for ( j = 0, k = start; j < len; j++, k++ ) + top[k] = 1; + } + } + for ( k = 0; k < l && top[k]; k++ ); + if ( k < l ) { + /* incomplete order specification; add @grlex */ + w_or_b[n].type = IS_BLOCK; + w_or_b[n].length = l; + w_or_b[n].body.block.start = 0; + w_or_b[n].body.block.order = 0; + spec->ord.composite.length = n+1; + } +} + +/* module order spec */ + +void create_modorder_spec(int id,LIST shift,struct modorder_spec **s) +{ + struct modorder_spec *spec; + NODE n,t; + LIST list; + int *ds; + int i,l; + Q q; + + *s = spec = (struct modorder_spec *)MALLOC(sizeof(struct modorder_spec)); + spec->id = id; + if ( shift ) { + n = BDY(shift); + spec->len = l = length(n); + spec->degree_shift = ds = (int *)MALLOC_ATOMIC(l*sizeof(int)); + for ( t = n, i = 0; t; t = NEXT(t), i++ ) + ds[i] = QTOS((Q)BDY(t)); + } else { + spec->len = 0; + spec->degree_shift = 0; + } + STOQ(id,q); + n = mknode(2,q,shift); + MKLIST(list,n); + spec->obj = (Obj)list; +} + /* * converters * */ -void dp_homo(p,rp) -DP p; -DP *rp; +void dp_homo(DP p,DP *rp) { MP m,mr,mr0; int i,n,nv,td; @@ -1397,9 +2156,7 @@ DP *rp; } } -void dp_dehomo(p,rp) -DP p; -DP *rp; +void dp_dehomo(DP p,DP *rp) { MP m,mr,mr0; int i,n,nv; @@ -1422,11 +2179,7 @@ DP *rp; } } -void dp_mod(p,mod,subst,rp) -DP p; -int mod; -NODE subst; -DP *rp; +void dp_mod(DP p,int mod,NODE subst,DP *rp) { MP m,mr,mr0; P t,s,s1; @@ -1453,9 +2206,7 @@ DP *rp; } } -void dp_rat(p,rp) -DP p; -DP *rp; +void dp_rat(DP p,DP *rp) { MP m,mr,mr0; @@ -1473,14 +2224,16 @@ DP *rp; } -void homogenize_order(old,n,new) -struct order_spec *old,*new; -int n; +void homogenize_order(struct order_spec *old,int n,struct order_spec **newp) { struct order_pair *l; int length,nv,row,i,j; int **newm,**oldm; + struct order_spec *new; + int onv,nnv,nlen,olen,owlen; + struct weight_or_block *owb,*nwb; + *newp = new = (struct order_spec *)MALLOC(sizeof(struct order_spec)); switch ( old->id ) { case 0: switch ( old->ord.simple ) { @@ -1506,17 +2259,18 @@ int n; error("homogenize_order : invalid input"); } break; - case 1: + case 1: case 257: length = old->ord.block.length; l = (struct order_pair *) MALLOC_ATOMIC((length+1)*sizeof(struct order_pair)); bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair)); l[length].order = 2; l[length].length = 1; - new->id = 1; new->nv = n+1; + new->id = old->id; new->nv = n+1; new->ord.block.order_pair = l; new->ord.block.length = length+1; + new->ispot = old->ispot; break; - case 2: + case 2: case 258: nv = old->nv; row = old->ord.matrix.row; oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1); for ( i = 0; i <= nv; i++ ) @@ -1526,17 +2280,77 @@ int n; newm[i+1][j] = oldm[i][j]; newm[i+1][j] = 0; } - new->id = 2; new->nv = nv+1; + new->id = old->id; new->nv = nv+1; new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm; + new->ispot = old->ispot; break; + case 3: case 259: + onv = old->nv; + nnv = onv+1; + olen = old->ord.composite.length; + nlen = olen+1; + owb = old->ord.composite.w_or_b; + nwb = (struct weight_or_block *) + MALLOC(nlen*sizeof(struct weight_or_block)); + for ( i = 0; i < olen; i++ ) { + nwb[i].type = owb[i].type; + switch ( owb[i].type ) { + case IS_DENSE_WEIGHT: + owlen = owb[i].length; + nwb[i].length = owlen+1; + nwb[i].body.dense_weight = (int *)MALLOC((owlen+1)*sizeof(int)); + for ( j = 0; j < owlen; j++ ) + nwb[i].body.dense_weight[j] = owb[i].body.dense_weight[j]; + nwb[i].body.dense_weight[owlen] = 0; + break; + case IS_SPARSE_WEIGHT: + nwb[i].length = owb[i].length; + nwb[i].body.sparse_weight = owb[i].body.sparse_weight; + break; + case IS_BLOCK: + nwb[i].length = owb[i].length; + nwb[i].body.block = owb[i].body.block; + break; + } + } + nwb[i].type = IS_SPARSE_WEIGHT; + nwb[i].body.sparse_weight = + (struct sparse_weight *)MALLOC(sizeof(struct sparse_weight)); + nwb[i].body.sparse_weight[0].pos = onv; + nwb[i].body.sparse_weight[0].value = 1; + new->id = old->id; + new->nv = nnv; + new->ord.composite.length = nlen; + new->ord.composite.w_or_b = nwb; + new->ispot = old->ispot; + print_composite_order_spec(new); + break; + case 256: /* simple module order */ + switch ( old->ord.simple ) { + case 0: + new->id = 256; new->ord.simple = 0; break; + case 1: + l = (struct order_pair *) + MALLOC_ATOMIC(2*sizeof(struct order_pair)); + l[0].length = n; l[0].order = old->ord.simple; + l[1].length = 1; l[1].order = 2; + new->id = 257; + new->ord.block.order_pair = l; + new->ord.block.length = 2; new->nv = n+1; + break; + case 2: + new->id = 256; new->ord.simple = 1; break; + default: + error("homogenize_order : invalid input"); + } + new->ispot = old->ispot; + break; default: error("homogenize_order : invalid input"); } } -void qltozl(w,n,dvr) -Q *w,*dvr; -int n; +void qltozl(Q *w,int n,Q *dvr) { N nm,dn; N g,l1,l2,l3; @@ -1551,9 +2365,12 @@ int n; v.id = O_VECT; v.len = n; v.body = (pointer *)w; igcdv(&v,dvr); return; } - c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c); - for ( i = 1; i < n; i++ ) { - c = w[i]; l1 = INT(c) ? ONEN : DN(c); + for ( i = 0; !w[i]; i++ ); + c = w[i]; nm = NM(c); dn = INT(c) ? ONEN : DN(c); + for ( i++; i < n; i++ ) { + c = w[i]; + if ( !c ) continue; + l1 = INT(c) ? ONEN : DN(c); gcdn(nm,NM(c),&g); nm = g; gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn); } @@ -1564,15 +2381,12 @@ int n; *dvr = d; } -int comp_nm(a,b) -Q *a,*b; +int comp_nm(Q *a,Q *b) { return cmpn((*a)?NM(*a):0,(*b)?NM(*b):0); } -void sortbynm(w,n) -Q *w; -int n; +void sortbynm(Q *w,int n) { qsort(w,n,sizeof(Q),(int (*)(const void *,const void *))comp_nm); } @@ -1583,8 +2397,7 @@ int n; * */ -int dp_redble(p1,p2) -DP p1,p2; +int dp_redble(DP p1,DP p2) { int i,n; DL d1,d2; @@ -1600,9 +2413,7 @@ DP p1,p2; } } -void dp_subd(p1,p2,rp) -DP p1,p2; -DP *rp; +void dp_subd(DP p1,DP p2,DP *rp) { int i,n; DL d1,d2,d; @@ -1618,10 +2429,7 @@ DP *rp; *rp = s; } -void dltod(d,n,rp) -DL d; -int n; -DP *rp; +void dltod(DL d,int n,DP *rp) { MP m; DP s; @@ -1631,9 +2439,7 @@ DP *rp; *rp = s; } -void dp_hm(p,rp) -DP p; -DP *rp; +void dp_hm(DP p,DP *rp) { MP m,mr; @@ -1646,9 +2452,21 @@ DP *rp; } } -void dp_rest(p,rp) -DP p,*rp; +void dp_ht(DP p,DP *rp) { + MP m,mr; + + if ( !p ) + *rp = 0; + else { + m = BDY(p); + NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0; + MKDP(p->nv,mr,*rp); (*rp)->sugar = mr->dl->td; /* XXX */ + } +} + +void dp_rest(DP p,DP *rp) +{ MP m; m = BDY(p); @@ -1661,24 +2479,21 @@ DP p,*rp; } } -DL lcm_of_DL(nv,dl1,dl2,dl) -int nv; -DL dl1,dl2; -register DL dl; +DL lcm_of_DL(int nv,DL dl1,DL dl2,DL dl) { - register int n, *d1, *d2, *d, td; + register int i, *d1, *d2, *d, td; if ( !dl ) NEWDL(dl,nv); d = dl->d, d1 = dl1->d, d2 = dl2->d; - for ( td = 0, n = nv; --n >= 0; d1++, d2++, d++ ) - td += (*d = *d1 > *d2 ? *d1 : *d2 ); + for ( td = 0, i = 0; i < nv; d1++, d2++, d++, i++ ) { + *d = *d1 > *d2 ? *d1 : *d2; + td += MUL_WEIGHT(*d,i); + } dl->td = td; return dl; } -int dl_equal(nv,dl1,dl2) -int nv; -DL dl1, dl2; +int dl_equal(int nv,DL dl1,DL dl2) { register int *d1, *d2, n; @@ -1688,8 +2503,7 @@ DL dl1, dl2; return 1; } -int dp_nt(p) -DP p; +int dp_nt(DP p) { int i; MP m; @@ -1702,8 +2516,7 @@ DP p; } } -int dp_homogeneous(p) -DP p; +int dp_homogeneous(DP p) { MP m; int d; @@ -1722,9 +2535,7 @@ DP p; } } -_print_mp(nv,m) -int nv; -MP m; +void _print_mp(int nv,MP m) { int i; @@ -1740,4 +2551,582 @@ MP m; fprintf(stderr,">",C(m)); } fprintf(stderr,"\n"); +} + +static int cmp_mp_nvar; + +int comp_mp(MP *a,MP *b) +{ + return -(*cmpdl)(cmp_mp_nvar,(*a)->dl,(*b)->dl); +} + +void dp_sort(DP p,DP *rp) +{ + MP t,mp,mp0; + int i,n; + DP r; + MP *w; + + if ( !p ) { + *rp = 0; + return; + } + for ( t = BDY(p), n = 0; t; t = NEXT(t), n++ ); + w = (MP *)ALLOCA(n*sizeof(MP)); + for ( t = BDY(p), i = 0; i < n; t = NEXT(t), i++ ) + w[i] = t; + cmp_mp_nvar = NV(p); + qsort(w,n,sizeof(MP),(int (*)(const void *,const void *))comp_mp); + mp0 = 0; + for ( i = n-1; i >= 0; i-- ) { + NEWMP(mp); mp->dl = w[i]->dl; C(mp) = C(w[i]); + NEXT(mp) = mp0; mp0 = mp; + } + MKDP(p->nv,mp0,r); + r->sugar = p->sugar; + *rp = r; +} + +DP extract_initial_term_from_dp(DP p,int *weight,int n); +LIST extract_initial_term(LIST f,int *weight,int n); + +DP extract_initial_term_from_dp(DP p,int *weight,int n) +{ + int w,t,i,top; + MP m,r0,r; + DP dp; + + if ( !p ) return 0; + top = 1; + for ( m = BDY(p); m; m = NEXT(m) ) { + for ( i = 0, t = 0; i < n; i++ ) + t += weight[i]*m->dl->d[i]; + if ( top || t > w ) { + r0 = 0; + w = t; + top = 0; + } + if ( t == w ) { + NEXTMP(r0,r); + r->dl = m->dl; + r->c = m->c; + } + } + NEXT(r) = 0; + MKDP(p->nv,r0,dp); + return dp; +} + +LIST extract_initial_term(LIST f,int *weight,int n) +{ + NODE nd,r0,r; + Obj p; + LIST l; + + nd = BDY(f); + for ( r0 = 0; nd; nd = NEXT(nd) ) { + NEXTNODE(r0,r); + p = (Obj)BDY(nd); + BDY(r) = (pointer)extract_initial_term_from_dp((DP)p,weight,n); + } + if ( r0 ) NEXT(r) = 0; + MKLIST(l,r0); + return l; +} + +LIST dp_initial_term(LIST f,struct order_spec *ord) +{ + int n,l,i; + struct weight_or_block *worb; + int *weight; + + switch ( ord->id ) { + case 2: /* matrix order */ + /* extract the first row */ + n = ord->nv; + weight = ord->ord.matrix.matrix[0]; + return extract_initial_term(f,weight,n); + case 3: /* composite order */ + /* the first w_or_b */ + worb = ord->ord.composite.w_or_b; + switch ( worb->type ) { + case IS_DENSE_WEIGHT: + n = worb->length; + weight = worb->body.dense_weight; + return extract_initial_term(f,weight,n); + case IS_SPARSE_WEIGHT: + n = ord->nv; + weight = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0; i < n; i++ ) weight[i] = 0; + l = worb->length; + for ( i = 0; i < l; i++ ) + weight[worb->body.sparse_weight[i].pos] + = worb->body.sparse_weight[i].value; + return extract_initial_term(f,weight,n); + default: + error("dp_initial_term : unsupported order"); + } + default: + error("dp_initial_term : unsupported order"); + } +} + +int highest_order_dp(DP p,int *weight,int n); +LIST highest_order(LIST f,int *weight,int n); + +int highest_order_dp(DP p,int *weight,int n) +{ + int w,t,i,top; + MP m; + + if ( !p ) return -1; + top = 1; + for ( m = BDY(p); m; m = NEXT(m) ) { + for ( i = 0, t = 0; i < n; i++ ) + t += weight[i]*m->dl->d[i]; + if ( top || t > w ) { + w = t; + top = 0; + } + } + return w; +} + +LIST highest_order(LIST f,int *weight,int n) +{ + int h; + NODE nd,r0,r; + Obj p; + LIST l; + Q q; + + nd = BDY(f); + for ( r0 = 0; nd; nd = NEXT(nd) ) { + NEXTNODE(r0,r); + p = (Obj)BDY(nd); + h = highest_order_dp((DP)p,weight,n); + STOQ(h,q); + BDY(r) = (pointer)q; + } + if ( r0 ) NEXT(r) = 0; + MKLIST(l,r0); + return l; +} + +LIST dp_order(LIST f,struct order_spec *ord) +{ + int n,l,i; + struct weight_or_block *worb; + int *weight; + + switch ( ord->id ) { + case 2: /* matrix order */ + /* extract the first row */ + n = ord->nv; + weight = ord->ord.matrix.matrix[0]; + return highest_order(f,weight,n); + case 3: /* composite order */ + /* the first w_or_b */ + worb = ord->ord.composite.w_or_b; + switch ( worb->type ) { + case IS_DENSE_WEIGHT: + n = worb->length; + weight = worb->body.dense_weight; + return highest_order(f,weight,n); + case IS_SPARSE_WEIGHT: + n = ord->nv; + weight = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0; i < n; i++ ) weight[i] = 0; + l = worb->length; + for ( i = 0; i < l; i++ ) + weight[worb->body.sparse_weight[i].pos] + = worb->body.sparse_weight[i].value; + return highest_order(f,weight,n); + default: + error("dp_initial_term : unsupported order"); + } + default: + error("dp_initial_term : unsupported order"); + } +} + +int dpv_ht(DPV p,DP *h) +{ + int len,max,maxi,i,t; + DP *e; + MP m,mr; + + len = p->len; + e = p->body; + max = -1; + maxi = -1; + for ( i = 0; i < len; i++ ) + if ( e[i] && (t = BDY(e[i])->dl->td) > max ) { + max = t; + maxi = i; + } + if ( max < 0 ) { + *h = 0; + return -1; + } else { + m = BDY(e[maxi]); + NEWMP(mr); mr->dl = m->dl; mr->c = (P)ONE; NEXT(mr) = 0; + MKDP(e[maxi]->nv,mr,*h); (*h)->sugar = mr->dl->td; /* XXX */ + return maxi; + } +} + +/* return 1 if 0 <_w1 v && v <_w2 0 */ + +int in_c12(int n,int *v,int row1,int **w1,int row2, int **w2) +{ + int t1,t2; + + t1 = compare_zero(n,v,row1,w1); + t2 = compare_zero(n,v,row2,w2); + if ( t1 > 0 && t2 < 0 ) return 1; + else return 0; +} + +/* 0 < u => 1, 0 > u => -1 */ + +int compare_zero(int n,int *u,int row,int **w) +{ + int i,j,t; + int *wi; + + for ( i = 0; i < row; i++ ) { + wi = w[i]; + for ( j = 0, t = 0; j < n; j++ ) t += u[j]*wi[j]; + if ( t > 0 ) return 1; + else if ( t < 0 ) return -1; + } + return 0; +} + +/* functions for generic groebner walk */ +/* u=0 means u=-infty */ + +int compare_facet_preorder(int n,int *u,int *v, + int row1,int **w1,int row2,int **w2) +{ + int i,j,s,t,tu,tv; + int *w2i,*uv; + + if ( !u ) return 1; + uv = W_ALLOC(n); + for ( i = 0; i < row2; i++ ) { + w2i = w2[i]; + for ( j = 0, tu = tv = 0; j < n; j++ ) + if ( s = w2i[j] ) { + tu += s*u[j]; tv += s*v[j]; + } + for ( j = 0; j < n; j++ ) uv[j] = u[j]*tv-v[j]*tu; + t = compare_zero(n,uv,row1,w1); + if ( t > 0 ) return 1; + else if ( t < 0 ) return 0; + } + return 1; +} + +Q inner_product_with_small_vector(VECT w,int *v) +{ + int n,i; + Q q,s,t,u; + + n = w->len; + s = 0; + for ( i = 0; i < n; i++ ) { + STOQ(v[i],q); mulq((Q)w->body[i],q,&t); addq(t,s,&u); s = u; + } + return s; +} + +Q compute_last_t(NODE g,NODE gh,Q t,VECT w1,VECT w2,NODE *homo,VECT *wp) +{ + int n,i; + int *wt; + Q last,d1,d2,dn,nm,s,t1; + VECT wd,wt1,wt2,w; + NODE tg,tgh; + MP f; + int *h; + NODE r0,r; + MP m0,m; + DP d; + + n = w1->len; + wt = W_ALLOC(n); + last = ONE; + /* t1 = 1-t */ + for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) { + f = BDY((DP)BDY(tg)); + h = BDY((DP)BDY(tgh))->dl->d; + for ( ; f; f = NEXT(f) ) { + for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i]; + for ( i = 0; i < n && !wt[i]; i++ ); + if ( i == n ) continue; + d1 = inner_product_with_small_vector(w1,wt); + d2 = inner_product_with_small_vector(w2,wt); + nm = d1; subq(d1,d2,&dn); + /* if d1=d2 then nothing happens */ + if ( !dn ) continue; + /* s satisfies ds = 0*/ + divq(nm,dn,&s); + + if ( cmpq(s,t) > 0 && cmpq(s,last) < 0 ) + last = s; + else if ( !cmpq(s,t) ) { + if ( cmpq(d2,0) < 0 ) { + last = t; + break; + } + } + } + } + if ( !last ) { + dn = ONE; nm = 0; + } else { + NTOQ(NM(last),1,nm); + if ( INT(last) ) dn = ONE; + else { + NTOQ(DN(last),1,dn); + } + } + /* (1-n/d)*w1+n/d*w2 -> w=(d-n)*w1+n*w2 */ + subq(dn,nm,&t1); mulvect(CO,(Obj)w1,(Obj)t1,(Obj *)&wt1); + mulvect(CO,(Obj)w2,(Obj)nm,(Obj *)&wt2); addvect(CO,wt1,wt2,&w); + + r0 = 0; + for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) { + f = BDY((DP)BDY(tg)); + h = BDY((DP)BDY(tgh))->dl->d; + for ( m0 = 0; f; f = NEXT(f) ) { + for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i]; + for ( i = 0; i < n && !wt[i]; i++ ); + if ( !inner_product_with_small_vector(w,wt) ) { + NEXTMP(m0,m); m->c = f->c; m->dl = f->dl; + } + } + NEXT(m) = 0; + MKDP(((DP)BDY(tg))->nv,m0,d); d->sugar = ((DP)BDY(tg))->sugar; + NEXTNODE(r0,r); BDY(r) = (pointer)d; + } + NEXT(r) = 0; + *homo = r0; + *wp = w; + return last; +} + +/* return 0 if last_w = infty */ + +NODE compute_last_w(NODE g,NODE gh,int n,int **w, + int row1,int **w1,int row2,int **w2) +{ + DP d; + MP f,m0,m; + int *wt,*v,*h; + NODE t,s,n0,tn,n1,r0,r; + int i; + + wt = W_ALLOC(n); + n0 = 0; + for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) { + f = BDY((DP)BDY(t)); + h = BDY((DP)BDY(s))->dl->d; + for ( ; f; f = NEXT(f) ) { + for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i]; + for ( i = 0; i < n && !wt[i]; i++ ); + if ( i == n ) continue; + + if ( in_c12(n,wt,row1,w1,row2,w2) && + compare_facet_preorder(n,*w,wt,row1,w1,row2,w2) ) { + v = (int *)MALLOC_ATOMIC(n*sizeof(int)); + for ( i = 0; i < n; i++ ) v[i] = wt[i]; + MKNODE(n1,v,n0); n0 = n1; + } + } + } + if ( !n0 ) return 0; + for ( t = n0; t; t = NEXT(t) ) { + v = (int *)BDY(t); + for ( s = n0; s; s = NEXT(s) ) + if ( !compare_facet_preorder(n,v,(int *)BDY(s),row1,w1,row2,w2) ) + break; + if ( !s ) { + *w = v; + break; + } + } + if ( !t ) + error("compute_last_w : cannot happen"); + r0 = 0; + for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) { + f = BDY((DP)BDY(t)); + h = BDY((DP)BDY(s))->dl->d; + for ( m0 = 0; f; f = NEXT(f) ) { + for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i]; + for ( i = 0; i < n && !wt[i]; i++ ); + if ( i == n || + (compare_facet_preorder(n,wt,*w,row1,w1,row2,w2) + && compare_facet_preorder(n,*w,wt,row1,w1,row2,w2)) ) { + NEXTMP(m0,m); m->c = f->c; m->dl = f->dl; + } + } + NEXT(m) = 0; + MKDP(((DP)BDY(t))->nv,m0,d); d->sugar = ((DP)BDY(t))->sugar; + NEXTNODE(r0,r); BDY(r) = (pointer)d; + } + NEXT(r) = 0; + return r0; +} + +/* compute a sufficient set of d(f)=u-v */ + +NODE compute_essential_df(DP *g,DP *gh,int ng) +{ + int nv,i,j,k,t,lj; + NODE r,r1,ri,rt,r0; + MP m; + MP *mj; + DL di,hj,dl,dlt; + int *d,*dt; + LIST l; + Q q; + + nv = g[0]->nv; + r = 0; + for ( j = 0; j < ng; j++ ) { + for ( m = BDY(g[j]), lj = 0; m; m = NEXT(m), lj++ ); + mj = (MP *)ALLOCA(lj*sizeof(MP)); + for ( m = BDY(g[j]), k = 0; m; m = NEXT(m), k++ ) + mj[k] = m; + for ( i = 0; i < lj; i++ ) { + for ( di = mj[i]->dl, k = i+1; k < lj; k++ ) + if ( _dl_redble(di,mj[k]->dl,nv) ) break; + if ( k < lj ) mj[i] = 0; + } + hj = BDY(gh[j])->dl; + _NEWDL(dl,nv); d = dl->d; + r0 = r; + for ( i = 0; i < lj; i++ ) { + if ( mj[i] && !dl_equal(nv,di=mj[i]->dl,hj) ) { + for ( k = 0, t = 0; k < nv; k++ ) { + d[k] = hj->d[k]-di->d[k]; + t += d[k]; + } + dl->td = t; +#if 1 + for ( rt = r0; rt; rt = NEXT(rt) ) { + dlt = (DL)BDY(rt); + if ( dlt->td != dl->td ) continue; + for ( dt = dlt->d, k = 0; k < nv; k++ ) + if ( d[k] != dt[k] ) break; + if ( k == nv ) break; + } +#else + rt = 0; +#endif + if ( !rt ) { + MKNODE(r1,dl,r); r = r1; + _NEWDL(dl,nv); d = dl->d; + } + } + } + } + for ( rt = r; rt; rt = NEXT(rt) ) { + dl = (DL)BDY(rt); d = dl->d; + ri = 0; + for ( k = nv-1; k >= 0; k-- ) { + STOQ(d[k],q); + MKNODE(r1,q,ri); ri = r1; + } + MKNODE(r1,0,ri); MKLIST(l,r1); + BDY(rt) = (pointer)l; + } + return r; +} + +int comp_bits_divisible(int *a,int *b,int n) +{ + int bpi,i,wi,bi; + + bpi = (sizeof(int)/sizeof(char))*8; + for ( i = 0; i < n; i++ ) { + wi = i/bpi; bi = i%bpi; + if ( !(a[wi]&(1< bb ) return 1; + else if ( ba < bb ) return -1; + } + return 0; +} + +NODE mono_raddec(NODE ideal) +{ + DP p; + int nv,w,i,bpi,di,c,len; + int *d,*s,*u,*new; + NODE t,t1,v,r,rem,prev; + + if( !ideal ) return 0; + p = (DP)BDY(ideal); + nv = NV(p); + bpi = (sizeof(int)/sizeof(char))*8; + w = (nv+(bpi-1))/bpi; + d = p->body->dl->d; + if ( !NEXT(ideal) ) { + for ( t = 0, i = nv-1; i >= 0; i-- ) { + if ( d[i] ) { + s = (int *)CALLOC(w,sizeof(int)); + s[i/bpi] |= 1<<(i%bpi); + MKNODE(t1,s,t); + t = t1; + } + } + return t; + } + rem = mono_raddec(NEXT(ideal)); + r = 0; + len = w*sizeof(int); + u = (int *)CALLOC(w,sizeof(int)); + for ( i = nv-1; i >= 0; i-- ) { + if ( d[i] ) { + for ( t = rem; t; t = NEXT(t) ) { + bcopy((char *)BDY(t),(char *)u,len); + u[i/bpi] |= 1<<(i%bpi); + for ( v = r; v; v = NEXT(v) ) { + if ( comp_bits_divisible(u,(int *)BDY(v),nv) ) break; + } + if ( v ) continue; + for ( v = r, prev = 0; v; v = NEXT(v) ) { + if ( comp_bits_divisible((int *)BDY(v),u,nv) ) { + if ( prev ) NEXT(prev) = NEXT(v); + else r = NEXT(r); + } else prev =v; + } + for ( v = r, prev = 0; v; prev = v, v = NEXT(v) ) { + if ( comp_bits_lex(u,(int *)BDY(v),nv) < 0 ) break; + } + new = (int *)CALLOC(w,sizeof(int)); + bcopy((char *)u,(char *)new,len); + MKNODE(t1,new,v); + if ( prev ) NEXT(prev) = t1; + else r = t1; + } + } + } + return r; }