=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v retrieving revision 1.25 retrieving revision 1.28 diff -u -p -r1.25 -r1.28 --- OpenXM_contrib2/asir2000/engine/nd.c 2003/08/07 08:46:50 1.25 +++ OpenXM_contrib2/asir2000/engine/nd.c 2003/08/11 06:58:01 1.28 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.24 2003/08/05 08:06:21 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.27 2003/08/10 01:31:24 noro Exp $ */ #include "ca.h" #include "inline.h" @@ -11,6 +11,8 @@ #define INLINE #endif +#define USE_GEOBUCKET 1 + #define REDTAB_LEN 32003 typedef struct oPGeoBucket { @@ -122,6 +124,7 @@ 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);} +int nd_check_candidate(NODE input,NODE cand); void nd_removecont(int mod,ND p); void nd_removecont2(ND p1,ND p2); void ndv_removecont(int mod,NDV p); @@ -130,13 +133,16 @@ void ndv_mul_c_q(NDV p,Q mul); void nd_mul_c_q(ND p,Q mul); ND normalize_pbucket(int mod,PGeoBucket g); int head_pbucket(int mod,PGeoBucket g); +int head_pbucket_q(PGeoBucket g); void add_pbucket(int mod,PGeoBucket g,ND d,int l); void free_pbucket(PGeoBucket b); +void mulq_pbucket(PGeoBucket g,Q c); PGeoBucket create_pbucket(); ND nd_remove_head(ND p); void GC_gcollect(); NODE append_one(NODE,int); +NODE nd_reducebase(NODE x); void removecont_array(Q *c,int n); ND_pairs crit_B( ND_pairs d, int s ); @@ -154,7 +160,7 @@ ND_pairs crit_F( ND_pairs d1 ); ND_pairs crit_M( ND_pairs d1 ); ND_pairs nd_newpairs( NODE g, int t ); ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t); -NODE nd_gb(int m); +NODE nd_gb(int m,int checkonly); NODE nd_gb_trace(int m); void nd_free_private_storage(); void _NM_alloc(); @@ -172,6 +178,8 @@ int nd_nf(int mod,ND g,int full,ND *nf); ND nd_reduce(ND p1,ND p2); ND nd_reduce_special(ND p1,ND p2); NODE nd_reduceall(int m,NODE f); +int nd_gbcheck(int m,NODE f); +int nd_membercheck(int m,NODE f); void nd_free(ND p); void ndv_free(NDV p); void ndl_print(unsigned int *dl); @@ -915,7 +923,7 @@ ND nd_add_q(ND p1,ND p2) } } -#if 1 +#if !USE_GEOBUCKET /* ret=1 : success, ret=0 : overflow */ int nd_nf(int mod,ND g,int full,ND *rp) { @@ -995,7 +1003,83 @@ afo: *rp = d; return 1; } + +int nd_nf_direct(int mod,ND g,NDV *ps,int len,int full,ND *rp) +{ + ND d; + NM m,mrd,tail; + NM mul; + int n,sugar,psugar,sugar0,stat,index; + int c,c1,c2; + RHist h; + NDV p,red; + Q cg,cred,gcd; + double hmag; + + if ( !g ) { + *rp = 0; + return 1; + } +#if 0 + if ( !mod ) + hmag = ((double)p_mag((P)HCQ(g)))*nd_scale; #else + /* XXX */ + hmag = 0; +#endif + + sugar0 = sugar = SG(g); + n = NV(g); + mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int)); + for ( d = 0; g; ) { + index = nd_find_reducer_direct(g,ps,len); + if ( index >= 0 ) { + p = ps[index]; + ndl_sub(HDL(g),HDL(p),DL(mul)); + TD(mul) = HTD(g)-HTD(p); + if ( ndl_check_bound2_direct(HDL(p),DL(mul)) ) { + nd_free(g); nd_free(d); + return 0; + } + if ( mod ) { + c1 = invm(HCM(p),mod); c2 = mod-HCM(g); + DMAR(c1,c2,0,mod,c); CM(mul) = c; + } else { + igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred); + chsgnq(cg,&CQ(mul)); + nd_mul_c_q(d,cred); nd_mul_c_q(g,cred); + } + g = nd_add(mod,g,ndv_mul_nm(mod,p,mul)); + sugar = MAX(sugar,SG(p)+TD(mul)); + if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) { + nd_removecont2(d,g); + hmag = ((double)p_mag((P)HCQ(g)))*nd_scale; + } + } else if ( !full ) { + *rp = g; + return 1; + } else { + m = BDY(g); + if ( NEXT(m) ) { + BDY(g) = NEXT(m); NEXT(m) = 0; + } else { + FREEND(g); g = 0; + } + if ( d ) { + NEXT(tail)=m; + tail=m; + } else { + MKND(n,m,d); + tail = BDY(d); + } + } + } + if ( d ) + SG(d) = sugar; + *rp = d; + return 1; +} +#else int nd_nf(int mod,ND g,int full,ND *rp) { int hindex,index; @@ -1006,8 +1090,9 @@ int nd_nf(int mod,ND g,int full,ND *rp) int sugar,psugar,n,h_reducible; PGeoBucket bucket; int c,c1,c2; - Q cg,cred,gcd; + Q cg,cred,gcd,zzz; RHist h; + double hmag,gmag; if ( !g ) { *rp = 0; @@ -1015,12 +1100,14 @@ int nd_nf(int mod,ND g,int full,ND *rp) } sugar = SG(g); n = NV(g); + if ( !mod ) + hmag = ((double)p_mag((P)HCQ(g)))*nd_scale; bucket = create_pbucket(); add_pbucket(mod,bucket,g,nd_length(g)); d = 0; mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int)); while ( 1 ) { - hindex = head_pbucket(mod,bucket); + hindex = mod?head_pbucket(mod,bucket):head_pbucket_q(bucket); if ( hindex < 0 ) { if ( d ) SG(d) = sugar; @@ -1034,8 +1121,8 @@ int nd_nf(int mod,ND g,int full,ND *rp) ndl_sub(HDL(g),DL(h),DL(mul)); TD(mul) = HTD(g)-TD(h); if ( ndl_check_bound2(index,DL(mul)) ) { - free_pbucket(bucket); nd_free(d); + free_pbucket(bucket); *rp = 0; return 0; } @@ -1047,13 +1134,28 @@ int nd_nf(int mod,ND g,int full,ND *rp) p = nd_psq[index]; igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred); chsgnq(cg,&CQ(mul)); - nd_mul_c_q(d,cred); nd_mul_c_q(g,cred); + nd_mul_c_q(d,cred); + mulq_pbucket(bucket,cred); + g = bucket->body[hindex]; + gmag = (double)p_mag((P)HCQ(g)); } red = ndv_mul_nm(mod,p,mul); bucket->body[hindex] = nd_remove_head(g); red = nd_remove_head(red); - add_pbucket(mod,bucket,red,LEN(p)); + add_pbucket(mod,bucket,red,LEN(p)-1); sugar = MAX(sugar,SG(p)+TD(mul)); + if ( !mod && hmag && (gmag > hmag) ) { + g = normalize_pbucket(mod,bucket); + if ( !g ) { + if ( d ) + SG(d) = sugar; + *rp = d; + return 1; + } + nd_removecont2(d,g); + hmag = ((double)p_mag((P)HCQ(g)))*nd_scale; + add_pbucket(mod,bucket,g,nd_length(g)-1); + } } else if ( !full ) { g = normalize_pbucket(mod,bucket); if ( g ) @@ -1078,24 +1180,27 @@ int nd_nf(int mod,ND g,int full,ND *rp) } } } -#endif int nd_nf_direct(int mod,ND g,NDV *ps,int len,int full,ND *rp) { - ND d; - NM m,mrd,tail; - NM mul; - int n,sugar,psugar,sugar0,stat,index; + int hindex,index; + NDV p; + ND u,d,red; + NODE l; + NM mul,m,mrd; + int sugar,psugar,n,h_reducible; + PGeoBucket bucket; int c,c1,c2; + Q cg,cred,gcd,zzz; RHist h; - NDV p,red; - Q cg,cred,gcd; - double hmag; + double hmag,gmag; if ( !g ) { *rp = 0; return 1; } + sugar = SG(g); + n = NV(g); #if 0 if ( !mod ) hmag = ((double)p_mag((P)HCQ(g)))*nd_scale; @@ -1103,35 +1208,62 @@ int nd_nf_direct(int mod,ND g,NDV *ps,int len,int full /* XXX */ hmag = 0; #endif - - sugar0 = sugar = SG(g); - n = NV(g); + bucket = create_pbucket(); + add_pbucket(mod,bucket,g,nd_length(g)); + d = 0; mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int)); - for ( d = 0; g; ) { + while ( 1 ) { + hindex = mod?head_pbucket(mod,bucket):head_pbucket_q(bucket); + if ( hindex < 0 ) { + if ( d ) + SG(d) = sugar; + *rp = d; + return 1; + } + g = bucket->body[hindex]; index = nd_find_reducer_direct(g,ps,len); if ( index >= 0 ) { p = ps[index]; ndl_sub(HDL(g),HDL(p),DL(mul)); TD(mul) = HTD(g)-HTD(p); if ( ndl_check_bound2_direct(HDL(p),DL(mul)) ) { - nd_free(g); nd_free(d); + nd_free(d); + free_pbucket(bucket); + *rp = 0; return 0; } - if ( mod ) { + if ( mod ) { c1 = invm(HCM(p),mod); c2 = mod-HCM(g); DMAR(c1,c2,0,mod,c); CM(mul) = c; } else { igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred); chsgnq(cg,&CQ(mul)); - nd_mul_c_q(d,cred); nd_mul_c_q(g,cred); + nd_mul_c_q(d,cred); + mulq_pbucket(bucket,cred); + g = bucket->body[hindex]; + gmag = (double)p_mag((P)HCQ(g)); } - g = nd_add(mod,g,ndv_mul_nm(mod,p,mul)); + red = ndv_mul_nm(mod,p,mul); + bucket->body[hindex] = nd_remove_head(g); + red = nd_remove_head(red); + add_pbucket(mod,bucket,red,LEN(p)-1); sugar = MAX(sugar,SG(p)+TD(mul)); - if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) { + if ( !mod && hmag && (gmag > hmag) ) { + g = normalize_pbucket(mod,bucket); + if ( !g ) { + if ( d ) + SG(d) = sugar; + *rp = d; + return 1; + } nd_removecont2(d,g); hmag = ((double)p_mag((P)HCQ(g)))*nd_scale; + add_pbucket(mod,bucket,g,nd_length(g)-1); } } else if ( !full ) { + g = normalize_pbucket(mod,bucket); + if ( g ) + SG(g) = sugar; *rp = g; return 1; } else { @@ -1141,18 +1273,42 @@ int nd_nf_direct(int mod,ND g,NDV *ps,int len,int full } else { FREEND(g); g = 0; } + bucket->body[hindex] = g; + NEXT(m) = 0; if ( d ) { - NEXT(tail)=m; - tail=m; + for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) ); + NEXT(mrd) = m; } else { MKND(n,m,d); - tail = BDY(d); } } } - if ( d ) - SG(d) = sugar; - *rp = d; +} +#endif + +/* input : list of DP, cand : list of DP */ + +int nd_check_candidate(NODE input,NODE cand) +{ + int n,i,stat; + ND nf,d; + NODE t; + + nd_setup(0,cand); + + /* membercheck : list is a subset of Id(cand) ? */ + for ( t = input; t; t = NEXT(t) ) { + d = dptond(0,(DP)BDY(t)); + stat = nd_nf_direct(0,d,nd_psq,n,0,&nf); + if ( !stat ) + nd_reconstruct_direct(0,nd_psq,n); + else if ( nf ) + return 0; + } + /* gbcheck : cand is a GB of Id(cand) ? */ + if ( !nd_gb(0,1) ) + return 0; + /* XXX */ return 1; } @@ -1182,7 +1338,7 @@ PGeoBucket create_pbucket() void free_pbucket(PGeoBucket b) { int i; - for ( i = 0; i < b->m; i++ ) + for ( i = 0; i <= b->m; i++ ) if ( b->body[i] ) { nd_free(b->body[i]); b->body[i] = 0; @@ -1205,6 +1361,14 @@ void add_pbucket(int mod,PGeoBucket g,ND d,int l) g->m = MAX(g->m,k); } +void mulq_pbucket(PGeoBucket g,Q c) +{ + int k; + + for ( k = 0; k <= g->m; k++ ) + nd_mul_c_q(g->body[k],c); +} + /* XXX not completed */ int head_pbucket(int mod,PGeoBucket g) { @@ -1259,18 +1423,76 @@ int head_pbucket(int mod,PGeoBucket g) } } +int head_pbucket_q(PGeoBucket g) +{ + int j,i,c,k,nv; + Q sum,t; + unsigned int *di,*dj; + 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]; + dj = HDL(gj); + sum = HCQ(gj); + } else { + di = HDL(gi); + nv = NV(gi); + if ( HTD(gi) > HTD(gj) ) + c = 1; + else if ( HTD(gi) < HTD(gj) ) + c = -1; + else + c = ndl_compare(di,dj); + if ( c > 0 ) { + if ( sum ) + HCQ(gj) = sum; + else + g->body[j] = nd_remove_head(gj); + j = i; + gj = g->body[j]; + dj = HDL(gj); + sum = HCQ(gj); + } else if ( c == 0 ) { + addq(sum,HCQ(gi),&t); + sum = t; + g->body[i] = nd_remove_head(gi); + } + } + } + if ( j < 0 ) + return -1; + else if ( sum ) { + HCQ(gj) = sum; + return j; + } else + g->body[j] = nd_remove_head(gj); + } +} + ND normalize_pbucket(int mod,PGeoBucket g) { int i; ND r,t; r = 0; - for ( i = 0; i <= g->m; i++ ) + for ( i = 0; i <= g->m; i++ ) { r = nd_add(mod,r,g->body[i]); + g->body[i] = 0; + } + g->m = -1; return r; } -NODE nd_gb(int m) +/* return value = 0 => input is not a GB */ + +NODE nd_gb(int m,int checkonly) { int i,nh,sugar,stat; NODE r,g,t; @@ -1303,6 +1525,7 @@ again: d = nd_reconstruct(m,0,d); goto again; } else if ( nf ) { + if ( checkonly ) return 0; printf("+"); fflush(stdout); nh = nd_newps(m,nf); d = update_pairs(d,g,nh); @@ -1361,6 +1584,8 @@ again: if ( nfq ) { printf("+"); fflush(stdout); nh = nd_newps_trace(m,nf,nfq); + /* failure; m|HC(nfq) */ + if ( nf < 0 ) return 0; d = update_pairs(d,g,nh); g = update_base(g,nh); } else { @@ -1761,6 +1986,7 @@ int nd_newps_trace(int mod,ND nf,ND nfq) nd_bound = (unsigned int **) REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *)); } + if ( !rem(NM(HCQ(nfq)),mod) ) return -1; nd_removecont(mod,nf); nd_ps[nd_psn] = ndtondv(mod,nf); @@ -1912,7 +2138,7 @@ void nd_gr(LIST f,LIST v,int m,struct order_spec *ord, } if ( fd0 ) NEXT(fd) = 0; nd_setup(m,fd0); - x = nd_gb(m); + x = nd_gb(m,0); fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n", nd_found,nd_notfirst,nd_create); x = nd_reduceall(m,x); @@ -1934,10 +2160,9 @@ void nd_gr_trace(LIST f,LIST v,int m,int homo,struct o { struct order_spec ord1; VL fv,vv,vc; - NODE fd,fd0,r,r0,t,x,s,xx,t0; + NODE fd,fd0,in0,in,r,r0,t,s,cand; DP a,b,c,h; - NDV *w; - int len,i,j; + P p; get_vars((Obj)f,&fv); pltovl(v,&vv); nd_nvar = length(vv); @@ -1951,38 +2176,17 @@ void nd_gr_trace(LIST f,LIST v,int m,int homo,struct o case 1: is_rlex = 0; break; default: error("nd_gr : unsupported order"); } - for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) { + for ( fd0 = 0, in0 = 0, t = BDY(f); t; t = NEXT(t) ) { ptod(CO,vv,(P)BDY(t),&c); if ( c ) { dp_homo(c,&h); NEXTNODE(fd0,fd); BDY(fd) = (pointer)h; + NEXTNODE(in0,in); BDY(in) = (pointer)c; } } if ( fd0 ) NEXT(fd) = 0; + if ( in0 ) NEXT(in) = 0; initd(&ord1); nd_nvar++; - nd_setup_trace(m,fd0); - x = nd_gb_trace(m); - - /* dehomogenization */ - for ( t = x; t; t = NEXT(t) ) - ndv_dehomogenize((NDV)BDY(t)); - nd_nvar--; - nd_setup_parameters(); - initd(ord); - len = length(x); - w = (NDV *)ALLOCA(len*sizeof(NDV)); - for ( i = 0, t = x; i < len; i++, t = NEXT(t) ) w[i] = BDY(t); - for ( i = 0; i < len; i++ ) { - for ( j = 0; j < i; j++ ) { - if ( w[i] && w[j] ) - if ( ndl_reducible(HDL(w[i]),HDL(w[j])) ) w[i] = 0; - else if ( ndl_reducible(HDL(w[j]),HDL(w[i])) ) w[j] = 0; - } - } - for ( i = len-1, t0 = 0; i >= 0; i-- ) { - if ( w[i] ) { NEXTNODE(t0,t); BDY(t) = (pointer)w[i]; } - } - NEXT(t) = 0; x = t0; } else { switch ( ord->ord.simple ) { case 0: is_rlex = 1; break; @@ -1996,20 +2200,38 @@ void nd_gr_trace(LIST f,LIST v,int m,int homo,struct o } } if ( fd0 ) NEXT(fd) = 0; - /* setup over GF(m) */ + in0 = fd0; + } + do { nd_setup_trace(m,fd0); - x = nd_gb_trace(m); + cand = nd_gb_trace(m); + if ( !cand ) continue; + if ( homo ) { + /* dehomogenization */ + for ( t = cand; t; t = NEXT(t) ) + ndv_dehomogenize((NDV)BDY(t)); + nd_nvar--; + nd_setup_parameters(); + initd(ord); + cand = nd_reducebase(cand); + } + fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n", + nd_found,nd_notfirst,nd_create); + cand = nd_reduceall(0,cand); + initd(ord); + for ( r0 = 0; cand; cand = NEXT(cand) ) { + NEXTNODE(r0,r); + BDY(r) = (pointer)ndvtodp(0,(NDV)BDY(cand)); + } + if ( r0 ) NEXT(r) = 0; + cand = r0; + } while ( !nd_check_candidate(in0,cand) ); + /* dp->p */ + for ( r = cand; r; r = NEXT(r) ) { + dtop(CO,vv,BDY(r),&p); + BDY(r) = (pointer)p; } - fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n", - nd_found,nd_notfirst,nd_create); - x = nd_reduceall(0,x); - for ( r0 = 0; x; x = NEXT(x) ) { - NEXTNODE(r0,r); - a = ndvtodp(0,(NDV)BDY(x)); - dtop(CO,vv,a,(P *)&BDY(r)); - } - if ( r0 ) NEXT(r) = 0; - MKLIST(*rp,r0); + MKLIST(*rp,cand); } void dltondl(int n,DL dl,unsigned int *r) @@ -2879,4 +3101,27 @@ int nd_equal(ND a,ND b) } if ( !ma && !mb ) return 1; else return 0; +} + +NODE nd_reducebase(NODE x) +{ + int len,i,j; + NDV *w; + NODE t,t0; + + len = length(x); + w = (NDV *)ALLOCA(len*sizeof(NDV)); + for ( i = 0, t = x; i < len; i++, t = NEXT(t) ) w[i] = BDY(t); + for ( i = 0; i < len; i++ ) { + for ( j = 0; j < i; j++ ) { + if ( w[i] && w[j] ) + if ( ndl_reducible(HDL(w[i]),HDL(w[j])) ) w[i] = 0; + else if ( ndl_reducible(HDL(w[j]),HDL(w[i])) ) w[j] = 0; + } + } + for ( i = len-1, t0 = 0; i >= 0; i-- ) { + if ( w[i] ) { NEXTNODE(t0,t); BDY(t) = (pointer)w[i]; } + } + NEXT(t) = 0; x = t0; + return x; }