=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/dp-supp.c,v retrieving revision 1.41 retrieving revision 1.63 diff -u -p -r1.41 -r1.63 --- OpenXM_contrib2/asir2000/builtin/dp-supp.c 2007/08/21 23:53:00 1.41 +++ OpenXM_contrib2/asir2000/builtin/dp-supp.c 2016/03/31 07:33:32 1.63 @@ -45,7 +45,7 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.40 2006/12/12 11:50:37 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.62 2016/03/31 01:40:10 noro Exp $ */ #include "ca.h" #include "base.h" @@ -71,6 +71,30 @@ void print_composite_order_spec(struct order_spec *spe * */ +static NODE RatDenomList; + +void init_denomlist() +{ + RatDenomList = 0; +} + +void add_denomlist(P f) +{ + NODE n; + + if ( OID(f)==O_P ) { + MKNODE(n,f,RatDenomList); RatDenomList = n; + } +} + +LIST get_denomlist() +{ + LIST l; + + MKLIST(l,RatDenomList); RatDenomList = 0; + return l; +} + void dp_ptozp(DP p,DP *rp) { MP m,mr,mr0; @@ -448,6 +472,7 @@ void dp_prim(DP p,DP *rp) P *w; Q *c; Q dvr; + NODE tn; if ( !p ) *rp = 0; @@ -502,6 +527,7 @@ void dp_prim(DP p,DP *rp) NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl; } NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar; + add_denomlist(g); } } } @@ -625,11 +651,11 @@ void dp_sp(DP p1,DP p2,DP *rp) LIST hist; NODE node; - node = mknode(4,ONE,0,s1,ONE); + node = mknode(4,ONE,NULLP,s1,ONE); MKLIST(hist,node); MKNODE(TraceList,hist,0); - node = mknode(4,ONE,0,0,ONE); + node = mknode(4,ONE,NULLP,NULLP,ONE); chsgnd(s2,(DP *)&ARG2(node)); MKLIST(hist,node); MKNODE(node,hist,TraceList); TraceList = node; @@ -678,11 +704,11 @@ void _dp_sp_dup(DP p1,DP p2,DP *rp) LIST hist; NODE node; - node = mknode(4,ONE,0,s1,ONE); + node = mknode(4,ONE,NULLP,s1,ONE); MKLIST(hist,node); MKNODE(TraceList,hist,0); - node = mknode(4,ONE,0,0,ONE); + node = mknode(4,ONE,NULLP,NULLP,ONE); chsgnd(s2,(DP *)&ARG2(node)); MKLIST(hist,node); MKNODE(node,hist,TraceList); TraceList = node; @@ -805,6 +831,7 @@ void dp_red(DP p0,DP p1,DP p2,DP *head,DP *rest,P *dnp } else { ezgcdpz(CO,(P)c1,(P)c2,&g); divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a; + add_denomlist(g); } NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td; *multp = s; @@ -853,13 +880,43 @@ void dp_red_marked(DP p0,DP p1,DP p2,DP hp2,DP *head,D 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; + NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td; *multp = s; - muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); addd(CO,s,t,&r); + muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); subd(CO,s,t,&r); muldc(CO,p0,(P)c2,&h); *head = h; *rest = r; *dnp = (P)c2; } +void dp_red_marked_mod(DP p0,DP p1,DP p2,DP hp2,int mod,DP *head,DP *rest,P *dnp,DP *multp) +{ + int i,n; + DL d1,d2,d; + MP m; + DP t,s,r,h; + P c1,c2,g,u; + + n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(hp2)->dl; + NEWDL(d,n); d->td = d1->td - d2->td; + for ( i = 0; i < n; i++ ) + d->d[i] = d1->d[i]-d2->d[i]; + c1 = (P)BDY(p1)->c; c2 = (P)BDY(hp2)->c; + gcdprsmp(CO,mod,c1,c2,&g); + divsmp(CO,mod,c1,g,&u); c1 = u; divsmp(CO,mod,c2,g,&u); c2 = u; + if ( NUM(c2) ) { + divsmp(CO,mod,c1,c2,&u); c1 = u; c2 = (P)ONEM; + } + NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0; + MKDP(n,m,s); s->sugar = d->td; + *multp = s; + mulmd(CO,mod,s,p2,&t); + if ( NUM(c2) ) { + submd(CO,mod,p1,t,&r); h = p0; + } else { + mulmdc(CO,mod,p1,c2,&s); submd(CO,mod,s,t,&r); mulmdc(CO,mod,p0,c2,&h); + } + *head = h; *rest = r; *dnp = c2; +} + /* m-reduction over a field */ void dp_red_f(DP p1,DP p2,DP *rest) @@ -1008,18 +1065,127 @@ void dp_true_nf(NODE b,DP g,DP *ps,int full,DP *rp,P * *rp = d; *dnp = dn; } +DP *dp_true_nf_and_quotient(NODE b,DP g,DP *ps,DP *rp,P *dnp) +{ + DP u,p,d,s,t,dmy,hp,mult; + DP *q; + NODE l; + MP m,mr; + int i,n,j; + int *wb; + int sugar,psugar,multiple; + P nm,tnm1,dn,tdn,tdn1; + Q cont; + + dn = (P)ONE; + if ( !g ) { + *rp = 0; *dnp = dn; return 0; + } + for ( n = 0, l = b; l; l = NEXT(l), n++ ); + wb = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0, l = b; i < n; l = NEXT(l), i++ ) + wb[i] = QTOS((Q)BDY(l)); + q = (DP *)MALLOC(n*sizeof(DP)); + for ( i = 0; i < n; i++ ) q[i] = 0; + sugar = g->sugar; + for ( d = 0; g; ) { + for ( u = 0, i = 0; i < n; i++ ) { + if ( dp_redble(g,ps[wb[i]]) ) { + p = ps[wb[i]]; + dp_red(d,g,p,&t,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + for ( j = 0; j < n; j++ ) { + muldc(CO,q[j],(P)tdn,&dmy); q[j] = dmy; + } + addd(CO,q[wb[i]],mult,&dmy); q[wb[i]] = dmy; + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + d = t; + if ( !u ) goto last; + break; + } + } + if ( u ) { + g = u; + } else { + m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c; + NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td; + addd(CO,d,t,&s); d = s; + dp_rest(g,&t); g = t; + } + } +last: + if ( d ) d->sugar = sugar; + *rp = d; *dnp = dn; + return q; +} + +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 *dnp) +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; + int i,n,hmag; int *wb; - int sugar,psugar; - P dn,tdn,tdn1; + 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; @@ -1037,12 +1203,73 @@ void dp_true_nf_marked(NODE b,DP g,DP *ps,DP *hps,DP * psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; sugar = MAX(sugar,psugar); if ( !u ) { + goto last; + } else { + d = t; + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + } + break; + } + } + if ( u ) { + g = u; + if ( multiple && ((d && HMAG(d)>hmag) || (HMAG(g)>hmag)) ) { + dp_removecont2(d,g,&t,&u,&cont); d = t; g = u; + mulp(CO,nm,(P)cont,&tnm1); nm = tnm1; + if ( d ) + hmag = multiple*HMAG(d); + else + hmag = multiple*HMAG(g); + } + } else { + m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c; + NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td; + addd(CO,d,t,&s); d = s; + dp_rest(g,&t); g = t; + } + } +last: + if ( d ) { + dp_removecont2(d,0,&t,&u,&cont); d = t; + mulp(CO,nm,(P)cont,&tnm1); nm = tnm1; + d->sugar = sugar; + } + *rp = d; *nmp = nm; *dnp = dn; +} + +void dp_true_nf_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp) +{ + DP hp,u,p,d,s,t,dmy; + NODE l; + MP m,mr; + int i,n; + int *wb; + int sugar,psugar; + P dn,tdn,tdn1; + + dn = (P)ONEM; + if ( !g ) { + *rp = 0; *dnp = dn; return; + } + for ( n = 0, l = b; l; l = NEXT(l), n++ ); + wb = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0, l = b; i < n; l = NEXT(l), i++ ) + wb[i] = QTOS((Q)BDY(l)); + sugar = g->sugar; + for ( d = 0; g; ) { + for ( u = 0, i = 0; i < n; i++ ) { + if ( dp_redble(g,hp = hps[wb[i]]) ) { + p = ps[wb[i]]; + dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn,&dmy); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + if ( !u ) { if ( d ) d->sugar = sugar; *rp = d; *dnp = dn; return; } else { d = t; - mulp(CO,dn,tdn,&tdn1); dn = tdn1; + mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1; } break; } @@ -1052,13 +1279,125 @@ void dp_true_nf_marked(NODE b,DP g,DP *ps,DP *hps,DP * else { m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c; NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td; + addmd(CO,mod,d,t,&s); d = s; + dp_rest(g,&t); g = t; + } + } + if ( d ) + d->sugar = sugar; + *rp = d; *dnp = dn; +} + +/* true nf by a marked GB and collect quotients */ + +DP *dp_true_nf_and_quotient_marked (NODE b,DP g,DP *ps,DP *hps,DP *rp,P *dnp) +{ + DP u,p,d,s,t,dmy,hp,mult; + DP *q; + NODE l; + MP m,mr; + int i,n,j; + int *wb; + int sugar,psugar,multiple; + P nm,tnm1,dn,tdn,tdn1; + Q cont; + + dn = (P)ONE; + if ( !g ) { + *rp = 0; *dnp = dn; return 0; + } + for ( n = 0, l = b; l; l = NEXT(l), n++ ); + wb = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0, l = b; i < n; l = NEXT(l), i++ ) + wb[i] = QTOS((Q)BDY(l)); + q = (DP *)MALLOC(n*sizeof(DP)); + for ( i = 0; i < n; i++ ) q[i] = 0; + sugar = g->sugar; + for ( d = 0; g; ) { + for ( u = 0, i = 0; i < n; i++ ) { + if ( dp_redble(g,hp = hps[wb[i]]) ) { + p = ps[wb[i]]; + dp_red_marked(d,g,p,hp,&t,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + for ( j = 0; j < n; j++ ) { + muldc(CO,q[j],(P)tdn,&dmy); q[j] = dmy; + } + addd(CO,q[wb[i]],mult,&dmy); q[wb[i]] = dmy; + mulp(CO,dn,tdn,&tdn1); dn = tdn1; + d = t; + if ( !u ) goto last; + break; + } + } + if ( u ) { + g = u; + } else { + m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c; + NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td; addd(CO,d,t,&s); d = s; dp_rest(g,&t); g = t; } } +last: + if ( d ) d->sugar = sugar; + *rp = d; *dnp = dn; + return q; +} + +DP *dp_true_nf_and_quotient_marked_mod(NODE b,DP g,DP *ps,DP *hps,int mod,DP *rp,P *dnp) +{ + DP u,p,d,s,t,dmy,hp,mult; + DP *q; + NODE l; + MP m,mr; + int i,n,j; + int *wb; + int sugar,psugar; + P dn,tdn,tdn1; + + for ( n = 0, l = b; l; l = NEXT(l), n++ ); + q = (DP *)MALLOC(n*sizeof(DP)); + for ( i = 0; i < n; i++ ) q[i] = 0; + dn = (P)ONEM; + if ( !g ) { + *rp = 0; *dnp = dn; return 0; + } + wb = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0, l = b; i < n; l = NEXT(l), i++ ) + wb[i] = QTOS((Q)BDY(l)); + sugar = g->sugar; + for ( d = 0; g; ) { + for ( u = 0, i = 0; i < n; i++ ) { + if ( dp_redble(g,hp = hps[wb[i]]) ) { + p = ps[wb[i]]; + dp_red_marked_mod(d,g,p,hp,mod,&t,&u,&tdn,&mult); + psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar; + sugar = MAX(sugar,psugar); + for ( j = 0; j < n; j++ ) { + mulmdc(CO,mod,q[j],(P)tdn,&dmy); q[j] = dmy; + } + addmd(CO,mod,q[wb[i]],mult,&dmy); q[wb[i]] = dmy; + mulmp(CO,mod,dn,tdn,&tdn1); dn = tdn1; + d = t; + if ( !u ) goto last; + break; + } + } + if ( u ) + g = u; + else { + m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c; + NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td; + addmd(CO,mod,d,t,&s); d = s; + dp_rest(g,&t); g = t; + } + } +last: if ( d ) d->sugar = sugar; *rp = d; *dnp = dn; + return q; } /* nf computation over Z */ @@ -1451,12 +1790,14 @@ void dp_nf_tab_f(DP p,LIST *tab,DP *rp) int create_order_spec(VL vl,Obj obj,struct order_spec **specp) { - int i,j,n,s,row,col,ret; + int i,j,n,s,row,col,ret,wlen; struct order_spec *spec; struct order_pair *l; - NODE node,t,tn; + Obj wp,wm; + NODE node,t,tn,wpair; MAT m; - pointer **b; + VECT v; + pointer **b,*bv; int **w; if ( vl && obj && OID(obj) == O_LIST ) { @@ -1472,18 +1813,71 @@ int create_order_spec(VL vl,Obj obj,struct order_spec spec->ord.simple = QTOS((Q)obj); return 1; } else if ( OID(obj) == O_LIST ) { + /* module order; obj = [0|1,w,ord] or [0|1,ord] */ node = BDY((LIST)obj); - 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)); - s += l[i].length; - } - spec->id = 1; spec->obj = obj; - spec->ord.block.order_pair = l; - spec->ord.block.length = n; spec->nv = s; - return 1; + if ( !BDY(node) || NUM(BDY(node)) ) { + switch ( length(node) ) { + case 2: + 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->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"); + + 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)); + + 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)); + break; + default: + error("create_order_spec : invalid arguments for module order"); + } + + *specp = spec; + return 1; + } else { + /* block order in polynomial ring */ + 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)); + s += l[i].length; + } + spec->id = 1; spec->obj = obj; + spec->ord.block.order_pair = l; + spec->ord.block.length = n; spec->nv = s; + return 1; + } } else if ( OID(obj) == O_MAT ) { m = (MAT)obj; row = m->row; col = m->col; b = BDY(m); w = almat(row,col); @@ -1957,17 +2351,18 @@ void homogenize_order(struct order_spec *old,int n,str error("homogenize_order : invalid input"); } break; - case 1: + case 1: case 257: length = old->ord.block.length; l = (struct order_pair *) MALLOC_ATOMIC((length+1)*sizeof(struct order_pair)); bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair)); l[length].order = 2; l[length].length = 1; - new->id = 1; new->nv = n+1; + new->id = old->id; new->nv = n+1; new->ord.block.order_pair = l; new->ord.block.length = length+1; + new->ispot = old->ispot; break; - case 2: + case 2: case 258: nv = old->nv; row = old->ord.matrix.row; oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1); for ( i = 0; i <= nv; i++ ) @@ -1977,10 +2372,11 @@ void homogenize_order(struct order_spec *old,int n,str newm[i+1][j] = oldm[i][j]; newm[i+1][j] = 0; } - new->id = 2; new->nv = nv+1; + new->id = old->id; new->nv = nv+1; new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm; + new->ispot = old->ispot; break; - case 3: + case 3: case 259: onv = old->nv; nnv = onv+1; olen = old->ord.composite.length; @@ -2014,12 +2410,33 @@ 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 = 3; + new->id = old->id; new->nv = nnv; new->ord.composite.length = nlen; new->ord.composite.w_or_b = nwb; + new->ispot = old->ispot; print_composite_order_spec(new); break; + case 256: /* simple module order */ + switch ( old->ord.simple ) { + case 0: + new->id = 256; new->ord.simple = 0; break; + case 1: + l = (struct order_pair *) + MALLOC_ATOMIC(2*sizeof(struct order_pair)); + l[0].length = n; l[0].order = old->ord.simple; + l[1].length = 1; l[1].order = 2; + new->id = 257; + new->ord.block.order_pair = l; + new->ord.block.length = 2; new->nv = n+1; + break; + case 2: + new->id = 256; new->ord.simple = 1; break; + default: + error("homogenize_order : invalid input"); + } + new->ispot = old->ispot; + break; default: error("homogenize_order : invalid input"); } @@ -2040,9 +2457,12 @@ void qltozl(Q *w,int n,Q *dvr) v.id = O_VECT; v.len = n; v.body = (pointer *)w; igcdv(&v,dvr); return; } - c = w[0]; nm = NM(c); dn = INT(c) ? ONEN : DN(c); - for ( i = 1; i < n; i++ ) { - c = w[i]; l1 = INT(c) ? ONEN : DN(c); + for ( i = 0; !w[i]; i++ ); + c = w[i]; nm = NM(c); dn = INT(c) ? ONEN : DN(c); + for ( i++; i < n; i++ ) { + c = w[i]; + if ( !c ) continue; + l1 = INT(c) ? ONEN : DN(c); gcdn(nm,NM(c),&g); nm = g; gcdn(dn,l1,&l2); muln(dn,l1,&l3); divsn(l3,l2,&dn); } @@ -2446,4 +2866,359 @@ int dpv_ht(DPV p,DP *h) MKDP(e[maxi]->nv,mr,*h); (*h)->sugar = mr->dl->td; /* XXX */ return maxi; } +} + +/* return 1 if 0 <_w1 v && v <_w2 0 */ + +int in_c12(int n,int *v,int row1,int **w1,int row2, int **w2) +{ + int t1,t2; + + t1 = compare_zero(n,v,row1,w1); + t2 = compare_zero(n,v,row2,w2); + if ( t1 > 0 && t2 < 0 ) return 1; + else return 0; +} + +/* 0 < u => 1, 0 > u => -1 */ + +int compare_zero(int n,int *u,int row,int **w) +{ + int i,j,t; + int *wi; + + for ( i = 0; i < row; i++ ) { + wi = w[i]; + for ( j = 0, t = 0; j < n; j++ ) t += u[j]*wi[j]; + if ( t > 0 ) return 1; + else if ( t < 0 ) return -1; + } + return 0; +} + +/* functions for generic groebner walk */ +/* u=0 means u=-infty */ + +int compare_facet_preorder(int n,int *u,int *v, + int row1,int **w1,int row2,int **w2) +{ + int i,j,s,t,tu,tv; + int *w2i,*uv; + + if ( !u ) return 1; + uv = W_ALLOC(n); + for ( i = 0; i < row2; i++ ) { + w2i = w2[i]; + for ( j = 0, tu = tv = 0; j < n; j++ ) + if ( s = w2i[j] ) { + tu += s*u[j]; tv += s*v[j]; + } + for ( j = 0; j < n; j++ ) uv[j] = u[j]*tv-v[j]*tu; + t = compare_zero(n,uv,row1,w1); + if ( t > 0 ) return 1; + else if ( t < 0 ) return 0; + } + return 1; +} + +Q inner_product_with_small_vector(VECT w,int *v) +{ + int n,i; + Q q,s,t,u; + + n = w->len; + s = 0; + for ( i = 0; i < n; i++ ) { + STOQ(v[i],q); mulq((Q)w->body[i],q,&t); addq(t,s,&u); s = u; + } + return s; +} + +Q compute_last_t(NODE g,NODE gh,Q t,VECT w1,VECT w2,NODE *homo,VECT *wp) +{ + int n,i; + int *wt; + Q last,d1,d2,dn,nm,s,t1; + VECT wd,wt1,wt2,w; + NODE tg,tgh; + MP f; + int *h; + NODE r0,r; + MP m0,m; + DP d; + + n = w1->len; + wt = W_ALLOC(n); + last = ONE; + /* t1 = 1-t */ + for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) { + f = BDY((DP)BDY(tg)); + h = BDY((DP)BDY(tgh))->dl->d; + for ( ; f; f = NEXT(f) ) { + for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i]; + for ( i = 0; i < n && !wt[i]; i++ ); + if ( i == n ) continue; + d1 = inner_product_with_small_vector(w1,wt); + d2 = inner_product_with_small_vector(w2,wt); + nm = d1; subq(d1,d2,&dn); + /* if d1=d2 then nothing happens */ + if ( !dn ) continue; + /* s satisfies ds = 0*/ + divq(nm,dn,&s); + + if ( cmpq(s,t) > 0 && cmpq(s,last) < 0 ) + last = s; + else if ( !cmpq(s,t) ) { + if ( cmpq(d2,0) < 0 ) { + last = t; + break; + } + } + } + } + if ( !last ) { + dn = ONE; nm = 0; + } else { + NTOQ(NM(last),1,nm); + if ( INT(last) ) dn = ONE; + else { + NTOQ(DN(last),1,dn); + } + } + /* (1-n/d)*w1+n/d*w2 -> w=(d-n)*w1+n*w2 */ + subq(dn,nm,&t1); mulvect(CO,(Obj)w1,(Obj)t1,(Obj *)&wt1); + mulvect(CO,(Obj)w2,(Obj)nm,(Obj *)&wt2); addvect(CO,wt1,wt2,&w); + + r0 = 0; + for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) { + f = BDY((DP)BDY(tg)); + h = BDY((DP)BDY(tgh))->dl->d; + for ( m0 = 0; f; f = NEXT(f) ) { + for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i]; + for ( i = 0; i < n && !wt[i]; i++ ); + if ( !inner_product_with_small_vector(w,wt) ) { + NEXTMP(m0,m); m->c = f->c; m->dl = f->dl; + } + } + NEXT(m) = 0; + MKDP(((DP)BDY(tg))->nv,m0,d); d->sugar = ((DP)BDY(tg))->sugar; + NEXTNODE(r0,r); BDY(r) = (pointer)d; + } + NEXT(r) = 0; + *homo = r0; + *wp = w; + return last; +} + +/* return 0 if last_w = infty */ + +NODE compute_last_w(NODE g,NODE gh,int n,int **w, + int row1,int **w1,int row2,int **w2) +{ + DP d; + MP f,m0,m; + int *wt,*v,*h; + NODE t,s,n0,tn,n1,r0,r; + int i; + + wt = W_ALLOC(n); + n0 = 0; + for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) { + f = BDY((DP)BDY(t)); + h = BDY((DP)BDY(s))->dl->d; + for ( ; f; f = NEXT(f) ) { + for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i]; + for ( i = 0; i < n && !wt[i]; i++ ); + if ( i == n ) continue; + + if ( in_c12(n,wt,row1,w1,row2,w2) && + compare_facet_preorder(n,*w,wt,row1,w1,row2,w2) ) { + v = (int *)MALLOC_ATOMIC(n*sizeof(int)); + for ( i = 0; i < n; i++ ) v[i] = wt[i]; + MKNODE(n1,v,n0); n0 = n1; + } + } + } + if ( !n0 ) return 0; + for ( t = n0; t; t = NEXT(t) ) { + v = (int *)BDY(t); + for ( s = n0; s; s = NEXT(s) ) + if ( !compare_facet_preorder(n,v,(int *)BDY(s),row1,w1,row2,w2) ) + break; + if ( !s ) { + *w = v; + break; + } + } + if ( !t ) + error("compute_last_w : cannot happen"); + r0 = 0; + for ( t = g, s = gh; t; t = NEXT(t), s = NEXT(s) ) { + f = BDY((DP)BDY(t)); + h = BDY((DP)BDY(s))->dl->d; + for ( m0 = 0; f; f = NEXT(f) ) { + for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i]; + for ( i = 0; i < n && !wt[i]; i++ ); + if ( i == n || + (compare_facet_preorder(n,wt,*w,row1,w1,row2,w2) + && compare_facet_preorder(n,*w,wt,row1,w1,row2,w2)) ) { + NEXTMP(m0,m); m->c = f->c; m->dl = f->dl; + } + } + NEXT(m) = 0; + MKDP(((DP)BDY(t))->nv,m0,d); d->sugar = ((DP)BDY(t))->sugar; + NEXTNODE(r0,r); BDY(r) = (pointer)d; + } + NEXT(r) = 0; + return r0; +} + +/* compute a sufficient set of d(f)=u-v */ + +NODE compute_essential_df(DP *g,DP *gh,int ng) +{ + int nv,i,j,k,t,lj; + NODE r,r1,ri,rt,r0; + MP m; + MP *mj; + DL di,hj,dl,dlt; + int *d,*dt; + LIST l; + Q q; + + nv = g[0]->nv; + r = 0; + for ( j = 0; j < ng; j++ ) { + for ( m = BDY(g[j]), lj = 0; m; m = NEXT(m), lj++ ); + mj = (MP *)ALLOCA(lj*sizeof(MP)); + for ( m = BDY(g[j]), k = 0; m; m = NEXT(m), k++ ) + mj[k] = m; + for ( i = 0; i < lj; i++ ) { + for ( di = mj[i]->dl, k = i+1; k < lj; k++ ) + if ( _dl_redble(di,mj[k]->dl,nv) ) break; + if ( k < lj ) mj[i] = 0; + } + hj = BDY(gh[j])->dl; + _NEWDL(dl,nv); d = dl->d; + r0 = r; + for ( i = 0; i < lj; i++ ) { + if ( mj[i] && !dl_equal(nv,di=mj[i]->dl,hj) ) { + for ( k = 0, t = 0; k < nv; k++ ) { + d[k] = hj->d[k]-di->d[k]; + t += d[k]; + } + dl->td = t; +#if 1 + for ( rt = r0; rt; rt = NEXT(rt) ) { + dlt = (DL)BDY(rt); + if ( dlt->td != dl->td ) continue; + for ( dt = dlt->d, k = 0; k < nv; k++ ) + if ( d[k] != dt[k] ) break; + if ( k == nv ) break; + } +#else + rt = 0; +#endif + if ( !rt ) { + MKNODE(r1,dl,r); r = r1; + _NEWDL(dl,nv); d = dl->d; + } + } + } + } + for ( rt = r; rt; rt = NEXT(rt) ) { + dl = (DL)BDY(rt); d = dl->d; + ri = 0; + for ( k = nv-1; k >= 0; k-- ) { + STOQ(d[k],q); + MKNODE(r1,q,ri); ri = r1; + } + MKNODE(r1,0,ri); MKLIST(l,r1); + BDY(rt) = (pointer)l; + } + return r; +} + +int comp_bits_divisible(int *a,int *b,int n) +{ + int bpi,i,wi,bi; + + bpi = (sizeof(int)/sizeof(char))*8; + for ( i = 0; i < n; i++ ) { + wi = i/bpi; bi = i%bpi; + if ( !(a[wi]&(1< bb ) return 1; + else if ( ba < bb ) return -1; + } + return 0; +} + +NODE mono_raddec(NODE ideal) +{ + DP p; + int nv,w,i,bpi,di,c,len; + int *d,*s,*u,*new; + NODE t,t1,v,r,rem,prev; + + if( !ideal ) return 0; + p = (DP)BDY(ideal); + nv = NV(p); + bpi = (sizeof(int)/sizeof(char))*8; + w = (nv+(bpi-1))/bpi; + d = p->body->dl->d; + if ( !NEXT(ideal) ) { + for ( t = 0, i = nv-1; i >= 0; i-- ) { + if ( d[i] ) { + s = (int *)CALLOC(w,sizeof(int)); + s[i/bpi] |= 1<<(i%bpi); + MKNODE(t1,s,t); + t = t1; + } + } + return t; + } + rem = mono_raddec(NEXT(ideal)); + r = 0; + len = w*sizeof(int); + u = (int *)CALLOC(w,sizeof(int)); + for ( i = nv-1; i >= 0; i-- ) { + if ( d[i] ) { + for ( t = rem; t; t = NEXT(t) ) { + bcopy((char *)BDY(t),(char *)u,len); + u[i/bpi] |= 1<<(i%bpi); + for ( v = r; v; v = NEXT(v) ) { + if ( comp_bits_divisible(u,(int *)BDY(v),nv) ) break; + } + if ( v ) continue; + for ( v = r, prev = 0; v; v = NEXT(v) ) { + if ( comp_bits_divisible((int *)BDY(v),u,nv) ) { + if ( prev ) NEXT(prev) = NEXT(v); + else r = NEXT(r); + } else prev =v; + } + for ( v = r, prev = 0; v; prev = v, v = NEXT(v) ) { + if ( comp_bits_lex(u,(int *)BDY(v),nv) < 0 ) break; + } + new = (int *)CALLOC(w,sizeof(int)); + bcopy((char *)u,(char *)new,len); + MKNODE(t1,new,v); + if ( prev ) NEXT(prev) = t1; + else r = t1; + } + } + } + return r; }