=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/dist.c,v retrieving revision 1.43 retrieving revision 1.50 diff -u -p -r1.43 -r1.50 --- OpenXM_contrib2/asir2000/engine/dist.c 2007/09/17 12:47:45 1.43 +++ OpenXM_contrib2/asir2000/engine/dist.c 2014/10/10 09:02:25 1.50 @@ -45,7 +45,7 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.42 2006/08/27 22:17:27 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.49 2014/09/12 06:28:46 noro Exp $ */ #include "ca.h" @@ -71,6 +71,9 @@ int (*cmpdl)()=cmpdl_revgradlex; int (*cmpdl_tie_breaker)(); int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex}; +Obj current_top_weight; +int current_top_weight_len; + int do_weyl; int dp_nelim,dp_fcoeffs; @@ -82,6 +85,8 @@ void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr); void comm_quod(VL vl,DP p1,DP p2,DP *pr); void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr); void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr); +int create_order_spec(VL vl,Obj obj,struct order_spec **specp); +void create_modorder_spec(int id,LIST shift,struct modorder_spec **s); void order_init() { @@ -92,6 +97,8 @@ void order_init() create_modorder_spec(0,0,&dp_current_modspec); } +int has_sfcoef_p(P f); + int has_sfcoef(DP f) { MP t; @@ -120,12 +127,21 @@ int has_sfcoef_p(P f) } } -extern N *current_top_weight_vector; -static int current_top_weight_len; +extern Obj nd_top_weight; +void reset_top_weight() +{ + cmpdl = cmpdl_tie_breaker; + cmpdl_tie_breaker = 0; + nd_top_weight = 0; + current_top_weight = 0; + current_top_weight_len = 0; +} + void initd(struct order_spec *spec) { - int len,i; + int len,i,k,row; + Q **mat; switch ( spec->id ) { case 3: @@ -172,12 +188,20 @@ void initd(struct order_spec *spec) } break; } - if ( current_top_weight_vector ) { + if ( current_top_weight ) { cmpdl_tie_breaker = cmpdl; cmpdl = cmpdl_top_weight; - for ( len = 0, i = 0; i < spec->nv; i++ ) - if ( current_top_weight_vector[i] ) - len = MAX(PL(current_top_weight_vector[i]),len); + if ( OID(current_top_weight) == O_VECT ) { + mat = (Q **)&BDY((VECT)current_top_weight); + row = 1; + } else { + mat = (Q **)BDY((MAT)current_top_weight); + row = ((MAT)current_top_weight)->row; + } + for ( k = 0, len = 0; k < row; k++ ) + for ( i = 0; i < spec->nv; i++ ) + if ( mat[k][i] ) + len = MAX(PL(NM(mat[k][i])),len); current_top_weight_len = len; } dp_current_spec = spec; @@ -185,7 +209,6 @@ void initd(struct order_spec *spec) void ptod(VL vl,VL dvl,P p,DP *pr) { - int isconst = 0; int n,i,j,k; VL tvl; V v; @@ -313,7 +336,7 @@ int sugard(MP m) void addd(VL vl,DP p1,DP p2,DP *pr) { int n; - MP m1,m2,mr,mr0; + MP m1,m2,mr=0,mr0; P t; DL d; @@ -373,7 +396,7 @@ void addd(VL vl,DP p1,DP p2,DP *pr) void symb_addd(DP p1,DP p2,DP *pr) { int n; - MP m1,m2,mr,mr0; + MP m1,m2,mr=0,mr0; if ( !p1 ) *pr = p2; @@ -424,8 +447,7 @@ void symb_addd(DP p1,DP p2,DP *pr) NODE symb_merge(NODE m1,NODE m2,int n) { - NODE top,prev,cur,m,t; - int c,i; + NODE top=0,prev,cur,m=0,t; DL d1,d2; if ( !m1 ) @@ -524,7 +546,7 @@ NODE symb_mul_merge(NODE m1,DL dl,DP f,int n) if ( !cur ) { MKDP(n,m,g); NEXT(prev) = mul_dllist(dl,g); - return; + return top; } m = NEXT(m); if ( m ) _adddl(n,m->dl,dl,t); @@ -534,7 +556,7 @@ NODE symb_mul_merge(NODE m1,DL dl,DP f,int n) if ( !cur ) { MKDP(n,m,g); NEXT(prev) = mul_dllist(dl,g); - return; + return top; } break; case -1: @@ -610,7 +632,7 @@ void subd(VL vl,DP p1,DP p2,DP *pr) void chsgnd(DP p,DP *pr) { - MP m,mr,mr0; + MP m,mr=0,mr0; Obj r; if ( !p ) @@ -657,7 +679,7 @@ void comm_muld(VL vl,DP p1,DP p2,DP *pr) l = l1; } if ( l > wlen ) { - if ( w ) GC_free(w); + if ( w ) GCFREE(w); w = (MP *)MALLOC(l*sizeof(MP)); wlen = l; } @@ -695,7 +717,7 @@ void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr) l = l1; } if ( l > wlen ) { - if ( w ) GC_free(w); + if ( w ) GCFREE(w); w = (MP *)MALLOC(l*sizeof(MP)); wlen = l; } @@ -711,7 +733,7 @@ void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr) void comm_quod(VL vl,DP p1,DP p2,DP *pr) { - MP m,m0; + MP m=0,m0; DP s,t; int i,n,sugar; DL d1,d2,d; @@ -752,7 +774,7 @@ void comm_quod(VL vl,DP p1,DP p2,DP *pr) void muldm(VL vl,DP p,MP m0,DP *pr) { - MP m,mr,mr0; + MP m,mr=0,mr0; P c; DL d; int n; @@ -777,7 +799,7 @@ void muldm(VL vl,DP p,MP m0,DP *pr) void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr) { - MP m,mr,mr0; + MP m,mr=0,mr0; P c; DL d,tdl; int n,i; @@ -829,7 +851,7 @@ void weyl_muld(VL vl,DP p1,DP p2,DP *pr) else { for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ ); if ( l > wlen ) { - if ( w ) GC_free(w); + if ( w ) GCFREE(w); w = (MP *)MALLOC(l*sizeof(MP)); wlen = l; } @@ -861,7 +883,7 @@ void weyl_muldm(VL vl,MP m0,DP p,DP *pr) else { for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ ); if ( l > wlen ) { - if ( w ) GC_free(w); + if ( w ) GCFREE(w); w = (MP *)MALLOC(l*sizeof(MP)); wlen = l; } @@ -873,8 +895,8 @@ void weyl_muldm(VL vl,MP m0,DP p,DP *pr) for ( i = 0, tlen = 1; i < n2; i++ ) tlen *= d0->d[n2+i]+1; if ( tlen > rtlen ) { - if ( tab ) GC_free(tab); - if ( psum ) GC_free(psum); + if ( tab ) GCFREE(tab); + if ( psum ) GCFREE(psum); rtlen = tlen; tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl)); psum = (MP *)MALLOC(rtlen*sizeof(MP)); @@ -936,7 +958,7 @@ void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rt rtab[0].d = d; if ( rtablen > tmptablen ) { - if ( tmptab ) GC_free(tmptab); + if ( tmptab ) GCFREE(tmptab); tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl)); tmptablen = rtablen; } @@ -962,8 +984,8 @@ void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rt continue; } if ( k+1 > tablen ) { - if ( tab ) GC_free(tab); - if ( ctab ) GC_free(ctab); + if ( tab ) GCFREE(tab); + if ( ctab ) GCFREE(ctab); tablen = k+1; tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl)); ctab = (Q *)MALLOC(tablen*sizeof(Q)); @@ -1029,7 +1051,7 @@ void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,st void muldc(VL vl,DP p,P c,DP *pr) { - MP m,mr,mr0; + MP m,mr=0,mr0; if ( !p || !c ) *pr = 0; @@ -1054,7 +1076,7 @@ void muldc(VL vl,DP p,P c,DP *pr) void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr) { - MP m,mr,mr0; + MP m,mr=0,mr0; DL mdl; int i,n; @@ -1083,7 +1105,7 @@ void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr) void divsdc(VL vl,DP p,P c,DP *pr) { - MP m,mr,mr0; + MP m,mr=0,mr0; if ( !c ) error("disvsdc : division by 0"); @@ -1104,16 +1126,10 @@ void adddl(int n,DL d1,DL d2,DL *dr) DL dt; int i; - if ( !d1->td ) - *dr = d2; - else if ( !d2->td ) - *dr = d1; - else { - *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int)); - dt->td = d1->td + d2->td; - for ( i = 0; i < n; i++ ) - dt->d[i] = d1->d[i]+d2->d[i]; - } + *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int)); + dt->td = d1->td + d2->td; + for ( i = 0; i < n; i++ ) + dt->d[i] = d1->d[i]+d2->d[i]; } /* d1 += d2 */ @@ -1136,9 +1152,10 @@ int compd(VL vl,DP p1,DP p2) return p2 ? -1 : 0; else if ( !p2 ) return 1; - else if ( NV(p1) != NV(p2) ) + else if ( NV(p1) != NV(p2) ) { error("compd : size mismatch"); - else { + return 0; /* XXX */ + } else { for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2); m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) ) if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) || @@ -1259,7 +1276,7 @@ int cmpdl_blex(int n,DL d1,DL d2) { int c; - if ( c = cmpdl_lex(n-1,d1,d2) ) + if ( (c = cmpdl_lex(n-1,d1,d2)) ) return c; else { c = d1->d[n-1] - d2->d[n-1]; @@ -1464,8 +1481,8 @@ int cmpdl_drl_zigzag(int n,DL d1,DL d2) else { m = n>>1; for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) { - if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1; - if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1; + if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1; + if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1; } return 0; } @@ -1499,8 +1516,8 @@ int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2) return -1; for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) { - if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1; - if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1; + if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1; + if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1; } return 0; } @@ -1659,9 +1676,10 @@ int cmpdl_top_weight(int n,DL d1,DL d2) { int *w; N sum,wm,wma,t; - N *a; + Q **mat; + Q *a; struct oN tn; - int len,i,sgn,tsgn; + int len,i,sgn,tsgn,row,k; int *t1,*t2; w = (int *)ALLOCA(n*sizeof(int)); @@ -1671,32 +1689,41 @@ int cmpdl_top_weight(int n,DL d1,DL d2) sum = (N)W_ALLOC(len); sgn = 0; wm = (N)W_ALLOC(len); wma = (N)W_ALLOC(len); - a = current_top_weight_vector; - for ( i = 0; i < n; i++ ) { - if ( !a[i] || !w[i] ) continue; - tn.p = 1; - if ( w[i] > 0 ) { - tn.b[0] = w[i]; tsgn = 1; - } else { - tn.b[0] = -w[i]; tsgn = -1; - } - _muln(a[i],&tn,wm); - if ( !sgn ) { - sgn = tsgn; - t = wm; wm = sum; sum = t; - } else if ( sgn == tsgn ) { - _addn(sum,wm,wma); - if ( !PL(wma) ) - sgn = 0; - t = wma; wma = sum; sum = t; - } else { - sgn *= _subn(sum,wm,wma); - t = wma; wma = sum; sum = t; - } + if ( OID(current_top_weight) == O_VECT ) { + mat = (Q **)&BDY((VECT)current_top_weight); + row = 1; + } else { + mat = (Q **)BDY((MAT)current_top_weight); + row = ((MAT)current_top_weight)->row; } - if ( sgn > 0 ) return 1; - else if ( sgn < 0 ) return -1; - else return (*cmpdl_tie_breaker)(n,d1,d2); + for ( k = 0; k < row; k++ ) { + a = mat[k]; + for ( i = 0; i < n; i++ ) { + if ( !a[i] || !w[i] ) continue; + tn.p = 1; + if ( w[i] > 0 ) { + tn.b[0] = w[i]; tsgn = 1; + } else { + tn.b[0] = -w[i]; tsgn = -1; + } + _muln(NM(a[i]),&tn,wm); + if ( !sgn ) { + sgn = tsgn; + t = wm; wm = sum; sum = t; + } else if ( sgn == tsgn ) { + _addn(sum,wm,wma); + if ( !PL(wma) ) + sgn = 0; + t = wma; wma = sum; sum = t; + } else { + sgn *= _subn(sum,wm,wma); + t = wma; wma = sum; sum = t; + } + } + if ( sgn > 0 ) return 1; + else if ( sgn < 0 ) return -1; + } + return (*cmpdl_tie_breaker)(n,d1,d2); } GeoBucket create_bucket() @@ -1708,6 +1735,8 @@ GeoBucket create_bucket() return g; } +int length(NODE d); + void add_bucket(GeoBucket g,NODE d,int nv) { int l,k,m; @@ -1831,12 +1860,13 @@ int compdv(VL vl,DPV p1,DPV p2) { int i,t,len; - if ( p1->len != p2->len ) + if ( p1->len != p2->len ) { error("compdv : size mismatch"); - else { + return 0; /* XXX */ + } else { len = p1->len; for ( i = 0; i < len; i++ ) - if ( t = compd(vl,p1->body[i],p2->body[i]) ) + if ( (t = compd(vl,p1->body[i],p2->body[i])) ) return t; return 0; } @@ -1862,7 +1892,7 @@ int ni_next(int *a,int n) int comp_nbm(NBM a,NBM b) { - int d,i,w,ai,bi; + int d,i,ai,bi; int *ab,*bb; if ( a->d > b->d ) return 1; @@ -1891,8 +1921,6 @@ NBM mul_nbm(NBM a,NBM b) int ad,bd,d,i,j; int *ab,*bb,*mb; NBM m; - NODE r; - NBP u; ad = a->d; bd = b->d; ab = a->b; bb = b->b; d = ad + bd; @@ -2085,7 +2113,7 @@ NBP harmonic_mul_nbm(NBM a,NBM b) void addnbp(VL vl,NBP p1,NBP p2, NBP *rp) { - NODE b1,b2,br,br0; + NODE b1,b2,br=0,br0; NBM m1,m2,m; P c; @@ -2142,7 +2170,7 @@ void subnbp(VL vl,NBP p1,NBP p2, NBP *rp) void chsgnnbp(NBP p,NBP *rp) { - NODE r0,r,b; + NODE r0,r=0,b; NBM m,m1; for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) { @@ -2193,7 +2221,7 @@ void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp) { - NODE b,r0,r; + NODE b,r0,r=0; if ( !p ) *rp = 0; else { @@ -2208,7 +2236,7 @@ void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp) void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp) { - NODE b,r0,r; + NODE b,r0,r=0; if ( !p ) *rp = 0; else {