=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/dist.c,v retrieving revision 1.38 retrieving revision 1.46 diff -u -p -r1.38 -r1.46 --- OpenXM_contrib2/asir2000/engine/dist.c 2005/11/27 00:07:05 1.38 +++ OpenXM_contrib2/asir2000/engine/dist.c 2013/11/05 02:55:03 1.46 @@ -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.37 2005/11/26 01:28:12 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.45 2012/12/17 07:20:44 noro Exp $ */ #include "ca.h" @@ -65,8 +65,10 @@ #define ORD_HOMO_WW_DRL_ZIGZAG 13 int cmpdl_drl_zigzag(), cmpdl_homo_ww_drl_zigzag(); +int cmpdl_top_weight(); int (*cmpdl)()=cmpdl_revgradlex; +int (*cmpdl_tie_breaker)(); int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex}; int do_weyl; @@ -118,8 +120,13 @@ int has_sfcoef_p(P f) } } +extern N *current_top_weight_vector; +static int current_top_weight_len; + void initd(struct order_spec *spec) { + int len,i; + switch ( spec->id ) { case 3: cmpdl = cmpdl_composite; @@ -165,6 +172,14 @@ void initd(struct order_spec *spec) } break; } + if ( current_top_weight_vector ) { + 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); + current_top_weight_len = len; + } dp_current_spec = spec; } @@ -509,7 +524,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); @@ -519,7 +534,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: @@ -642,7 +657,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; } @@ -680,7 +695,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; } @@ -814,7 +829,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; } @@ -846,7 +861,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; } @@ -858,8 +873,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)); @@ -921,7 +936,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; } @@ -947,8 +962,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)); @@ -1089,16 +1104,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 */ @@ -1121,6 +1130,8 @@ int compd(VL vl,DP p1,DP p2) return p2 ? -1 : 0; else if ( !p2 ) return 1; + else if ( NV(p1) != NV(p2) ) + error("compd : size mismatch"); else { for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2); m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) ) @@ -1496,6 +1507,8 @@ int cmpdl_order_pair(int n,DL d1,DL d2) struct order_pair *pair; len = dp_current_spec->ord.block.length; + if ( n != dp_current_spec->nv ) + error("cmpdl_order_pair : incompatible order specification"); pair = dp_current_spec->ord.block.order_pair; head = 0; @@ -1636,6 +1649,50 @@ int cmpdl_matrix(int n,DL d1,DL d2) return 0; } +int cmpdl_top_weight(int n,DL d1,DL d2) +{ + int *w; + N sum,wm,wma,t; + N *a; + struct oN tn; + int len,i,sgn,tsgn; + int *t1,*t2; + + w = (int *)ALLOCA(n*sizeof(int)); + len = current_top_weight_len+3; + t1 = d1->d; t2 = d2->d; + for ( i = 0; i < n; i++ ) w[i] = t1[i]-t2[i]; + 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 ( sgn > 0 ) return 1; + else if ( sgn < 0 ) return -1; + else return (*cmpdl_tie_breaker)(n,d1,d2); +} + GeoBucket create_bucket() { GeoBucket g; @@ -1799,17 +1856,26 @@ int ni_next(int *a,int n) int comp_nbm(NBM a,NBM b) { - int d,i,w; + int d,i,w,ai,bi; int *ab,*bb; if ( a->d > b->d ) return 1; else if ( a->d < b->d ) return -1; else { d = a->d; ab = a->b; bb = b->b; +#if 0 w = (d+31)/32; for ( i = 0; i < w; i++ ) if ( ab[i] > bb[i] ) return 1; else if ( ab[i] < bb[i] ) return -1; +#else + for ( i = 0; i < d; i++ ) { + ai = NBM_GET(ab,i); + bi = NBM_GET(bb,i); + if ( ai > bi ) return 1; + else if ( ai < bi ) return -1; + } +#endif return 0; } } @@ -1819,14 +1885,13 @@ NBM mul_nbm(NBM a,NBM b) int ad,bd,d,i,j; int *ab,*bb,*mb; NBM m; - Q c,c1; NODE r; NBP u; ad = a->d; bd = b->d; ab = a->b; bb = b->b; d = ad + bd; NEWNBM(m); NEWNBMBDY(m,d); - m->d = d; mulq(a->c,b->c,&m->c); mb = m->b; + m->d = d; mulp(CO,a->c,b->c,&m->c); mb = m->b; j = 0; for ( i = 0; i < ad; i++, j++ ) if ( NBM_GET(ab,i) ) NBM_SET(mb,j); @@ -1849,7 +1914,7 @@ NBP nbmtonbp(NBM m) /* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */ -Q separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) +P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) { int i,d1; NBM t; @@ -1857,12 +1922,12 @@ Q separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) if ( !a->d ) error("separate_nbm : invalid argument"); if ( a0 ) { - NEWNBM(t); t->d = a->d; t->b = a->b; t->c = ONE; + NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE; *a0 = nbmtonbp(t); } if ( ah ) { - NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = ONE; + NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE; if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0); else NBM_CLR(t->b,0); *ah = nbmtonbp(t); @@ -1870,7 +1935,7 @@ Q separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) if ( ar ) { d1 = a->d-1; - NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = ONE; + NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE; for ( i = 0; i < d1; i++ ) { if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i); else NBM_CLR(t->b,i); @@ -1881,12 +1946,46 @@ Q separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) return a->c; } +/* a=c*rest*x -> a0= rest*x, ar=rest, at=x */ + +P separate_tail_nbm(NBM a,NBP *a0,NBP *ar,NBP *at) +{ + int i,d,d1; + NBM t; + + if ( !(d=a->d) ) error("separate_tail_nbm : invalid argument"); + + if ( a0 ) { + NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE; + *a0 = nbmtonbp(t); + } + + d1 = a->d-1; + if ( at ) { + NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE; + if ( NBM_GET(a->b,d1) ) NBM_SET(t->b,0); + else NBM_CLR(t->b,0); + *at = nbmtonbp(t); + } + + if ( ar ) { + NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE; + for ( i = 0; i < d1; i++ ) { + if ( NBM_GET(a->b,i) ) NBM_SET(t->b,i); + else NBM_CLR(t->b,i); + } + *ar = nbmtonbp(t); + } + + return a->c; +} + NBP make_xky(int k) { int k1,i; NBM t; - NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = ONE; + NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE; k1 = k-1; for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i); NBM_CLR(t->b,i); @@ -1895,7 +1994,7 @@ NBP make_xky(int k) /* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */ -Q separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) +P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) { int i,d1,k,k1; NBM t; @@ -1909,12 +2008,12 @@ Q separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) k = i+1; if ( a0 ) { - NEWNBM(t); t->d = a->d; t->b = a->b; t->c = ONE; + NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE; *a0 = nbmtonbp(t); } if ( ah ) { - NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = ONE; + NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE; for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i); NBM_CLR(t->b,i); *ah = nbmtonbp(t); @@ -1922,7 +2021,7 @@ Q separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) if ( ar ) { d1 = a->d-k; - NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = ONE; + NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE; for ( i = 0; i < d1; i++ ) { if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i); else NBM_CLR(t->b,i); @@ -1941,14 +2040,14 @@ void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp); NBP shuffle_mul_nbm(NBM a,NBM b) { NBP u,a0,ah,ar,b0,bh,br,a1,b1,t; - Q ac,bc,c; + P ac,bc,c; if ( !a->d || !b->d ) u = nbmtonbp(mul_nbm(a,b)); else { ac = separate_nbm(a,&a0,&ah,&ar); bc = separate_nbm(b,&b0,&bh,&br); - mulq(ac,bc,&c); + mulp(CO,ac,bc,&c); shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1); shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1); addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u); @@ -1959,15 +2058,15 @@ NBP shuffle_mul_nbm(NBM a,NBM b) NBP harmonic_mul_nbm(NBM a,NBM b) { NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1; - Q ac,bc,c; + P ac,bc,c; if ( !a->d || !b->d ) u = nbmtonbp(mul_nbm(a,b)); else { - mulq(a->c,b->c,&c); + mulp(CO,a->c,b->c,&c); ac = separate_xky_nbm(a,&a0,&ah,&ar); bc = separate_xky_nbm(b,&b0,&bh,&br); - mulq(ac,bc,&c); + mulp(CO,ac,bc,&c); harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1); harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1); abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d); @@ -1982,7 +2081,7 @@ void addnbp(VL vl,NBP p1,NBP p2, NBP *rp) { NODE b1,b2,br,br0; NBM m1,m2,m; - Q c; + P c; if ( !p1 ) *rp = p2; @@ -1993,7 +2092,7 @@ void addnbp(VL vl,NBP p1,NBP p2, NBP *rp) m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2); switch ( comp_nbm(m1,m2) ) { case 0: - addq(m1->c,m2->c,&c); + addp(CO,m1->c,m2->c,&c); if ( c ) { NEXTNODE(br0,br); NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b; @@ -2043,7 +2142,7 @@ void chsgnnbp(NBP p,NBP *rp) for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) { NEXTNODE(r0,r); m = (NBM)BDY(b); - NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnq(m->c,&m1->c); + NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnp(m->c,&m1->c); BDY(r) = m1; } if ( r0 ) NEXT(r) = 0; @@ -2060,15 +2159,15 @@ void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) *rp = 0; return; } if ( OID(p1) != O_NBP ) { - if ( !NUM(p1) || !RATN(p1) ) + if ( !POLY(p1) ) error("mulnbp : invalid argument"); - NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p1; + NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1; MKNODE(n,m,0); MKNBP(p1,n); } if ( OID(p2) != O_NBP ) { - if ( !NUM(p2) || !RATN(p2) ) + if ( !POLY(p2) ) error("mulnbp : invalid argument"); - NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p2; + NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2; MKNODE(n,m,0); MKNBP(p2,n); } if ( length(BDY(p1)) < length(BDY(p2)) ) { @@ -2126,7 +2225,7 @@ void pwrnbp(VL vl,NBP a,Q q,NBP *c) NODE r; if ( !q ) { - NEWNBM(m); m->d = 0; m->c = ONE; m->b = 0; + NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0; MKNODE(r,m,0); MKNBP(*c,r); } else if ( !a ) *c = 0; @@ -2157,7 +2256,7 @@ int compnbp(VL vl,NBP p1,NBP p2) for ( n1 = BDY(p1), n2 = BDY(p2); n1 && n2; n1 = NEXT(n1), n2 = NEXT(n2) ) { m1 = (NBM)BDY(n1); m2 = (NBM)BDY(n2); - if ( (t = comp_nbm(m1,m2)) || (t = cmpq(m1->c,m2->c) ) ) + if ( (t = comp_nbm(m1,m2)) || (t = compp(CO,m1->c,m2->c) ) ) return t; } if ( n1 ) @@ -2179,15 +2278,15 @@ void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) *rp = 0; return; } if ( OID(p1) != O_NBP ) { - if ( !NUM(p1) || !RATN(p1) ) + if ( !POLY(p1) ) error("shuffle_mulnbp : invalid argument"); - NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p1; + NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1; MKNODE(n,m,0); MKNBP(p1,n); } if ( OID(p2) != O_NBP ) { - if ( !NUM(p2) || !RATN(p2) ) + if ( !POLY(p2) ) error("shuffle_mulnbp : invalid argument"); - NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p2; + NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2; MKNODE(n,m,0); MKNBP(p2,n); } for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) ) @@ -2208,15 +2307,15 @@ void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) *rp = 0; return; } if ( OID(p1) != O_NBP ) { - if ( !NUM(p1) || !RATN(p1) ) + if ( !POLY(p1) ) error("harmonic_mulnbp : invalid argument"); - NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p1; + NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1; MKNODE(n,m,0); MKNBP(p1,n); } if ( OID(p2) != O_NBP ) { - if ( !NUM(p2) || !RATN(p2) ) + if ( !POLY(p2) ) error("harmonic_mulnbp : invalid argument"); - NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p2; + NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2; MKNODE(n,m,0); MKNBP(p2,n); } for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) ) @@ -2233,7 +2332,7 @@ NBP shuffle_mul_nbm(NBM a,NBM b) int ad,bd,d,i,ai,bi,bit,s; int *ab,*bb,*wmb,*w; NBM wm,tm; - Q c,c1; + P c,c1; NODE r,t,t1,p; NBP u; @@ -2243,7 +2342,7 @@ NBP shuffle_mul_nbm(NBM a,NBM b) NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b; for ( i = 0; i < ad; i++ ) w[i] = 1; for ( ; i < d; i++ ) w[i] = 0; - mulq(a->c,b->c,&c); + mulp(CO,a->c,b->c,&c); r = 0; do { wm->d = d; wm->c = c; @@ -2266,7 +2365,7 @@ NBP shuffle_mul_nbm(NBM a,NBM b) break; } else if ( s == 0 ) { /* add coefs */ - addq(tm->c,c,&c1); + addp(CO,tm->c,c,&c1); if ( c1 ) tm->c = c1; else NEXT(p) = NEXT(t); break; @@ -2304,7 +2403,7 @@ NBP harmonic_mul_nbm(NBM a,NBM b) int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1; int i,j,k,ia,ib,s; int *wa,*wb,*w,*wab,*wa1,*wmb; - Q c,c1; + P c,c1; NBM wm,tm; NODE r,t1,t,p; NBP u; @@ -2314,7 +2413,7 @@ NBP harmonic_mul_nbm(NBM a,NBM b) wb = (int *)ALLOCA(db*sizeof(int)); la = nbmtoxky(a,wa); lb = nbmtoxky(b,wb); - mulq(a->c,b->c,&c); + mulp(CO,a->c,b->c,&c); /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */ /* lmax : total length */ lmax = la+lb; @@ -2361,7 +2460,7 @@ NBP harmonic_mul_nbm(NBM a,NBM b) break; } else if ( s == 0 ) { /* add coefs */ - addq(tm->c,c,&c1); + addp(CO,tm->c,c,&c1); if ( c1 ) tm->c = c1; else NEXT(p) = NEXT(t); break;