=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2018/builtin/dp-supp.c,v retrieving revision 1.1 retrieving revision 1.8 diff -u -p -r1.1 -r1.8 --- OpenXM_contrib2/asir2018/builtin/dp-supp.c 2018/09/19 05:45:05 1.1 +++ OpenXM_contrib2/asir2018/builtin/dp-supp.c 2019/11/01 04:28:52 1.8 @@ -45,7 +45,7 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM$ + * $OpenXM: OpenXM_contrib2/asir2018/builtin/dp-supp.c,v 1.7 2019/10/11 03:45:56 noro Exp $ */ #include "ca.h" #include "base.h" @@ -147,7 +147,7 @@ void dp_ptozp2(DP p0,DP p1,DP *hp,DP *rp) *hp = h; *rp = r; } -void dpm_ptozp(DPM p,DPM *rp) +void dpm_ptozp(DPM p,Z *cont,DPM *rp) { DMM m,mr,mr0; int i,n; @@ -155,9 +155,9 @@ void dpm_ptozp(DPM p,DPM *rp) Z dvr; P t; - if ( !p ) - *rp = 0; - else { + if ( !p ) { + *rp = 0; *cont = ONE; + } 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++ ) @@ -171,6 +171,7 @@ void dpm_ptozp(DPM p,DPM *rp) NEXTDMM(mr0,mr); divsp(CO,(P)m->c,(P)dvr,(P *)&mr->c); mr->dl = m->dl; mr->pos = m->pos; } NEXT(mr) = 0; MKDPM(p->nv,mr0,*rp); (*rp)->sugar = p->sugar; + *cont = dvr; } } @@ -178,8 +179,9 @@ void dpm_ptozp2(DPM p0,DPM p1,DPM *hp,DPM *rp) { DPM t,s,h,r; DMM m,mr,mr0,m0; + Z cont; - adddpm(CO,p0,p1,&t); dpm_ptozp(t,&s); + adddpm(CO,p0,p1,&t); dpm_ptozp(t,&cont,&s); if ( !p0 ) { h = 0; r = s; } else if ( !p1 ) { @@ -238,7 +240,7 @@ void dp_idiv(DP p,Z c,DP *rp) else { for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { NEXTMP(mr0,mr); - divz((Z)(m->c),c,(Z *)&mr->c); + divsz((Z)(m->c),c,(Z *)&mr->c); mr->dl = m->dl; } NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); @@ -475,11 +477,11 @@ void heu_nezgcdnpz(VL vl,P *pl,int m,P *pr) while ( 1 ) { for ( i = 0, s1 = 0; i < m; i++ ) { - r = random(); UTOQ(r,rq); + r = random(); UTOZ(r,rq); mulp(vl,pl[i],(P)rq,&t); addp(vl,s1,t,&u); s1 = u; } for ( i = 0, s2 = 0; i < m; i++ ) { - r = random(); UTOQ(r,rq); + r = random(); UTOZ(r,rq); mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u; } ezgcdp(vl,s1,s2,&gcd); @@ -565,8 +567,8 @@ void dp_sp(DP p1,DP p2,DP *rp) if ( INT(c1) && INT(c2) ) { gcdz(c1,c2,&gn); if ( !UNIQ(gn) ) { - divz(c1,gn,&c); c1 = c; - divz(c2,gn,&c);c2 = c; + divsz(c1,gn,&c); c1 = c; + divsz(c2,gn,&c);c2 = c; } } @@ -595,7 +597,7 @@ void dp_sp(DP p1,DP p2,DP *rp) } } -void dpm_sp(DPM p1,DPM p2,DPM *rp) +void dpm_sp(DPM p1,DPM p2,DPM *rp,DP *mul1,DP *mul2) { int i,n,td; int *w; @@ -608,7 +610,7 @@ void dpm_sp(DPM p1,DPM p2,DPM *rp) n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; if ( BDY(p1)->pos != BDY(p2)->pos ) { - *rp = 0; + *mul1 = 0; *mul2 = 0; *rp = 0; return; } w = (int *)ALLOCA(n*sizeof(int)); @@ -623,19 +625,21 @@ void dpm_sp(DPM p1,DPM p2,DPM *rp) if ( INT(c1) && INT(c2) ) { gcdz(c1,c2,&gn); if ( !UNIQ(gn) ) { - divz(c1,gn,&c); c1 = c; - divz(c2,gn,&c);c2 = c; + divsz(c1,gn,&c); c1 = c; + divsz(c2,gn,&c);c2 = c; } } NEWMP(m); m->dl = d; m->c = (Obj)c2; NEXT(m) = 0; MKDP(n,m,s1); s1->sugar = d->td; mulobjdpm(CO,(Obj)s1,p1,&t); + *mul1 = s1; NEWDL(d,n); d->td = td - d2->td; for ( i = 0; i < n; i++ ) d->d[i] = w[i] - d2->d[i]; NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0; MKDP(n,m,s2); s2->sugar = d->td; mulobjdpm(CO,(Obj)s2,p2,&u); + *mul2 = s2; subdpm(CO,t,u,rp); if ( GenTrace ) { @@ -653,6 +657,32 @@ void dpm_sp(DPM p1,DPM p2,DPM *rp) } } +DP dpm_sp_hm(DPM p1,DPM p2) +{ + int i,n,td; + int *w; + DL d1,d2,d; + MP m; + DP s1; + + n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; + if ( BDY(p1)->pos != BDY(p2)->pos ) { + return 0; + } + w = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0, td = 0; i < n; 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++ ) + d->d[i] = w[i] - d1->d[i]; + + NEWMP(m); m->dl = d; m->c = (Obj)ONE; NEXT(m) = 0; + MKDP(n,m,s1); s1->sugar = d->td; + return s1; +} + void _dp_sp_dup(DP p1,DP p2,DP *rp) { int i,n,td; @@ -676,8 +706,8 @@ void _dp_sp_dup(DP p1,DP p2,DP *rp) if ( INT(c1) && INT(c2) ) { gcdz(c1,c2,&gn); if ( !UNIQ(gn) ) { - divz(c1,gn,&c); c1 = c; - divz(c2,gn,&c);c2 = c; + divsz(c1,gn,&c); c1 = c; + divsz(c2,gn,&c);c2 = c; } } @@ -789,6 +819,7 @@ void _dp_sp_mod(DP p1,DP p2,int mod,DP *rp) * do content reduction over Z or Q(x,...) * do nothing over finite fields * + * head+rest = dn*(p0+p1)+mult*p2 */ void dp_red(DP p0,DP p1,DP p2,DP *head,DP *rest,P *dnp,DP *multp) @@ -815,8 +846,8 @@ void dp_red(DP p0,DP p1,DP p2,DP *head,DP *rest,P *dnp } else if ( INT(c1) && INT(c2) ) { gcdz(c1,c2,&gn); if ( !UNIQ(gn) ) { - divz(c1,gn,&c); c1 = c; - divz(c2,gn,&c); c2 = c; + divsz(c1,gn,&c); c1 = c; + divsz(c2,gn,&c); c2 = c; } } else { ezgcdpz(CO,(P)c1,(P)c2,&g); @@ -830,12 +861,13 @@ void dp_red(DP p0,DP p1,DP p2,DP *head,DP *rest,P *dnp *head = h; *rest = r; *dnp = (P)c2; } +// head+rest = dn*(p0+p1)-mult*p2 void dpm_red(DPM p0,DPM p1,DPM p2,DPM *head,DPM *rest,P *dnp,DP *multp) { int i,n,pos; DL d1,d2,d; MP m; - DP s; + DP s,ms; DPM t,r,h,u,w; Z c,c1,c2,gn; P g,a; @@ -857,22 +889,62 @@ void dpm_red(DPM p0,DPM p1,DPM p2,DPM *head,DPM *rest, } else if ( INT(c1) && INT(c2) ) { gcdz(c1,c2,&gn); if ( !UNIQ(gn) ) { - divz(c1,gn,&c); c1 = c; - divz(c2,gn,&c); c2 = c; + 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; chsgnp((P)c1,(P *)&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td; + NEWMP(m); m->dl = d; m->c = (Obj)c1; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td; *multp = s; - mulobjdpm(CO,(Obj)s,p2,&u); mulobjdpm(CO,(Obj)c2,p1,&w); adddpm(CO,u,w,&r); + chsgnd(s,&ms); mulobjdpm(CO,(Obj)ms,p2,&u); mulobjdpm(CO,(Obj)c2,p1,&w); adddpm(CO,w,u,&r); mulobjdpm(CO,(Obj)c2,p0,&h); *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,...) @@ -905,8 +977,8 @@ void dp_red_marked(DP p0,DP p1,DP p2,DP hp2,DP *head,D } else if ( INT(c1) && INT(c2) ) { gcdz(c1,c2,&gn); if ( !UNIQ(gn) ) { - divz(c1,gn,&c); c1 = c; - divz(c2,gn,&c); c2 = c; + divsz(c1,gn,&c); c1 = c; + divsz(c2,gn,&c); c2 = c; } } else { ezgcdpz(CO,(P)c1,(P)c2,&g); @@ -1084,7 +1156,7 @@ void dp_true_nf(NODE b,DP g,DP *ps,int full,DP *rp,P * 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)); + wb[i] = ZTOS((Q)BDY(l)); sugar = g->sugar; for ( d = 0; g; ) { for ( u = 0, i = 0; i < n; i++ ) { @@ -1152,7 +1224,7 @@ void dp_removecont2(DP p1,DP p2,DP *r1p,DP *r2p,Z *con 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); divz(h,w[0],contp); + 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) ) { @@ -1172,6 +1244,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) @@ -1195,7 +1317,7 @@ void dp_true_nf_marked(NODE b,DP g,DP *ps,DP *hps,DP * 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((Z)BDY(l)); + wb[i] = ZTOS((Z)BDY(l)); sugar = g->sugar; for ( d = 0; g; ) { for ( u = 0, i = 0; i < n; i++ ) { @@ -1253,10 +1375,11 @@ void dp_true_nf_marked_mod(NODE b,DP g,DP *ps,DP *hps, if ( !g ) { *rp = 0; *dnp = dn; return; } - for ( n = 0, l = b; l; l = NEXT(l), n++ ); - wb = (int *)ALLOCA(n*sizeof(int)); + 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)); + wb[i] = ZTOS((Q)BDY(l)); sugar = g->sugar; for ( d = 0; g; ) { for ( u = 0, i = 0; i < n; i++ ) { @@ -1311,7 +1434,7 @@ DP *dp_true_nf_and_quotient_marked (NODE b,DP g,DP *ps 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)); + wb[i] = ZTOS((Q)BDY(l)); q = (DP *)MALLOC(n*sizeof(DP)); for ( i = 0; i < n; i++ ) q[i] = 0; sugar = g->sugar; @@ -1347,6 +1470,102 @@ last: return q; } +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,s,t,d; + DP dmy,mult,zzz; + DPM *ps; + DP *q; + NODE l; + DMM m,mr; + MP mp; + int i,n,j,len,nv; + int *wb; + int sugar,psugar,multiple; + P nm,tnm1,dn,tdn,tdn1; + Q cont; + 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; + } + q = (DP *)MALLOC(len*sizeof(DP)); + for ( i = 0; i < len; i++ ) q[i] = 0; + sugar = g->sugar; + for ( d = 0; g; ) { + for ( u = 0, i = 0; i < n; i++ ) { + if ( dpm_redble(g,p = ps[wb[i]]) ) { +// 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); + for ( j = 0; j < len; j++ ) { + if ( q[j] ) { mulcmp((Obj)tdn,BDY(q[j])); } + } + q[wb[i]] = appendd(q[wb[i]],mult); + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + if ( d ) mulcdmm((Obj)tdn,BDY(d)); + if ( !u ) goto last; + break; + } + } + if ( u ) { + g = u; + } 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; + d = appenddpm(d,t); + dpm_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; @@ -1367,7 +1586,7 @@ DP *dp_true_nf_and_quotient_marked_mod(NODE b,DP g,DP } wb = (int *)ALLOCA(n*sizeof(int)); for ( i = 0, l = b; i < n; l = NEXT(l), i++ ) - wb[i] = QTOS((Q)BDY(l)); + wb[i] = ZTOS((Q)BDY(l)); sugar = g->sugar; for ( d = 0; g; ) { for ( u = 0, i = 0; i < n; i++ ) { @@ -1421,7 +1640,7 @@ void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple, 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)); + wb[i] = ZTOS((Q)BDY(l)); hmag = multiple*HMAG(g); sugar = g->sugar; @@ -1473,11 +1692,13 @@ void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple, *rp = d; } -void dpm_nf_z(NODE b,DPM g,DPM *ps,int full,int multiple,DPM *rp) +void dpm_nf_z(NODE b,DPM g,VECT psv,int full,int multiple,DPM *rp) { + DPM *ps; DPM u,p,d,s,t; DP dmy1; P dmy; + Z cont; NODE l; DMM m,mr; int i,n; @@ -1488,26 +1709,34 @@ void dpm_nf_z(NODE b,DPM g,DPM *ps,int full,int multip if ( !g ) { *rp = 0; 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)); + 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)); + ps = (DPM *)BDY(psv); + } else { + n = psv->len; + wb = (int *)MALLOC(n*sizeof(int)); + for ( i = 0; i < n; i++ ) wb[i] = i; + ps = (DPM *)BDY(psv); + } hmag = multiple*HMAG(g); sugar = g->sugar; 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,&dmy,&dmy1); + if ( (p=ps[wb[i]])!=0 && dpm_redble(g,p) ) { + 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; } } @@ -1520,7 +1749,7 @@ void dpm_nf_z(NODE b,DPM g,DPM *ps,int full,int multip } } else { if ( multiple && HMAG(g) > hmag ) { - dpm_ptozp(g,&t); g = t; + dpm_ptozp(g,&cont,&t); g = t; hmag = multiple*HMAG(g); } } @@ -1533,7 +1762,7 @@ void dpm_nf_z(NODE b,DPM g,DPM *ps,int full,int multip } 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; } } @@ -1542,6 +1771,57 @@ void dpm_nf_z(NODE b,DPM g,DPM *ps,int full,int multip *rp = d; } +void dpm_shift(DPM p,int s,DPM *r) +{ + DMM m,mr0,mr; + DPM t; + + if ( !p ) *r = 0; + else { + for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) { + NEXTDMM(mr0,mr); + mr->dl = m->dl; mr->c = m->c; mr->pos = m->pos-s; + if ( mr->pos <= 0 ) + error("dpm_shift : too large shift value"); + } + NEXT(mr) = 0; + MKDPM(p->nv,mr0,t); t->sugar = p->sugar; + *r = t; + } +} + +// up=sum{c*<<...:i>>|i<=s}, lo=sum{c*<<...:i>>|i>s} + +void dpm_split(DPM p,int s,DPM *up,DPM *lo) +{ + DMM m,mu0,mu,ml0,ml; + DPM t; + + if ( !p ) { + *up = 0; *lo = 0; + } else { + for ( m = BDY(p), mu0 = ml0 = 0; m; m = NEXT(m) ) { + if ( m->pos <= s ) { + NEXTDMM(mu0,mu); + mu->dl = m->dl; mu->c = m->c; mu->pos = m->pos; + } else { + NEXTDMM(ml0,ml); + ml->dl = m->dl; ml->c = m->c; ml->pos = m->pos; + } + } + if ( mu0 ) { + NEXT(mu) = 0; MKDPM(p->nv,mu0,t); t->sugar = p->sugar; + *up = t; + } else + *up = 0; + if ( ml0 ) { + NEXT(ml) = 0; MKDPM(p->nv,ml0,t); t->sugar = p->sugar; + *lo = t; + } else + *lo = 0; + } +} + /* nf computation over a field */ void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp) @@ -1559,7 +1839,7 @@ void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp) 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)); + wb[i] = ZTOS((Q)BDY(l)); sugar = g->sugar; for ( d = 0; g; ) { @@ -1595,8 +1875,9 @@ void dp_nf_f(NODE b,DP g,DP *ps,int full,DP *rp) *rp = d; } -void dpm_nf_f(NODE b,DPM g,DPM *ps,int full,DPM *rp) +void dpm_nf_f(NODE b,DPM g,VECT psv,int full,DPM *rp) { + DPM *ps; DPM u,p,d,s,t; NODE l; DMM m,mr; @@ -1607,15 +1888,23 @@ void dpm_nf_f(NODE b,DPM g,DPM *ps,int full,DPM *rp) if ( !g ) { *rp = 0; 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)); + 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)); + ps = (DPM *)BDY(psv); + } else { + n = psv->len; + wb = (int *)MALLOC(n*sizeof(int)); + for ( i = 0; i < n; i++ ) wb[i] = i; + ps = (DPM *)BDY(psv); + } sugar = g->sugar; for ( d = 0; g; ) { for ( u = 0, i = 0; i < n; i++ ) { - if ( dpm_redble(g,p = ps[wb[i]]) ) { + if ( ( (p=ps[wb[i]]) != 0 ) && dpm_redble(g,p) ) { dpm_red_f(g,p,&u); psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; sugar = MAX(sugar,psugar); @@ -1708,10 +1997,11 @@ void dp_true_nf_mod(NODE b,DP g,DP *ps,int mod,int ful if ( !g ) { *rp = 0; *dnp = dn; return; } - for ( n = 0, l = b; l; l = NEXT(l), n++ ); - wb = (int *)ALLOCA(n*sizeof(int)); + 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)); + wb[i] = ZTOS((Q)BDY(l)); sugar = g->sugar; for ( d = 0; g; ) { for ( u = 0, i = 0; i < n; i++ ) { @@ -1932,54 +2222,91 @@ int create_order_spec(VL vl,Obj obj,struct order_spec *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); + 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 = QTOS((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,[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] = QTOS((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] = QTOS((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 { /* weighted order */ + 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"); } @@ -1991,8 +2318,8 @@ int create_order_spec(VL vl,Obj obj,struct order_spec 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++ ) { - tn = BDY((LIST)BDY(t)); l[i].order = QTOS((Q)BDY(tn)); - tn = NEXT(tn); l[i].length = QTOS((Q)BDY(tn)); + tn = BDY((LIST)BDY(t)); l[i].order = ZTOS((Q)BDY(tn)); + tn = NEXT(tn); l[i].length = ZTOS((Q)BDY(tn)); s += l[i].length; } spec->id = 1; spec->obj = obj; @@ -2005,7 +2332,7 @@ int create_order_spec(VL vl,Obj obj,struct order_spec w = almat(row,col); for ( i = 0; i < row; i++ ) for ( j = 0; j < col; j++ ) - w[i][j] = QTOS((Q)b[i][j]); + w[i][j] = ZTOS((Q)b[i][j]); spec->id = 2; spec->obj = obj; spec->nv = col; spec->ord.matrix.row = row; spec->ord.matrix.matrix = w; @@ -2077,9 +2404,9 @@ struct order_spec *append_block(struct order_spec *spe r = (struct order_spec *)MALLOC(sizeof(struct order_spec)); switch ( spec->id ) { case 0: - STOQ(spec->ord.simple,oq); STOQ(nv,nq); + STOZ(spec->ord.simple,oq); STOZ(nv,nq); t = mknode(2,oq,nq); MKLIST(list0,t); - STOQ(ord,oq); STOQ(nalg,nq); + STOZ(ord,oq); STOZ(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)); @@ -2105,7 +2432,7 @@ struct order_spec *append_block(struct order_spec *spe for ( t = BDY(list0), s0 = 0; t; t = NEXT(t) ) { NEXTNODE(s0,s); BDY(s) = BDY(t); } - STOQ(ord,oq); STOQ(nalg,nq); + STOZ(ord,oq); STOZ(nalg,nq); t = mknode(2,oq,nq); MKLIST(list,t); NEXTNODE(s0,s); BDY(s) = (pointer)list; NEXT(s) = 0; MKLIST(list,s0); @@ -2123,7 +2450,7 @@ struct order_spec *append_block(struct order_spec *spe MKMAT(mat,row+nalg,col+nalg); wp = (Z **)BDY(mat); for ( i = 0; i < row; i++ ) for ( j = 0; j < col; j++ ) { - w[i][j] = QTOS(b[i][j]); + w[i][j] = ZTOS(b[i][j]); wp[i][j] = b[i][j]; } for ( i = 0; i < nalg; i++ ) { @@ -2198,7 +2525,7 @@ int create_composite_order_spec(VL vl,LIST order,struc 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)); + dw[j] = ZTOS((Q)BDY(p)); } w_or_b[i].type = IS_DENSE_WEIGHT; w_or_b[i].length = len; @@ -2224,7 +2551,7 @@ int create_composite_order_spec(VL vl,LIST order,struc 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); + sw[j].value = ZTOS((Q)BDY(p)); p = NEXT(p); } qsort(sw,len,sizeof(struct sparse_weight), (int (*)(const void *,const void *))comp_sw); @@ -2248,7 +2575,7 @@ int create_composite_order_spec(VL vl,LIST order,struc len = end-start+1; sw = (struct sparse_weight *) MALLOC(sizeof(struct sparse_weight)*len); - w = QTOS((Q)BDY(NEXT(a))); + w = ZTOS((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; @@ -2310,6 +2637,7 @@ int create_composite_order_spec(VL vl,LIST order,struc w_or_b[n].body.block.order = 0; spec->ord.composite.length = n+1; } + return 1; } /* module order spec */ @@ -2330,12 +2658,12 @@ void create_modorder_spec(int id,LIST shift,struct mod 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)); + ds[i] = ZTOS((Q)BDY(t)); } else { spec->len = 0; spec->degree_shift = 0; } - STOQ(id,q); + STOZ(id,q); n = mknode(2,q,shift); MKLIST(list,n); spec->obj = (Obj)list; @@ -2346,6 +2674,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; @@ -2393,6 +2771,7 @@ void dp_dehomo(DP p,DP *rp) } } + void dp_mod(DP p,int mod,NODE subst,DP *rp) { MP m,mr,mr0; @@ -2448,11 +2827,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)); @@ -2463,12 +2843,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"); } @@ -2479,10 +2859,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; @@ -2494,9 +2873,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; @@ -2532,17 +2910,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)); @@ -2553,11 +2929,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"); @@ -2725,6 +3100,26 @@ void dpm_rest(DPM p,DPM *rp) } } +int dpm_getdeg(DPM p,int *r) +{ + int max,n,i,rank; + DMM m; + int *d; + + if ( !p ) return 0; + n = p->nv; + max = 0; + rank = 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]; + rank = MAX(rank,m->pos); + } + *r = rank; + return max; +} + DL lcm_of_DL(int nv,DL dl1,DL dl2,DL dl) { register int i, *d1, *d2, *d, td; @@ -2788,13 +3183,13 @@ void _print_mp(int nv,MP m) if ( !m ) return; for ( ; m; m = NEXT(m) ) { - fprintf(stderr,"%d<",ITOS(C(m))); + fprintf(stderr,"%ld<",ITOS(C(m))); for ( i = 0; i < nv; i++ ) { fprintf(stderr,"%d",m->dl->d[i]); if ( i != nv-1 ) fprintf(stderr," "); } - fprintf(stderr,">",C(m)); + fprintf(stderr,">"); } fprintf(stderr,"\n"); } @@ -2915,6 +3310,7 @@ LIST dp_initial_term(LIST f,struct order_spec *ord) default: error("dp_initial_term : unsupported order"); } + return 0; } int highest_order_dp(DP p,int *weight,int n); @@ -2951,7 +3347,7 @@ LIST highest_order(LIST f,int *weight,int n) NEXTNODE(r0,r); p = (Obj)BDY(nd); h = highest_order_dp((DP)p,weight,n); - STOQ(h,q); + STOZ(h,q); BDY(r) = (pointer)q; } if ( r0 ) NEXT(r) = 0; @@ -2994,6 +3390,7 @@ LIST dp_order(LIST f,struct order_spec *ord) default: error("dp_initial_term : unsupported order"); } + return 0; } int dpv_ht(DPV p,DP *h) @@ -3064,7 +3461,7 @@ int compare_facet_preorder(int n,int *u,int *v, for ( i = 0; i < row2; i++ ) { w2i = w2[i]; for ( j = 0, tu = tv = 0; j < n; j++ ) - if ( s = w2i[j] ) { + if ( (s = w2i[j]) != 0 ) { tu += s*u[j]; tv += s*v[j]; } for ( j = 0; j < n; j++ ) uv[j] = u[j]*tv-v[j]*tu; @@ -3084,7 +3481,7 @@ Q inner_product_with_small_vector(VECT w,int *v) n = w->len; s = 0; for ( i = 0; i < n; i++ ) { - STOQ(v[i],q); mulq((Q)w->body[i],(Q)q,&t); addq(t,s,&u); s = u; + STOZ(v[i],q); mulq((Q)w->body[i],(Q)q,&t); addq(t,s,&u); s = u; } return s; } @@ -3278,7 +3675,7 @@ NODE compute_essential_df(DP *g,DP *gh,int ng) dl = (DL)BDY(rt); d = dl->d; ri = 0; for ( k = nv-1; k >= 0; k-- ) { - STOQ(d[k],q); + STOZ(d[k],q); MKNODE(r1,q,ri); ri = r1; } MKNODE(r1,0,ri); MKLIST(l,r1);