=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2018/builtin/dp-supp.c,v retrieving revision 1.4 retrieving revision 1.17 diff -u -p -r1.4 -r1.17 --- OpenXM_contrib2/asir2018/builtin/dp-supp.c 2019/09/04 01:12:02 1.4 +++ OpenXM_contrib2/asir2018/builtin/dp-supp.c 2020/12/15 07:40:09 1.17 @@ -45,7 +45,7 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2018/builtin/dp-supp.c,v 1.3 2019/08/21 00:37:47 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2018/builtin/dp-supp.c,v 1.16 2020/12/05 03:27:20 noro Exp $ */ #include "ca.h" #include "base.h" @@ -315,6 +315,21 @@ int _dl_redble(DL d1,DL d2,int nvar) return 1; } +int _dl_redble_ext(DL d1,DL d2,DL quo,int nvar) +{ + int i; + + if ( (quo->td = d2->td-d1->td) < 0 ) + return 0; + for ( i = 0; i < nvar; i++ ) + if ( (quo->d[i] = d2->d[i]-d1->d[i]) < 0 ) + break; + if ( i < nvar ) + return 0; + else + return 1; +} + void insert_to_node(DL d,NODE *n,int nvar) { DL d1; @@ -904,7 +919,47 @@ void dpm_red(DPM p0,DPM p1,DPM p2,DPM *head,DPM *rest, *head = h; *rest = r; *dnp = (P)c2; } +void dpm_red2(DPM p1,DPM p2,DPM *rest,P *dnp,DP *multp) +{ + int i,n,pos; + DL d1,d2,d; + MP m; + DP s,ms; + DPM t,r,h,u,w; + Z c,c1,c2,gn; + P g,a; + P p[2]; + n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; pos = BDY(p1)->pos; + if ( pos != BDY(p2)->pos ) + error("dpm_red : cannot happen"); + 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 = (Z)BDY(p1)->c; c2 = (Z)BDY(p2)->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 = (Z)a; divsp(CO,(P)c2,g,&a); c2 = (Z)a; + } else if ( dp_fcoeffs ) { + /* do nothing */ + } else if ( INT(c1) && INT(c2) ) { + gcdz(c1,c2,&gn); + if ( !UNIQ(gn) ) { + divsz(c1,gn,&c); c1 = c; + divsz(c2,gn,&c); c2 = c; + } + } else { + ezgcdpz(CO,(P)c1,(P)c2,&g); + divsp(CO,(P)c1,g,&a); c1 = (Z)a; divsp(CO,(P)c2,g,&a); c2 = (Z)a; + add_denomlist(g); + } + NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td; + *multp = s; + chsgnd(s,&ms); mulobjdpm(CO,(Obj)ms,p2,&u); mulobjdpm(CO,(Obj)c2,p1,&w); adddpm(CO,w,u,&r); + *rest = r; *dnp = (P)c2; +} + /* * m-reduction by a marked poly * do content reduction over Z or Q(x,...) @@ -1204,6 +1259,56 @@ void dp_removecont2(DP p1,DP p2,DP *r1p,DP *r2p,Z *con *r2p = 0; } +void dpm_removecont2(DPM p1,DPM p2,DPM *r1p,DPM *r2p,Z *contp) +{ + struct oVECT v; + int i,n1,n2,n; + DMM m,m0,t; + Z *w; + Z 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 = (Z *)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] = (Z)m->c; + if ( p2 ) + for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) w[i] = (Z)m->c; + h = w[0]; removecont_array((P *)w,n,1); divsz(h,w[0],contp); + i = 0; + if ( p1 ) { + for ( m0 = 0, t = BDY(p1); i < n1; i++, t = NEXT(t) ) { + NEXTDMM(m0,m); m->c = (Obj)w[i]; m->dl = t->dl; m->pos = t->pos; + } + NEXT(m) = 0; + MKDPM(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) ) { + NEXTDMM(m0,m); m->c = (Obj)w[i]; m->dl = t->dl; m->pos = t->pos; + } + NEXT(m) = 0; + MKDPM(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) @@ -1380,17 +1485,41 @@ last: return q; } -// struct oEGT egc; +struct oEGT egred; +void mulcmp(Obj c,MP m); +void mulcdmm(Obj c,DMM m); + +DP appendd(DP d,DP m) +{ + MP t; + + if ( !d ) return m; + for ( t = BDY(d); NEXT(t); t = NEXT(t) ); + NEXT(t) = BDY(m); + return d; +} + +DPM appenddpm(DPM d,DPM m) +{ + DMM t; + + if ( !d ) return m; + for ( t = BDY(d); NEXT(t); t = NEXT(t) ); + NEXT(t) = BDY(m); + return d; +} + DP *dpm_nf_and_quotient(NODE b,DPM g,VECT psv,DPM *rp,P *dnp) { - DPM u,p,d,s,t; - DP dmy,mult; + DPM u,p,s,t,d; + DP dmy,mult,zzz; DPM *ps; DP *q; NODE l; DMM m,mr; - int i,n,j,len; + MP mp; + int i,n,j,len,nv; int *wb; int sugar,psugar,multiple; P nm,tnm1,dn,tdn,tdn1; @@ -1401,6 +1530,7 @@ DP *dpm_nf_and_quotient(NODE b,DPM g,VECT psv,DPM *rp, if ( !g ) { *rp = 0; *dnp = dn; return 0; } + nv = NV(g); ps = (DPM *)BDY(psv); len = psv->len; if ( b ) { @@ -1421,17 +1551,17 @@ DP *dpm_nf_and_quotient(NODE b,DPM g,VECT psv,DPM *rp, for ( d = 0; g; ) { for ( u = 0, i = 0; i < n; i++ ) { if ( dpm_redble(g,p = ps[wb[i]]) ) { - dpm_red(d,g,p,&t,&u,&tdn,&mult); +// get_eg(&eg0); + dpm_red2(g,p,&u,&tdn,&mult); +// get_eg(&eg1); add_eg(&egred,&eg0,&eg1); psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; sugar = MAX(sugar,psugar); -// get_eg(&eg0); for ( j = 0; j < len; j++ ) { - muldc(CO,q[j],(Obj)tdn,&dmy); q[j] = dmy; + if ( q[j] ) { mulcmp((Obj)tdn,BDY(q[j])); } } - addd(CO,q[wb[i]],mult,&dmy); q[wb[i]] = dmy; -// get_eg(&eg1); add_eg(&egc,&eg0,&eg1); + q[wb[i]] = appendd(q[wb[i]],mult); mulp(CO,dn,tdn,&tdn1); dn = tdn1; - d = t; + if ( d ) mulcdmm((Obj)tdn,BDY(d)); if ( !u ) goto last; break; } @@ -1441,7 +1571,7 @@ DP *dpm_nf_and_quotient(NODE b,DPM g,VECT psv,DPM *rp, } else { m = BDY(g); NEWDMM(mr); mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos; NEXT(mr) = 0; MKDPM(g->nv,mr,t); t->sugar = mr->dl->td; - adddpm(CO,d,t,&s); d = s; + d = appenddpm(d,t); dpm_rest(g,&t); g = t; } } @@ -1451,6 +1581,487 @@ last: return q; } +DPM dpm_nf_and_quotient2(NODE b,DPM g,VECT psv,DPM *rp,P *dnp) +{ + DPM u,p,s,t,d,q; + DP dmy,mult,zzz; + DPM *ps; + NODE l; + DMM mr0,mq0,mr,mq,m; + MP mp; + int i,n,j,len,nv; + int *wb; + int sugar,psugar,multiple; + P nm,tnm1,dn,tdn,tdn1; + Q cont; + Obj c1; + struct oEGT eg0,eg1; + + dn = (P)ONE; + if ( !g ) { + *rp = 0; *dnp = dn; return 0; + } + nv = NV(g); + ps = (DPM *)BDY(psv); + len = psv->len; + if ( b ) { + 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] = ZTOS((Q)BDY(l)); + } else { + wb = (int *)ALLOCA(len*sizeof(int)); + for ( i = j = 0; i < len; i++ ) + if ( ps[i] ) wb[j++] = i; + n = j; + } + sugar = g->sugar; + mq0 = 0; + mr0 = 0; + for ( ; g; ) { + for ( u = 0, i = 0; i < n; i++ ) { + if ( dpm_redble(g,p = ps[wb[i]]) ) { + dpm_red2(g,p,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + for ( m = mq0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + for ( m = mr0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + NEXTDMM(mq0,mq); + mq->c = BDY(mult)->c; mq->dl = BDY(mult)->dl; mq->pos = wb[i]+1; + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + if ( !u ) goto last; + break; + } + } + if ( u ) { + g = u; + } else { + m = BDY(g); + NEXTDMM(mr0,mr); + mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos; + dpm_rest(g,&t); g = t; + } + } +last: + if ( mr0 ) { + NEXT(mr) = 0; + MKDPM(nv,mr0,d); d->sugar = sugar; + } else + d = 0; + if ( mq0 ) { + NEXT(mq) = 0; + MKDPM(nv,mq0,q); q->sugar = sugar; + } else + q = 0; + *rp = d; *dnp = dn; + return q; +} + +DPM dpm_nf_and_quotient3(DPM g,VECT psv,DPM *rp,P *dnp) +{ + DPM u,p,s,t,d,q; + DP dmy,mult,zzz; + DPM *ps; + NODE2 nd; + DMM mr0,mq0,mr,mq,m; + MP mp; + int i,n,j,len,nv,pos,max; + int *wb; + int sugar,psugar,multiple; + P nm,tnm1,dn,tdn,tdn1; + Q cont; + Obj c1; + struct oEGT eg0,eg1; + + dn = (P)ONE; + if ( !g ) { + *rp = 0; *dnp = dn; return 0; + } + nv = NV(g); + sugar = g->sugar; + mq0 = 0; + mr0 = 0; + max = psv->len; + for ( ; g; ) { + pos = BDY(g)->pos; + u = 0; + if ( pos < max ) { + nd = (NODE2)BDY(psv)[pos]; + for ( u = 0; nd; nd = NEXT(nd) ) { + if ( dpm_redble(g,p = (DPM)(nd->body1)) ) { + dpm_red2(g,p,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + if ( !UNIZ(tdn) ) { + for ( m = mq0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + for ( m = mr0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + } + NEXTDMM(mq0,mq); + mq->c = BDY(mult)->c; mq->dl = BDY(mult)->dl; mq->pos = (long)nd->body2; + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + if ( !u ) goto last; + break; + } + } + } + if ( u ) { + g = u; + } else { + m = BDY(g); + NEXTDMM(mr0,mr); + mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos; + dpm_rest(g,&t); g = t; + } + } +last: + if ( mr0 ) { + NEXT(mr) = 0; + MKDPM(nv,mr0,d); d->sugar = sugar; + } else + d = 0; + if ( mq0 ) { + NEXT(mq) = 0; + MKDPM(nv,mq0,q); q->sugar = sugar; + } else + q = 0; + *rp = d; *dnp = dn; + return q; +} + +DPM dpm_nf_and_quotient4(DPM g,DPM *ps,VECT psiv,DPM head,DPM *rp,P *dnp) +{ + DPM u,p,s,t,d,q; + DP dmy,mult,zzz; + NODE nd; + DMM mr0,mq0,mr,mq,m; + MP mp; + int i,n,j,len,nv,pos,max; + int *wb; + int sugar,psugar,multiple; + P nm,tnm1,dn,tdn,tdn1,c; + Q cont; + Obj c1; + struct oEGT eg0,eg1; + + dn = (P)ONE; + if ( !g ) { + *rp = 0; *dnp = dn; return 0; + } + nv = NV(g); + sugar = g->sugar; + mq0 = 0; + if ( head ) { + for ( m = BDY(head); m; m = NEXT(m) ) { + NEXTDMM(mq0,mq); + mq->c = m->c; mq->dl = m->dl; mq->pos = m->pos; + } + } + mr0 = 0; + max = psiv->len; + for ( ; g; ) { + pos = BDY(g)->pos; + u = 0; + if ( pos < max ) { + nd = (NODE)BDY(psiv)[pos]; + for ( u = 0; nd; nd = NEXT(nd) ) { + if ( dpm_redble(g,p = ps[(long)(BDY(nd))-1]) ) { + dpm_red2(g,p,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + if ( !UNIZ(tdn) ) { + for ( m = mq0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + for ( m = mr0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + } + NEXTDMM(mq0,mq); + mq->c = BDY(mult)->c; + mq->dl = BDY(mult)->dl; mq->pos = (long)BDY(nd); + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + if ( !u ) goto last; + break; + } + } + } + if ( u ) { + g = u; + } else { + m = BDY(g); + NEXTDMM(mr0,mr); + mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos; + dpm_rest(g,&t); g = t; + } + } +last: + if ( mr0 ) { + NEXT(mr) = 0; + MKDPM(nv,mr0,d); d->sugar = sugar; + } else + d = 0; + if ( mq0 ) { + NEXT(mq) = 0; + MKDPM(nv,mq0,q); q->sugar = sugar; + } else + q = 0; + *rp = d; *dnp = dn; + return q; +} + +/* an intermediate version for calling from the user language */ +/* XXX : i, j must be positive */ + +DPM dpm_sp_nf_asir(VECT psv,int i,int j,DPM *nf) +{ + DPM *ps; + int n,k,nv,s1,s2,sugar,max,pos,psugar; + DPM g,u,p,d,q,t; + DMM mq0,mq,mr0,mr,m; + DP mult,t1,t2; + P dn,tdn,tdn1; + NODE nd; + Obj c1; + + ps = (DPM *)BDY(psv); + n = psv->len; + nv = ps[1]->nv; + dpm_sp(ps[i],ps[j],&g,&t1,&t2); + mq0 = 0; + NEXTDMM(mq0,mq); mq->c = BDY(t1)->c; mq->pos = i; mq->dl = BDY(t1)->dl; + NEXTDMM(mq0,mq); chsgnp((P)BDY(t2)->c,(P *)&mq->c); mq->pos = j; mq->dl = BDY(t2)->dl; + + if ( !g ) { + NEXT(mq) = 0; + MKDPM(nv,mq0,d); + s1 = BDY(t1)->dl->td + ps[i]->sugar; + s2 = BDY(t2)->dl->td + ps[j]->sugar; + d->sugar = MAX(s1,s2); + *nf = 0; + return d; + } + + dn = (P)ONE; + sugar = g->sugar; + mr0 = 0; + while ( g ) { + pos = BDY(g)->pos; + for ( u = 0, k = 1; k < n; k++ ) { + if ( (p=ps[k])!=0 && pos == BDY(p)->pos && dpm_redble(g,p) ) { + dpm_red2(g,p,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + if ( !UNIZ(tdn) ) { + for ( m = mq0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + for ( m = mr0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + } + NEXTDMM(mq0,mq); + chsgnp((P)BDY(mult)->c,(P *)&mq->c); + mq->dl = BDY(mult)->dl; mq->pos = k; + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + if ( !u ) goto last; + break; + } + } + if ( u ) { + g = u; + } else { + m = BDY(g); + NEXTDMM(mr0,mr); + mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos; + dpm_rest(g,&t); g = t; + } + } +last: + if ( mr0 ) { + NEXT(mr) = 0; MKDPM(nv,mr0,d); d->sugar = sugar; + } else + d = 0; + NEXT(mq) = 0; MKDPM(nv,mq0,q); q->sugar = sugar; + *nf = d; + return q; +} + +DPM dpm_sp_nf(VECT psv,VECT psiv,int i,int j,DPM *nf) +{ + DPM *ps; + int n,nv,s1,s2,sugar,max,pos,psugar; + DPM g,u,p,d,q,t; + DMM mq0,mq,mr0,mr,m; + DP mult,t1,t2; + P dn,tdn,tdn1; + NODE nd; + Obj c1; + + ps = (DPM *)BDY(psv); + n = psv->len; + nv = ps[i]->nv; + dpm_sp(ps[i],ps[j],&g,&t1,&t2); + mq0 = 0; + NEXTDMM(mq0,mq); mq->c = BDY(t1)->c; mq->pos = i; mq->dl = BDY(t1)->dl; + NEXTDMM(mq0,mq); chsgnp((P)BDY(t2)->c,(P *)&mq->c); mq->pos = j; mq->dl = BDY(t2)->dl; + + if ( !g ) { + NEXT(mq) = 0; + MKDPM(nv,mq0,d); + s1 = BDY(t1)->dl->td + ps[i]->sugar; + s2 = BDY(t2)->dl->td + ps[j]->sugar; + d->sugar = MAX(s1,s2); + *nf = 0; + return d; + } + + dn = (P)ONE; + sugar = g->sugar; + mr0 = 0; + max = psiv->len; + while ( g ) { + pos = BDY(g)->pos; + u = 0; + if ( pos < max ) { + nd = (NODE)BDY(psiv)[pos]; + for ( u = 0; nd; nd = NEXT(nd) ) { + if ( dpm_redble(g,p = ps[(long)(BDY(nd))]) ) { + dpm_red2(g,p,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + if ( !UNIZ(tdn) ) { + for ( m = mq0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + for ( m = mr0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + } + NEXTDMM(mq0,mq); + chsgnp((P)BDY(mult)->c,(P *)&mq->c); + mq->dl = BDY(mult)->dl; mq->pos = (long)BDY(nd); + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + if ( !u ) goto last; + break; + } + } + } + if ( u ) { + g = u; + } else { + m = BDY(g); + NEXTDMM(mr0,mr); + mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos; + dpm_rest(g,&t); g = t; + } + } +last: + if ( mr0 ) { + NEXT(mr) = 0; MKDPM(nv,mr0,d); d->sugar = sugar; + } else + d = 0; + NEXT(mq) = 0; MKDPM(nv,mq0,q); q->sugar = sugar; + *nf = d; + return q; +} + +/* psiv is a vector of lists of Z */ + +DPM dpm_sp_nf_zlist(VECT psv,VECT psiv,int i,int j,int top,DPM *nf) +{ + DPM *ps; + int n,nv,s1,s2,sugar,max,pos,psugar; + DPM g,u,p,d,q,t; + DMM mq0,mq,mr0,mr,m; + DP mult,t1,t2; + P dn,tdn,tdn1; + NODE nd; + Obj c1; + + ps = (DPM *)BDY(psv); + n = psv->len; + nv = ps[i]->nv; + dpm_sp(ps[i],ps[j],&g,&t1,&t2); + mq0 = 0; + NEXTDMM(mq0,mq); mq->c = BDY(t1)->c; mq->pos = i; mq->dl = BDY(t1)->dl; + NEXTDMM(mq0,mq); chsgnp((P)BDY(t2)->c,(P *)&mq->c); mq->pos = j; mq->dl = BDY(t2)->dl; + + if ( !g ) { + NEXT(mq) = 0; + MKDPM(nv,mq0,d); + s1 = BDY(t1)->dl->td + ps[i]->sugar; + s2 = BDY(t2)->dl->td + ps[j]->sugar; + d->sugar = MAX(s1,s2); + *nf = 0; + return d; + } + + dn = (P)ONE; + sugar = g->sugar; + mr0 = 0; + max = psiv->len; + while ( g ) { + pos = BDY(g)->pos; + u = 0; + if ( pos < max ) { + nd = BDY((LIST)BDY(psiv)[pos]); + for ( u = 0; nd; nd = NEXT(nd) ) { + if ( dpm_redble(g,p = ps[ZTOS((Q)BDY(nd))]) ) { + dpm_red2(g,p,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + if ( !UNIZ(tdn) ) { + for ( m = mq0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + for ( m = mr0; m; m = NEXT(m) ) { + arf_mul(CO,(Obj)tdn,m->c,&c1); m->c = c1; + } + } + NEXTDMM(mq0,mq); + chsgnp((P)BDY(mult)->c,(P *)&mq->c); + mq->dl = BDY(mult)->dl; mq->pos = ZTOS((Q)BDY(nd)); + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + if ( !u ) goto last; + break; + } + } + } + if ( u ) { + g = u; + } else if ( !top ) { + m = BDY(g); + NEXTDMM(mr0,mr); + mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos; + dpm_rest(g,&t); g = t; + } else { + *nf = g; + if ( mq0 ) { + NEXT(mq) = 0; MKDPM(nv,mq0,q); q->sugar = sugar; + } else + q = 0; + return q; + } + } +last: + if ( mr0 ) { + NEXT(mr) = 0; MKDPM(nv,mr0,d); d->sugar = sugar; + } else + d = 0; + NEXT(mq) = 0; MKDPM(nv,mq0,q); q->sugar = sugar; + *nf = d; + 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; @@ -1613,15 +2224,15 @@ void dpm_nf_z(NODE b,DPM g,VECT psv,int full,int multi for ( d = 0; g; ) { for ( u = 0, i = 0; i < n; i++ ) { if ( (p=ps[wb[i]])!=0 && dpm_redble(g,p) ) { - dpm_red(d,g,p,&t,&u,&dmy,&dmy1); + dpm_red2(g,p,&u,&dmy,&dmy1); psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; sugar = MAX(sugar,psugar); + if ( d ) mulcdmm((Obj)dmy,BDY(d)); if ( !u ) { if ( d ) d->sugar = sugar; *rp = d; return; } - d = t; break; } } @@ -1647,7 +2258,7 @@ void dpm_nf_z(NODE b,DPM g,VECT psv,int full,int multi } else { m = BDY(g); NEWDMM(mr); mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos; NEXT(mr) = 0; MKDPM(g->nv,mr,t); t->sugar = mr->dl->td; - adddpm(CO,d,t,&s); d = s; + d = appenddpm(d,t); dpm_rest(g,&t); g = t; } } @@ -1707,6 +2318,29 @@ void dpm_split(DPM p,int s,DPM *up,DPM *lo) } } +/* extract the component in DP of position s */ +void dpm_extract(DPM p,int s,DP *r) +{ + DMM m; + MP mu0,mu; + DP t; + + if ( !p ) { + *r = 0; return; + } + for ( m = BDY(p), mu0 = 0; m; m = NEXT(m) ) { + if ( m->pos == s ) { + NEXTMP(mu0,mu); + mu->dl = m->dl; mu->c = m->c; + } + } + if ( mu0 ) { + NEXT(mu) = 0; MKDP(p->nv,mu0,t); t->sugar = p->sugar; + *r = t; + } else + *r = 0; +} + /* nf computation over a field */ void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp) @@ -2110,51 +2744,108 @@ int create_order_spec(VL vl,Obj obj,struct order_spec spec->ord.simple = ZTOS((Q)obj); return 1; } else if ( OID(obj) == O_LIST ) { - /* module order; obj = [0|1,w,ord] or [0|1,ord] */ + /* module order */ node = BDY((LIST)obj); if ( !BDY(node) || NUM(BDY(node)) ) { switch ( length(node) ) { - case 2: + case 2: /* [n,ord] */ create_order_spec(0,(Obj)BDY(NEXT(node)),&spec); spec->id += 256; spec->obj = obj; spec->top_weight = 0; spec->module_rank = 0; spec->module_top_weight = 0; - spec->ispot = (BDY(node)!=0); - if ( spec->ispot ) { - n = ZTOS((Q)BDY(node)); - if ( n < 0 ) - spec->pot_nelim = -n; - else + spec->module_ordtype = ZTOS((Z)BDY(node)); + if ( spec->module_ordtype < 0 ) { + spec->pot_nelim = -spec->module_ordtype; + spec->module_ordtype = 1; + } else spec->pot_nelim = 0; - } break; - case 3: - create_order_spec(0,(Obj)BDY(NEXT(NEXT(node))),&spec); - spec->id += 256; spec->obj = obj; - spec->ispot = (BDY(node)!=0); - node = NEXT(node); - if ( !BDY(node) || OID(BDY(node)) != O_LIST ) - error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight"); - wpair = BDY((LIST)BDY(node)); - if ( length(wpair) != 2 ) - error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight"); + case 3: /* [n,[mlist1,mlist2,...],ord] or [n,[wv,wm],ord] */ + spec->module_ordtype = ZTOS((Z)BDY(node)); + if ( spec->module_ordtype < 0 ) { + spec->pot_nelim = -spec->module_ordtype; + spec->module_ordtype = 1; + } else + spec->pot_nelim = 0; - wp = BDY(wpair); - wm = BDY(NEXT(wpair)); - if ( !wp || OID(wp) != O_LIST || !wm || OID(wm) != O_LIST ) - error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight"); - spec->nv = length(BDY((LIST)wp)); - spec->top_weight = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int)); - for ( i = 0, t = BDY((LIST)wp); i < spec->nv; t = NEXT(t), i++ ) - spec->top_weight[i] = ZTOS((Q)BDY(t)); + if ( spec->module_ordtype == 3 ) { /* schreyer order */ + Obj baseobj; + struct order_spec *basespec; + int len; + NODE in; + LIST *la; + DMMstack stack; + DMMstack push_schreyer_order(LIST l,DMMstack s); - spec->module_rank = length(BDY((LIST)wm)); - spec->module_top_weight = (int *)MALLOC_ATOMIC(spec->module_rank*sizeof(int)); - for ( i = 0, t = BDY((LIST)wm); i < spec->module_rank; t = NEXT(t), i++ ) - spec->module_top_weight[i] = ZTOS((Q)BDY(t)); + spec->id = 300; spec->obj = obj; + node = NEXT(node); + if ( !BDY(node) || OID(BDY(node)) != O_LIST ) + error("create_order_spec : [mlist1,mlist,...] must be specified for defining a schreyer order"); + stack = 0; + in = BDY((LIST)BDY(node)); + len = length(in); + la = (LIST *)MALLOC(len*sizeof(LIST)); + for ( i = 0; i < len; i++, in = NEXT(in) ) la[i] = (LIST)(BDY(in)); + for ( i = len-1; i >= 0; i-- ) stack = push_schreyer_order(la[i],stack); + spec->dmmstack = stack; + + node = NEXT(node); + baseobj = (Obj)BDY(node); + create_order_spec(0,baseobj,&basespec); + basespec->obj = baseobj; + spec->base = basespec; + } else if ( spec->module_ordtype == 4 ) { /* POT with base order [n,bord,ord] */ + NODE base_ord; + int rank; + + create_order_spec(0,(Obj)BDY(NEXT(NEXT(node))),&spec); + spec->id += 256; spec->obj = obj; + spec->top_weight = 0; + spec->module_rank = 0; + spec->module_top_weight = 0; + spec->pot_nelim = 0; + spec->module_ordtype = 4; + node = NEXT(node); + if ( !BDY(node) || OID(BDY(node)) != O_LIST ) + error("create_order_spec : a permitation list must be specified"); + base_ord = BDY((LIST)BDY(node)); + spec->module_rank = rank = length(base_ord); + spec->module_base_ord = (int *)MALLOC_ATOMIC((rank+1)*sizeof(int)); + for ( i = 1, t = base_ord; i <= rank; i++, t = NEXT(t) ) + spec->module_base_ord[ZTOS((Q)BDY(t))] = i; + break; + } else { /* weighted order * [n,[wv,wm],ord] */ + int ordtype; + + ordtype = spec->module_ordtype; + create_order_spec(0,(Obj)BDY(NEXT(NEXT(node))),&spec); + spec->module_ordtype = ordtype; + spec->id += 256; spec->obj = obj; + node = NEXT(node); + if ( !BDY(node) || OID(BDY(node)) != O_LIST ) + error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight"); + wpair = BDY((LIST)BDY(node)); + if ( length(wpair) != 2 ) + error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight"); + + wp = BDY(wpair); + wm = BDY(NEXT(wpair)); + if ( !wp || OID(wp) != O_LIST || !wm || OID(wm) != O_LIST ) + error("create_order_spec : [weight_for_poly,weight_for_modlue] must be specified as a module topweight"); + spec->nv = length(BDY((LIST)wp)); + spec->top_weight = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int)); + for ( i = 0, t = BDY((LIST)wp); i < spec->nv; t = NEXT(t), i++ ) + spec->top_weight[i] = ZTOS((Q)BDY(t)); + + spec->module_rank = length(BDY((LIST)wm)); + spec->module_top_weight = (int *)MALLOC_ATOMIC(spec->module_rank*sizeof(int)); + for ( i = 0, t = BDY((LIST)wm); i < spec->module_rank; t = NEXT(t), i++ ) + spec->module_top_weight[i] = ZTOS((Q)BDY(t)); + } break; + default: error("create_order_spec : invalid arguments for module order"); } @@ -2522,6 +3213,56 @@ void create_modorder_spec(int id,LIST shift,struct mod * */ +void dpm_homo(DPM p,DPM *rp) +{ + DMM m,mr,mr0,t; + int i,n,nv,td; + DL dl,dlh; + + if ( !p ) + *rp = 0; + else { + n = p->nv; nv = n + 1; + m = BDY(p); + td = 0; + for ( t = m; t; t = NEXT(t) ) + if ( m->dl->td > td ) td = m->dl->td; + for ( mr0 = 0; m; m = NEXT(m) ) { + NEXTDMM(mr0,mr); mr->c = m->c; mr->pos = m->pos; + dl = m->dl; + mr->dl = dlh = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int)); + dlh->td = td; + for ( i = 0; i < n; i++ ) + dlh->d[i] = dl->d[i]; + dlh->d[n] = td - dl->td; + } + NEXT(mr) = 0; MKDPM(nv,mr0,*rp); (*rp)->sugar = p->sugar; + } +} + +void dpm_dehomo(DPM p,DPM *rp) +{ + DMM m,mr,mr0; + int i,n,nv; + DL dl,dlh; + + if ( !p ) + *rp = 0; + else { + n = p->nv; nv = n - 1; + m = BDY(p); + for ( mr0 = 0; m; m = NEXT(m) ) { + NEXTDMM(mr0,mr); mr->c = m->c; mr->pos = m->pos; + dlh = m->dl; + mr->dl = dl = (DL)MALLOC_ATOMIC((nv+1)*sizeof(int)); + dl->td = dlh->td - dlh->d[nv]; + for ( i = 0; i < nv; i++ ) + dl->d[i] = dlh->d[i]; + } + NEXT(mr) = 0; MKDPM(nv,mr0,*rp); (*rp)->sugar = p->sugar; + } +} + void dp_homo(DP p,DP *rp) { MP m,mr,mr0; @@ -2569,6 +3310,7 @@ void dp_dehomo(DP p,DP *rp) } } + void dp_mod(DP p,int mod,NODE subst,DP *rp) { MP m,mr,mr0; @@ -2596,6 +3338,29 @@ void dp_mod(DP p,int mod,NODE subst,DP *rp) } } +void dpm_mod(DPM p,int mod,DPM *rp) +{ + DMM m,mr,mr0; + P t; + V v; + NODE tn; + + if ( !p ) + *rp = 0; + else { + for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { + ptomp(mod,(P)m->c,&t); + if ( t ) { + NEXTDMM(mr0,mr); mr->c = (Obj)t; mr->dl = m->dl; + } + } + if ( mr0 ) { + NEXT(mr) = 0; MKDPM(p->nv,mr0,*rp); (*rp)->sugar = p->sugar; + } else + *rp = 0; + } +} + void dp_rat(DP p,DP *rp) { MP m,mr,mr0; @@ -2624,11 +3389,12 @@ void homogenize_order(struct order_spec *old,int n,str struct weight_or_block *owb,*nwb; *newp = new = (struct order_spec *)MALLOC(sizeof(struct order_spec)); + bcopy((char *)old,(char *)new,sizeof(struct order_spec)); switch ( old->id ) { case 0: switch ( old->ord.simple ) { case 0: - new->id = 0; new->ord.simple = 0; break; + break; case 1: l = (struct order_pair *) MALLOC_ATOMIC(2*sizeof(struct order_pair)); @@ -2639,12 +3405,12 @@ void homogenize_order(struct order_spec *old,int n,str new->ord.block.length = 2; new->nv = n+1; break; case 2: - new->id = 0; new->ord.simple = 1; break; + new->ord.simple = 1; break; case 3: case 4: case 5: - new->id = 0; new->ord.simple = old->ord.simple+3; + new->ord.simple = old->ord.simple+3; dp_nelim = n-1; break; case 6: case 7: case 8: case 9: - new->id = 0; new->ord.simple = old->ord.simple; break; + break; default: error("homogenize_order : invalid input"); } @@ -2655,10 +3421,9 @@ void homogenize_order(struct order_spec *old,int n,str 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 = old->id; new->nv = n+1; + new->nv = n+1; new->ord.block.order_pair = l; new->ord.block.length = length+1; - new->ispot = old->ispot; break; case 2: case 258: nv = old->nv; row = old->ord.matrix.row; @@ -2670,9 +3435,8 @@ void homogenize_order(struct order_spec *old,int n,str newm[i+1][j] = oldm[i][j]; newm[i+1][j] = 0; } - new->id = old->id; new->nv = nv+1; + 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; @@ -2708,17 +3472,15 @@ void homogenize_order(struct order_spec *old,int n,str (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; + break; case 1: l = (struct order_pair *) MALLOC_ATOMIC(2*sizeof(struct order_pair)); @@ -2729,11 +3491,10 @@ void homogenize_order(struct order_spec *old,int n,str new->ord.block.length = 2; new->nv = n+1; break; case 2: - new->id = 256; new->ord.simple = 1; break; + new->ord.simple = 1; break; default: error("homogenize_order : invalid input"); } - new->ispot = old->ispot; break; default: error("homogenize_order : invalid input"); @@ -2899,6 +3660,23 @@ void dpm_rest(DPM p,DPM *rp) if ( *rp ) (*rp)->sugar = p->sugar; } +} + +int dp_getdeg(DP p) +{ + int max,n,i; + MP m; + int *d; + + if ( !p ) return 0; + n = p->nv; + max = 0; + for ( m = BDY(p); m; m = NEXT(m) ) { + d = m->dl->d; + for ( i = 0; i < n; i++ ) + if ( d[i] > max ) max = d[i]; + } + return max; } int dpm_getdeg(DPM p,int *r)