=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2018/engine/nd.c,v retrieving revision 1.11 retrieving revision 1.20 diff -u -p -r1.11 -r1.20 --- OpenXM_contrib2/asir2018/engine/nd.c 2018/10/23 04:53:37 1.11 +++ OpenXM_contrib2/asir2018/engine/nd.c 2019/09/15 08:46:12 1.20 @@ -1,13 +1,15 @@ -/* $OpenXM: OpenXM_contrib2/asir2018/engine/nd.c,v 1.10 2018/10/19 23:27:38 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2018/engine/nd.c,v 1.19 2019/09/04 05:32:10 noro Exp $ */ #include "nd.h" int Nnd_add,Nf4_red; -struct oEGT eg_search; +struct oEGT eg_search,f4_symb,f4_conv,f4_elim1,f4_elim2; int diag_period = 6; int weight_check = 1; int (*ndl_compare_function)(UINT *a1,UINT *a2); +/* for schreyer order */ +int (*ndl_base_compare_function)(UINT *a1,UINT *a2); int nd_dcomp; int nd_rref2; NM _nm_free_list; @@ -78,6 +80,7 @@ NDV pltondv(VL vl,VL dvl,LIST p); void pltozpl(LIST l,Q *cont,LIST *pp); void ndl_max(UINT *d1,unsigned *d2,UINT *d); void nmtodp(int mod,NM m,DP *r); +void ndltodp(UINT *d,DP *r); NODE reverse_node(NODE n); P ndc_div(int mod,union oNDC a,union oNDC b); P ndctop(int mod,union oNDC c); @@ -87,6 +90,11 @@ void parse_nd_option(NODE opt); void dltondl(int n,DL dl,UINT *r); DP ndvtodp(int mod,NDV p); DP ndtodp(int mod,ND p); +DPM ndvtodpm(int mod,NDV p); +NDV dpmtondv(int mod,DPM p); +int dpm_getdeg(DPM p,int *rank); +void dpm_ptozp(DPM p,Z *cont,DPM *r); +int compdmm(int nv,DMM a,DMM b); void Pdp_set_weight(NODE,VECT *); void Pox_cmo_rpc(NODE,Obj *); @@ -470,8 +478,11 @@ int ndl_weight(UINT *d) for ( j = 0; j < nd_epw; j++, u>>=nd_bpe ) t += (u&nd_mask0); } - if ( nd_module && current_module_weight_vector && MPOS(d) ) - t += current_module_weight_vector[MPOS(d)]; + if ( nd_module && nd_module_rank && MPOS(d) ) + t += nd_module_weight[MPOS(d)-1]; + for ( i = nd_exporigin; i < nd_wpd; i++ ) + if ( d[i] && !t ) + printf("afo\n"); return t; } @@ -486,8 +497,8 @@ int ndl_weight2(UINT *d) u = GET_EXP(d,i); t += nd_sugarweight[i]*u; } - if ( nd_module && current_module_weight_vector && MPOS(d) ) - t += current_module_weight_vector[MPOS(d)]; + if ( nd_module && nd_module_rank && MPOS(d) ) + t += nd_module_weight[MPOS(d)-1]; return t; } @@ -705,19 +716,19 @@ int ndl_module_grlex_compare(UINT *d1,UINT *d2) { int i,c; - if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; +// if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; if ( nd_ispot ) { - if ( nd_pot_nelim && MPOS(d1)>=nd_pot_nelim+1 && MPOS(d2) >= nd_pot_nelim+1 ) { - if ( TD(d1) > TD(d2) ) return 1; - else if ( TD(d1) < TD(d2) ) return -1; - if ( (c = ndl_lex_compare(d1,d2)) != 0 ) return c; - if ( MPOS(d1) < MPOS(d2) ) return 1; - else if ( MPOS(d1) > MPOS(d2) ) return -1; - return 0; + if ( nd_pot_nelim && MPOS(d1)>=nd_pot_nelim+1 && MPOS(d2) >= nd_pot_nelim+1 ) { + if ( TD(d1) > TD(d2) ) return 1; + else if ( TD(d1) < TD(d2) ) return -1; + if ( (c = ndl_lex_compare(d1,d2)) != 0 ) return c; + if ( MPOS(d1) < MPOS(d2) ) return 1; + else if ( MPOS(d1) > MPOS(d2) ) return -1; + return 0; + } + if ( MPOS(d1) < MPOS(d2) ) return 1; + else if ( MPOS(d1) > MPOS(d2) ) return -1; } - if ( MPOS(d1) < MPOS(d2) ) return 1; - else if ( MPOS(d1) > MPOS(d2) ) return -1; - } if ( TD(d1) > TD(d2) ) return 1; else if ( TD(d1) < TD(d2) ) return -1; if ( (c = ndl_lex_compare(d1,d2)) != 0 ) return c; @@ -732,7 +743,7 @@ int ndl_module_glex_compare(UINT *d1,UINT *d2) { int i,c; - if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; +// if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; if ( nd_ispot ) { if ( MPOS(d1) < MPOS(d2) ) return 1; else if ( MPOS(d1) > MPOS(d2) ) return -1; @@ -751,7 +762,7 @@ int ndl_module_lex_compare(UINT *d1,UINT *d2) { int i,c; - if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; +// if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; if ( nd_ispot ) { if ( MPOS(d1) < MPOS(d2) ) return 1; else if ( MPOS(d1) > MPOS(d2) ) return -1; @@ -768,7 +779,7 @@ int ndl_module_block_compare(UINT *d1,UINT *d2) { int i,c; - if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; +// if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; if ( nd_ispot ) { if ( MPOS(d1) < MPOS(d2) ) return 1; else if ( MPOS(d1) > MPOS(d2) ) return -1; @@ -785,7 +796,7 @@ int ndl_module_matrix_compare(UINT *d1,UINT *d2) { int i,c; - if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; +// if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; if ( nd_ispot ) { if ( MPOS(d1) < MPOS(d2) ) return 1; else if ( MPOS(d1) > MPOS(d2) ) return -1; @@ -802,7 +813,7 @@ int ndl_module_composite_compare(UINT *d1,UINT *d2) { int i,c; - if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; +// if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c; if ( nd_ispot ) { if ( MPOS(d1) > MPOS(d2) ) return 1; else if ( MPOS(d1) < MPOS(d2) ) return -1; @@ -1130,7 +1141,7 @@ INLINE int ndl_hash_value(UINT *d) r = 0; for ( i = 0; i < nd_wpd; i++ ) - r = (r*10007+d[i]); + r = (r*1511+d[i]); r %= REDTAB_LEN; return r; } @@ -2605,7 +2616,7 @@ ND_pairs nd_newpairs( NODE g, int t ) dl = DL(nd_psh[t]); ts = SG(nd_psh[t]) - TD(dl); - if ( nd_module && nd_intersect && (MPOS(dl) > 1) ) return 0; + if ( nd_module && nd_intersect && (MPOS(dl) > nd_intersect) ) return 0; for ( r0 = 0, h = g; h; h = NEXT(h) ) { if ( nd_module && (MPOS(DL(nd_psh[(long)BDY(h)])) != MPOS(dl)) ) continue; @@ -3224,6 +3235,7 @@ void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int int *perm; EPOS oepos; int obpe,oadv,ompos,cbpe; + VECT hvect; nd_module = 0; if ( !m && Demand ) nd_demand = 1; @@ -3267,12 +3279,18 @@ void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int for ( t = BDY(f), max = 1; t; t = NEXT(t) ) for ( tv = vv; tv; tv = NEXT(tv) ) { if ( nd_module ) { - s = BDY((LIST)BDY(t)); - trank = length(s); - mrank = MAX(mrank,trank); - for ( ; s; s = NEXT(s) ) { - e = getdeg(tv->v,(P)BDY(s)); - max = MAX(e,max); + if ( OID(BDY(t)) == O_DPM ) { + e = dpm_getdeg((DPM)BDY(t),&trank); + max = MAX(e,max); + mrank = MAX(mrank,trank); + } else { + s = BDY((LIST)BDY(t)); + trank = length(s); + mrank = MAX(mrank,trank); + for ( ; s; s = NEXT(s) ) { + e = getdeg(tv->v,(P)BDY(s)); + max = MAX(e,max); + } } } else { e = getdeg(tv->v,(P)BDY(t)); @@ -3284,9 +3302,18 @@ void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int ishomo = 1; for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) { if ( nd_module ) { - if ( !m && !nd_gentrace ) pltozpl((LIST)BDY(t),&dmy,&zpl); - else zpl = (LIST)BDY(t); + if ( OID(BDY(t)) == O_DPM ) { + Z cont; + DPM zdpm; + + if ( !m && !nd_gentrace ) dpm_ptozp((DPM)BDY(t),&cont,&zdpm); + else zdpm = (DPM)BDY(t); + b = (pointer)dpmtondv(m,zdpm); + } else { + if ( !m && !nd_gentrace ) pltozpl((LIST)BDY(t),&dmy,&zpl); + else zpl = (LIST)BDY(t); b = (pointer)pltondv(CO,vv,zpl); + } } else { if ( !m && !nd_gentrace ) ptozp((P)BDY(t),1,&dmy,&zp); else zp = (P)BDY(t); @@ -3334,6 +3361,11 @@ void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int if ( !x ) { *rp = 0; return; } + if ( nd_gentrace ) { + MKVECT(hvect,nd_psn); + for ( i = 0; i < nd_psn; i++ ) + ndltodp(nd_psh[i]->dl,(DP *)&BDY(hvect)[i]); + } if ( !ishomo && homo ) { /* dehomogenization */ for ( t = x; t; t = NEXT(t) ) ndv_dehomogenize((NDV)BDY(t),ord); @@ -3343,7 +3375,7 @@ void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int nd_demand = 0; if ( nd_module && nd_intersect ) { for ( j = nd_psn-1, x = 0; j >= 0; j-- ) - if ( MPOS(DL(nd_psh[j])) > 1 ) { + if ( MPOS(DL(nd_psh[j])) > nd_intersect ) { MKNODE(xx,(pointer)((unsigned long)j),x); x = xx; } conv_ilist(nd_demand,0,x,0); @@ -3367,10 +3399,12 @@ void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int nd_setup_parameters(nd_nvar,0); FINAL: for ( r0 = 0, t = x; t; t = NEXT(t) ) { - NEXTNODE(r0,r); - if ( nd_module ) BDY(r) = ndvtopl(m,CO,vv,BDY(t),mrank); - else if ( retdp ) BDY(r) = ndvtodp(m,BDY(t)); - else BDY(r) = ndvtop(m,CO,vv,BDY(t)); + NEXTNODE(r0,r); + if ( nd_module ) { + if ( retdp ) BDY(r) = ndvtodpm(m,BDY(t)); + else BDY(r) = ndvtopl(m,CO,vv,BDY(t),mrank); + } else if ( retdp ) BDY(r) = ndvtodp(m,BDY(t)); + else BDY(r) = ndvtop(m,CO,vv,BDY(t)); } if ( r0 ) NEXT(r) = 0; if ( !m && nd_nalg ) @@ -3380,8 +3414,7 @@ FINAL: if ( f4 ) { STOZ(16,bpe); STOZ(nd_last_nonzero,last_nonzero); - tr = mknode(5,*rp,(!ishomo&&homo)?ONE:0,BDY(nzlist),bpe,last_nonzero); MKLIST(*rp,tr); - + tr = mknode(6,*rp,(!ishomo&&homo)?ONE:0,BDY(nzlist),bpe,last_nonzero,hvect); MKLIST(*rp,tr); } else { tl1 = reverse_node(tl1); tl2 = reverse_node(tl2); tl3 = reverse_node(tl3); @@ -3401,7 +3434,7 @@ FINAL: MKLIST(l1,tl1); MKLIST(l2,tl2); MKLIST(l3,t); MKLIST(l4,tl3); MKLIST(l5,tl4); STOZ(nd_bpe,bpe); - tr = mknode(8,*rp,(!ishomo&&homo)?ONE:0,l1,l2,l3,l4,l5,bpe); MKLIST(*rp,tr); + tr = mknode(9,*rp,(!ishomo&&homo)?ONE:0,l1,l2,l3,l4,l5,bpe,hvect); MKLIST(*rp,tr); } } #if 0 @@ -3645,7 +3678,7 @@ void nd_gr_recompute_trace(LIST f,LIST v,int m,struct if ( DP_Print ) fprintf(asir_out,"\n"); } -void nd_gr_trace(LIST f,LIST v,int trace,int homo,int f4,struct order_spec *ord,LIST *rp) +void nd_gr_trace(LIST f,LIST v,int trace,int homo,int retdp,int f4,struct order_spec *ord,LIST *rp) { VL tv,fv,vv,vc,av; NODE fd,fd0,in0,in,r,r0,t,s,cand,alist; @@ -3668,6 +3701,7 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int int *perm; int j,ret; Z jq,bpe; + VECT hvect; nd_module = 0; nd_lf = 0; @@ -3722,6 +3756,11 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int for ( t = BDY(f), max = 1; t; t = NEXT(t) ) for ( tv = vv; tv; tv = NEXT(tv) ) { if ( nd_module ) { + if ( OID(BDY(t)) == O_DPM ) { + e = dpm_getdeg((DPM)BDY(t),&trank); + max = MAX(e,max); + mrank = MAX(mrank,trank); + } else { s = BDY((LIST)BDY(t)); trank = length(s); mrank = MAX(mrank,trank); @@ -3729,6 +3768,7 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int e = getdeg(tv->v,(P)BDY(s)); max = MAX(e,max); } + } } else { e = getdeg(tv->v,(P)BDY(t)); max = MAX(e,max); @@ -3739,13 +3779,22 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int ishomo = 1; for ( in0 = 0, fd0 = 0, t = BDY(f); t; t = NEXT(t) ) { if ( nd_module ) { - if ( !nd_gentrace ) pltozpl((LIST)BDY(t),&dmy,&zpl); - else zpl = (LIST)BDY(t); + if ( OID(BDY(t)) == O_DPM ) { + Z cont; + DPM zdpm; + + if ( !nd_gentrace ) dpm_ptozp((DPM)BDY(t),&cont,&zdpm); + else zdpm = (DPM)BDY(t); + c = (pointer)dpmtondv(m,zdpm); + } else { + if ( !nd_gentrace ) pltozpl((LIST)BDY(t),&dmy,&zpl); + else zpl = (LIST)BDY(t); c = (pointer)pltondv(CO,vv,zpl); + } } else { - if ( !nd_gentrace ) ptozp((P)BDY(t),1,&dmy,&zp); - else zp = (P)BDY(t); - c = (pointer)ptondv(CO,vv,zp); + if ( !nd_gentrace ) ptozp((P)BDY(t),1,&dmy,&zp); + else zp = (P)BDY(t); + c = (pointer)ptondv(CO,vv,zp); } if ( ishomo ) ishomo = ishomo && ndv_ishomo(c); @@ -3785,6 +3834,11 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int else m = get_lprime(++mindex); continue; } + if ( nd_gentrace ) { + MKVECT(hvect,nd_psn); + for ( i = 0; i < nd_psn; i++ ) + ndltodp(nd_psh[i]->dl,(DP *)&BDY(hvect)[i]); + } if ( !ishomo && homo ) { /* dehomogenization */ for ( t = cand; t; t = NEXT(t) ) ndv_dehomogenize((NDV)BDY(t),ord); @@ -3837,8 +3891,11 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int nd_bpe = cbpe; nd_setup_parameters(nd_nvar,0); for ( r = cand; r; r = NEXT(r) ) { - if ( nd_module ) BDY(r) = ndvtopl(0,CO,vv,BDY(r),mrank); - else BDY(r) = (pointer)ndvtop(0,CO,vv,BDY(r)); + if ( nd_module ) { + if ( retdp ) BDY(r) = ndvtodpm(0,BDY(r)); + else BDY(r) = ndvtopl(0,CO,vv,BDY(r),mrank); + } else if ( retdp ) BDY(r) = ndvtodp(0,BDY(r)); + else BDY(r) = (pointer)ndvtop(0,CO,vv,BDY(r)); } if ( nd_nalg ) cand = postprocess_algcoef(av,alist,cand); @@ -3862,7 +3919,7 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int MKLIST(l1,tl1); MKLIST(l2,tl2); MKLIST(l3,t); MKLIST(l4,tl3); MKLIST(l5,tl4); STOZ(nd_bpe,bpe); - tr = mknode(8,*rp,(!ishomo&&homo)?ONE:0,l1,l2,l3,l4,l5,bpe); MKLIST(*rp,tr); + tr = mknode(9,*rp,(!ishomo&&homo)?ONE:0,l1,l2,l3,l4,l5,bpe,hvect); MKLIST(*rp,tr); } } @@ -3926,6 +3983,18 @@ void nmtodp(int mod,NM m,DP *r) *r = dp; } +void ndltodp(UINT *d,DP *r) +{ + DP dp; + MP mr; + + NEWMP(mr); + mr->dl = ndltodl(nd_nvar,d); + mr->c = (Obj)ONE; + NEXT(mr) = 0; MKDP(nd_nvar,mr,dp); dp->sugar = mr->dl->td; + *r = dp; +} + void ndl_print(UINT *dl) { int n; @@ -4090,7 +4159,7 @@ void ndv_homogenize(NDV p,int obpe,int oadv,EPOS oepos NMV m,mr0,mr,t; len = p->len; - for ( m = BDY(p), i = 0, max = 1; i < len; NMV_OADV(m), i++ ) + for ( m = BDY(p), i = 0, max = 0; i < len; NMV_OADV(m), i++ ) max = MAX(max,TD(DL(m))); mr0 = nmv_adv>oadv?(NMV)REALLOC(BDY(p),len*nmv_adv):BDY(p); m = (NMV)((char *)mr0+(len-1)*oadv); @@ -4235,14 +4304,18 @@ void mpz_removecont_array(mpz_t *c,int n) { mpz_t d0,a,u,u1,gcd; int i,j; - mpz_t *q,*r; + static mpz_t *q,*r; + static int c_len = 0; for ( i = 0; i < n; i++ ) if ( mpz_sgn(c[i]) ) break; if ( i == n ) return; gcdv_mpz_estimate(d0,c,n); - q = (mpz_t *)MALLOC(n*sizeof(mpz_t)); - r = (mpz_t *)MALLOC(n*sizeof(mpz_t)); + if ( n > c_len ) { + q = (mpz_t *)MALLOC(n*sizeof(mpz_t)); + r = (mpz_t *)MALLOC(n*sizeof(mpz_t)); + c_len = n; + } for ( i = 0; i < n; i++ ) { mpz_init(q[i]); mpz_init(r[i]); mpz_fdiv_qr(q[i],r[i],c[i],d0); @@ -5192,31 +5265,32 @@ NDV ptondv(VL vl,VL dvl,P p) void pltozpl(LIST l,Q *cont,LIST *pp) { - NODE nd,nd1; - int n; - P *pl; - Q *cl; - int i; - P dmy; - Z dvr; - LIST r; + NODE nd,nd1; + int n; + P *pl; + Q *cl; + int i; + P dmy; + Z dvr,inv; + LIST r; - nd = BDY(l); n = length(nd); - pl = (P *)MALLOC(n*sizeof(P)); - cl = (Q *)MALLOC(n*sizeof(P)); - for ( i = 0; i < n; i++, nd = NEXT(nd) ) - ptozp((P)BDY(nd),1,&cl[i],&dmy); - qltozl(cl,n,&dvr); - nd = BDY(l); - for ( i = 0; i < n; i++, nd = NEXT(nd) ) { - divsp(CO,(P)BDY(nd),(P)dvr,&pl[i]); - } - nd = 0; - for ( i = n-1; i >= 0; i-- ) { - MKNODE(nd1,pl[i],nd); nd = nd1; - } - MKLIST(r,nd); - *pp = r; + nd = BDY(l); n = length(nd); + pl = (P *)MALLOC(n*sizeof(P)); + cl = (Q *)MALLOC(n*sizeof(Q)); + for ( i = 0; i < n; i++, nd = NEXT(nd) ) { + ptozp((P)BDY(nd),1,&cl[i],&dmy); + } + qltozl(cl,n,&dvr); + divz(ONE,dvr,&inv); + nd = BDY(l); + for ( i = 0; i < n; i++, nd = NEXT(nd) ) + divsp(CO,(P)BDY(nd),(P)dvr,&pl[i]); + nd = 0; + for ( i = n-1; i >= 0; i-- ) { + MKNODE(nd1,pl[i],nd); nd = nd1; + } + MKLIST(r,nd); + *pp = r; } /* (a1,a2,...,an) -> a1*e(1)+...+an*e(n) */ @@ -5405,6 +5479,106 @@ NDV ndtondv(int mod,ND p) return d; } +static int dmm_comp_nv; + +int dmm_comp(DMM *a,DMM *b) +{ + return -compdmm(dmm_comp_nv,*a,*b); +} + +void dmm_sort_by_ord(DMM *a,int len,int nv) +{ + dmm_comp_nv = nv; + qsort(a,len,sizeof(DMM),(int (*)(const void *,const void *))dmm_comp); +} + +void dpm_sort(DPM p,DPM *rp) +{ + DMM t,t1; + int len,i,n; + DMM *a; + DPM d; + + if ( !p ) *rp = 0; + for ( t = BDY(p), len = 0; t; t = NEXT(t), len++ ); + a = (DMM *)MALLOC(len*sizeof(DMM)); + for ( i = 0, t = BDY(p); i < len; i++, t = NEXT(t) ) a[i] = t; + n = p->nv; + dmm_sort_by_ord(a,len,n); + t = 0; + for ( i = len-1; i >= 0; i-- ) { + NEWDMM(t1); + t1->c = a[i]->c; + t1->dl = a[i]->dl; + t1->pos = a[i]->pos; + t1->next = t; + t = t1; + } + MKDPM(n,t,d); + SG(d) = SG(p); + *rp = d; +} + +int dpm_comp(DPM *a,DPM *b) +{ + return compdpm(CO,*a,*b); +} + +NODE dpm_sort_list(NODE l) +{ + int i,len; + NODE t,t1; + DPM *a; + + len = length(l); + a = (DPM *)MALLOC(len*sizeof(DPM)); + for ( t = l, i = 0; i < len; i++, t = NEXT(t) ) a[i] = (DPM)BDY(t); + qsort(a,len,sizeof(DPM),(int (*)(const void *,const void *))dpm_comp); + t = 0; + for ( i = len-1; i >= 0; i-- ) { + MKNODE(t1,(pointer)a[i],t); t = t1; + } + return t; +} + +int nmv_comp(NMV a,NMV b) +{ + return -DL_COMPARE(a->dl,b->dl); +} + +NDV dpmtondv(int mod,DPM p) +{ + NDV d; + NMV m,m0; + DMM t; + DMM *a; + int i,len,n; + + if ( !p ) return 0; + for ( t = BDY(p), len = 0; t; t = NEXT(t), len++ ); + a = (DMM *)MALLOC(len*sizeof(DMM)); + for ( i = 0, t = BDY(p); i < len; i++, t = NEXT(t) ) a[i] = t; + n = p->nv; + dmm_sort_by_ord(a,len,n); + if ( mod > 0 || mod == -1 ) + m0 = m = (NMV)MALLOC_ATOMIC_IGNORE_OFF_PAGE(len*nmv_adv); + else + m0 = m = MALLOC(len*nmv_adv); +#if 0 + ndv_alloc += nmv_adv*len; +#endif + for ( i = 0; i < len; i++, NMV_ADV(m) ) { + dltondl(n,a[i]->dl,DL(m)); + MPOS(DL(m)) = a[i]->pos; + TD(DL(m)) = ndl_weight(DL(m)); + CZ(m) = (Z)a[i]->c; + } + qsort(m0,len,nmv_adv,(int (*)(const void *,const void *))nmv_comp); + MKNDV(NV(p),m0,len,d); + SG(d) = SG(p); + return d; +} + ND ndvtond(int mod,NDV p) { ND d; @@ -5447,6 +5621,29 @@ DP ndvtodp(int mod,NDV p) return d; } +DPM ndvtodpm(int mod,NDV p) +{ + DMM m,m0; + DPM d; + NMV t; + int i,len; + + if ( !p ) return 0; + m0 = 0; + len = p->len; + for ( t = BDY(p), i = 0; i < len; NMV_ADV(t), i++ ) { + NEXTDMM(m0,m); + m->dl = ndltodl(nd_nvar,DL(t)); + m->c = (Obj)ndctop(mod,t->c); + m->pos = MPOS(DL(t)); + } + NEXT(m) = 0; + MKDPM(nd_nvar,m0,d); + SG(d) = SG(p); + return d; +} + + DP ndtodp(int mod,ND p) { MP m,m0; @@ -5534,6 +5731,8 @@ NODE ndv_reducebase(NODE x,int *perm) /* XXX incomplete */ +extern int dpm_ordtype; + void nd_init_ord(struct order_spec *ord) { nd_module = (ord->id >= 256); @@ -5545,6 +5744,7 @@ void nd_init_ord(struct order_spec *ord) nd_poly_weight = ord->top_weight; nd_module_rank = ord->module_rank; nd_module_weight = ord->module_top_weight; + dpm_ordtype = ord->ispot; } nd_matrix = 0; nd_matrix_len = 0; @@ -6056,7 +6256,6 @@ int ndv_reduce_vect_q(Z *svect,int trace,int col,IndAr int ndv_reduce_vect_q(Z *svect0,int trace,int col,IndArray *imat,NM_ind_pair *rp0,int nred) { int i,j,k,len,pos,prev; - mpz_t *svect; mpz_t cs,cr,gcd; IndArray ivect; unsigned char *ivc; @@ -6068,12 +6267,17 @@ int ndv_reduce_vect_q(Z *svect0,int trace,int col,IndA int maxrs; double hmag; int l; + static mpz_t *svect; + static int svect_len=0; maxrs = 0; for ( i = 0; i < col && !svect0[i]; i++ ); if ( i == col ) return maxrs; hmag = p_mag((P)svect0[i])*nd_scale; - svect = (mpz_t *)MALLOC(col*sizeof(mpz_t)); + if ( col > svect_len ) { + svect = (mpz_t *)MALLOC(col*sizeof(mpz_t)); + svect_len = col; + } for ( i = 0; i < col; i++ ) { mpz_init(svect[i]); if ( svect0[i] ) @@ -6679,6 +6883,8 @@ NODE nd_f4(int m,int checkonly,int **indp) PGeoBucket bucket; struct oEGT eg0,eg1,eg_f4; Z i1,i2,sugarq; + + init_eg(&f4_symb); init_eg(&f4_conv); init_eg(&f4_conv); init_eg(&f4_elim1); init_eg(&f4_elim2); #if 0 ndv_alloc = 0; #endif @@ -6724,7 +6930,7 @@ NODE nd_f4(int m,int checkonly,int **indp) d = nd_reconstruct(0,d); continue; } - get_eg(&eg1); init_eg(&eg_f4); add_eg(&eg_f4,&eg0,&eg1); + get_eg(&eg1); init_eg(&eg_f4); add_eg(&eg_f4,&eg0,&eg1); add_eg(&f4_symb,&eg0,&eg1); if ( DP_Print ) fprintf(asir_out,"sugar=%d,symb=%.3fsec,", sugar,eg_f4.exectime); @@ -6773,8 +6979,11 @@ NODE nd_f4(int m,int checkonly,int **indp) #if 0 fprintf(asir_out,"ndv_alloc=%d\n",ndv_alloc); #endif - if ( DP_Print ) - fprintf(asir_out,"number of red=%d\n",Nf4_red); + if ( DP_Print ) { + fprintf(asir_out,"number of red=%d,",Nf4_red); + fprintf(asir_out,"symb=%.3fsec,conv=%.3fsec,elim1=%.3fsec,elim2=%.3fsec\n", + f4_symb.exectime,f4_conv.exectime,f4_elim1.exectime,f4_elim2.exectime); + } conv_ilist(nd_demand,0,g,indp); return g; } @@ -7083,7 +7292,7 @@ NODE nd_f4_red(int m,ND_pairs sp0,int trace,UINT *s0ve rhead[imat[i]->head] = 1; start = imat[i]->head; } - get_eg(&eg1); init_eg(&eg_conv); add_eg(&eg_conv,&eg0,&eg1); + get_eg(&eg1); init_eg(&eg_conv); add_eg(&eg_conv,&eg0,&eg1); add_eg(&f4_conv,&eg0,&eg1); if ( DP_Print ) { fprintf(asir_out,"conv=%.3fsec,",eg_conv.exectime); fflush(asir_out); @@ -7666,7 +7875,9 @@ int ndv_ishomo(NDV p) h = TD(DL(m)); NMV_ADV(m); for ( len--; len; len--, NMV_ADV(m) ) - if ( TD(DL(m)) != h ) return 0; + if ( TD(DL(m)) != h ) { + return 0; + } return 1; } @@ -8316,6 +8527,8 @@ void parse_nd_option(NODE opt) nd_newelim = value?1:0; else if ( !strcmp(key,"intersect") ) nd_intersect = value?1:0; + else if ( !strcmp(key,"syzgen") ) + nd_intersect = ZTOS((Q)value); else if ( !strcmp(key,"lf") ) nd_lf = value?1:0; else if ( !strcmp(key,"trace") ) { @@ -9088,33 +9301,27 @@ int ndv_reduce_vect64(int m,mp_limb_t *svect,mp_limb_t ivc = ivect->index.c; for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { pos = prev+ivc[j]; c1 = CM(mr); prev = pos; - if ( c1 ) { - c2 = svect[pos]+c1*c; - if ( c2 < svect[pos] ) cvect[pos]++; - svect[pos] = c2; - } + c2 = svect[pos]+c1*c; + if ( c2 < svect[pos] ) cvect[pos]++; + svect[pos] = c2; } break; case 2: ivs = ivect->index.s; for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { pos = prev+ivs[j]; c1 = CM(mr); prev = pos; - if ( c1 ) { - c2 = svect[pos]+c1*c; - if ( c2 < svect[pos] ) cvect[pos]++; - svect[pos] = c2; - } + c2 = svect[pos]+c1*c; + if ( c2 < svect[pos] ) cvect[pos]++; + svect[pos] = c2; } break; case 4: ivi = ivect->index.i; for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { pos = prev+ivi[j]; c1 = CM(mr); prev = pos; - if ( c1 ) { - c2 = svect[pos]+c1*c; - if ( c2 < svect[pos] ) cvect[pos]++; - svect[pos] = c2; - } + c2 = svect[pos]+c1*c; + if ( c2 < svect[pos] ) cvect[pos]++; + svect[pos] = c2; } break; } @@ -9170,7 +9377,7 @@ NODE nd_f4_red_mod64_main(int m,ND_pairs sp0,int nsp,U } nd_free(spol); } - get_eg(&eg1); init_eg(&eg_f4_1); add_eg(&eg_f4_1,&eg0,&eg1); + get_eg(&eg1); init_eg(&eg_f4_1); add_eg(&eg_f4_1,&eg0,&eg1); add_eg(&f4_elim1,&eg0,&eg1); if ( DP_Print ) { fprintf(asir_out,"elim1=%.3fsec,",eg_f4_1.exectime); fflush(asir_out); @@ -9191,7 +9398,7 @@ NODE nd_f4_red_mod64_main(int m,ND_pairs sp0,int nsp,U if ( r0 ) NEXT(r) = 0; for ( ; i < sprow; i++ ) GCFREE(spmat[i]); - get_eg(&eg2); init_eg(&eg_f4_2); add_eg(&eg_f4_2,&eg1,&eg2); + get_eg(&eg2); init_eg(&eg_f4_2); add_eg(&eg_f4_2,&eg1,&eg2); add_eg(&f4_elim2,&eg1,&eg2); init_eg(&eg_f4); add_eg(&eg_f4,&eg0,&eg2); if ( DP_Print ) { fprintf(asir_out,"elim2=%.3fsec,",eg_f4_2.exectime);