=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v retrieving revision 1.223 retrieving revision 1.233 diff -u -p -r1.223 -r1.233 --- OpenXM_contrib2/asir2000/engine/nd.c 2015/08/14 13:51:54 1.223 +++ OpenXM_contrib2/asir2000/engine/nd.c 2017/02/07 08:30:30 1.233 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.222 2015/08/08 14:19:41 fujimoto Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.232 2017/01/08 03:05:39 noro Exp $ */ #include "nd.h" @@ -56,12 +56,15 @@ static int nd_found,nd_create,nd_notfirst; static int nmv_adv; static int nd_demand; static int nd_module,nd_ispot,nd_mpos,nd_pot_nelim; +static int nd_module_rank,nd_poly_weight_len; +static int *nd_poly_weight,*nd_module_weight; static NODE nd_tracelist; static NODE nd_alltracelist; static int nd_gentrace,nd_gensyz,nd_nora,nd_newelim,nd_intersect; static int *nd_gbblock; static NODE nd_nzlist,nd_check_splist; static int nd_splist; +static int *nd_sugarweight; NumberField get_numberfield(); UINT *nd_det_compute_bound(NDV **dm,int n,int j); @@ -91,8 +94,17 @@ ND ndgztond(ND p); void Pdp_set_weight(NODE,VECT *); void Pox_cmo_rpc(NODE,Obj *); -extern int Denominator,DP_Multiple; +ND nd_add_lf(ND p1,ND p2); +void nd_mul_c_lf(ND p,GZ mul); +void ndv_mul_c_lf(NDV p,GZ mul); +NODE nd_f4_red_lf_main(int m,ND_pairs sp0,int nsp,UINT *s0vect,int col, + NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred); +int nd_gauss_elim_lf(mpz_t **mat0,int *sugar,int row,int col,int *colstat); +extern int lf_lazy; +extern GZ current_mod_lf; +extern int Denominator,DP_Multiple,MaxDeg; + #define BLEN (8*sizeof(unsigned long)) typedef struct matrix { @@ -171,6 +183,8 @@ INLINE int nd_length(ND p) } } +extern int dp_negative_weight; + INLINE int ndl_reducible(UINT *d1,UINT *d2) { UINT u1,u2; @@ -178,7 +192,7 @@ INLINE int ndl_reducible(UINT *d1,UINT *d2) if ( nd_module && (MPOS(d1) != MPOS(d2)) ) return 0; - if ( TD(d1) < TD(d2) ) return 0; + if ( !dp_negative_weight && TD(d1) < TD(d2) ) return 0; #if USE_UNROLL switch ( nd_bpe ) { case 3: @@ -457,6 +471,22 @@ int ndl_weight(UINT *d) return t; } +/* for sugarweight */ + +int ndl_weight2(UINT *d) +{ + int t,u; + int i,j; + + for ( i = 0, t = 0; i < nd_nvar; i++ ) { + 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)]; + return t; +} + void ndl_weight_mask(UINT *d) { UINT t,u; @@ -647,10 +677,28 @@ int ndl_ww_lex_compare(UINT *d1,UINT *d2) return ndl_lex_compare(d1,d2); } +int ndl_module_weight_compare(UINT *d1,UINT *d2) +{ + int s,j; + + if ( nd_nvar != nd_poly_weight_len ) + error("invalid module weight : the length of polynomial weight != the number of variables"); + s = 0; + for ( j = 0; j < nd_nvar; j++ ) + s += (GET_EXP(d1,j)-GET_EXP(d2,j))*nd_poly_weight[j]; + if ( MPOS(d1) >= 1 && MPOS(d2) >= 1 ) { + s += nd_module_weight[MPOS(d1)-1]-nd_module_weight[MPOS(d2)-1]; + } + if ( s > 0 ) return 1; + else if ( s < 0 ) return -1; + else return 0; +} + 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_ispot ) { if ( nd_pot_nelim && MPOS(d1)>=nd_pot_nelim+1 && MPOS(d2) >= nd_pot_nelim+1 ) { if ( TD(d1) > TD(d2) ) return 1; @@ -677,6 +725,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_ispot ) { if ( MPOS(d1) < MPOS(d2) ) return 1; else if ( MPOS(d1) > MPOS(d2) ) return -1; @@ -695,6 +744,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_ispot ) { if ( MPOS(d1) < MPOS(d2) ) return 1; else if ( MPOS(d1) > MPOS(d2) ) return -1; @@ -711,6 +761,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_ispot ) { if ( MPOS(d1) < MPOS(d2) ) return 1; else if ( MPOS(d1) > MPOS(d2) ) return -1; @@ -727,6 +778,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_ispot ) { if ( MPOS(d1) < MPOS(d2) ) return 1; else if ( MPOS(d1) > MPOS(d2) ) return -1; @@ -743,6 +795,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_ispot ) { if ( MPOS(d1) > MPOS(d2) ) return 1; else if ( MPOS(d1) < MPOS(d2) ) return -1; @@ -1160,6 +1213,7 @@ ND nd_add(int mod,ND p1,ND p2) if ( !p1 ) return p2; else if ( !p2 ) return p1; else if ( mod == -1 ) return nd_add_sf(p1,p2); + else if ( mod == -2 ) return nd_add_lf(p1,p2); else if ( !mod ) return nd_add_q(p1,p2); else { can = 0; @@ -1320,6 +1374,7 @@ ND nd_add_sf(ND p1,ND p2) } } + ND nd_reduce2(int mod,ND d,ND g,NDV p,NM mul,NDC dn,Obj *divp) { int c,c1,c2; @@ -1330,6 +1385,12 @@ ND nd_reduce2(int mod,ND d,ND g,NDV p,NM mul,NDC dn,Ob if ( mod == -1 ) { CM(mul) = _mulsf(_invsf(HCM(p)),_chsgnsf(HCM(g))); *divp = (Obj)ONE; + } else if ( mod == -2 ) { + GZ inv,t; + divlf(ONEGZ,HCZ(p),&inv); + chsgnlf(HCZ(g),&t); + mullf(inv,t,&CZ(mul)); + *divp = (Obj)ONE; } else if ( mod ) { c1 = invm(HCM(p),mod); c2 = mod-HCM(g); DMAR(c1,c2,0,mod,c); CM(mul) = c; @@ -1481,7 +1542,12 @@ int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp d = 0; mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(UINT)); while ( 1 ) { - hindex = mod?head_pbucket(mod,bucket):head_pbucket_q(bucket); + if ( mod > 0 || mod == -1 ) + hindex = head_pbucket(mod,bucket); + else if ( mod == -2 ) + hindex = head_pbucket_lf(bucket); + else + hindex = head_pbucket_q(bucket); if ( hindex < 0 ) { if ( DP_Print > 3 ) printf("(%d %d)",count,hcount); if ( d ) SG(d) = sugar; @@ -1504,7 +1570,12 @@ int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp p = ps[index]; if ( mod == -1 ) CM(mul) = _mulsf(_invsf(HCM(p)),_chsgnsf(HCM(g))); - else if ( mod ) { + else if ( mod == -2 ) { + GZ inv,t; + divlf(ONEGZ,HCZ(p),&inv); + chsgnlf(HCZ(g),&t); + mullf(inv,t,&CZ(mul)); + } else if ( mod ) { c1 = invm(HCM(p),mod); c2 = mod-HCM(g); DMAR(c1,c2,0,mod,c); CM(mul) = c; } else { @@ -1812,6 +1883,46 @@ int head_pbucket_q(PGeoBucket g) } } +int head_pbucket_lf(PGeoBucket g) +{ + int j,i,c,k,nv; + GZ sum,t; + ND gi,gj; + + k = g->m; + while ( 1 ) { + j = -1; + for ( i = 0; i <= k; i++ ) { + if ( !(gi = g->body[i]) ) continue; + if ( j < 0 ) { + j = i; + gj = g->body[j]; + sum = HCZ(gj); + } else { + nv = NV(gi); + c = DL_COMPARE(HDL(gi),HDL(gj)); + if ( c > 0 ) { + if ( sum ) HCZ(gj) = sum; + else g->body[j] = nd_remove_head(gj); + j = i; + gj = g->body[j]; + sum = HCZ(gj); + } else if ( c == 0 ) { + addlf(sum,HCZ(gi),&t); + sum = t; + g->body[i] = nd_remove_head(gi); + } + } + } + if ( j < 0 ) return -1; + else if ( sum ) { + HCZ(gj) = sum; + return j; + } else + g->body[j] = nd_remove_head(gj); + } +} + ND normalize_pbucket(int mod,PGeoBucket g) { int i; @@ -1980,6 +2091,7 @@ NODE nd_gb(int m,int ishomo,int checkonly,int gensyz,i while ( d ) { again: l = nd_minp(d,&d); + if ( MaxDeg > 0 && SG(l) > MaxDeg ) break; if ( SG(l) != sugar ) { if ( ishomo ) { diag_count = 0; @@ -2238,6 +2350,7 @@ NODE nd_gb_trace(int m,int ishomo,int **indp) while ( d ) { again: l = nd_minp(d,&d); + if ( MaxDeg > 0 && SG(l) > MaxDeg ) break; if ( SG(l) != sugar ) { #if 1 if ( ishomo ) { @@ -2335,9 +2448,11 @@ again: FREENDP(l); } if ( nd_nalg ) { - print_eg("monic",&eg_monic); - print_eg("invdalg",&eg_invdalg); - print_eg("le",&eg_le); + if ( DP_Print ) { + print_eg("monic",&eg_monic); + print_eg("invdalg",&eg_invdalg); + print_eg("le",&eg_le); + } } conv_ilist(nd_demand,1,g,indp); if ( DP_Print ) { printf("nd_gb_trace done.\n"); fflush(stdout); } @@ -2503,11 +2618,38 @@ ND_pairs nd_newpairs( NODE g, int t ) ndl_lcm(DL(nd_psh[r->i1]),dl,r->lcm); s = SG(nd_psh[r->i1])-TD(DL(nd_psh[r->i1])); SG(r) = MAX(s,ts) + TD(LCM(r)); + /* experimental */ + if ( nd_sugarweight ) + r->sugar2 = ndl_weight2(r->lcm); } if ( r0 ) NEXT(r) = 0; return r0; } +/* ipair = [i1,i2],[i1,i2],... */ +ND_pairs nd_ipairtospair(NODE ipair) +{ + int s1,s2; + NODE tn,t; + ND_pairs r,r0; + + for ( r0 = 0, t = ipair; t; t = NEXT(t) ) { + NEXTND_pairs(r0,r); + tn = BDY((LIST)BDY(t)); + r->i1 = QTOS((Q)ARG0(tn)); + r->i2 = QTOS((Q)ARG1(tn)); + ndl_lcm(DL(nd_psh[r->i1]),DL(nd_psh[r->i2]),r->lcm); + s1 = SG(nd_psh[r->i1])-TD(DL(nd_psh[r->i1])); + s2 = SG(nd_psh[r->i2])-TD(DL(nd_psh[r->i2])); + SG(r) = MAX(s1,s2) + TD(LCM(r)); + /* experimental */ + if ( nd_sugarweight ) + r->sugar2 = ndl_weight2(r->lcm); + } + if ( r0 ) NEXT(r) = 0; + return r0; +} + /* kokokara */ ND_pairs crit_B( ND_pairs d, int s ) @@ -2688,13 +2830,22 @@ ND_pairs nd_minp( ND_pairs d, ND_pairs *prest ) NEXT(m) = 0; return m; } - s = SG(m); if ( !NoSugar ) { - for ( ml = 0, l = m; p; p = NEXT(l = p) ) - if ( (SG(p) < s) - || ((SG(p) == s) && (DL_COMPARE(LCM(p),LCM(m)) < 0)) ) { - ml = l; m = p; s = SG(m); - } + if ( nd_sugarweight ) { + s = m->sugar2; + for ( ml = 0, l = m; p; p = NEXT(l = p) ) + if ( (p->sugar2 < s) + || ((p->sugar2 == s) && (DL_COMPARE(LCM(p),LCM(m)) < 0)) ) { + ml = l; m = p; s = m->sugar2; + } + } else { + s = SG(m); + for ( ml = 0, l = m; p; p = NEXT(l = p) ) + if ( (SG(p) < s) + || ((SG(p) == s) && (DL_COMPARE(LCM(p),LCM(m)) < 0)) ) { + ml = l; m = p; s = SG(m); + } + } } else { for ( ml = 0, l = m; p; p = NEXT(l = p) ) if ( DL_COMPARE(LCM(p),LCM(m)) < 0 ) { @@ -2715,20 +2866,37 @@ ND_pairs nd_minsugarp( ND_pairs d, ND_pairs *prest ) int msugar,i; ND_pairs t,dm0,dm,dr0,dr; - for ( msugar = SG(d), t = NEXT(d); t; t = NEXT(t) ) - if ( SG(t) < msugar ) msugar = SG(t); - dm0 = 0; dr0 = 0; - for ( i = 0, t = d; t; t = NEXT(t) ) - if ( i < nd_f4_nsp && SG(t) == msugar ) { - if ( dm0 ) NEXT(dm) = t; - else dm0 = t; - dm = t; - i++; - } else { - if ( dr0 ) NEXT(dr) = t; - else dr0 = t; - dr = t; - } + if ( nd_sugarweight ) { + for ( msugar = d->sugar2, t = NEXT(d); t; t = NEXT(t) ) + if ( t->sugar2 < msugar ) msugar = t->sugar2; + dm0 = 0; dr0 = 0; + for ( i = 0, t = d; t; t = NEXT(t) ) + if ( i < nd_f4_nsp && t->sugar2 == msugar ) { + if ( dm0 ) NEXT(dm) = t; + else dm0 = t; + dm = t; + i++; + } else { + if ( dr0 ) NEXT(dr) = t; + else dr0 = t; + dr = t; + } + } else { + for ( msugar = SG(d), t = NEXT(d); t; t = NEXT(t) ) + if ( SG(t) < msugar ) msugar = SG(t); + dm0 = 0; dr0 = 0; + for ( i = 0, t = d; t; t = NEXT(t) ) + if ( i < nd_f4_nsp && SG(t) == msugar ) { + if ( dm0 ) NEXT(dm) = t; + else dm0 = t; + dm = t; + i++; + } else { + if ( dr0 ) NEXT(dr) = t; + else dr0 = t; + dr = t; + } + } NEXT(dm) = 0; if ( dr0 ) NEXT(dr) = 0; *prest = dr0; @@ -2905,7 +3073,7 @@ void preprocess_algcoef(VL vv,VL av,struct order_spec LIST f1,f2; struct order_spec *current_spec; VECT obj,obj0; - Obj tmp; + VECT tmp; for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++); for ( nalg = 0, tv = av; tv; tv = NEXT(tv), nalg++); @@ -2919,7 +3087,7 @@ void preprocess_algcoef(VL vv,VL av,struct order_spec /* simplification, making polynomials monic */ setfield_dalg(alist); - obj_algtodalg(f,&f1); + obj_algtodalg((Obj)f,(Obj *)&f1); for ( t = BDY(f); t; t = NEXT(t) ) { initd(ord); ptod(vv,vv,(P)BDY(t),&d); hc = (DAlg)BDY(d)->c; @@ -2931,7 +3099,7 @@ void preprocess_algcoef(VL vv,VL av,struct order_spec } initd(ord); dtop(vv,vv,d,&poly); BDY(f) = (pointer)poly; } - obj_dalgtoalg(f1,&f); + obj_dalgtoalg((Obj)f1,(Obj *)&f); /* append alg vars to the var list */ for ( tv = vv; NEXT(tv); tv = NEXT(tv) ); @@ -2948,7 +3116,7 @@ void preprocess_algcoef(VL vv,VL av,struct order_spec } MKLIST(f1,t); *alistp = alist; - algobjtorat(f1,f1p); + algobjtorat((Obj)f1,(Obj *)f1p); /* creating a new weight vector */ prev_weight_vector_obj = obj0 = current_dl_weight_vector_obj; @@ -2969,7 +3137,8 @@ NODE postprocess_algcoef(VL av,NODE alist,NODE r) NODE s,t,u0,u; P p; VL tv; - Obj obj,tmp; + Obj obj; + VECT tmp; NODE arg; u0 = 0; @@ -3023,6 +3192,8 @@ void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int ndv_alloc = 0; #endif get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + if ( m && nd_vc ) + error("nd_{gr,f4} : computation over Fp(X) is unsupported. Use dp_gr_mod_main()."); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); switch ( ord->id ) { case 1: @@ -3316,15 +3487,23 @@ NDV recompute_trace(NODE ti,NDV *p,int mod) } } if ( !m ) return 0; /* failure */ - else { + else { BDY(r) = m; - c1 = invm(HCM(rm),mod); c2 = mod-HCM(r); - DMAR(c1,c2,0,mod,c); - nd_mul_c(mod,rm,c); + if ( mod > 0 || mod == -1 ) { + c1 = invm(HCM(rm),mod); c2 = mod-HCM(r); + DMAR(c1,c2,0,mod,c); + nd_mul_c(mod,rm,c); + } else { + GZ t,u; + + chsgnlf(HCZ(r),&t); + divlf(t,HCZ(rm),&u); + nd_mul_c_lf(rm,u); + } r = nd_add(mod,r,rm); - } - } - } + } + } + } } if ( tail ) NEXT(tail) = 0; d = nd_add(mod,d,r); @@ -3537,6 +3716,7 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int for ( t = fd0; t; t = NEXT(t) ) ndv_homogenize((NDV)BDY(t),obpe,oadv,oepos,ompos); } + if ( MaxDeg > 0 ) nocheck = 1; while ( 1 ) { tl1 = tl2 = tl3 = tl4 = 0; if ( Demand ) @@ -3766,7 +3946,11 @@ void nd_removecont(int mod,ND p) N q,r; if ( mod == -1 ) nd_mul_c(mod,p,_invsf(HCM(p))); - else if ( mod ) nd_mul_c(mod,p,invm(HCM(p),mod)); + else if ( mod == -2 ) { + GZ inv; + divlf(ONEGZ,HCZ(p),&inv); + nd_mul_c_lf(p,inv); + } else if ( mod ) nd_mul_c(mod,p,invm(HCM(p),mod)); else { for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ ); w = (Q *)ALLOCA(n*sizeof(Q)); @@ -3817,7 +4001,11 @@ void ndv_removecont(int mod,NDV p) if ( mod == -1 ) ndv_mul_c(mod,p,_invsf(HCM(p))); - else if ( mod ) + else if ( mod == -2 ) { + GZ inv; + divlf(ONEGZ,HCZ(p),&inv); + ndv_mul_c_lf(p,inv); + } else if ( mod ) ndv_mul_c(mod,p,invm(HCM(p),mod)); else { len = p->len; @@ -4005,6 +4193,18 @@ void nd_mul_c(int mod,ND p,int mul) } } +void nd_mul_c_lf(ND p,GZ mul) +{ + NM m; + GZ c; + + if ( !p ) return; + if ( UNIGZ(mul) ) return; + for ( m = BDY(p); m; m = NEXT(m) ) { + mullf(CZ(m),mul,&c); CZ(m) = c; + } +} + void nd_mul_c_q(ND p,P mul) { NM m; @@ -4375,8 +4575,10 @@ int nd_sp(int mod,int trace,ND_pairs p,ND *rp) if ( mod == -1 ) { CM(m1) = HCM(p2); CM(m2) = _chsgnsf(HCM(p1)); - } else if ( mod ) { + } else if ( mod > 0 ) { CM(m1) = HCM(p2); CM(m2) = mod-HCM(p1); + } else if ( mod == -2 ) { + CZ(m1) = HCZ(p2); chsgnlf(HCZ(p1),&CZ(m2)); } else if ( nd_vc ) { ezgcdpz(nd_vc,HCP(p1),HCP(p2),&gp); divsp(nd_vc,HCP(p2),gp,&CP(m1)); @@ -4414,6 +4616,19 @@ void ndv_mul_c(int mod,NDV p,int mul) } } +void ndv_mul_c_lf(NDV p,GZ mul) +{ + NMV m; + GZ c; + int len,i; + + if ( !p ) return; + len = LEN(p); + for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) { + mullf(CZ(m),mul,&c); CZ(m) = c; + } +} + void ndv_mul_c_q(NDV p,Q mul) { NMV m; @@ -4620,7 +4835,7 @@ ND ndv_mul_nm(int mod,NM m0,NDV p) if ( !p ) return 0; else if ( do_weyl ) - if ( mod == -1 ) + if ( mod < 0 ) error("ndv_mul_nm : not implemented (weyl)"); else return weyl_ndv_mul_nm(mod,m0,p); @@ -4637,6 +4852,14 @@ ND ndv_mul_nm(int mod,NM m0,NDV p) CM(mr) = _mulsf(CM(m),c); ndl_add(DL(m),d,DL(mr)); } + } else if ( mod == -2 ) { + GZ cl; + cl = CZ(m0); + for ( i = 0; i < len; i++, NMV_ADV(m) ) { + NEXTNM(mr0,mr); + mullf(CZ(m),cl,&CZ(mr)); + ndl_add(DL(m),d,DL(mr)); + } } else if ( mod ) { c = CM(m0); for ( i = 0; i < len; i++, NMV_ADV(m) ) { @@ -4676,7 +4899,12 @@ ND nd_quo(int mod,PGeoBucket bucket,NDV d) mq0 = 0; tm = (NMV)ALLOCA(nmv_adv); while ( 1 ) { - hindex = mod?head_pbucket(mod,bucket):head_pbucket_q(bucket); + if ( mod > 0 || mod == -1 ) + hindex = head_pbucket(mod,bucket); + else if ( mod == -2 ) + hindex = head_pbucket_lf(bucket); + else + hindex = head_pbucket_q(bucket); if ( hindex < 0 ) break; p = bucket->body[hindex]; NEXTNM(mq0,mq); @@ -4759,7 +4987,7 @@ NDV ndv_dup(int mod,NDV p) if ( !p ) return 0; len = LEN(p); - m0 = m = (NMV)(mod?MALLOC_ATOMIC(len*nmv_adv):MALLOC(len*nmv_adv)); + m0 = m = (NMV)((mod>0 || mod==-1)?MALLOC_ATOMIC(len*nmv_adv):MALLOC(len*nmv_adv)); for ( t = BDY(p), i = 0; i < len; i++, NMV_ADV(t), NMV_ADV(m) ) { ndl_copy(DL(t),DL(m)); CQ(m) = CQ(t); @@ -4801,7 +5029,7 @@ NDV ndv_symbolic(int mod,NDV p) if ( !p ) return 0; len = LEN(p); - m0 = m = (NMV)(mod?MALLOC_ATOMIC(len*nmv_adv):MALLOC(len*nmv_adv)); + m0 = m = (NMV)((mod>0||mod==-1)?MALLOC_ATOMIC(len*nmv_adv):MALLOC(len*nmv_adv)); for ( t = BDY(p), i = 0; i < len; i++, NMV_ADV(t), NMV_ADV(m) ) { ndl_copy(DL(t),DL(m)); CQ(m) = ONE; @@ -4873,6 +5101,13 @@ void ndv_mod(int mod,NDV p) NMV_ADV(d); dlen++; } + else if ( mod == -2 ) + for ( t = d = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) { + simp_ff((Obj)CP(t),&gfs); lmtolf(gfs,&CZ(d)); + ndl_copy(DL(t),DL(d)); + NMV_ADV(d); + dlen++; + } else for ( t = d = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) { if ( nd_vc ) { @@ -5034,8 +5269,10 @@ P ndvtop(int mod,VL vl,VL dvl,NDV p) m = (NMV)(((char *)BDY(p))+nmv_adv*(len-1)); for ( j = len-1, s = 0; j >= 0; j--, NMV_PREV(m) ) { if ( mod == -1 ) { - e = IFTOF(CM(m)); MKGFS(e,gfs); c = (P)gfs; - } else if ( mod ) { + e = IFTOF(CM(m)); MKGFS(e,gfs); c = (P)gfs; + } else if ( mod == -2 ) { + c = gztoz(CZ(m)); + } else if ( mod > 0 ) { STOQ(CM(m),q); c = (P)q; } else c = CP(m); @@ -5103,7 +5340,7 @@ NDV ndtondv(int mod,ND p) if ( !p ) return 0; len = LEN(p); - if ( mod ) + if ( mod > 0 || mod == -1 ) m0 = m = (NMV)MALLOC_ATOMIC_IGNORE_OFF_PAGE(len*nmv_adv); else m0 = m = MALLOC(len*nmv_adv); @@ -5249,7 +5486,16 @@ NODE ndv_reducebase(NODE x,int *perm) void nd_init_ord(struct order_spec *ord) { - nd_module = (ord->id >= 256); + nd_module = (ord->id >= 256); + if ( nd_module ) { + nd_dcomp = -1; + nd_ispot = ord->ispot; + nd_pot_nelim = ord->pot_nelim; + nd_poly_weight_len = ord->nv; + nd_poly_weight = ord->top_weight; + nd_module_rank = ord->module_rank; + nd_module_weight = ord->module_top_weight; + } nd_matrix = 0; nd_matrix_len = 0; switch ( ord->id ) { @@ -5305,9 +5551,6 @@ void nd_init_ord(struct order_spec *ord) /* module order */ case 256: - nd_ispot = ord->ispot; - nd_pot_nelim = ord->pot_nelim; - nd_dcomp = -1; switch ( ord->ord.simple ) { case 0: nd_isrlex = 1; @@ -5327,17 +5570,11 @@ void nd_init_ord(struct order_spec *ord) break; case 257: /* block order */ - nd_ispot = ord->ispot; - nd_pot_nelim = ord->pot_nelim; - nd_dcomp = -1; nd_isrlex = 0; ndl_compare_function = ndl_module_block_compare; break; case 258: /* matrix order */ - nd_ispot = ord->ispot; - nd_pot_nelim = ord->pot_nelim; - nd_dcomp = -1; nd_isrlex = 0; nd_matrix_len = ord->ord.matrix.row; nd_matrix = ord->ord.matrix.matrix; @@ -5345,9 +5582,6 @@ void nd_init_ord(struct order_spec *ord) break; case 259: /* composite order */ - nd_ispot = ord->ispot; - nd_pot_nelim = ord->pot_nelim; - nd_dcomp = -1; nd_isrlex = 0; nd_worb_len = ord->ord.composite.length; nd_worb = ord->ord.composite.w_or_b; @@ -5532,6 +5766,22 @@ int nd_to_vect_q(UINT *s0,int n,ND d,Q *r) return i; } +int nd_to_vect_lf(UINT *s0,int n,ND d,mpz_t *r) +{ + NM m; + UINT *t,*s; + int i; + + for ( i = 0; i < n; i++ ) { mpz_init(r[i]); mpz_set_ui(r[i],0); } + for ( i = 0, s = s0, m = BDY(d); m; m = NEXT(m) ) { + t = DL(m); + for ( ; !ndl_equal(t,s); s += nd_wpd, i++ ); + mpz_set(r[i],BDY(CZ(m))); + } + for ( i = 0; !mpz_sgn(r[i]); i++ ); + return i; +} + unsigned long *nd_to_vect_2(UINT *s0,int n,int *s0hash,ND p) { NM m; @@ -5907,7 +6157,7 @@ int ndv_reduce_vect_sf(int m,UINT *svect,int col,IndAr maxrs = 0; for ( i = 0; i < nred; i++ ) { ivect = imat[i]; - k = ivect->head; svect[k] %= m; + k = ivect->head; if ( c = svect[k] ) { maxrs = MAX(maxrs,rp0[i]->sugar); c = _chsgnsf(c); redv = nd_ps[rp0[i]->index]; @@ -5941,6 +6191,183 @@ int ndv_reduce_vect_sf(int m,UINT *svect,int col,IndAr return maxrs; } +ND nd_add_lf(ND p1,ND p2) +{ + int n,c,can; + ND r; + NM m1,m2,mr0,mr,s; + GZ t; + + if ( !p1 ) return p2; + else if ( !p2 ) return p1; + else { + can = 0; + for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) { + c = DL_COMPARE(DL(m1),DL(m2)); + switch ( c ) { + case 0: + addlf(CZ(m1),CZ(m2),&t); + s = m1; m1 = NEXT(m1); + if ( t ) { + can++; NEXTNM2(mr0,mr,s); CZ(mr) = (t); + } else { + can += 2; FREENM(s); + } + s = m2; m2 = NEXT(m2); FREENM(s); + break; + case 1: + s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s); + break; + case -1: + s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s); + break; + } + } + if ( !mr0 ) + if ( m1 ) mr0 = m1; + else if ( m2 ) mr0 = m2; + else return 0; + else if ( m1 ) NEXT(mr) = m1; + else if ( m2 ) NEXT(mr) = m2; + else NEXT(mr) = 0; + BDY(p1) = mr0; + SG(p1) = MAX(SG(p1),SG(p2)); + LEN(p1) = LEN(p1)+LEN(p2)-can; + FREEND(p2); + return p1; + } +} + +int ndv_reduce_vect_lf(mpz_t *svect,int col,IndArray *imat,NM_ind_pair *rp0,int nred) +{ + int i,j,k,len,pos,prev; + mpz_t c,mc,c1; + IndArray ivect; + unsigned char *ivc; + unsigned short *ivs; + unsigned int *ivi; + NDV redv; + NMV mr; + NODE rp; + int maxrs; + + maxrs = 0; + lf_lazy = 1; + for ( i = 0; i < nred; i++ ) { + ivect = imat[i]; + k = ivect->head; + mpz_mod(svect[k],svect[k],BDY(current_mod_lf)); + if ( mpz_sgn(svect[k]) ) { + maxrs = MAX(maxrs,rp0[i]->sugar); + mpz_neg(svect[k],svect[k]); + redv = nd_ps[rp0[i]->index]; + len = LEN(redv); mr = BDY(redv); + prev = k; + switch ( ivect->width ) { + case 1: + ivc = ivect->index.c; + for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { + pos = prev+ivc[j]; prev = pos; + mpz_addmul(svect[pos],svect[k],BDY(CZ(mr))); + } + break; + case 2: + ivs = ivect->index.s; + for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { + pos = prev+ivs[j]; prev = pos; + mpz_addmul(svect[pos],svect[k],BDY(CZ(mr))); + } + break; + case 4: + ivi = ivect->index.i; + for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) { + pos = prev+ivi[j]; prev = pos; + mpz_addmul(svect[pos],svect[k],BDY(CZ(mr))); + } + break; + } + mpz_set_ui(svect[k],0); + } + } + lf_lazy=0; + for ( i = 0; i < col; i++ ) { + mpz_mod(svect[i],svect[i],BDY(current_mod_lf)); + } + return maxrs; +} + +int nd_gauss_elim_lf(mpz_t **mat0,int *sugar,int row,int col,int *colstat) +{ + int i,j,k,l,rank,s; + mpz_t a,a1,inv; + mpz_t *t,*pivot,*pk; + mpz_t **mat; + struct oEGT eg0,eg1,eg_forward,eg_mod,eg_back; + + mpz_init(inv); + mpz_init(a); + mat = (mpz_t **)mat0; + for ( rank = 0, j = 0; j < col; j++ ) { + for ( i = rank; i < row; i++ ) { + mpz_mod(mat[i][j],mat[i][j],BDY(current_mod_lf)); + } + for ( i = rank; i < row; i++ ) + if ( mpz_sgn(mat[i][j]) ) + break; + if ( i == row ) { + colstat[j] = 0; + continue; + } else + colstat[j] = 1; + if ( i != rank ) { + t = mat[i]; mat[i] = mat[rank]; mat[rank] = t; + s = sugar[i]; sugar[i] = sugar[rank]; sugar[rank] = s; + } + pivot = mat[rank]; + s = sugar[rank]; + mpz_invert(inv,pivot[j],BDY(current_mod_lf)); + for ( k = j, pk = pivot+k; k < col; k++, pk++ ) + if ( mpz_sgn(*pk) ) { + mpz_mul(a,*pk,inv); mpz_mod(*pk,a,BDY(current_mod_lf)); + } + for ( i = rank+1; i < row; i++ ) { + t = mat[i]; + if ( mpz_sgn(t[j]) ) { + sugar[i] = MAX(sugar[i],s); + mpz_neg(a,t[j]); + red_by_vect_lf(t+j,pivot+j,a,col-j); + } + } + rank++; + } + for ( j = col-1, l = rank-1; j >= 0; j-- ) + if ( colstat[j] ) { + pivot = mat[l]; + s = sugar[l]; + for ( k = j; k < col; k++ ) + mpz_mod(pivot[k],pivot[k],BDY(current_mod_lf)); + for ( i = 0; i < l; i++ ) { + t = mat[i]; + if ( mpz_sgn(t[j]) ) { + sugar[i] = MAX(sugar[i],s); + mpz_neg(a,t[j]); + red_by_vect_lf(t+j,pivot+j,a,col-j); + } + } + l--; + } + for ( j = 0, l = 0; l < rank; j++ ) + if ( colstat[j] ) { + t = mat[l]; + for ( k = j; k < col; k++ ) { + mpz_mod(t[k],t[k],BDY(current_mod_lf)); + } + l++; + } + return rank; +} + + NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhead,UINT *s0vect) { int j,k,len; @@ -6027,7 +6454,7 @@ NDV vect_to_ndv_gz(GZ *vect,int spcol,int col,int *rhe { int j,k,len; UINT *p; - Q c; + GZ c; NDV r; NMV mr0,mr; @@ -6051,6 +6478,35 @@ NDV vect_to_ndv_gz(GZ *vect,int spcol,int col,int *rhe } } +NDV vect_to_ndv_lf(mpz_t *vect,int spcol,int col,int *rhead,UINT *s0vect) +{ + int j,k,len; + UINT *p; + mpz_t c; + NDV r; + NMV mr0,mr; + + for ( j = 0, len = 0; j < spcol; j++ ) if ( mpz_sgn(vect[j]) ) len++; + if ( !len ) return 0; + else { + mr0 = (NMV)MALLOC(nmv_adv*len); +#if 0 + ndv_alloc += nmv_adv*len; +#endif + mr = mr0; + p = s0vect; + for ( j = k = 0; j < col; j++, p += nd_wpd ) + if ( !rhead[j] ) { + c[0] = vect[k++][0]; + if ( mpz_sgn(c) ) { + ndl_copy(p,DL(mr)); MPZTOGZ(c,CZ(mr)); NMV_ADV(mr); + } + } + MKNDV(nd_nvar,mr0,len,r); + return r; + } +} + /* for plain vector */ NDV plain_vect_to_ndv_q(Q *vect,int col,UINT *s0vect) @@ -6153,7 +6609,7 @@ NODE nd_f4(int m,int **indp) ND spol,red; NDV nf,redv; NM s0,s; - NODE rp0,srp0,nflist,nzlist; + NODE rp0,srp0,nflist,nzlist,nzlist_t; int nsp,nred,col,rank,len,k,j,a,i1s,i2s; UINT c; UINT **spmat; @@ -6171,42 +6627,36 @@ NODE nd_f4(int m,int **indp) #endif g = 0; d = 0; for ( i = 0; i < nd_psn; i++ ) { - d = update_pairs(d,g,i,0); + if ( !nd_nzlist ) d = update_pairs(d,g,i,0); g = update_base(g,i); } nzlist = 0; - while ( d ) { + nzlist_t = nd_nzlist; + while ( d || nzlist_t ) { get_eg(&eg0); - l = nd_minsugarp(d,&d); - sugar = SG(l); if ( nd_nzlist ) { - for ( tn = nd_nzlist; tn; tn = NEXT(tn) ) { - node = BDY((LIST)BDY(tn)); - if ( QTOS((Q)ARG0(node)) == sugar ) break; + node = BDY((LIST)BDY(nzlist_t)); + sugar = (int)ARG0(node); + tn = BDY((LIST)ARG1(node)); + if ( !tn ) { + nzlist_t = NEXT(nzlist_t); + continue; } - if ( tn ) { - for ( t = l, ll0 = 0; t; t = NEXT(t) ) { - for ( tn = BDY((LIST)ARG1(node)); tn; tn = NEXT(tn) ) { - i1s = QTOS((Q)ARG0(BDY((LIST)BDY(tn)))); - i2s = QTOS((Q)ARG1(BDY((LIST)BDY(tn)))); - if ( t->i1 == i1s && t->i2 == i2s ) break; - } - if ( tn ) { - if ( !ll0 ) ll0 = t; - else NEXT(ll) = t; - ll = t; - } - } - if ( ll0 ) NEXT(ll) = 0; - l = ll0; - } else l = 0; + /* tn = [[i1,i2],...] */ + l = nd_ipairtospair(tn); + } else { + l = nd_minsugarp(d,&d); + sugar = nd_sugarweight?l->sugar2:SG(l); + if ( MaxDeg > 0 && sugar > MaxDeg ) break; } bucket = create_pbucket(); stat = nd_sp_f4(m,0,l,bucket); if ( !stat ) { - for ( t = l; NEXT(t); t = NEXT(t) ); - NEXT(t) = d; d = l; - d = nd_reconstruct(0,d); + if ( !nd_nzlist ) { + for ( t = l; NEXT(t); t = NEXT(t) ); + NEXT(t) = d; d = l; + d = nd_reconstruct(0,d); + } continue; } if ( bucket->m < 0 ) continue; @@ -6235,7 +6685,7 @@ NODE nd_f4(int m,int **indp) nf = ndtondv(m,nf1); } nh = ndv_newps(m,nf,0,1); - d = update_pairs(d,g,nh,0); + if ( !nd_nzlist ) d = update_pairs(d,g,nh,0); g = update_base(g,nh); } if ( nd_gentrace ) { @@ -6249,6 +6699,7 @@ NODE nd_f4(int m,int **indp) STOQ(sugar,sugarq); node = mknode(2,sugarq,l0); MKLIST(l1,node); MKNODE(node,l1,nzlist); nzlist = node; } + if ( nd_nzlist ) nzlist_t = NEXT(nzlist_t); } if ( nd_gentrace ) { MKLIST(l0,reverse_node(nzlist)); @@ -6291,6 +6742,7 @@ NODE nd_f4_trace(int m,int **indp) get_eg(&eg0); l = nd_minsugarp(d,&d); sugar = SG(l); + if ( MaxDeg > 0 && sugar > MaxDeg ) break; bucket = create_pbucket(); stat = nd_sp_f4(m,0,l,bucket); if ( !stat ) { @@ -6600,7 +7052,9 @@ init_eg(&eg_search); } } mat->row = i; - fprintf(asir_out,"%dx%d,",mat->row,mat->col); fflush(asir_out); + if ( DP_Print ) { + fprintf(asir_out,"%dx%d,",mat->row,mat->col); fflush(asir_out); + } size = ((col+BLEN-1)/BLEN)*sizeof(unsigned long); v = CALLOC((col+BLEN-1)/BLEN,sizeof(unsigned long)); for ( rp = rp0, i = 0; rp; rp = NEXT(rp), i++ ) { @@ -6654,7 +7108,9 @@ init_eg(&eg_search); for ( i = 0; i < col; i++ ) rhead[i] = 0; /* construction of index arrays */ - fprintf(stderr,"%dx%d,",nsp+nred,col); + if ( DP_Print ) { + fprintf(stderr,"%dx%d,",nsp+nred,col); + } rvect = (NM_ind_pair *)ALLOCA(nred*sizeof(NM_ind_pair)); s0hash = (int *)ALLOCA(col*sizeof(int)); for ( i = 0, s = s0vect; i < col; i++, s += nd_wpd ) @@ -6664,14 +7120,17 @@ init_eg(&eg_search); imat[i] = nm_ind_pair_to_vect_compress(m,s0vect,col,s0hash,rvect[i]); rhead[imat[i]->head] = 1; } - if ( m ) + if ( m > 0 || m == -1 ) r0 = nd_f4_red_main(m,sp0,nsp,s0vect,col,rvect,rhead,imat,nred,nz); + else if ( m == -2 ) + r0 = nd_f4_red_lf_main(m,sp0,nsp,s0vect,col,rvect,rhead,imat,nred); else r0 = nd_f4_red_gz_main(sp0,nsp,trace,s0vect,col,rvect,rhead,imat,nred); -print_eg("search",&eg_search); + if ( DP_Print ) print_eg("search",&eg_search); return r0; } +/* for small finite fields */ NODE nd_f4_red_main(int m,ND_pairs sp0,int nsp,UINT *s0vect,int col, NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred,ND_pairs *nz) { @@ -6693,7 +7152,7 @@ NODE nd_f4_red_main(int m,ND_pairs sp0,int nsp,UINT *s /* elimination (1st step) */ spmat = (int **)ALLOCA(nsp*sizeof(UINT *)); svect = (UINT *)ALLOCA(col*sizeof(UINT)); - spsugar = (int *)ALLOCA(nsp*sizeof(UINT)); + spsugar = (int *)ALLOCA(nsp*sizeof(int)); spactive = !nz?0:(ND_pairs *)ALLOCA(nsp*sizeof(ND_pairs)); for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) { nd_sp(m,0,sp,&spol); @@ -6758,6 +7217,88 @@ NODE nd_f4_red_main(int m,ND_pairs sp0,int nsp,UINT *s return r0; } +NODE nd_f4_red_lf_main(int m,ND_pairs sp0,int nsp,UINT *s0vect,int col, + NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred) +{ + int spcol,sprow,a; + int i,j,k,l,rank; + NODE r0,r; + ND_pairs sp; + ND spol; + mpz_t **spmat; + mpz_t *svect,*v; + int *colstat; + struct oEGT eg0,eg1,eg2,eg_f4,eg_f4_1,eg_f4_2; + int maxrs; + int *spsugar; + pointer *w; + + spcol = col-nred; + get_eg(&eg0); + /* elimination (1st step) */ + spmat = (mpz_t **)ALLOCA(nsp*sizeof(mpz_t *)); + svect = (mpz_t *)ALLOCA(col*sizeof(mpz_t)); + spsugar = (int *)ALLOCA(nsp*sizeof(int)); + for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) { + nd_sp(m,0,sp,&spol); + if ( !spol ) continue; + nd_to_vect_lf(s0vect,col,spol,svect); + maxrs = ndv_reduce_vect_lf(svect,col,imat,rvect,nred); + for ( i = 0; i < col; i++ ) if ( mpz_sgn(svect[i]) ) break; + if ( i < col ) { + spmat[sprow] = v = (mpz_t *)MALLOC(spcol*sizeof(mpz_t)); + for ( j = k = 0; j < col; j++ ) + if ( !rhead[j] ) v[k++][0] = svect[j][0]; + spsugar[sprow] = MAX(maxrs,SG(spol)); + sprow++; + } +/* nd_free(spol); */ + } + get_eg(&eg1); init_eg(&eg_f4_1); add_eg(&eg_f4_1,&eg0,&eg1); + if ( DP_Print ) { + fprintf(asir_out,"elim1=%fsec,",eg_f4_1.exectime+eg_f4_1.gctime); + fflush(asir_out); + } + /* free index arrays */ +/* for ( i = 0; i < nred; i++ ) GCFREE(imat[i]->index.c); */ + + /* elimination (2nd step) */ + colstat = (int *)ALLOCA(spcol*sizeof(int)); + rank = nd_gauss_elim_lf(spmat,spsugar,sprow,spcol,colstat); + w = (pointer *)ALLOCA(rank*sizeof(pointer)); + for ( i = 0; i < rank; i++ ) { +#if 0 + w[rank-i-1] = (pointer)vect_to_ndv_lf(spmat[i],spcol,col,rhead,s0vect); + SG((NDV)w[rank-i-1]) = spsugar[i]; +#else + w[i] = (pointer)vect_to_ndv_lf(spmat[i],spcol,col,rhead,s0vect); + SG((NDV)w[i]) = spsugar[i]; +#endif +/* GCFREE(spmat[i]); */ + + } +#if 0 + qsort(w,rank,sizeof(NDV), + (int (*)(const void *,const void *))ndv_compare); +#endif + r0 = 0; + for ( i = 0; i < rank; i++ ) { + NEXTNODE(r0,r); BDY(r) = w[i]; + } + 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); + init_eg(&eg_f4); add_eg(&eg_f4,&eg0,&eg2); + if ( DP_Print ) { + fprintf(asir_out,"elim2=%fsec\n",eg_f4_2.exectime+eg_f4_2.gctime); + fprintf(asir_out,"nsp=%d,nred=%d,spmat=(%d,%d),rank=%d ", + nsp,nred,sprow,spcol,rank); + fprintf(asir_out,"%fsec\n",eg_f4.exectime+eg_f4.gctime); + } + return r0; +} + #if 1 NODE nd_f4_red_q_main(ND_pairs sp0,int nsp,int trace,UINT *s0vect,int col, NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred) @@ -6780,7 +7321,7 @@ NODE nd_f4_red_q_main(ND_pairs sp0,int nsp,int trace,U /* elimination (1st step) */ spmat = (Q **)ALLOCA(nsp*sizeof(Q *)); svect = (Q *)ALLOCA(col*sizeof(Q)); - spsugar = (int *)ALLOCA(nsp*sizeof(Q)); + spsugar = (int *)ALLOCA(nsp*sizeof(int)); for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) { nd_sp(0,trace,sp,&spol); if ( !spol ) continue; @@ -6861,7 +7402,7 @@ NODE nd_f4_red_gz_main(ND_pairs sp0,int nsp,int trace, /* elimination (1st step) */ spmat = (GZ **)ALLOCA(nsp*sizeof(GZ *)); svect = (GZ *)ALLOCA(col*sizeof(GZ)); - spsugar = (int *)ALLOCA(nsp*sizeof(Q)); + spsugar = (int *)ALLOCA(nsp*sizeof(int)); for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) { nd_sp(0,trace,sp,&spol); if ( !spol ) continue; @@ -6896,7 +7437,7 @@ NODE nd_f4_red_gz_main(ND_pairs sp0,int nsp,int trace, w[rank-i-1] = ndvgztondv(w[rank-i-1]); SG((NDV)w[rank-i-1]) = spsugar[i]; #else - w[i] = (pointer)vect_to_ndv_gz((Q *)spmat[i],spcol,col,rhead,s0vect); + w[i] = (pointer)vect_to_ndv_gz(spmat[i],spcol,col,rhead,s0vect); w[i] = ndvgztondv(w[i]); SG((NDV)w[i]) = spsugar[i]; #endif @@ -7207,7 +7748,7 @@ void nd_exec_f4_red_dist() /* elimination (1st step) */ spmat = (int **)MALLOC(nsp*sizeof(UINT *)); svect = (UINT *)MALLOC(col*sizeof(UINT)); - spsugar = (int *)ALLOCA(nsp*sizeof(UINT)); + spsugar = (int *)ALLOCA(nsp*sizeof(int)); for ( a = sprow = 0; a < nsp; a++ ) { nd_sp(m,0,sp0[a],&spol); if ( !spol ) continue; @@ -7962,6 +8503,7 @@ P ndc_div(int mod,union oNDC a,union oNDC b) int inv,t; if ( mod == -1 ) c.m = _mulsf(a.m,_invsf(b.m)); + else if ( mod == -2 ) divlf(a.z,b.z,&c.z); else if ( mod ) { inv = invm(b.m,mod); DMAR(a.m,inv,0,mod,t); c.m = t; @@ -7980,7 +8522,7 @@ P ndctop(int mod,union oNDC c) if ( mod == -1 ) { e = IFTOF(c.m); MKGFS(e,gfs); return (P)gfs; - } else if ( mod ) { + } else if ( mod > 0 ) { STOQ(c.m,q); return (P)q; } else return (P)c.p; @@ -8025,13 +8567,15 @@ void conv_ilist(int demand,int trace,NODE g,int **indp void parse_nd_option(NODE opt) { NODE t,p,u; - int i,s; + int i,s,n; char *key; Obj value; nd_gentrace = 0; nd_gensyz = 0; nd_nora = 0; nd_gbblock = 0; nd_newelim = 0; nd_intersect = 0; nd_nzlist = 0; nd_splist = 0; nd_check_splist = 0; + nd_sugarweight = 0; + for ( t = opt; t; t = NEXT(t) ) { p = BDY((LIST)BDY(t)); key = BDY((STRING)BDY(p)); @@ -8043,16 +8587,17 @@ void parse_nd_option(NODE opt) else if ( !strcmp(key,"nora") ) nd_nora = value?1:0; else if ( !strcmp(key,"gbblock") ) { - if ( !value || OID(value) != O_LIST ) - error("nd_* : invalid value for gbblock option"); - u = BDY((LIST)value); - nd_gbblock = MALLOC((2*length(u)+1)*sizeof(int)); - for ( i = 0; u; u = NEXT(u) ) { - p = BDY((LIST)BDY(u)); - s = nd_gbblock[i++] = QTOS((Q)BDY(p)); - nd_gbblock[i++] = s+QTOS((Q)BDY(NEXT(p)))-1; - } - nd_gbblock[i] = -1; + if ( value && OID(value) == O_LIST ) { + u = BDY((LIST)value); + nd_gbblock = MALLOC((2*length(u)+1)*sizeof(int)); + for ( i = 0; u; u = NEXT(u) ) { + p = BDY((LIST)BDY(u)); + s = nd_gbblock[i++] = QTOS((Q)BDY(p)); + nd_gbblock[i++] = s+QTOS((Q)BDY(NEXT(p)))-1; + } + nd_gbblock[i] = -1; + } else + nd_gbblock = 0; } else if ( !strcmp(key,"newelim") ) nd_newelim = value?1:0; else if ( !strcmp(key,"intersect") ) @@ -8065,6 +8610,12 @@ void parse_nd_option(NODE opt) nd_splist = value?1:0; else if ( !strcmp(key,"check_splist") ) { nd_check_splist = BDY((LIST)value); + } else if ( !strcmp(key,"sugarweight") ) { + u = BDY((LIST)value); + n = length(u); + nd_sugarweight = MALLOC(n*sizeof(int)); + for ( i = 0; i < n; i++, u = NEXT(u) ) + nd_sugarweight[i] = QTOS((Q)BDY(u)); } } } @@ -8346,4 +8897,22 @@ VECT nd_btog_one(LIST f,LIST v,int mod,struct order_sp BDY(vect)[j] = ndtodp(mod,u); } return vect; +} + +void ndv_print_lf(NDV p) +{ + NMV m; + int i,len; + + if ( !p ) printf("0\n"); + else { + len = LEN(p); + for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) { + printf("+"); + mpz_out_str(asir_out,10,BDY(CZ(m))); + printf("*"); + ndl_print(DL(m)); + } + printf("\n"); + } }