=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v retrieving revision 1.199 retrieving revision 1.207 diff -u -p -r1.199 -r1.207 --- OpenXM_contrib2/asir2000/engine/nd.c 2012/08/27 05:38:00 1.199 +++ OpenXM_contrib2/asir2000/engine/nd.c 2013/09/12 06:46:16 1.207 @@ -1,10 +1,11 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.198 2012/04/10 07:15:07 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.206 2013/09/10 02:10:00 noro Exp $ */ #include "nd.h" struct oEGT eg_search; int diag_period = 6; +int weight_check = 1; int (*ndl_compare_function)(UINT *a1,UINT *a2); int nd_dcomp; NM _nm_free_list; @@ -64,7 +65,6 @@ LIST ndvtopl(int mod,VL vl,VL dvl,NDV p,int rank); 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); -pointer GC_malloc_atomic_ignore_off_page(int); void nmtodp(int mod,NM m,DP *r); NODE reverse_node(NODE n); P ndc_div(int mod,union oNDC a,union oNDC b); @@ -74,6 +74,7 @@ void conv_ilist(int demand,int trace,NODE g,int **indp 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); extern int Denominator,DP_Multiple; @@ -92,7 +93,7 @@ void _NM_alloc() int i; for ( i = 0; i < 1024; i++ ) { - p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(UINT)); + p = (NM)MALLOC(sizeof(struct oNM)+(nd_wpd-1)*sizeof(UINT)); p->next = _nm_free_list; _nm_free_list = p; } } @@ -103,7 +104,7 @@ void _ND_alloc() int i; for ( i = 0; i < 1024; i++ ) { - p = (ND)GC_malloc(sizeof(struct oND)); + p = (ND)MALLOC(sizeof(struct oND)); p->body = (NM)_nd_free_list; _nd_free_list = p; } } @@ -114,7 +115,7 @@ void _NDP_alloc() int i; for ( i = 0; i < 1024; i++ ) { - p = (ND_pairs)GC_malloc(sizeof(struct oND_pairs) + p = (ND_pairs)MALLOC(sizeof(struct oND_pairs) +(nd_wpd-1)*sizeof(UINT)); p->next = _ndp_free_list; _ndp_free_list = p; } @@ -1581,7 +1582,7 @@ void free_pbucket(PGeoBucket b) { nd_free(b->body[i]); b->body[i] = 0; } - GC_free(b); + GCFREE(b); } void add_pbucket_symbolic(PGeoBucket g,ND d) @@ -3842,7 +3843,7 @@ void nd_free(ND p) void ndv_free(NDV p) { - GC_free(BDY(p)); + GCFREE(BDY(p)); } void nd_append_red(UINT *d,int i) @@ -3953,6 +3954,31 @@ void nd_setup_parameters(int nvar,int max) { else if ( max < 65536 ) nd_bpe = 16; else nd_bpe = 32; } + if ( !do_weyl && weight_check && (current_dl_weight_vector || nd_matrix) ) { + UINT t; + int st; + int *v; + /* t = max(weights) */ + t = 0; + if ( current_dl_weight_vector ) + for ( i = 0, t = 0; i < nd_nvar; i++ ) { + if ( (st=current_dl_weight_vector[i]) < 0 ) st = -st; + if ( t < st ) t = st; + } + if ( nd_matrix ) + for ( i = 0; i < nd_matrix_len; i++ ) + for ( j = 0, v = nd_matrix[i]; j < nd_nvar; j++ ) { + if ( (st=v[j]) < 0 ) st = -st; + if ( t < st ) t = st; + } + /* i = bitsize of t */ + for ( i = 0; t; t >>=1, i++ ); + /* i += bitsize of nd_nvar */ + for ( t = nd_nvar; t; t >>=1, i++); + /* nd_bpe+i = bitsize of max(weights)*max(exp)*nd_nvar */ + if ( (nd_bpe+i) >= 31 ) + error("nd_setup_parameters : too large weight"); + } nd_epw = (sizeof(UINT)*8)/nd_bpe; elen = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0); nd_exporigin = nd_get_exporigin(nd_ord); @@ -4807,7 +4833,7 @@ NDV ndtondv(int mod,ND p) if ( !p ) return 0; len = LEN(p); if ( mod ) - m0 = m = (NMV)GC_malloc_atomic_ignore_off_page(len*nmv_adv); + m0 = m = (NMV)MALLOC_ATOMIC_IGNORE_OFF_PAGE(len*nmv_adv); else m0 = m = MALLOC(len*nmv_adv); #if 0 @@ -4864,6 +4890,27 @@ DP ndvtodp(int mod,NDV p) return d; } +DP ndtodp(int mod,ND p) +{ + MP m,m0; + DP d; + NM t; + int i,len; + + if ( !p ) return 0; + m0 = 0; + len = p->len; + for ( t = BDY(p); t; t = NEXT(t) ) { + NEXTMP(m0,m); + m->dl = ndltodl(nd_nvar,DL(t)); + m->c = ndctop(mod,t->c); + } + NEXT(m) = 0; + MKDP(nd_nvar,m0,d); + SG(d) = SG(p); + return d; +} + void ndv_print(NDV p) { NMV m; @@ -4932,6 +4979,8 @@ NODE ndv_reducebase(NODE x,int *perm) void nd_init_ord(struct order_spec *ord) { nd_module = (ord->id >= 256); + nd_matrix = 0; + nd_matrix_len = 0; switch ( ord->id ) { case 0: switch ( ord->ord.simple ) { @@ -5502,7 +5551,7 @@ NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea for ( j = 0, len = 0; j < spcol; j++ ) if ( vect[j] ) len++; if ( !len ) return 0; else { - mr0 = (NMV)GC_malloc_atomic_ignore_off_page(nmv_adv*len); + mr0 = (NMV)MALLOC_ATOMIC_IGNORE_OFF_PAGE(nmv_adv*len); #if 0 ndv_alloc += nmv_adv*len; #endif @@ -5532,7 +5581,7 @@ NDV vect_to_ndv_q(Q *vect,int spcol,int col,int *rhead for ( j = 0, len = 0; j < spcol; j++ ) if ( vect[j] ) len++; if ( !len ) return 0; else { - mr0 = (NMV)GC_malloc(nmv_adv*len); + mr0 = (NMV)MALLOC(nmv_adv*len); #if 0 ndv_alloc += nmv_adv*len; #endif @@ -5564,7 +5613,7 @@ NDV plain_vect_to_ndv_q(Q *vect,int col,UINT *s0vect) for ( j = 0, len = 0; j < col; j++ ) if ( vect[j] ) len++; if ( !len ) return 0; else { - mr0 = (NMV)GC_malloc(nmv_adv*len); + mr0 = (NMV)MALLOC(nmv_adv*len); #if 0 ndv_alloc += nmv_adv*len; #endif @@ -6005,7 +6054,7 @@ NODE nd_f4_red_main(int m,ND_pairs sp0,int nsp,UINT *s fflush(asir_out); } /* free index arrays */ - for ( i = 0; i < nred; i++ ) GC_free(imat[i]->index.c); + for ( i = 0; i < nred; i++ ) GCFREE(imat[i]->index.c); /* elimination (2nd step) */ colstat = (int *)ALLOCA(spcol*sizeof(int)); @@ -6018,11 +6067,11 @@ NODE nd_f4_red_main(int m,ND_pairs sp0,int nsp,UINT *s NEXTNODE(r0,r); BDY(r) = (pointer)vect_to_ndv(spmat[i],spcol,col,rhead,s0vect); SG((NDV)BDY(r)) = spsugar[i]; - GC_free(spmat[i]); + GCFREE(spmat[i]); } if ( r0 ) NEXT(r) = 0; - for ( ; i < sprow; i++ ) GC_free(spmat[i]); + 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 ) { @@ -6086,7 +6135,7 @@ NODE nd_f4_red_q_main(ND_pairs sp0,int nsp,int trace,U fflush(asir_out); } /* free index arrays */ -/* for ( i = 0; i < nred; i++ ) GC_free(imat[i]->index.c); */ +/* for ( i = 0; i < nred; i++ ) GCFREE(imat[i]->index.c); */ /* elimination (2nd step) */ colstat = (int *)ALLOCA(spcol*sizeof(int)); @@ -6095,7 +6144,7 @@ NODE nd_f4_red_q_main(ND_pairs sp0,int nsp,int trace,U for ( i = 0; i < rank; i++ ) { w[rank-i-1] = (pointer)vect_to_ndv_q(spmat[i],spcol,col,rhead,s0vect); SG((NDV)w[rank-i-1]) = spsugar[i]; -/* GC_free(spmat[i]); */ +/* GCFREE(spmat[i]); */ } #if 0 qsort(w,rank,sizeof(NDV), @@ -6107,7 +6156,7 @@ NODE nd_f4_red_q_main(ND_pairs sp0,int nsp,int trace,U } if ( r0 ) NEXT(r) = 0; -/* for ( ; i < sprow; i++ ) GC_free(spmat[i]); */ +/* 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 ) { @@ -6239,7 +6288,7 @@ NDV nd_recv_ndv() len = nd_recv_int(); if ( !len ) return 0; else { - m0 = m = (NMV)GC_malloc_atomic_ignore_off_page(nmv_adv*len); + m0 = m = (NMV)MALLOC_ATOMIC_IGNORE_OFF_PAGE(nmv_adv*len); #if 0 ndv_alloc += len*nmv_adv; #endif @@ -6630,6 +6679,32 @@ void ndv_save(NDV p,int index) fclose(s); } +void nd_save_mod(ND p,int index) +{ + FILE *s; + char name[BUFSIZ]; + int nv,sugar,len,c; + NM m; + + sprintf(name,"%s/%d",Demand,index); + s = fopen(name,"w"); + if ( !p ) { + len = 0; + write_int(s,&len); + fclose(s); + return; + } + nv = NV(p); + sugar = SG(p); + len = LEN(p); + write_int(s,&nv); write_int(s,&sugar); write_int(s,&len); + for ( m = BDY(p); m; m = NEXT(m) ) { + c = CM(m); write_int(s,&c); + write_intarray(s,DL(m),nd_wpd); + } + fclose(s); +} + NDV ndv_load(int index) { FILE *s; @@ -6674,6 +6749,36 @@ NDV ndv_load(int index) return d; } +ND nd_load_mod(int index) +{ + FILE *s; + char name[BUFSIZ]; + int nv,sugar,len,i,c; + ND d; + NM m0,m; + + sprintf(name,"%s/%d",Demand,index); + s = fopen(name,"r"); + /* if the file does not exist, it means p[index]=0 */ + if ( !s ) return 0; + + read_int(s,&nv); + if ( !nv ) { fclose(s); return 0; } + + read_int(s,&sugar); + read_int(s,&len); + for ( m0 = 0, i = 0; i < len; i++ ) { + NEXTNM(m0,m); + read_int(s,&c); CM(m) = c; + read_intarray(s,DL(m),nd_wpd); + } + NEXT(m) = 0; + MKND(nv,m0,len,d); + SG(d) = sugar; + fclose(s); + return d; +} + void nd_det(int mod,MAT f,P *rp) { VL fv,tv; @@ -7168,4 +7273,281 @@ void parse_nd_option(NODE opt) else if ( !strcmp(key,"intersect") ) nd_intersect = value?1:0; } +} + +ND mdptond(DP d); +ND nd_mul_nm(int mod,NM m0,ND p); +ND *btog(NODE ti,ND **p,int nb,int mod); +ND btog_one(NODE ti,ND *p,int nb,int mod); +MAT nd_btog(LIST f,LIST v,int m,struct order_spec *ord,LIST tlist,MAT *rp); +VECT nd_btog_one(LIST f,LIST v,int m,struct order_spec *ord,LIST tlist,int pos,MAT *rp); + +/* d:monomial */ +ND mdptond(DP d) +{ + NM m; + ND r; + + if ( OID(d) == 1 ) + r = ptond(CO,CO,(P)d); + else { + NEWNM(m); + dltondl(NV(d),BDY(d)->dl,DL(m)); + CQ(m) = (Q)BDY(d)->c; + NEXT(m) = 0; + MKND(NV(d),m,1,r); + } + return r; +} + +ND nd_mul_nm(int mod,NM m0,ND p) +{ + UINT *d0; + int c0,c1,c; + NM tm,mr,mr0; + ND r; + + if ( !p ) return 0; + d0 = DL(m0); + c0 = CM(m0); + mr0 = 0; + for ( tm = BDY(p); tm; tm = NEXT(tm) ) { + NEXTNM(mr0,mr); + c = CM(tm); DMAR(c0,c,0,mod,c1); CM(mr) = c1; + ndl_add(d0,DL(tm),DL(mr)); + } + NEXT(mr) = 0; + MKND(NV(p),mr0,LEN(p),r); + return r; +} + +ND *btog(NODE ti,ND **p,int nb,int mod) +{ + PGeoBucket *r; + int i,ci; + NODE t,s; + ND m,tp; + ND *pi,*rd; + P c; + + r = (PGeoBucket *)MALLOC(nb*sizeof(PGeoBucket)); + for ( i = 0; i < nb; i++ ) + r[i] = create_pbucket(); + for ( t = ti; t; t = NEXT(t) ) { + s = BDY((LIST)BDY(t)); + if ( ARG0(s) ) { + m = mdptond((DP)ARG2(s)); + ptomp(mod,(P)HCQ(m),&c); + if ( ci = ((MQ)c)->cont ) { + HCM(m) = ci; + pi = p[QTOS((Q)ARG1(s))]; + for ( i = 0; i < nb; i++ ) { + tp = nd_mul_nm(mod,BDY(m),pi[i]); + add_pbucket(mod,r[i],tp); + } + } + ci = 1; + } else { + ptomp(mod,(P)ARG3(s),&c); ci = ((MQ)c)->cont; + ci = invm(ci,mod); + } + } + rd = (ND *)MALLOC(nb*sizeof(ND)); + for ( i = 0; i < nb; i++ ) + rd[i] = normalize_pbucket(mod,r[i]); + if ( ci != 1 ) + for ( i = 0; i < nb; i++ ) nd_mul_c(mod,rd[i],ci); + return rd; +} + +ND btog_one(NODE ti,ND *p,int nb,int mod) +{ + PGeoBucket r; + int i,ci,j; + NODE t,s; + ND m,tp; + ND pi,rd; + P c; + + r = create_pbucket(); + for ( t = ti; t; t = NEXT(t) ) { + s = BDY((LIST)BDY(t)); + if ( ARG0(s) ) { + m = mdptond((DP)ARG2(s)); + ptomp(mod,(P)HCQ(m),&c); + if ( ci = ((MQ)c)->cont ) { + HCM(m) = ci; + pi = p[j=QTOS((Q)ARG1(s))]; + if ( !pi ) { + pi = nd_load_mod(j); + tp = nd_mul_nm(mod,BDY(m),pi); + nd_free(pi); + add_pbucket(mod,r,tp); + } else { + tp = nd_mul_nm(mod,BDY(m),pi); + add_pbucket(mod,r,tp); + } + } + ci = 1; + } else { + ptomp(mod,(P)ARG3(s),&c); ci = ((MQ)c)->cont; + ci = invm(ci,mod); + } + } + rd = normalize_pbucket(mod,r); + free_pbucket(r); + if ( ci != 1 ) nd_mul_c(mod,rd,ci); + return rd; +} + +MAT nd_btog(LIST f,LIST v,int mod,struct order_spec *ord,LIST tlist,MAT *rp) +{ + int i,j,n,m,nb,pi0,pi1,nvar; + VL fv,tv,vv; + NODE permtrace,perm,trace,intred,ind,t,pi,ti; + ND **p; + ND *c; + ND u; + P inv; + MAT mat; + + parse_nd_option(current_option); + get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); + switch ( ord->id ) { + case 1: + if ( ord->nv != nvar ) + error("nd_check : invalid order specification"); + break; + default: + break; + } + nd_init_ord(ord); +#if 0 + nd_bpe = QTOS((Q)ARG7(BDY(tlist))); +#else + nd_bpe = 32; +#endif + nd_setup_parameters(nvar,0); + permtrace = BDY((LIST)ARG2(BDY(tlist))); + intred = BDY((LIST)ARG3(BDY(tlist))); + ind = BDY((LIST)ARG4(BDY(tlist))); + perm = BDY((LIST)BDY(permtrace)); trace =NEXT(permtrace); + for ( i = length(perm)-1, t = trace; t; t = NEXT(t) ) { + j = QTOS((Q)BDY(BDY((LIST)BDY(t)))); + if ( j > i ) i = j; + } + n = i+1; + nb = length(BDY(f)); + p = (ND **)MALLOC(n*sizeof(ND *)); + for ( t = perm, i = 0; t; t = NEXT(t), i++ ) { + pi = BDY((LIST)BDY(t)); + pi0 = QTOS((Q)ARG0(pi)); pi1 = QTOS((Q)ARG1(pi)); + p[pi0] = c = (ND *)MALLOC(nb*sizeof(ND)); + ptomp(mod,(P)ARG2(pi),&inv); + u = ptond(CO,vv,(P)ONE); + HCM(u) = ((MQ)inv)->cont; + c[pi1] = u; + } + for ( t = trace,i=0; t; t = NEXT(t), i++ ) { + printf("%d ",i); fflush(stdout); + ti = BDY((LIST)BDY(t)); + p[j=QTOS((Q)ARG0(ti))] = btog(BDY((LIST)ARG1(ti)),p,nb,mod); + if ( j == 441 ) + printf("afo"); + } + for ( t = intred, i=0; t; t = NEXT(t), i++ ) { + printf("%d ",i); fflush(stdout); + ti = BDY((LIST)BDY(t)); + p[j=QTOS((Q)ARG0(ti))] = btog(BDY((LIST)ARG1(ti)),p,nb,mod); + if ( j == 441 ) + printf("afo"); + } + m = length(ind); + MKMAT(mat,nb,m); + for ( j = 0, t = ind; j < m; j++, t = NEXT(t) ) + for ( i = 0, c = p[QTOS((Q)BDY(t))]; i < nb; i++ ) + BDY(mat)[i][j] = ndtodp(mod,c[i]); + return mat; +} + +VECT nd_btog_one(LIST f,LIST v,int mod,struct order_spec *ord, + LIST tlist,int pos,MAT *rp) +{ + int i,j,n,m,nb,pi0,pi1,nvar; + VL fv,tv,vv; + NODE permtrace,perm,trace,intred,ind,t,pi,ti; + ND *p; + ND *c; + ND u; + P inv; + VECT vect; + + parse_nd_option(current_option); + get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc); + for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); + switch ( ord->id ) { + case 1: + if ( ord->nv != nvar ) + error("nd_check : invalid order specification"); + break; + default: + break; + } + nd_init_ord(ord); +#if 0 + nd_bpe = QTOS((Q)ARG7(BDY(tlist))); +#else + nd_bpe = 32; +#endif + nd_setup_parameters(nvar,0); + permtrace = BDY((LIST)ARG2(BDY(tlist))); + intred = BDY((LIST)ARG3(BDY(tlist))); + ind = BDY((LIST)ARG4(BDY(tlist))); + perm = BDY((LIST)BDY(permtrace)); trace =NEXT(permtrace); + for ( i = length(perm)-1, t = trace; t; t = NEXT(t) ) { + j = QTOS((Q)BDY(BDY((LIST)BDY(t)))); + if ( j > i ) i = j; + } + n = i+1; + nb = length(BDY(f)); + p = (ND *)MALLOC(n*sizeof(ND *)); + for ( t = perm, i = 0; t; t = NEXT(t), i++ ) { + pi = BDY((LIST)BDY(t)); + pi0 = QTOS((Q)ARG0(pi)); pi1 = QTOS((Q)ARG1(pi)); + if ( pi1 == pos ) { + ptomp(mod,(P)ARG2(pi),&inv); + u = ptond(CO,vv,(P)ONE); + HCM(u) = ((MQ)inv)->cont; + p[pi0] = u; + } + } + for ( t = trace,i=0; t; t = NEXT(t), i++ ) { + printf("%d ",i); fflush(stdout); + ti = BDY((LIST)BDY(t)); + p[j=QTOS((Q)ARG0(ti))] = btog_one(BDY((LIST)ARG1(ti)),p,nb,mod); + if ( Demand ) { + nd_save_mod(p[j],j); nd_free(p[j]); p[j] = 0; + } + } + for ( t = intred, i=0; t; t = NEXT(t), i++ ) { + printf("%d ",i); fflush(stdout); + ti = BDY((LIST)BDY(t)); + p[j=QTOS((Q)ARG0(ti))] = btog_one(BDY((LIST)ARG1(ti)),p,nb,mod); + if ( Demand ) { + nd_save_mod(p[j],j); nd_free(p[j]); p[j] = 0; + } + } + m = length(ind); + MKVECT(vect,m); + for ( j = 0, t = ind; j < m; j++, t = NEXT(t) ) { + u = p[QTOS((Q)BDY(t))]; + if ( !u ) { + u = nd_load_mod(QTOS((Q)BDY(t))); + BDY(vect)[j] = ndtodp(mod,u); + nd_free(u); + } else + BDY(vect)[j] = ndtodp(mod,u); + } + return vect; }