Return to dp-supp.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / builtin |
version 1.23, 2003/01/04 09:06:15 | version 1.44, 2007/09/15 10:17:08 | ||
---|---|---|---|
|
|
||
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, | * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, | ||
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. | * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. | ||
* | * | ||
* $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.22 2002/12/27 07:37:57 noro Exp $ | * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.43 2007/09/07 00:45:50 noro Exp $ | ||
*/ | */ | ||
#include "ca.h" | #include "ca.h" | ||
#include "base.h" | #include "base.h" | ||
|
|
||
extern int GenTrace; | extern int GenTrace; | ||
extern NODE TraceList; | extern NODE TraceList; | ||
int show_orderspec; | |||
void print_composite_order_spec(struct order_spec *spec); | |||
/* | /* | ||
* content reduction | * content reduction | ||
* | * | ||
|
|
||
*hp = h; *rp = r; | *hp = h; *rp = r; | ||
} | } | ||
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) | void dp_idiv(DP p,Q c,DP *rp) | ||
{ | { | ||
Q t; | Q t; | ||
|
|
||
} | } | ||
} | } | ||
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) | void dp_prim(DP p,DP *rp) | ||
{ | { | ||
P t,g; | P t,g; | ||
|
|
||
for ( m = BDY(p); m; m = NEXT(m) ) | for ( m = BDY(p); m; m = NEXT(m) ) | ||
if ( OID(m->c) == O_N ) { | if ( OID(m->c) == O_N ) { | ||
/* GCD of coeffs = 1 */ | /* GCD of coeffs = 1 */ | ||
*rp = p; | dp_monic_sf(p,rp); | ||
return; | return; | ||
} else break; | } else break; | ||
/* compute GCD over the finite fieid */ | /* compute GCD over the finite fieid */ | ||
|
|
||
w[i] = m->c; | w[i] = m->c; | ||
gcdsf(CO,w,n,&g); | gcdsf(CO,w,n,&g); | ||
if ( NUM(g) ) | if ( NUM(g) ) | ||
*rp = p; | dp_monic_sf(p,rp); | ||
else { | else { | ||
for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { | for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { | ||
NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl; | 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; | NEXT(mr) = 0; MKDP(p->nv,mr0,p1); p1->sugar = p->sugar; | ||
dp_monic_sf(p1,rp); | |||
} | } | ||
return; | return; | ||
} else if ( dp_fcoeffs ) | } else if ( dp_fcoeffs ) | ||
|
|
||
int i,r; | int i,r; | ||
P gcd,t,s1,s2,u; | P gcd,t,s1,s2,u; | ||
Q rq; | Q rq; | ||
DCP dc; | |||
extern int DP_Print; | |||
while ( 1 ) { | while ( 1 ) { | ||
for ( i = 0, s1 = 0; i < m; i++ ) { | for ( i = 0, s1 = 0; i < m; i++ ) { | ||
r = random(); UTOQ(r,rq); | r = random(); UTOQ(r,rq); | ||
|
|
||
mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u; | mulp(vl,pl[i],(P)rq,&t); addp(vl,s2,t,&u); s2 = u; | ||
} | } | ||
ezgcdp(vl,s1,s2,&gcd); | ezgcdp(vl,s1,s2,&gcd); | ||
if ( DP_Print > 2 ) | |||
{ fprintf(asir_out,"(%d)",nmonop(gcd)); fflush(asir_out); } | |||
for ( i = 0; i < m; i++ ) { | for ( i = 0; i < m; i++ ) { | ||
if ( !divtpz(vl,pl[i],gcd,&t) ) | if ( !divtpz(vl,pl[i],gcd,&t) ) | ||
break; | break; | ||
|
|
||
*head = h; *rest = r; *dnp = (P)c2; | *head = h; *rest = r; *dnp = (P)c2; | ||
} | } | ||
/* | |||
* m-reduction by a marked poly | |||
* do content reduction over Z or Q(x,...) | |||
* do nothing over finite fields | |||
* | |||
*/ | |||
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; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td; | |||
*multp = s; | |||
muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); addd(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) | |||
{ | |||
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; chsgnmp(mod,(P)c1,&m->c); NEXT(m) = 0; | |||
MKDP(n,m,s); s->sugar = d->td; mulmd(CO,mod,s,p2,&t); | |||
if ( NUM(c2) ) { | |||
addmd(CO,mod,p1,t,&r); h = p0; | |||
} else { | |||
mulmdc(CO,mod,p1,c2,&s); addmd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h); | |||
} | |||
*head = h; *rest = r; *dnp = c2; | |||
} | |||
/* m-reduction over a field */ | /* m-reduction over a field */ | ||
void dp_red_f(DP p1,DP p2,DP *rest) | void dp_red_f(DP p1,DP p2,DP *rest) | ||
|
|
||
*rp = d; *dnp = dn; | *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; | |||
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); | |||
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; | |||
} | |||
/* nf computation over Z */ | /* nf computation over Z */ | ||
void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp) | void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp) | ||
|
|
||
*rp = s; | *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 | * setting flags | ||
* call create_order_spec with vl=0 to set old type order. | |||
* | * | ||
*/ | */ | ||
int create_order_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; | struct order_pair *l; | ||
NODE node,t,tn; | NODE node,t,tn; | ||
MAT m; | MAT m; | ||
pointer **b; | pointer **b; | ||
int **w; | 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) ) { | if ( !obj || NUM(obj) ) { | ||
spec->id = 0; spec->obj = obj; | spec->id = 0; spec->obj = obj; | ||
spec->ord.simple = QTOS((Q)obj); | spec->ord.simple = QTOS((Q)obj); | ||
|
|
||
return 0; | 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 | * converters | ||
* | * | ||
|
|
||
} | } | ||
void homogenize_order(struct order_spec *old,int n,struct order_spec *new) | void homogenize_order(struct order_spec *old,int n,struct order_spec **newp) | ||
{ | { | ||
struct order_pair *l; | struct order_pair *l; | ||
int length,nv,row,i,j; | int length,nv,row,i,j; | ||
int **newm,**oldm; | 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 ) { | switch ( old->id ) { | ||
case 0: | case 0: | ||
switch ( old->ord.simple ) { | switch ( old->ord.simple ) { | ||
|
|
||
new->id = 2; new->nv = nv+1; | new->id = 2; new->nv = nv+1; | ||
new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm; | new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm; | ||
break; | break; | ||
case 3: | |||
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 = 3; | |||
new->nv = nnv; | |||
new->ord.composite.length = nlen; | |||
new->ord.composite.w_or_b = nwb; | |||
print_composite_order_spec(new); | |||
break; | |||
default: | default: | ||
error("homogenize_order : invalid input"); | error("homogenize_order : invalid input"); | ||
} | } | ||
|
|
||
} | } | ||
} | } | ||
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) | void dp_rest(DP p,DP *rp) | ||
{ | { | ||
MP m; | MP m; | ||
|
|
||
fprintf(stderr,">",C(m)); | fprintf(stderr,">",C(m)); | ||
} | } | ||
fprintf(stderr,"\n"); | 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; | |||
} | |||
/* 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 */ | |||
static int comp_vector_lex_nv; | |||
int comp_vector_lex(int **a,int **b) | |||
{ | |||
int i; | |||
int *pa,*pb; | |||
pa = *a; pb = *b; | |||
for ( i = 0; i < comp_vector_lex_nv; i++ ) | |||
if ( pa[i] < pb[i] ) return -1; | |||
else if ( pa[i] > pb[i] ) return 1; | |||
return 0; | |||
} | |||
NODE compute_essential_df(DP *g,DP *gh,int ng) | |||
{ | |||
VECT v; | |||
Q q; | |||
MP m; | |||
NODE r,r1; | |||
int nv,len,i,j,k; | |||
int *p,*dm,*mi,*mj,*h; | |||
int **mat; | |||
nv = comp_vector_lex_nv = g[0]->nv; | |||
for ( len = 0, j = 0; j < ng; j++ ) { | |||
for ( m = BDY(g[j]); m; m = NEXT(m), len++ ); | |||
} | |||
mat = almat(len,nv); | |||
for ( i = 0, j = 0; j < ng; j++ ) { | |||
h = BDY(gh[j])->dl->d; | |||
for ( m = BDY(g[j]); m; m = NEXT(m) ) { | |||
dm = m->dl->d; | |||
for ( k = 0; k < nv; k++ ) | |||
if ( dm[k] ) break; | |||
if ( k == nv ) continue; | |||
else { | |||
p = mat[i]; | |||
for ( k = 0; k < nv; k++ ) | |||
p[k] = h[k]-dm[k]; | |||
i++; | |||
} | |||
} | |||
} | |||
len = i; | |||
qsort(mat,len,sizeof(int *), | |||
(int (*)(const void *,const void *))comp_vector_lex); | |||
for ( i = 0; i < len; i++ ) { | |||
for ( j = 0; j < nv; j++ ) | |||
printf("%d ",mat[i][j]); | |||
printf("\n"); | |||
} | |||
for ( i = 0; i < len; i++ ) { | |||
mi = mat[i]; | |||
if ( !mi ) continue; | |||
for ( j = i+1; j < len; j++ ) { | |||
mj = mat[j]; | |||
if ( !mj ) continue; | |||
for ( k = 0; k < nv; k++ ) | |||
if ( mi[k] > mj[k] ) break; | |||
if ( k == nv ) mat[j] = 0; | |||
} | |||
} | |||
for ( i = 0; i < len; i++ ) { | |||
if ( mat[i] ) { | |||
for ( j = 0; j < nv; j++ ) | |||
printf("%d ",mat[i][j]); | |||
printf("\n"); | |||
} | |||
} | |||
r = 0; | |||
for ( i = 0; i < len; i++ ) { | |||
if ( mi = mat[i] ) { | |||
MKVECT(v,nv); | |||
for ( k = 0; k < nv; k++ ) { | |||
STOQ(mi[k],q); | |||
v->body[k] = (pointer)q; | |||
} | |||
MKNODE(r1,v,r); r = r1; | |||
} | |||
} | |||
return r; | |||
} | } |