=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/dist.c,v retrieving revision 1.28 retrieving revision 1.48 diff -u -p -r1.28 -r1.48 --- OpenXM_contrib2/asir2000/engine/dist.c 2004/02/05 08:28:53 1.28 +++ OpenXM_contrib2/asir2000/engine/dist.c 2014/08/19 06:35:01 1.48 @@ -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.27 2004/02/03 23:31:57 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.47 2013/12/20 02:02:24 noro Exp $ */ #include "ca.h" @@ -65,21 +65,37 @@ #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; int dp_nelim,dp_fcoeffs; struct order_spec *dp_current_spec; +struct modorder_spec *dp_current_modspec; int *dp_dl_work; 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() +{ + struct order_spec *spec; + + create_order_spec(0,0,&spec); + initd(spec); + create_modorder_spec(0,0,&dp_current_modspec); +} + +int has_sfcoef_p(P f); + int has_sfcoef(DP f) { MP t; @@ -108,8 +124,14 @@ int has_sfcoef_p(P f) } } +extern Obj current_top_weight; +int current_top_weight_len; + void initd(struct order_spec *spec) { + int len,i,k,row; + Q **mat; + switch ( spec->id ) { case 3: cmpdl = cmpdl_composite; @@ -155,12 +177,27 @@ void initd(struct order_spec *spec) } break; } + if ( current_top_weight ) { + cmpdl_tie_breaker = cmpdl; + cmpdl = cmpdl_top_weight; + if ( OID(current_top_weight) == O_VECT ) { + mat = (Q **)&BDY((MAT)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; } void ptod(VL vl,VL dvl,P p,DP *pr) { - int isconst = 0; int n,i,j,k; VL tvl; V v; @@ -288,14 +325,25 @@ 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; if ( !p1 ) *pr = p2; else if ( !p2 ) *pr = p1; else { + if ( OID(p1) <= O_R ) { + n = NV(p2); NEWDL(d,n); + NEWMP(m1); m1->dl = d; C(m1) = (P)p1; NEXT(m1) = 0; + MKDP(n,m1,p1); (p1)->sugar = 0; + } + if ( OID(p2) <= O_R ) { + n = NV(p1); NEWDL(d,n); + NEWMP(m2); m2->dl = d; C(m2) = (P)p2; NEXT(m2) = 0; + MKDP(n,m2,p2); (p2)->sugar = 0; + } for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) switch ( (*cmpdl)(n,m1->dl,m2->dl) ) { case 0: @@ -337,7 +385,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; @@ -388,8 +436,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 ) @@ -488,7 +535,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); @@ -498,7 +545,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: @@ -574,11 +621,14 @@ 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 ) *pr = 0; - else { + else if ( OID(p) <= O_R ) { + chsgnr((Obj)p,&r); *pr = (DP)r; + } else { for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl; } @@ -618,7 +668,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; } @@ -656,7 +706,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; } @@ -672,7 +722,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; @@ -713,7 +763,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; @@ -738,7 +788,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; @@ -790,7 +840,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; } @@ -822,7 +872,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; } @@ -834,8 +884,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)); @@ -897,7 +947,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; } @@ -923,8 +973,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)); @@ -990,7 +1040,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; @@ -1015,7 +1065,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; @@ -1044,7 +1094,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"); @@ -1065,16 +1115,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 */ @@ -1097,7 +1141,10 @@ int compd(VL vl,DP p1,DP p2) return p2 ? -1 : 0; else if ( !p2 ) return 1; - else { + else if ( NV(p1) != NV(p2) ) { + error("compd : size mismatch"); + 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)) || @@ -1218,7 +1265,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]; @@ -1423,8 +1470,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; } @@ -1458,8 +1505,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; } @@ -1472,6 +1519,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; @@ -1612,6 +1661,60 @@ 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; + Q **mat; + Q *a; + struct oN tn; + int len,i,sgn,tsgn,row,k; + 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); + 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; 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() { GeoBucket g; @@ -1621,6 +1724,8 @@ GeoBucket create_bucket() return g; } +int length(NODE d); + void add_bucket(GeoBucket g,NODE d,int nv) { int l,k,m; @@ -1664,3 +1769,737 @@ DL remove_head_bucket(GeoBucket g,int nv) return d; } } + +/* DPV functions */ + +void adddv(VL vl,DPV p1,DPV p2,DPV *pr) +{ + int i,len; + DP *e; + + if ( !p1 || !p2 ) + error("adddv : invalid argument"); + else if ( p1->len != p2->len ) + error("adddv : size mismatch"); + else { + len = p1->len; + e = (DP *)MALLOC(p1->len*sizeof(DP)); + for ( i = 0; i < len; i++ ) + addd(vl,p1->body[i],p2->body[i],&e[i]); + MKDPV(len,e,*pr); + (*pr)->sugar = MAX(p1->sugar,p2->sugar); + } +} + +void subdv(VL vl,DPV p1,DPV p2,DPV *pr) +{ + int i,len; + DP *e; + + if ( !p1 || !p2 ) + error("subdv : invalid argument"); + else if ( p1->len != p2->len ) + error("subdv : size mismatch"); + else { + len = p1->len; + e = (DP *)MALLOC(p1->len*sizeof(DP)); + for ( i = 0; i < len; i++ ) + subd(vl,p1->body[i],p2->body[i],&e[i]); + MKDPV(len,e,*pr); + (*pr)->sugar = MAX(p1->sugar,p2->sugar); + } +} + +void chsgndv(DPV p1,DPV *pr) +{ + int i,len; + DP *e; + + if ( !p1 ) + error("subdv : invalid argument"); + else { + len = p1->len; + e = (DP *)MALLOC(p1->len*sizeof(DP)); + for ( i = 0; i < len; i++ ) + chsgnd(p1->body[i],&e[i]); + MKDPV(len,e,*pr); + (*pr)->sugar = p1->sugar; + } +} + +void muldv(VL vl,DP p1,DPV p2,DPV *pr) +{ + int i,len; + DP *e; + + len = p2->len; + e = (DP *)MALLOC(p2->len*sizeof(DP)); + if ( !p1 ) { + MKDPV(len,e,*pr); + (*pr)->sugar = 0; + } else { + for ( i = 0; i < len; i++ ) + muld(vl,p1,p2->body[i],&e[i]); + MKDPV(len,e,*pr); + (*pr)->sugar = p1->sugar + p2->sugar; + } +} + +int compdv(VL vl,DPV p1,DPV p2) +{ + int i,t,len; + + if ( p1->len != p2->len ) { + error("compdv : size mismatch"); + return 0; /* XXX */ + } else { + len = p1->len; + for ( i = 0; i < len; i++ ) + if ( (t = compd(vl,p1->body[i],p2->body[i])) ) + return t; + return 0; + } +} + +int ni_next(int *a,int n) +{ + int i,j,k,kj; + + /* find the first nonzero a[j] */ + for ( j = 0; j < n && a[j] == 0; j++ ); + /* find the first zero a[k] after a[j] */ + for ( k = j; k < n && a[k] == 1; k++ ); + if ( k == n ) return 0; + /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */ + /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */ + kj = k-j-1; + for ( i = 0; i < kj; i++ ) a[i] = 1; + for ( ; i < k; i++ ) a[i] = 0; + a[k] = 1; + return 1; +} + +int comp_nbm(NBM a,NBM b) +{ + int d,i,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; + } +} + +NBM mul_nbm(NBM a,NBM b) +{ + int ad,bd,d,i,j; + int *ab,*bb,*mb; + NBM m; + + ad = a->d; bd = b->d; ab = a->b; bb = b->b; + d = ad + bd; + NEWNBM(m); NEWNBMBDY(m,d); + 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); + else NBM_CLR(mb,j); + for ( i = 0; i < bd; i++, j++ ) + if ( NBM_GET(bb,i) ) NBM_SET(mb,j); + else NBM_CLR(mb,j); + return m; +} + +NBP nbmtonbp(NBM m) +{ + NODE n; + NBP u; + + MKNODE(n,m,0); + MKNBP(u,n); + return u; +} + +/* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */ + +P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) +{ + int i,d1; + NBM t; + + if ( !a->d ) error("separate_nbm : invalid argument"); + + if ( a0 ) { + 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 = (P)ONE; + if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0); + else NBM_CLR(t->b,0); + *ah = nbmtonbp(t); + } + + if ( ar ) { + d1 = a->d-1; + 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); + } + *ar = nbmtonbp(t); + } + + 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 = (P)ONE; + k1 = k-1; + for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i); + NBM_CLR(t->b,i); + return nbmtonbp(t); +} + +/* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */ + +P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) +{ + int i,d1,k,k1; + NBM t; + + if ( !a->d ) + error("separate_nbm : invalid argument"); + for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ ); + if ( i == a->d ) + error("separate_nbm : invalid argument"); + k1 = i; + k = i+1; + + if ( a0 ) { + 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 = (P)ONE; + for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i); + NBM_CLR(t->b,i); + *ah = nbmtonbp(t); + } + + if ( ar ) { + d1 = a->d-k; + 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); + } + *ar = nbmtonbp(t); + } + + return a->c; +} + +void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp); +void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp); +void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp); +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; + 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); + 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); + } + return u; +} + +NBP harmonic_mul_nbm(NBM a,NBM b) +{ + NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1; + P ac,bc,c; + + if ( !a->d || !b->d ) + u = nbmtonbp(mul_nbm(a,b)); + else { + mulp(CO,a->c,b->c,&c); + ac = separate_xky_nbm(a,&a0,&ah,&ar); + bc = separate_xky_nbm(b,&b0,&bh,&br); + 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); + harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1); + addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u); + } + return u; + +} + +void addnbp(VL vl,NBP p1,NBP p2, NBP *rp) +{ + NODE b1,b2,br=0,br0; + NBM m1,m2,m; + P c; + + if ( !p1 ) + *rp = p2; + else if ( !p2 ) + *rp = p1; + else { + for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) { + m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2); + switch ( comp_nbm(m1,m2) ) { + case 0: + 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; + BDY(br) = (pointer)m; + } + b1 = NEXT(b1); b2 = NEXT(b2); break; + case 1: + NEXTNODE(br0,br); BDY(br) = BDY(b1); + b1 = NEXT(b1); break; + case -1: + NEXTNODE(br0,br); BDY(br) = BDY(b2); + b2 = NEXT(b2); break; + } + } + if ( !br0 ) + if ( b1 ) + br0 = b1; + else if ( b2 ) + br0 = b2; + else { + *rp = 0; + return; + } + else if ( b1 ) + NEXT(br) = b1; + else if ( b2 ) + NEXT(br) = b2; + else + NEXT(br) = 0; + MKNBP(*rp,br0); + } +} + +void subnbp(VL vl,NBP p1,NBP p2, NBP *rp) +{ + NBP t; + + chsgnnbp(p2,&t); + addnbp(vl,p1,t,rp); +} + +void chsgnnbp(NBP p,NBP *rp) +{ + NODE r0,r=0,b; + NBM m,m1; + + 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; chsgnp(m->c,&m1->c); + BDY(r) = m1; + } + if ( r0 ) NEXT(r) = 0; + MKNBP(*rp,r0); +} + +void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) +{ + NODE b,n; + NBP r,t,s; + NBM m; + + if ( !p1 || !p2 ) { + *rp = 0; return; + } + if ( OID(p1) != O_NBP ) { + if ( !POLY(p1) ) + error("mulnbp : invalid argument"); + 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 ( !POLY(p2) ) + error("mulnbp : invalid argument"); + 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)) ) { + for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) { + mulnbmnbp(vl,(NBM)BDY(b),p2,&t); + addnbp(vl,r,t,&s); r = s; + } + *rp = r; + } else { + for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) { + mulnbpnbm(vl,p1,(NBM)BDY(b),&t); + addnbp(vl,r,t,&s); r = s; + } + *rp = r; + } +} + +void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp) +{ + NODE b,r0,r=0; + + if ( !p ) *rp = 0; + else { + for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) { + NEXTNODE(r0,r); + BDY(r) = mul_nbm(m,(NBM)BDY(b)); + } + if ( r0 ) NEXT(r) = 0; + MKNBP(*rp,r0); + } +} + +void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp) +{ + NODE b,r0,r=0; + + if ( !p ) *rp = 0; + else { + for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) { + NEXTNODE(r0,r); + BDY(r) = mul_nbm((NBM)BDY(b),m); + } + if ( r0 ) NEXT(r) = 0; + MKNBP(*rp,r0); + } +} + +void pwrnbp(VL vl,NBP a,Q q,NBP *c) +{ + int t; + NBP a1,a2; + N n1; + Q q1; + NBM m; + NODE r; + + if ( !q ) { + NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0; + MKNODE(r,m,0); MKNBP(*c,r); + } else if ( !a ) + *c = 0; + else if ( UNIQ(q) ) + *c = a; + else { + t = divin(NM(q),2,&n1); NTOQ(n1,1,q1); + pwrnbp(vl,a,q1,&a1); + mulnbp(vl,a1,a1,&a2); + if ( t ) + mulnbp(vl,a2,a,c); + else + *c = a2; + } +} + +int compnbp(VL vl,NBP p1,NBP p2) +{ + NODE n1,n2; + NBM m1,m2; + int t; + + if ( !p1 ) + return p2 ? -1 : 0; + else if ( !p2 ) + return 1; + else { + 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 = compp(CO,m1->c,m2->c) ) ) + return t; + } + if ( n1 ) + return 1; + else if ( n2 ) + return -1; + else + return 0; + } +} + +void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) +{ + NODE b1,b2,n; + NBP r,t,s; + NBM m; + + if ( !p1 || !p2 ) { + *rp = 0; return; + } + if ( OID(p1) != O_NBP ) { + if ( !POLY(p1) ) + error("shuffle_mulnbp : invalid argument"); + 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 ( !POLY(p2) ) + error("shuffle_mulnbp : invalid argument"); + 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) ) + for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) { + t = shuffle_mul_nbm(m,(NBM)BDY(b2)); + addnbp(vl,r,t,&s); r = s; + } + *rp = r; +} + +void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) +{ + NODE b1,b2,n; + NBP r,t,s; + NBM m; + + if ( !p1 || !p2 ) { + *rp = 0; return; + } + if ( OID(p1) != O_NBP ) { + if ( !POLY(p1) ) + error("harmonic_mulnbp : invalid argument"); + 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 ( !POLY(p2) ) + error("harmonic_mulnbp : invalid argument"); + 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) ) + for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) { + t = harmonic_mul_nbm(m,(NBM)BDY(b2)); + addnbp(vl,r,t,&s); r = s; + } + *rp = r; +} + +#if 0 +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; + P c,c1; + NODE r,t,t1,p; + NBP u; + + ad = a->d; bd = b->d; ab = a->b; bb = b->b; + d = ad + bd; + w = (int *)ALLOCA(d*sizeof(int)); + NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b; + for ( i = 0; i < ad; i++ ) w[i] = 1; + for ( ; i < d; i++ ) w[i] = 0; + mulp(CO,a->c,b->c,&c); + r = 0; + do { + wm->d = d; wm->c = c; + ai = 0; bi = 0; + for ( i = 0; i < d; i++ ) { + if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; } + else { bit = NBM_GET(bb,bi); bi++; } + if ( bit ) NBM_SET(wmb,i); + else NBM_CLR(wmb,i); + } + for ( p = 0, t = r; t; p = t, t = NEXT(t) ) { + tm = (NBM)BDY(t); + s = comp_nbm(tm,wm); + if ( s < 0 ) { + /* insert */ + MKNODE(t1,wm,t); + if ( !p ) r = t1; + else NEXT(p) = t1; + NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b; + break; + } else if ( s == 0 ) { + /* add coefs */ + addp(CO,tm->c,c,&c1); + if ( c1 ) tm->c = c1; + else NEXT(p) = NEXT(t); + break; + } + } + if ( !t ) { + /* append */ + MKNODE(t1,wm,t); + if ( !p ) r = t1; + else NEXT(p) = t1; + NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b; + } + } while ( ni_next(w,d) ); + MKNBP(u,r); + return u; +} + +int nbmtoxky(NBM a,int *b) +{ + int d,i,j,k; + int *p; + + d = a->d; p = a->b; + for ( i = j = 0, k = 1; i < d; i++ ) { + if ( !NBM_GET(p,i) ) { + b[j++] = k; + k = 1; + } else k++; + } + return j; +} + +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; + P c,c1; + NBM wm,tm; + NODE r,t1,t,p; + NBP u; + + da = a->d; db = b->d; d = da+db; + wa = (int *)ALLOCA(da*sizeof(int)); + wb = (int *)ALLOCA(db*sizeof(int)); + la = nbmtoxky(a,wa); + lb = nbmtoxky(b,wb); + 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; + lmin = la>lb?la:lb; + w = (int *)ALLOCA(lmax*sizeof(int)); + /* position of a+b */ + wab = (int *)ALLOCA(lmax*sizeof(int)); + /* position of a */ + wa1 = (int *)ALLOCA(lmax*sizeof(int)); + NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b; + for ( l = lmin, r = 0; l <= lmax; l++ ) { + lab = lmax - l; + la1 = la - lab; + lb1 = lb - lab; + lab1 = l-lab; + /* partion l into three parts: a, b, a+b */ + /* initialize wab */ + for ( i = 0; i < lab; i++ ) wab[i] = 1; + for ( ; i < l; i++ ) wab[i] = 0; + do { + /* initialize wa1 */ + for ( i = 0; i < la1; i++ ) wa1[i] = 1; + for ( ; i < lab1; i++ ) wa1[i] = 0; + do { + ia = 0; ib = 0; + for ( i = j = 0; i < l; i++ ) + if ( wab[i] ) w[i] = wa[ia++]+wb[ib++]; + else if ( wa1[j++] ) w[i] = wa[ia++]; + else w[i] = wb[ib++]; + for ( i = j = 0; i < l; i++ ) { + for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j); + NBM_CLR(wmb,j); j++; + } + wm->d = j; wm->c = c; + for ( p = 0, t = r; t; p = t, t = NEXT(t) ) { + tm = (NBM)BDY(t); + s = comp_nbm(tm,wm); + if ( s < 0 ) { + /* insert */ + MKNODE(t1,wm,t); + if ( !p ) r = t1; + else NEXT(p) = t1; + NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b; + break; + } else if ( s == 0 ) { + /* add coefs */ + addp(CO,tm->c,c,&c1); + if ( c1 ) tm->c = c1; + else NEXT(p) = NEXT(t); + break; + } + } + if ( !t ) { + /* append */ + MKNODE(t1,wm,t); + if ( !p ) r = t1; + else NEXT(p) = t1; + NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b; + } + } while ( ni_next(wa1,lab1) ); + } while ( ni_next(wab,l) ); + } + MKNBP(u,r); + return u; +} +#endif