=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v retrieving revision 1.62 retrieving revision 1.65 diff -u -p -r1.62 -r1.65 --- OpenXM_contrib2/asir2000/engine/nd.c 2003/09/11 01:52:25 1.62 +++ OpenXM_contrib2/asir2000/engine/nd.c 2003/09/12 08:01:51 1.65 @@ -1,7 +1,8 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.61 2003/09/10 05:14:32 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.64 2003/09/12 01:12:40 noro Exp $ */ #include "ca.h" #include "inline.h" +#include #if defined(__GNUC__) #define INLINE inline @@ -14,6 +15,7 @@ typedef unsigned int UINT; #define USE_GEOBUCKET 1 +#define USE_UNROLL 1 #define REDTAB_LEN 32003 @@ -95,6 +97,13 @@ typedef struct oBaseSet { UINT **bound; } *BaseSet; +typedef struct oNM_ind_pair +{ + NM mul; + int index; +} *NM_ind_pair; + + int (*ndl_compare_function)(UINT *a1,UINT *a2); static double nd_scale=2; @@ -179,6 +188,8 @@ if(!_nd_free_list)_ND_alloc();\ NV(d)=(n); LEN(d)=(len); BDY(d)=(m) #define NEWNDV(d) ((d)=(NDV)MALLOC(sizeof(struct oNDV))) #define MKNDV(n,m,l,d) NEWNDV(d); NV(d)=(n); BDY(d)=(m); LEN(d) = l; +#define NEWNM_ind_pair(p)\ +((p)=(NM_ind_pair)MALLOC(sizeof(struct oNM_ind_pair))) /* allocate and link a new object */ #define NEXTRHist(r,c) \ @@ -189,6 +200,7 @@ if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEX if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);} #define NEXTND_pairs(r,c) \ if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);} +#define MKNM_ind_pair(p,m,i) (NEWNM_ind_pair(p),(p)->mul=(m),(p)->index=(i)) /* deallocators */ #define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m) @@ -202,7 +214,6 @@ if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT #define NMV_PREV(m) (m = (NMV)(((char *)m)-nmv_adv)) #define NMV_OPREV(m) (m = (NMV)(((char *)m)-oadv)) - /* external functions */ void GC_gcollect(); NODE append_one(NODE,int); @@ -216,15 +227,18 @@ void removecont_array(Q *c,int n); ND normalize_pbucket(int mod,PGeoBucket g); int head_pbucket(int mod,PGeoBucket g); int head_pbucket_q(PGeoBucket g); +void add_pbucket_symbolic(PGeoBucket g,ND d); void add_pbucket(int mod,PGeoBucket g,ND d); void free_pbucket(PGeoBucket b); void mulq_pbucket(PGeoBucket g,Q c); +NM remove_head_pbucket_symbolic(PGeoBucket g); PGeoBucket create_pbucket(); /* manipulation of pairs and bases */ int nd_newps(int mod,ND a,ND aq); ND_pairs nd_newpairs( NODE g, int t ); ND_pairs nd_minp( ND_pairs d, ND_pairs *prest ); +ND_pairs nd_minsugarp( ND_pairs d, ND_pairs *prest ); NODE update_base(NODE nd,int ndp); ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t); ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest ); @@ -232,10 +246,12 @@ ND_pairs crit_B( ND_pairs d, int s ); ND_pairs crit_M( ND_pairs d1 ); ND_pairs crit_F( ND_pairs d1 ); int crit_2( int dp1, int dp2 ); +int ndv_newps(NDV a,NDV aq); /* top level functions */ -void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp); +void nd_gr(LIST f,LIST v,int m,int f4,struct order_spec *ord,LIST *rp); void nd_gr_trace(LIST f,LIST v,int trace,int homo,struct order_spec *ord,LIST *rp); +NODE nd_f4(int m); NODE nd_gb(int m,int checkonly); NODE nd_gb_trace(int m); @@ -257,8 +273,8 @@ INLINE void ndl_sub(UINT *d1,UINT *d2,UINT *d); INLINE int ndl_hash_value(UINT *d); /* normal forms */ -INLINE int nd_find_reducer(ND g); -INLINE int nd_find_reducer_direct(ND g,NDV *ps,int len); +INLINE int ndl_find_reducer(UINT *g); +INLINE int ndl_find_reducer_direct(UINT *g,NDV *ps,int len); int nd_sp(int mod,int trace,ND_pairs p,ND *nf); int nd_nf(int mod,ND g,NDV *ps,int full,ND *nf); int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *nf); @@ -293,6 +309,7 @@ EPOS nd_create_epos(struct order_spec *ord); int nd_get_exporigin(struct order_spec *ord); void ndv_mod(int mod,NDV p); NDV ndv_dup(int mod,NDV p); +ND nd_dup(ND p); /* ND functions */ int ndv_check_candidate(NODE input,int obpe,int oadv,EPOS oepos,NODE cand); @@ -304,6 +321,7 @@ int nd_length(ND p); void nd_append_red(UINT *d,int i); UINT *ndv_compute_bound(NDV p); ND nd_copy(ND p); +ND nd_merge(ND p1,ND p2); ND nd_add(int mod,ND p1,ND p2); ND nd_add_q(ND p1,ND p2); INLINE int nd_length(ND p); @@ -313,6 +331,7 @@ ND weyl_ndv_mul_nm(int mod,NM m0,NDV p); void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *tab,int tlen); void ndv_mul_c(int mod,NDV p,int mul); void ndv_mul_c_q(NDV p,Q mul); +ND ndv_mul_nm_symbolic(NM m0,NDV p); ND ndv_mul_nm(int mod,NM m0,NDV p); void ndv_realloc(NDV p,int obpe,int oadv,EPOS oepos); NDV ndv_dup_realloc(NDV p,int obpe,int oadv,EPOS oepos); @@ -329,6 +348,9 @@ NDV ptondv(VL vl,VL dvl,P p); P ndvtop(int mod,VL vl,VL dvl,NDV p); NDV ndtondv(int mod,ND p); ND ndvtond(int mod,NDV p); +int nm_ind_pair_to_vect(int m,UINT *s0,int n,NM_ind_pair pair,UINT *r); +void nm_ind_pair_to_vect_compress(int m,UINT *s0,int n,NM_ind_pair pair,int *ind); +int nd_to_vect(int mod,UINT *s0,int n,ND d,UINT *r); void nd_free_private_storage() { @@ -392,6 +414,7 @@ INLINE int ndl_reducible(UINT *d1,UINT *d2) int i,j; if ( TD(d1) < TD(d2) ) return 0; +#if USE_UNROLL switch ( nd_bpe ) { case 3: for ( i = nd_exporigin; i < nd_wpd; i++ ) { @@ -465,6 +488,14 @@ INLINE int ndl_reducible(UINT *d1,UINT *d2) } return 1; } +#else + for ( i = nd_exporigin; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + for ( j = 0; j < nd_epw; j++ ) + if ( (u1&nd_mask[j]) < (u2&nd_mask[j]) ) return 0; + } + return 1; +#endif } /* @@ -532,6 +563,7 @@ void ndl_lcm(UINT *d1,unsigned *d2,UINT *d) UINT t1,t2,u,u1,u2; int i,j,l; +#if USE_UNROLL switch ( nd_bpe ) { case 3: for ( i = nd_exporigin; i < nd_wpd; i++ ) { @@ -608,6 +640,15 @@ void ndl_lcm(UINT *d1,unsigned *d2,UINT *d) } break; } +#else + for ( i = nd_exporigin; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + for ( j = 0, u = 0; j < nd_epw; j++ ) { + t1 = (u1&nd_mask[j]); t2 = (u2&nd_mask[j]); u |= t1>t2?t1:t2; + } + d[i] = u; + } +#endif TD(d) = ndl_weight(d); if ( nd_blockmask ) ndl_weight_mask(d); } @@ -813,6 +854,7 @@ int ndl_disjoint(UINT *d1,UINT *d2) UINT t1,t2,u,u1,u2; int i,j; +#if USE_UNROLL switch ( nd_bpe ) { case 3: for ( i = nd_exporigin; i < nd_wpd; i++ ) { @@ -889,6 +931,16 @@ int ndl_disjoint(UINT *d1,UINT *d2) return 1; break; } +#else + for ( i = nd_exporigin; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + for ( j = 0; j < nd_epw; j++ ) { + if ( (u1&nd_mask0) && (u2&nd_mask0) ) return 0; + u1 >>= nd_bpe; u2 >>= nd_bpe; + } + } + return 1; +#endif } int ndl_check_bound2(int index,UINT *d2) @@ -899,6 +951,7 @@ int ndl_check_bound2(int index,UINT *d2) d1 = nd_bound[index]; ind = 0; +#if USE_UNROLL switch ( nd_bpe ) { case 3: for ( i = nd_exporigin; i < nd_wpd; i++ ) { @@ -974,6 +1027,15 @@ int ndl_check_bound2(int index,UINT *d2) return 0; break; } +#else + for ( i = nd_exporigin; i < nd_wpd; i++ ) { + u2 = d2[i]; + k = (nd_epw-1)*nd_bpe; + for ( j = 0; j < nd_epw; j++, k -= nd_bpe ) + if ( d1[ind++]+((u2>>k)&nd_mask0) > nd_mask0 ) return 1; + } + return 0; +#endif } int ndl_check_bound2_direct(UINT *d1,UINT *d2) @@ -982,6 +1044,7 @@ int ndl_check_bound2_direct(UINT *d1,UINT *d2) int i,j,ind,k; ind = 0; +#if USE_UNROLL switch ( nd_bpe ) { case 3: for ( i = nd_exporigin; i < nd_wpd; i++ ) { @@ -1057,6 +1120,15 @@ int ndl_check_bound2_direct(UINT *d1,UINT *d2) return 0; break; } +#else + for ( i = nd_exporigin; i < nd_wpd; i++ ) { + u2 = d2[i]; + k = (nd_epw-1)*nd_bpe; + for ( j = 0; j < nd_epw; j++, k -= nd_bpe ) + if ( d1[ind++]+((u2>>k)&nd_mask0) > nd_mask0 ) return 1; + } + return 0; +#endif } INLINE int ndl_hash_value(UINT *d) @@ -1070,15 +1142,12 @@ INLINE int ndl_hash_value(UINT *d) return r; } -INLINE int nd_find_reducer(ND g) +INLINE int ndl_find_reducer(UINT *dg) { RHist r; - UINT *dg; int d,k,i; - dg = HDL(g); -#if 1 - d = ndl_hash_value(HDL(g)); + d = ndl_hash_value(dg); for ( r = nd_red[d], k = 0; r; r = NEXT(r), k++ ) { if ( ndl_equal(dg,DL(r)) ) { if ( k > 0 ) nd_notfirst++; @@ -1086,7 +1155,6 @@ INLINE int nd_find_reducer(ND g) return r->index; } } -#endif if ( Reverse ) for ( i = nd_psn-1; i >= 0; i-- ) { r = nd_psh[i]; @@ -1108,7 +1176,7 @@ INLINE int nd_find_reducer(ND g) return -1; } -INLINE int nd_find_reducer_direct(ND g,NDV *ps,int len) +INLINE int ndl_find_reducer_direct(UINT *dg,NDV *ps,int len) { NDV r; RHist s; @@ -1117,18 +1185,60 @@ INLINE int nd_find_reducer_direct(ND g,NDV *ps,int len if ( Reverse ) for ( i = len-1; i >= 0; i-- ) { r = ps[i]; - if ( ndl_reducible(HDL(g),HDL(r)) ) + if ( ndl_reducible(dg,HDL(r)) ) return i; } else for ( i = 0; i < len; i++ ) { r = ps[i]; - if ( ndl_reducible(HDL(g),HDL(r)) ) + if ( ndl_reducible(dg,HDL(r)) ) return i; } return -1; } +ND nd_merge(ND p1,ND p2) +{ + int n,c; + int t,can,td1,td2; + ND r; + NM m1,m2,mr0,mr,s; + + 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: + s = m1; m1 = NEXT(m1); + can++; NEXTNM2(mr0,mr,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; + } +} + ND nd_add(int mod,ND p1,ND p2) { int n,c; @@ -1248,7 +1358,7 @@ int nd_nf(int mod,ND g,NDV *ps,int full,ND *rp) n = NV(g); mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(UINT)); for ( d = 0; g; ) { - index = nd_find_reducer(g); + index = ndl_find_reducer(HDL(g)); if ( index >= 0 ) { h = nd_psh[index]; ndl_sub(HDL(g),DL(h),DL(mul)); @@ -1326,7 +1436,7 @@ int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp return 1; } g = bucket->body[hindex]; - index = nd_find_reducer(g); + index = ndl_find_reducer(HDL(g)); if ( index >= 0 ) { h = nd_psh[index]; ndl_sub(HDL(g),DL(h),DL(mul)); @@ -1421,7 +1531,7 @@ int nd_nf_direct(int mod,ND g,BaseSet base,int full,ND n = NV(g); mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(UINT)); for ( d = 0; g; ) { - index = nd_find_reducer_direct(g,ps,len); + index = ndl_find_reducer_direct(HDL(g),ps,len); if ( index >= 0 ) { p = ps[index]; ndl_sub(HDL(g),HDL(p),DL(mul)); @@ -1533,6 +1643,24 @@ void free_pbucket(PGeoBucket b) { GC_free(b); } +void add_pbucket_symbolic(PGeoBucket g,ND d) +{ + int l,i,k,m; + + if ( !d ) + return; + l = LEN(d); + for ( k = 0, m = 1; l > m; k++, m <<= 1 ); + /* 2^(k-1) < l <= 2^k (=m) */ + d = nd_merge(g->body[k],d); + for ( ; d && LEN(d) > m; k++, m <<= 1 ) { + g->body[k] = 0; + d = nd_merge(g->body[k+1],d); + } + g->body[k] = d; + g->m = MAX(g->m,k); +} + void add_pbucket(int mod,PGeoBucket g,ND d) { int l,i,k,m; @@ -1559,6 +1687,38 @@ void mulq_pbucket(PGeoBucket g,Q c) nd_mul_c_q(g->body[k],c); } +NM remove_head_pbucket_symbolic(PGeoBucket g) +{ + int j,i,k,c; + NM head; + + k = g->m; + j = -1; + for ( i = 0; i <= k; i++ ) { + if ( !g->body[i] ) continue; + if ( j < 0 ) j = i; + else { + c = DL_COMPARE(HDL(g->body[i]),HDL(g->body[j])); + if ( c > 0 ) + j = i; + else if ( c == 0 ) + g->body[i] = nd_remove_head(g->body[i]); + } + } + if ( j < 0 ) return 0; + else { + head = BDY(g->body[j]); + if ( !NEXT(head) ) { + FREEND(g->body[j]); + g->body[j] = 0; + } else { + BDY(g->body[j]) = NEXT(head); + LEN(g->body[j])--; + } + return head; + } +} + int head_pbucket(int mod,PGeoBucket g) { int j,i,c,k,nv,sum; @@ -1664,6 +1824,7 @@ NODE nd_gb(int m,int checkonly) ND_pairs d; ND_pairs l; ND h,nf; + NDV nfv; g = 0; d = 0; for ( i = 0; i < nd_psn; i++ ) { @@ -1696,7 +1857,9 @@ again: } else if ( nf ) { if ( checkonly ) return 0; printf("+"); fflush(stdout); - nh = nd_newps(m,nf,0); + nd_removecont(m,nf); + nfv = ndtondv(m,nf); nd_free(nf); + nh = ndv_newps(nfv,0); d = update_pairs(d,g,nh); g = update_base(g,nh); FREENDP(l); @@ -1716,6 +1879,7 @@ NODE nd_gb_trace(int m) ND_pairs d; ND_pairs l; ND h,nf,nfq; + NDV nfv,nfqv; g = 0; d = 0; for ( i = 0; i < nd_psn; i++ ) { @@ -1750,10 +1914,13 @@ again: nd_sp(0,1,l,&h); nd_nf(0,h,nd_ps_trace,!Top,&nfq); if ( nfq ) { - printf("+"); fflush(stdout); - nh = nd_newps(m,nf,nfq); /* failure; m|HC(nfq) */ - if ( nh < 0 ) return 0; + if ( !rem(NM(HCQ(nfq)),m) ) return 0; + + printf("+"); fflush(stdout); + nd_removecont(m,nf); nfv = ndtondv(m,nf); nd_free(nf); + nd_removecont(0,nfq); nfqv = ndtondv(0,nfq); nd_free(nfq); + nh = ndv_newps(nfv,nfqv); d = update_pairs(d,g,nh); g = update_base(g,nh); } else { @@ -2074,17 +2241,36 @@ ND_pairs nd_minp( ND_pairs d, ND_pairs *prest ) return m; } -int nd_newps(int mod,ND a,ND aq) +ND_pairs nd_minsugarp( ND_pairs d, ND_pairs *prest ) { + int msugar; + 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 ( t = d; t; t = NEXT(t) ) + if ( SG(t) == msugar ) { + if ( dm0 ) NEXT(dm) = t; + else dm0 = t; + dm = t; + } else { + if ( dr0 ) NEXT(dr) = t; + else dr0 = t; + dr = t; + } + NEXT(dm) = 0; + if ( dr0 ) NEXT(dr) = 0; + *prest = dr0; + return dm0; +} + +int ndv_newps(NDV a,NDV aq) +{ int len; RHist r; NDV b; - if ( aq ) { - /* trace lifting */ - if ( !rem(NM(HCQ(aq)),mod) ) - return -1; - } if ( nd_psn == nd_pslen ) { nd_pslen *= 2; nd_ps = (NDV *)REALLOC((char *)nd_ps,nd_pslen*sizeof(NDV)); @@ -2094,16 +2280,14 @@ int nd_newps(int mod,ND a,ND aq) REALLOC((char *)nd_bound,nd_pslen*sizeof(UINT *)); } NEWRHist(r); nd_psh[nd_psn] = r; - nd_removecont(mod,a); nd_ps[nd_psn] = ndtondv(mod,a); + nd_ps[nd_psn] = a; if ( aq ) { - nd_removecont(0,aq); nd_ps_trace[nd_psn] = ndtondv(0,aq); - nd_bound[nd_psn] = ndv_compute_bound(nd_ps_trace[nd_psn]); + nd_ps_trace[nd_psn] = aq; + nd_bound[nd_psn] = ndv_compute_bound(aq); SG(r) = SG(aq); ndl_copy(HDL(aq),DL(r)); - nd_free(a); nd_free(aq); } else { - nd_bound[nd_psn] = ndv_compute_bound(nd_ps[nd_psn]); + nd_bound[nd_psn] = ndv_compute_bound(a); SG(r) = SG(a); ndl_copy(HDL(a),DL(r)); - nd_free(a); } return nd_psn++; } @@ -2146,7 +2330,7 @@ void ndv_setup(int mod,int trace,NODE f) } } -void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp) +void nd_gr(LIST f,LIST v,int m,int f4,struct order_spec *ord,LIST *rp) { VL tv,fv,vv,vc; NODE fd,fd0,r,r0,t,x,s,xx; @@ -2168,7 +2352,7 @@ void nd_gr(LIST f,LIST v,int m,struct order_spec *ord, } if ( fd0 ) NEXT(fd) = 0; ndv_setup(m,0,fd0); - x = nd_gb(m,0); + x = f4?nd_f4(m):nd_gb(m,0); x = ndv_reducebase(x); x = ndv_reduceall(m,x); for ( r0 = 0, t = x; t; t = NEXT(t) ) { @@ -2651,11 +2835,14 @@ void nd_setup_parameters(int nvar,int max) { /* if max == 0, don't touch nd_bpe */ if ( max > 0 ) { + if ( max < 2 ) nd_bpe = 1; if ( max < 4 ) nd_bpe = 2; else if ( max < 8 ) nd_bpe = 3; else if ( max < 16 ) nd_bpe = 4; + else if ( max < 32 ) nd_bpe = 5; else if ( max < 64 ) nd_bpe = 6; else if ( max < 256 ) nd_bpe = 8; + else if ( max < 1024 ) nd_bpe = 10; else if ( max < 65536 ) nd_bpe = 16; else nd_bpe = 32; } @@ -2698,10 +2885,13 @@ ND_pairs nd_reconstruct(int mod,int trace,ND_pairs d) obpe = nd_bpe; oadv = nmv_adv; oepos = nd_epos; - if ( obpe < 3 ) nd_bpe = 3; + if ( obpe < 2 ) nd_bpe = 2; + else if ( obpe < 3 ) nd_bpe = 3; else if ( obpe < 4 ) nd_bpe = 4; + else if ( obpe < 5 ) nd_bpe = 5; else if ( obpe < 6 ) nd_bpe = 6; else if ( obpe < 8 ) nd_bpe = 8; + else if ( obpe < 10 ) nd_bpe = 10; else if ( obpe < 16 ) nd_bpe = 16; else if ( obpe < 32 ) nd_bpe = 32; else error("nd_reconstruct : exponent too large"); @@ -2767,10 +2957,13 @@ void nd_reconstruct_direct(int mod,NDV *ps,int len) obpe = nd_bpe; oadv = nmv_adv; oepos = nd_epos; - if ( obpe < 3 ) nd_bpe = 3; + if ( obpe < 2 ) nd_bpe = 2; + else if ( obpe < 3 ) nd_bpe = 3; else if ( obpe < 4 ) nd_bpe = 4; + else if ( obpe < 5 ) nd_bpe = 5; else if ( obpe < 6 ) nd_bpe = 6; else if ( obpe < 8 ) nd_bpe = 8; + else if ( obpe < 10 ) nd_bpe = 10; else if ( obpe < 16 ) nd_bpe = 16; else if ( obpe < 32 ) nd_bpe = 32; else error("nd_reconstruct_direct : exponent too large"); @@ -3028,6 +3221,35 @@ void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *ta FREENM(m); } +ND ndv_mul_nm_symbolic(NM m0,NDV p) +{ + NM mr,mr0; + NMV m; + UINT *d,*dt,*dm; + int c,n,td,i,c1,c2,len; + Q q; + ND r; + + if ( !p ) return 0; + else { + n = NV(p); m = BDY(p); + d = DL(m0); + len = LEN(p); + mr0 = 0; + td = TD(d); + c = CM(m0); + for ( i = 0; i < len; i++, NMV_ADV(m) ) { + NEXTNM(mr0,mr); + CM(mr) = 1; + ndl_add(DL(m),d,DL(mr)); + } + NEXT(mr) = 0; + MKND(NV(p),mr0,len,r); + SG(r) = SG(p) + TD(d); + return r; + } +} + ND ndv_mul_nm(int mod,NM m0,NDV p) { NM mr,mr0; @@ -3130,6 +3352,23 @@ NDV ndv_dup(int mod,NDV p) return d; } +ND nd_dup(ND p) +{ + ND d; + NM t,m,m0; + + if ( !p ) return 0; + for ( m0 = 0, t = BDY(p); t; t = NEXT(t) ) { + NEXTNM(m0,m); + ndl_copy(DL(t),DL(m)); + CQ(m) = CQ(t); + } + if ( m0 ) NEXT(m) = 0; + MKND(NV(p),m0,LEN(p),d); + SG(d) = SG(p); + return d; +} + /* XXX if p->len == 0 then it represents 0 */ void ndv_mod(int mod,NDV p) @@ -3507,4 +3746,267 @@ void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec break; } *rp = ndvtop(m,CO,vv,ndtondv(m,nf)); +} + +int nd_to_vect(int mod,UINT *s0,int n,ND d,UINT *r) +{ + NM m; + UINT *t,*s; + int i; + + for ( i = 0; i < n; i++ ) 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++ ); + r[i] = CM(m); + } + for ( i = 0; !r[i]; i++ ); + return i; +} + +int nm_ind_pair_to_vect(int mod,UINT *s0,int n,NM_ind_pair pair,UINT *r) +{ + NM m; + NMV mr; + UINT *d,*t,*s; + NDV p; + int i,j,len; + + m = pair->mul; + d = DL(m); + p = nd_ps[pair->index]; + t = (UINT *)ALLOCA(nd_wpd*sizeof(UINT)); + for ( i = 0; i < n; i++ ) r[i] = 0; + len = LEN(p); + for ( i = j = 0, s = s0, mr = BDY(p); j < len; j++, NMV_ADV(mr) ) { + ndl_add(d,DL(mr),t); + for ( ; !ndl_equal(t,s); s += nd_wpd, i++ ); + r[i] = CM(mr); + } + for ( i = 0; !r[i]; i++ ); + return i; +} + +void nm_ind_pair_to_vect_compress(int mod,UINT *s0,int n,NM_ind_pair pair,int *ind) +{ + NM m; + NMV mr; + UINT *d,*t,*s; + NDV p; + int i,j,len; + + m = pair->mul; + d = DL(m); + p = nd_ps[pair->index]; + len = LEN(p); + t = (UINT *)ALLOCA(nd_wpd*sizeof(UINT)); + for ( i = j = 0, s = s0, mr = BDY(p); j < len; j++, NMV_ADV(mr) ) { + ndl_add(d,DL(mr),t); + for ( ; !ndl_equal(t,s); s += nd_wpd, i++ ); + ind[j] = i; + } +} + + +void ndv_reduce_vect(int m,UINT *svect,int **imat,NODE rp0) +{ + int i,j,k,len,pos; + UINT c,c1,c2,c3; + UINT *ivect; + NDV redv; + NMV mr0,mr; + NODE rp; + + for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) { + ivect = imat[i]; + k = ivect[0]; + if ( c = svect[k] ) { + c = m-c; + redv = nd_ps[((NM_ind_pair)BDY(rp))->index]; + len = LEN(redv); + mr0 = BDY(redv); + for ( j = 0, mr = mr0; j < len; j++, NMV_ADV(mr) ) { + pos = ivect[j]; + c1 = CM(mr); + c2 = svect[pos]; DMAR(c1,c,c2,m,c3); + svect[pos] = c3; + } + } + } +} + +NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhead,UINT *s0vect) +{ + int j,k,len; + UINT *p; + UINT c; + NDV r; + NMV mr0,mr; + + for ( j = 0, len = 0; j < spcol; j++ ) if ( vect[j] ) len++; + if ( !len ) return 0; + else { + mr0 = (NMV)MALLOC_ATOMIC(nmv_adv*len); + mr = mr0; + p = s0vect; + for ( j = k = 0; j < col; j++, p += nd_wpd ) + if ( !rhead[j] ) { + if ( c = vect[k++] ) { + ndl_copy(p,DL(mr)); CM(mr) = c; NMV_ADV(mr); + } + } + MKNDV(nd_nvar,mr0,len,r); + return r; + } +} + +NODE nd_sp_f4(int m,ND_pairs l,PGeoBucket bucket) +{ + ND_pairs t; + NODE sp0,sp; + int stat; + ND spol; + + sp0 = 0; + for ( t = l; t; t = NEXT(t) ) { + stat = nd_sp(m,0,t,&spol); + if ( !stat ) return 0; + if ( spol ) { + NEXTNODE(sp0,sp); BDY(sp) = (pointer)nd_dup(spol); + add_pbucket_symbolic(bucket,spol); + } + } + return sp0; +} + +int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vect,NODE *r) +{ + NODE rp0,rp; + NM mul,head,s0,s; + int index,col,i; + RHist h; + UINT *s0v,*p; + NM_ind_pair pair; + ND red; + + s0 = 0; rp0 = 0; col = 0; + while ( 1 ) { + head = remove_head_pbucket_symbolic(bucket); + if ( !head ) break; + if ( !s0 ) s0 = head; + else NEXT(s) = head; + s = head; + index = ndl_find_reducer(DL(head)); + if ( index >= 0 ) { + h = nd_psh[index]; + NEWNM(mul); + ndl_sub(DL(head),DL(h),DL(mul)); + if ( ndl_check_bound2(index,DL(mul)) ) return 0; + MKNM_ind_pair(pair,mul,index); + red = ndv_mul_nm_symbolic(mul,nd_ps[index]); + add_pbucket_symbolic(bucket,nd_remove_head(red)); + NEXTNODE(rp0,rp); BDY(rp) = (pointer)pair; + } + col++; + } + NEXT(rp) = 0; NEXT(s) = 0; + s0v = (UINT *)MALLOC_ATOMIC(col*nd_wpd*sizeof(UINT)); + for ( i = 0, p = s0v, s = s0; i < col; + i++, p += nd_wpd, s = NEXT(s) ) ndl_copy(DL(s),p); + *s0vect = s0v; + *r = rp0; + return col; +} + +NODE nd_f4(int m) +{ + int i,nh,stat,index; + NODE r,g; + ND_pairs d,l,t; + ND spol,red; + NDV nf,redv; + NM s0,s; + NODE sp0,sp,rp0,rp; + int nsp,nred,col,rank,len,k,j,a; + UINT c; + UINT **spmat; + UINT *s0vect,*svect,*p,*v; + int *colstat; + int **imat; + int *rhead; + int spcol,sprow; + int sugar; + PGeoBucket bucket; + + if ( !m ) + error("nd_f4 : not implemented"); + + g = 0; d = 0; + for ( i = 0; i < nd_psn; i++ ) { + d = update_pairs(d,g,i); + g = update_base(g,i); + } + while ( d ) { + l = nd_minsugarp(d,&d); + sugar = SG(l); + bucket = create_pbucket(); + if ( !(sp0 = nd_sp_f4(m,l,bucket)) + || !(col = nd_symbolic_preproc(bucket,&s0vect,&rp0)) ) { + for ( t = l; NEXT(t); t = NEXT(t) ); + NEXT(t) = d; d = l; + d = nd_reconstruct(m,0,d); + continue; + } + + nsp = length(sp0); nred = length(rp0); spcol = col-nred; + imat = (int **)MALLOC(nred*sizeof(int *)); + rhead = (int *)MALLOC_ATOMIC(col*sizeof(int)); + for ( i = 0; i < col; i++ ) rhead[i] = 0; + + /* construction of index arrays */ + for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) { + redv = nd_ps[((NM_ind_pair)BDY(rp))->index]; + imat[i] = (int *)MALLOC_ATOMIC(LEN(redv)*sizeof(int)); + nm_ind_pair_to_vect_compress(m,s0vect,col, + (NM_ind_pair)BDY(rp),imat[i]); + rhead[imat[i][0]] = 1; + } + + /* elimination (1st step) */ + spmat = (UINT **)MALLOC(nsp*sizeof(UINT *)); + svect = (UINT *)MALLOC_ATOMIC(col*sizeof(UINT)); + for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) { + nd_to_vect(m,s0vect,col,BDY(sp),svect); + ndv_reduce_vect(m,svect,imat,rp0); + for ( i = 0; i < col; i++ ) if ( svect[i] ) break; + if ( i < col ) { + spmat[sprow] = v = (UINT *)MALLOC_ATOMIC(spcol*sizeof(UINT)); + for ( j = k = 0; j < col; j++ ) + if ( !rhead[j] ) v[k++] = svect[j]; + sprow++; + } + } + /* free index arrays */ + for ( i = 0; i < nred; i++ ) GC_free(imat[i]); + + /* elimination (2nd step) */ + colstat = (int *)ALLOCA(spcol*sizeof(int)); + rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat); + + fprintf(asir_out,"sugar=%d,rank=%d\n",sugar,rank); fflush(asir_out); + + /* adding new bases */ + for ( i = sprow-1; i >= rank; i-- ) GC_free(spmat[i]); + for ( ; i >= 0; i-- ) { + nf = vect_to_ndv(spmat[i],spcol,col,rhead,s0vect); + SG(nf) = sugar; + ndv_removecont(m,nf); + nh = ndv_newps(nf,0); + d = update_pairs(d,g,nh); + g = update_base(g,nh); + GC_free(spmat[i]); + } + } + for ( r = g; r; r = NEXT(r) ) BDY(r) = (pointer)nd_ps[(int)BDY(r)]; + return g; }