=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/distm.c,v retrieving revision 1.13 retrieving revision 1.20 diff -u -p -r1.13 -r1.20 --- OpenXM_contrib2/asir2000/engine/distm.c 2003/07/21 01:57:51 1.13 +++ OpenXM_contrib2/asir2000/engine/distm.c 2017/08/31 02:36:21 1.20 @@ -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/distm.c,v 1.12 2003/07/20 08:55:23 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.19 2012/12/17 07:20:44 noro Exp $ */ #include "ca.h" #include "inline.h" @@ -78,7 +78,7 @@ void mptomd(VL vl,int mod,VL dvl,P p,DP *pr) for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ ); if ( NUM(p) ) { NEWDL(d,n); - NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr); + NEWMP(m); m->dl = d; C(m) = (Obj)p; NEXT(m) = 0; MKDP(n,m,*pr); } else { for ( i = 0, tvl = dvl, v = VR(p); tvl && tvl->v != v; tvl = NEXT(tvl), i++ ); @@ -93,7 +93,7 @@ void mptomd(VL vl,int mod,VL dvl,P p,DP *pr) mptomd(vl,mod,dvl,COEF(dc),&t); NEWDL(d,n); d->d[i] = QTOS(DEG(dc)); d->td = MUL_WEIGHT(d->d[i],i); - NEWMP(m); m->dl = d; C(m) = (P)ONEM; NEXT(m) = 0; MKDP(n,m,u); + NEWMP(m); m->dl = d; C(m) = (Obj)ONEM; NEXT(m) = 0; MKDP(n,m,u); comm_mulmd(vl,mod,t,u,&r); addmd(vl,mod,r,s,&t); s = t; } *pr = s; @@ -102,6 +102,42 @@ void mptomd(VL vl,int mod,VL dvl,P p,DP *pr) } } +void mdtodp(DP p,DP *pr) +{ + MP m,mr0,mr; + + if ( !p ) + *pr = 0; + else { + for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) { + NEXTMP(mr0,mr); mr->dl = m->dl; + mptop((P)C(m),(P *)&C(mr)); + } + NEXT(mr) = 0; + MKDP(NV(p),mr0,*pr); + (*pr)->sugar = p->sugar; + } +} + +void _mdtodp(DP p,DP *pr) +{ + MP m,mr0,mr; + int i; + Q q; + + if ( !p ) + *pr = 0; + else { + for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) { + NEXTMP(mr0,mr); mr->dl = m->dl; + i = ITOS(m->c); STOQ(i,q); C(mr) = (Obj)q; + } + NEXT(mr) = 0; + MKDP(NV(p),mr0,*pr); + (*pr)->sugar = p->sugar; + } +} + void mdtop(VL vl,int mod,VL dvl,DP p,P *pr) { int n,i; @@ -115,7 +151,7 @@ void mdtop(VL vl,int mod,VL dvl,DP p,P *pr) *pr = 0; else { for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) { - for ( i = 0, t = C(m), d = m->dl, tvl = dvl; + for ( i = 0, t = (P)C(m), d = m->dl, tvl = dvl; i < n; tvl = NEXT(tvl), i++ ) { MKMV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u); mulmp(vl,mod,t,u,&w); t = w; @@ -151,9 +187,9 @@ void addmd(VL vl,int mod,DP p1,DP p2,DP *pr) } else t = 0; } else - addmp(vl,mod,C(m1),C(m2),&t); + addmp(vl,mod,(P)C(m1),(P)C(m2),&t); if ( t ) { - NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t; + NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = (Obj)t; } m1 = NEXT(m1); m2 = NEXT(m2); break; case 1: @@ -203,7 +239,7 @@ void chsgnmd(int mod,DP p,DP *pr) *pr = 0; else { for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { - NEXTMP(mr0,mr); chsgnmp(mod,C(m),&C(mr)); mr->dl = m->dl; + NEXTMP(mr0,mr); chsgnmp(mod,(P)C(m),(P *)&C(mr)); mr->dl = m->dl; } NEXT(mr) = 0; MKDP(NV(p),mr0,*pr); if ( *pr ) @@ -241,7 +277,7 @@ void comm_mulmd(VL vl,int mod,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; } @@ -272,7 +308,7 @@ void weyl_mulmd(VL vl,int mod,DP p1,DP p2,DP *pr) else { for ( m = BDY(p2), 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; } @@ -297,14 +333,14 @@ void mulmdm(VL vl,int mod,DP p,MP m0,DP *pr) if ( !p ) *pr = 0; else { - for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p); + for ( mr0 = 0, m = BDY(p), c = (P)C(m0), d = m0->dl, n = NV(p); m; m = NEXT(m) ) { NEXTMP(mr0,mr); if ( NUM(C(m)) && NUM(c) ) { t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod); - STOMQ(t,q); C(mr) = (P)q; + STOMQ(t,q); C(mr) = (Obj)q; } else - mulmp(vl,mod,C(m),c,&C(mr)); + mulmp(vl,mod,(P)C(m),c,(P *)&C(mr)); adddl(n,m->dl,d,&mr->dl); } NEXT(mr) = 0; MKDP(NV(p),mr0,*pr); @@ -326,7 +362,7 @@ void weyl_mulmdm(VL vl,int mod,DP p,MP m0,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; } @@ -359,7 +395,7 @@ void weyl_mulmmm(VL vl,int mod,MP m0,MP m1,int n,DP *p if ( !m0 || !m1 ) *pr = 0; else { - c0 = C(m0); c1 = C(m1); + c0 = (P)C(m0); c1 = (P)C(m1); mulmp(vl,mod,c0,c1,&c); d0 = m0->dl; d1 = m1->dl; n2 = n>>1; @@ -368,7 +404,7 @@ void weyl_mulmmm(VL vl,int mod,MP m0,MP m1,int n,DP *p /* offset of h-degree */ NEWDL(d,n); d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1]; - NEWMP(mr); mr->c = (P)ONEM; mr->dl = d; NEXT(mr) = 0; + NEWMP(mr); mr->c = (Obj)ONEM; mr->dl = d; NEXT(mr) = 0; MKDP(n,mr,r); r->sugar = d->td; } else r = (DP)ONEM; @@ -385,7 +421,7 @@ void weyl_mulmmm(VL vl,int mod,MP m0,MP m1,int n,DP *p min = MIN(k,l); if ( min+1 > tablen ) { - if ( tab ) GC_free(tab); + if ( tab ) GCFREE(tab); tab = (int *)MALLOC((min+1)*sizeof(int)); tablen = min+1; } @@ -396,7 +432,7 @@ void weyl_mulmmm(VL vl,int mod,MP m0,MP m1,int n,DP *p d->d[i] = a-j; d->d[n2+i] = b-j; d->td = s; d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i)); - STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d; + STOMQ(tab[j],mq); mr->c = (Obj)mq; mr->dl = d; } else for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) { @@ -404,7 +440,7 @@ void weyl_mulmmm(VL vl,int mod,MP m0,MP m1,int n,DP *p d->d[i] = a-j; d->d[n2+i] = b-j; d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */ s = MAX(s,d->td); /* XXX */ - STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d; + STOMQ(tab[j],mq); mr->c = (Obj)mq; mr->dl = d; } bzero(tab,(min+1)*sizeof(int)); if ( mr0 ) @@ -431,9 +467,9 @@ void mulmdc(VL vl,int mod,DP p,P c,DP *pr) NEXTMP(mr0,mr); if ( NUM(C(m)) && NUM(c) ) { t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod); - STOMQ(t,q); C(mr) = (P)q; + STOMQ(t,q); C(mr) = (Obj)q; } else - mulmp(vl,mod,C(m),c,&C(mr)); + mulmp(vl,mod,(P)C(m),c,(P *)&C(mr)); mr->dl = m->dl; } NEXT(mr) = 0; MKDP(NV(p),mr0,*pr); @@ -452,7 +488,7 @@ void divsmdc(VL vl,int mod,DP p,P c,DP *pr) *pr = 0; else { for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { - NEXTMP(mr0,mr); divsmp(vl,mod,C(m),c,&C(mr)); mr->dl = m->dl; + NEXTMP(mr0,mr); divsmp(vl,mod,(P)C(m),c,(P *)&C(mr)); mr->dl = m->dl; } NEXT(mr) = 0; MKDP(NV(p),mr0,*pr); if ( *pr ) @@ -495,13 +531,13 @@ void _dp_mod(DP p,int mod,NODE subst,DP *rp) *rp = 0; else { for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { - for ( tn = subst, s = m->c; tn; tn = NEXT(NEXT(tn)) ) { + for ( tn = subst, s = (P)m->c; tn; tn = NEXT(NEXT(tn)) ) { substp(CO,s,BDY(tn),BDY(NEXT(tn)),&t); s = t; } ptomp(mod,s,&t); if ( t ) { - NEXTMP(mr0,mr); mr->c = STOI(CONT((MQ)t)); mr->dl = m->dl; + NEXTMP(mr0,mr); mr->c = (Obj)STOI(CONT((MQ)t)); mr->dl = m->dl; } } if ( mr0 ) { @@ -522,7 +558,7 @@ void _dp_monic(DP p,int mod,DP *rp) c = invm(ITOS(BDY(p)->c),mod); for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { c1 = dmar(ITOS(m->c),c,0,mod); - NEXTMP(mr0,mr); mr->c = STOI(c1); mr->dl = m->dl; + NEXTMP(mr0,mr); mr->c = (Obj)STOI(c1); mr->dl = m->dl; } NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar; } @@ -569,7 +605,7 @@ void addmd_destructive(int mod,DP p1,DP p2,DP *pr) t += mod; s = m1; m1 = NEXT(m1); if ( t ) { - NEXTMP2(mr0,mr,s); C(mr) = STOI(t); + NEXTMP2(mr0,mr,s); C(mr) = (Obj)STOI(t); } s = m2; m2 = NEXT(m2); break; @@ -627,7 +663,7 @@ void comm_mulmd_dup(int mod,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; } @@ -654,7 +690,7 @@ void weyl_mulmd_dup(int mod,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; } @@ -681,7 +717,7 @@ void mulmdm_dup(int mod,DP p,MP m0,DP *pr) for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p); m; m = NEXT(m) ) { NEXTMP(mr0,mr); - C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod)); + C(mr) = (Obj)STOI(dmar(ITOS(C(m)),c,0,mod)); NEWDL_NOINIT(dt,n); mr->dl = dt; dm = m->dl; dt->td = d->td + dm->td; @@ -710,7 +746,7 @@ void weyl_mulmdm_dup(int mod,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; } @@ -722,8 +758,8 @@ void weyl_mulmdm_dup(int mod,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 cdlm *)MALLOC(rtlen*sizeof(struct cdlm)); psum = (MP *)MALLOC(rtlen*sizeof(MP)); @@ -735,7 +771,7 @@ void weyl_mulmdm_dup(int mod,MP m0,DP p,DP *pr) for ( j = 0; j < tlen; j++ ) { if ( tab[j].c ) { NEWMP(m); m->dl = tab[j].d; - C(m) = STOI(tab[j].c); NEXT(m) = psum[j]; + C(m) = (Obj)STOI(tab[j].c); NEXT(m) = psum[j]; psum[j] = m; } } @@ -785,7 +821,7 @@ void weyl_mulmmm_dup(int mod,MP m0,MP m1,int n,struct rtab[0].d = d; if ( rtablen > tmptablen ) { - if ( tmptab ) GC_free(tmptab); + if ( tmptab ) GCFREE(tmptab); tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm)); tmptablen = rtablen; } @@ -812,8 +848,8 @@ void weyl_mulmmm_dup(int mod,MP m0,MP m1,int n,struct 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 cdlm *)MALLOC(tablen*sizeof(struct cdlm)); ctab = (int *)MALLOC(tablen*sizeof(int)); @@ -884,1429 +920,4 @@ void adddl_dup(int n,DL d1,DL d2,DL *dr) dt->td = d1->td + d2->td; for ( i = 0; i < n; i++ ) dt->d[i] = d1->d[i]+d2->d[i]; -} - -/* ------------------------------------------------------ */ - -#if defined(__GNUC__) -#define INLINE inline -#elif defined(VISUAL) -#define INLINE __inline -#else -#define INLINE -#endif - -typedef struct oPGeoBucket { - int m; - struct oND *body[32]; -} *PGeoBucket; - -typedef struct oND { - struct oNM *body; - int nv; - int sugar; -} *ND; - -typedef struct oNM { - struct oNM *next; - int td; - int c; - unsigned int dl[1]; -} *NM; - -typedef struct oND_pairs { - struct oND_pairs *next; - int i1,i2; - int td,sugar; - unsigned int lcm[1]; -} *ND_pairs; - -static ND *nps; -int nd_mod,nd_nvar; -int is_rlex; -int nd_epw,nd_bpe,nd_wpd; -unsigned int nd_mask[32]; -unsigned int nd_mask0; -NM _nm_free_list; -ND _nd_free_list; -ND_pairs _ndp_free_list; - -extern int Top,Reverse; -int nd_psn,nd_pslen; - -void GC_gcollect(); -NODE append_one(NODE,int); - -#define HTD(d) ((d)->body->td) -#define HDL(d) ((d)->body->dl) -#define HC(d) ((d)->body->c) - -#define NEWND_pairs(m) if(!_ndp_free_list)_NDP_alloc(); (m)=_ndp_free_list; _ndp_free_list = NEXT(_ndp_free_list) -#define NEWNM(m) if(!_nm_free_list)_NM_alloc(); (m)=_nm_free_list; _nm_free_list = NEXT(_nm_free_list) -#define MKND(n,m,d) if(!_nd_free_list)_ND_alloc(); (d)=_nd_free_list; _nd_free_list = (ND)BDY(_nd_free_list); (d)->nv=(n); BDY(d)=(m) - -#define NEXTNM(r,c) \ -if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEXT(c);} -#define NEXTNM2(r,c,s) \ -if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);} -#define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m) -#define FREEND(m) BDY(m)=(NM)_nd_free_list; _nd_free_list=(m) - -#define NEXTND_pairs(r,c) \ -if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);} - -ND_pairs crit_B( ND_pairs d, int s ); -void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp); -NODE nd_setup(NODE f); -int nd_newps(ND a); -ND_pairs nd_minp( ND_pairs d, ND_pairs *prest ); -NODE update_base(NODE nd,int ndp); -static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest ); -int crit_2( int dp1, int dp2 ); -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(NODE f); -void nd_free_private_storage(); -void _NM_alloc(); -void _ND_alloc(); -int ndl_td(unsigned int *d); -ND nd_add(ND p1,ND p2); -ND nd_mul_nm(ND p,NM m0); -ND nd_mul_term(ND p,int td,unsigned int *d); -ND nd_sp(ND_pairs p); -ND nd_find_reducer(ND g); -ND nd_nf(ND g,int full); -void ndl_print(unsigned int *dl); -void nd_print(ND p); -void ndp_print(ND_pairs d); -int nd_length(ND p); -void nd_monic(ND p); -void nd_mul_c(ND p,int mul); - -void nd_free_private_storage() -{ - _nd_free_list = 0; - _nm_free_list = 0; - GC_gcollect(); -} - -void _NM_alloc() -{ - NM p; - int i; - - for ( i = 0; i < 10240; i++ ) { - p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int)); - p->next = _nm_free_list; _nm_free_list = p; - } -} - -void _ND_alloc() -{ - ND p; - int i; - - for ( i = 0; i < 1024; i++ ) { - p = (ND)GC_malloc(sizeof(struct oND)); - p->body = (NM)_nd_free_list; _nd_free_list = p; - } -} - -void _NDP_alloc() -{ - ND_pairs p; - int i; - - for ( i = 0; i < 10240; i++ ) { - p = (ND_pairs)GC_malloc(sizeof(struct oND_pairs) - +(nd_wpd-1)*sizeof(unsigned int)); - p->next = _ndp_free_list; _ndp_free_list = p; - } -} - -INLINE nd_length(ND p) -{ - NM m; - int i; - - if ( !p ) - return 0; - else { - for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ ); - return i; - } -} - -INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2) -{ - unsigned int u1,u2; - int i,j; - - switch ( nd_bpe ) { - case 4: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - if ( (u1&0xf0000000) < (u2&0xf0000000) ) return 0; - if ( (u1&0xf000000) < (u2&0xf000000) ) return 0; - if ( (u1&0xf00000) < (u2&0xf00000) ) return 0; - if ( (u1&0xf0000) < (u2&0xf0000) ) return 0; - if ( (u1&0xf000) < (u2&0xf000) ) return 0; - if ( (u1&0xf00) < (u2&0xf00) ) return 0; - if ( (u1&0xf0) < (u2&0xf0) ) return 0; - if ( (u1&0xf) < (u2&0xf) ) return 0; - } - return 1; - break; - case 6: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - if ( (u1&0x3f000000) < (u2&0x3f000000) ) return 0; - if ( (u1&0xfc0000) < (u2&0xfc0000) ) return 0; - if ( (u1&0x3f000) < (u2&0x3f000) ) return 0; - if ( (u1&0xfc0) < (u2&0xfc0) ) return 0; - if ( (u1&0x3f) < (u2&0x3f) ) return 0; - } - return 1; - break; - case 8: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - if ( (u1&0xff000000) < (u2&0xff000000) ) return 0; - if ( (u1&0xff0000) < (u2&0xff0000) ) return 0; - if ( (u1&0xff00) < (u2&0xff00) ) return 0; - if ( (u1&0xff) < (u2&0xff) ) return 0; - } - return 1; - break; - case 16: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - if ( (u1&0xffff0000) < (u2&0xffff0000) ) return 0; - if ( (u1&0xffff) < (u2&0xffff) ) return 0; - } - return 1; - break; - case 32: - for ( i = 0; i < nd_wpd; i++ ) - if ( d1[i] < d2[i] ) return 0; - return 1; - break; - default: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - for ( j = 0; j < nd_epw; j++ ) - if ( (u1&nd_mask[j]) < (u2&nd_mask[j]) ) return 0; - } - return 1; - } -} - -void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d) -{ - unsigned int t1,t2,u,u1,u2; - int i,j; - - switch ( nd_bpe ) { - case 4: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - t1 = (u1&0xf0000000); t2 = (u2&0xf0000000); u = t1>t2?t1:t2; - t1 = (u1&0xf000000); t2 = (u2&0xf000000); u |= t1>t2?t1:t2; - t1 = (u1&0xf00000); t2 = (u2&0xf00000); u |= t1>t2?t1:t2; - t1 = (u1&0xf0000); t2 = (u2&0xf0000); u |= t1>t2?t1:t2; - t1 = (u1&0xf000); t2 = (u2&0xf000); u |= t1>t2?t1:t2; - t1 = (u1&0xf00); t2 = (u2&0xf00); u |= t1>t2?t1:t2; - t1 = (u1&0xf0); t2 = (u2&0xf0); u |= t1>t2?t1:t2; - t1 = (u1&0xf); t2 = (u2&0xf); u |= t1>t2?t1:t2; - d[i] = u; - } - break; - case 6: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - t1 = (u1&0x3f000000); t2 = (u2&0x3f000000); u = t1>t2?t1:t2; - t1 = (u1&0xfc0000); t2 = (u2&0xfc0000); u |= t1>t2?t1:t2; - t1 = (u1&0x3f000); t2 = (u2&0x3f000); u |= t1>t2?t1:t2; - t1 = (u1&0xfc0); t2 = (u2&0xfc0); u |= t1>t2?t1:t2; - t1 = (u1&0x3f); t2 = (u2&0x3f); u |= t1>t2?t1:t2; - d[i] = u; - } - break; - case 8: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - t1 = (u1&0xff000000); t2 = (u2&0xff000000); u = t1>t2?t1:t2; - t1 = (u1&0xff0000); t2 = (u2&0xff0000); u |= t1>t2?t1:t2; - t1 = (u1&0xff00); t2 = (u2&0xff00); u |= t1>t2?t1:t2; - t1 = (u1&0xff); t2 = (u2&0xff); u |= t1>t2?t1:t2; - d[i] = u; - } - break; - case 16: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - t1 = (u1&0xffff0000); t2 = (u2&0xffff0000); u = t1>t2?t1:t2; - t1 = (u1&0xffff); t2 = (u2&0xffff); u |= t1>t2?t1:t2; - d[i] = u; - } - break; - case 32: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - d[i] = u1>u2?u1:u2; - } - break; - default: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - for ( j = 0, u = 0; j < nd_epw; j++ ) { - t1 = (u1&nd_mask[j]); t2 = (u2&nd_mask[j]); u |= t1>t2?t1:t2; - } - d[i] = u; - } - break; - } -} - -int ndl_td(unsigned int *d) -{ - unsigned int t,u; - int i,j; - - for ( t = 0, i = 0; i < nd_wpd; i++ ) { - u = d[i]; - for ( j = 0; j < nd_epw; j++, u>>=nd_bpe ) - t += (u&nd_mask0); - } - return t; -} - -INLINE int ndl_compare(unsigned int *d1,unsigned int *d2) -{ - int i; - - for ( i = 0; i < nd_wpd; i++, d1++, d2++ ) - if ( *d1 > *d2 ) - return is_rlex ? -1 : 1; - else if ( *d1 < *d2 ) - return is_rlex ? 1 : -1; - return 0; -} - -INLINE int ndl_equal(unsigned int *d1,unsigned int *d2) -{ - int i; - - for ( i = 0; i < nd_wpd; i++ ) - if ( d1[i] != d2[i] ) - return 0; - return 1; -} - -INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d) -{ - int i; - - for ( i = 0; i < nd_wpd; i++ ) - d[i] = d1[i]+d2[i]; -} - -void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d) -{ - int i; - - for ( i = 0; i < nd_wpd; i++ ) - d[i] = d1[i]-d2[i]; -} - -int ndl_disjoint(unsigned int *d1,unsigned int *d2) -{ - unsigned int t1,t2,u,u1,u2; - int i,j; - - switch ( nd_bpe ) { - case 4: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - t1 = u1&0xf0000000; t2 = u2&0xf0000000; if ( t1&&t2 ) return 0; - t1 = u1&0xf000000; t2 = u2&0xf000000; if ( t1&&t2 ) return 0; - t1 = u1&0xf00000; t2 = u2&0xf00000; if ( t1&&t2 ) return 0; - t1 = u1&0xf0000; t2 = u2&0xf0000; if ( t1&&t2 ) return 0; - t1 = u1&0xf000; t2 = u2&0xf000; if ( t1&&t2 ) return 0; - t1 = u1&0xf00; t2 = u2&0xf00; if ( t1&&t2 ) return 0; - t1 = u1&0xf0; t2 = u2&0xf0; if ( t1&&t2 ) return 0; - t1 = u1&0xf; t2 = u2&0xf; if ( t1&&t2 ) return 0; - } - return 1; - break; - case 6: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - t1 = u1&0x3f000000; t2 = u2&0x3f000000; if ( t1&&t2 ) return 0; - t1 = u1&0xfc0000; t2 = u2&0xfc0000; if ( t1&&t2 ) return 0; - t1 = u1&0x3f000; t2 = u2&0x3f000; if ( t1&&t2 ) return 0; - t1 = u1&0xfc0; t2 = u2&0xfc0; if ( t1&&t2 ) return 0; - t1 = u1&0x3f; t2 = u2&0x3f; if ( t1&&t2 ) return 0; - } - return 1; - break; - case 8: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - t1 = u1&0xff000000; t2 = u2&0xff000000; if ( t1&&t2 ) return 0; - t1 = u1&0xff0000; t2 = u2&0xff0000; if ( t1&&t2 ) return 0; - t1 = u1&0xff00; t2 = u2&0xff00; if ( t1&&t2 ) return 0; - t1 = u1&0xff; t2 = u2&0xff; if ( t1&&t2 ) return 0; - } - return 1; - break; - case 16: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - t1 = u1&0xffff0000; t2 = u2&0xffff0000; if ( t1&&t2 ) return 0; - t1 = u1&0xffff; t2 = u2&0xffff; if ( t1&&t2 ) return 0; - } - return 1; - break; - case 32: - for ( i = 0; i < nd_wpd; i++ ) - if ( d1[i] && d2[i] ) return 0; - return 1; - break; - default: - for ( i = 0; i < nd_wpd; i++ ) { - u1 = d1[i]; u2 = d2[i]; - for ( j = 0; j < nd_epw; j++ ) { - if ( (u1&nd_mask0) && (u2&nd_mask0) ) return 0; - u1 >>= nd_bpe; u2 >>= nd_bpe; - } - } - return 1; - break; - } -} - -ND nd_reduce(ND p1,ND p2) -{ - int c,c1,c2,t,td,td2,mul; - NM m2,prev,head,cur,new; - unsigned int *d; - - if ( !p1 ) - return 0; - else { - c2 = invm(HC(p2),nd_mod); - c1 = nd_mod-HC(p1); - DMAR(c1,c2,0,nd_mod,mul); - td = HTD(p1)-HTD(p2); - d = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int)); - ndl_sub(HDL(p1),HDL(p2),d); - prev = 0; head = cur = BDY(p1); - NEWNM(new); - for ( m2 = BDY(p2); m2; ) { - td2 = new->td = m2->td+td; - ndl_add(m2->dl,d,new->dl); - if ( !cur ) { - c1 = C(m2); - DMAR(c1,mul,0,nd_mod,c2); - C(new) = c2; - if ( !prev ) { - prev = new; - NEXT(prev) = 0; - head = prev; - } else { - NEXT(prev) = new; - NEXT(new) = 0; - prev = new; - } - m2 = NEXT(m2); - NEWNM(new); - continue; - } - if ( cur->td > td2 ) - c = 1; - else if ( cur->td < td2 ) - c = -1; - else - c = ndl_compare(cur->dl,new->dl); - switch ( c ) { - case 0: - c2 = C(m2); - c1 = C(cur); - DMAR(c2,mul,c1,nd_mod,t); - if ( t ) - C(cur) = t; - else if ( !prev ) { - head = NEXT(cur); - FREENM(cur); - cur = head; - } else { - NEXT(prev) = NEXT(cur); - FREENM(cur); - cur = NEXT(prev); - } - m2 = NEXT(m2); - break; - case 1: - prev = cur; - cur = NEXT(cur); - break; - case -1: - if ( !prev ) { - /* cur = head */ - prev = new; - c2 = C(m2); - DMAR(c2,mul,0,nd_mod,c1); - C(prev) = c1; - NEXT(prev) = head; - head = prev; - } else { - c2 = C(m2); - DMAR(c2,mul,0,nd_mod,c1); - C(new) = c1; - NEXT(prev) = new; - NEXT(new) = cur; - prev = new; - } - NEWNM(new); - m2 = NEXT(m2); - break; - } - } - FREENM(new); - if ( head ) { - BDY(p1) = head; - p1->sugar = MAX(p1->sugar,p2->sugar+td); - return p1; - } else { - FREEND(p1); - return 0; - } - - } -} - -ND nd_sp(ND_pairs p) -{ - NM m; - ND p1,p2,t1,t2; - unsigned int *lcm; - int td; - - p1 = nps[p->i1]; - p2 = nps[p->i2]; - lcm = p->lcm; - td = p->td; - NEWNM(m); - C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl); NEXT(m) = 0; - t1 = nd_mul_nm(p1,m); - C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl); - t2 = nd_mul_nm(p2,m); - FREENM(m); - return nd_add(t1,t2); -} - -ND nd_find_reducer(ND g) -{ - NM m; - ND r,p; - int i,c1,c2,c; - - for ( i = 0; i < nd_psn; i++ ) { - p = nps[i]; - if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) { -#if 1 - NEWNM(m); - c1 = invm(HC(p),nd_mod); - c2 = nd_mod-HC(g); - DMAR(c1,c2,0,nd_mod,c); - C(m) = c; - m->td = HTD(g)-HTD(p); - ndl_sub(HDL(g),HDL(p),m->dl); - NEXT(m) = 0; - r = nd_mul_nm(p,m); - FREENM(m); - r->sugar = m->td + p->sugar; - return r; -#else - return p; -#endif - } - } - return 0; -} - -ND nd_find_monic_reducer(ND g) -{ - int *d; - ND p,r; - int i; - - for ( i = 0; i < nd_psn; i++ ) { - p = nps[i]; - if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) { - d = (int *)ALLOCA(nd_wpd*sizeof(int)); - ndl_sub(HDL(g),HDL(p),d); - r = nd_mul_term(p,HTD(g)-HTD(p),d); - return r; - } - } - return 0; -} - -ND nd_add(ND p1,ND p2) -{ - int n,c; - int t; - ND r; - NM m1,m2,mr0,mr,s; - - if ( !p1 ) - return p2; - else if ( !p2 ) - return p1; - else { - for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) { - if ( m1->td > m2->td ) - c = 1; - else if ( m1->td < m2->td ) - c = -1; - else - c = ndl_compare(m1->dl,m2->dl); - switch ( c ) { - case 0: - t = ((C(m1))+(C(m2))) - nd_mod; - if ( t < 0 ) - t += nd_mod; - s = m1; m1 = NEXT(m1); - if ( t ) { - NEXTNM2(mr0,mr,s); C(mr) = (t); - } else { - FREENM(s); - } - s = m2; m2 = NEXT(m2); FREENM(s); - break; - case 1: - s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s); - break; - case -1: - s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s); - break; - } - } - if ( !mr0 ) - if ( m1 ) - mr0 = m1; - else if ( m2 ) - mr0 = m2; - else - return 0; - else if ( m1 ) - NEXT(mr) = m1; - else if ( m2 ) - NEXT(mr) = m2; - else - NEXT(mr) = 0; - BDY(p1) = mr0; - p1->sugar = MAX(p1->sugar,p2->sugar); - FREEND(p2); - return p1; - } -} - -ND nd_mul_nm(ND p,NM m0) -{ - NM m,mr,mr0; - unsigned int *d,*dt,*dm; - int c,n,td,i,c1,c2; - int *pt,*p1,*p2; - ND r; - - if ( !p ) - return 0; - else { - n = NV(p); m = BDY(p); - d = m0->dl; td = m0->td; c = C(m0); - mr0 = 0; - for ( ; m; m = NEXT(m) ) { - NEXTNM(mr0,mr); - c1 = C(m); - DMAR(c1,c,0,nd_mod,c2); - C(mr) = c2; - mr->td = m->td+td; - ndl_add(m->dl,d,mr->dl); - } - NEXT(mr) = 0; - MKND(NV(p),mr0,r); - r->sugar = p->sugar + td; - return r; - } -} - -ND nd_mul_term(ND p,int td,unsigned int *d) -{ - NM m,mr,mr0; - int c,n; - ND r; - - if ( !p ) - return 0; - else { - n = NV(p); m = BDY(p); - mr0 = 0; - for ( ; m; m = NEXT(m) ) { - NEXTNM(mr0,mr); - C(mr) = C(m); - mr->td = m->td+td; - ndl_add(m->dl,d,mr->dl); - } - NEXT(mr) = 0; - MKND(NV(p),mr0,r); - r->sugar = p->sugar + td; - return r; - } -} - -#if 1 -ND nd_nf(ND g,int full) -{ - ND p,d,red; - NM m,mrd,tail; - int n,sugar,psugar; - - if ( !g ) - return 0; - sugar = g->sugar; - n = NV(g); - for ( d = 0; g; ) { - red = nd_find_reducer(g); - if ( red ) { -#if 1 - g = nd_add(g,red); - sugar = MAX(sugar,red->sugar); -#else - psugar = (HTD(g)-HTD(red))+red->sugar; - g = nd_reduce(g,red); - sugar = MAX(sugar,psugar); -#endif - } else if ( !full ) - return g; - 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 ) - d->sugar = sugar; - return d; -} -#else - -ND nd_remove_head(ND p) -{ - NM m; - - m = BDY(p); - if ( !NEXT(m) ) { - FREEND(p); - p = 0; - } else - BDY(p) = NEXT(m); - FREENM(m); - return p; -} - -PGeoBucket create_pbucket() -{ - PGeoBucket g; - - g = CALLOC(1,sizeof(struct oPGeoBucket)); - g->m = -1; - return g; -} - -void add_pbucket(PGeoBucket g,ND d) -{ - int l,k,m; - - l = nd_length(d); - for ( k = 0, m = 1; l > m; k++, m <<= 2 ); - /* 4^(k-1) < l <= 4^k */ - d = nd_add(g->body[k],d); - for ( ; d && nd_length(d) > 1<<(2*k); k++ ) { - g->body[k] = 0; - d = nd_add(g->body[k+1],d); - } - g->body[k] = d; - g->m = MAX(g->m,k); -} - -int head_pbucket(PGeoBucket g) -{ - int j,i,c,k,nv,sum; - 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 = HC(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 ) - HC(gj) = sum; - else - g->body[j] = nd_remove_head(gj); - j = i; - gj = g->body[j]; - dj = HDL(gj); - sum = HC(gj); - } else if ( c == 0 ) { - sum = sum+HC(gi)-nd_mod; - if ( sum < 0 ) - sum += nd_mod; - g->body[i] = nd_remove_head(gi); - } - } - } - if ( j < 0 ) - return -1; - else if ( sum ) { - HC(gj) = sum; - return j; - } else - g->body[j] = nd_remove_head(gj); - } -} - -ND normalize_pbucket(PGeoBucket g) -{ - int i; - ND r,t; - - r = 0; - for ( i = 0; i <= g->m; i++ ) - r = nd_add(r,g->body[i]); - return r; -} - -ND nd_nf(ND g,int full) -{ - ND u,p,d,red; - NODE l; - NM m,mrd; - int sugar,psugar,n,h_reducible,h; - PGeoBucket bucket; - - if ( !g ) { - return 0; - } - sugar = g->sugar; - n = g->nv; - bucket = create_pbucket(); - add_pbucket(bucket,g); - d = 0; - while ( 1 ) { - h = head_pbucket(bucket); - if ( h < 0 ) { - if ( d ) - d->sugar = sugar; - return d; - } - g = bucket->body[h]; - red = nd_find_reducer(g); - if ( red ) { - bucket->body[h] = nd_remove_head(g); - red = nd_remove_head(red); - add_pbucket(bucket,red); - sugar = MAX(sugar,red->sugar); - } else if ( !full ) { - g = normalize_pbucket(bucket); - if ( g ) - g->sugar = sugar; - return g; - } else { - m = BDY(g); - if ( NEXT(m) ) { - BDY(g) = NEXT(m); NEXT(m) = 0; - } else { - FREEND(g); g = 0; - } - bucket->body[h] = g; - NEXT(m) = 0; - if ( d ) { - for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) ); - NEXT(mrd) = m; - } else { - MKND(n,m,d); - } - } - } -} -#endif - -NODE nd_gb(NODE f) -{ - int i,nh; - NODE r,g,gall; - ND_pairs d; - ND_pairs l; - ND h,nf; - - for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) { - i = (int)BDY(r); - d = update_pairs(d,g,i); - g = update_base(g,i); - gall = append_one(gall,i); - } - while ( d ) { -#if 0 - ndp_print(d); -#endif - l = nd_minp(d,&d); - h = nd_sp(l); - nf = nd_nf(h,!Top); - if ( nf ) { - printf("+"); fflush(stdout); -#if 0 - ndl_print(HDL(nf)); fflush(stdout); -#endif - nh = nd_newps(nf); - d = update_pairs(d,g,nh); - g = update_base(g,nh); - gall = append_one(gall,nh); - } else { - printf("."); fflush(stdout); - } - } - return g; -} - -ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t) -{ - ND_pairs d1,nd,cur,head,prev; - - if ( !g ) return d; - d = crit_B(d,t); - d1 = nd_newpairs(g,t); - d1 = crit_M(d1); - d1 = crit_F(d1); - prev = 0; cur = head = d1; - while ( cur ) { - if ( crit_2( cur->i1,cur->i2 ) ) { - if ( !prev ) { - head = cur = NEXT(cur); - } else { - cur = NEXT(prev) = NEXT(cur); - } - } else { - prev = cur; - cur = NEXT(cur); - } - } - if ( !d ) - return head; - else { - nd = d; - while ( NEXT(nd) ) - nd = NEXT(nd); - NEXT(nd) = head; - return d; - } -} - -ND_pairs nd_newpairs( NODE g, int t ) -{ - NODE h; - unsigned int *dl; - int td,ts,s; - ND_pairs r,r0; - - dl = HDL(nps[t]); - td = HTD(nps[t]); - ts = nps[t]->sugar - td; - for ( r0 = 0, h = g; h; h = NEXT(h) ) { - NEXTND_pairs(r0,r); - r->i1 = (int)BDY(h); - r->i2 = t; - ndl_lcm(HDL(nps[r->i1]),dl,r->lcm); - r->td = ndl_td(r->lcm); - s = nps[r->i1]->sugar-HTD(nps[r->i1]); - r->sugar = MAX(s,ts) + r->td; - } - NEXT(r) = 0; - return r0; -} - -ND_pairs crit_B( ND_pairs d, int s ) -{ - ND_pairs cur,head,prev; - unsigned int *t,*tl,*lcm; - int td,tdl; - - if ( !d ) return 0; - t = HDL(nps[s]); - prev = 0; - head = cur = d; - lcm = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int)); - while ( cur ) { - tl = cur->lcm; - if ( ndl_reducible(tl,t) - && (ndl_lcm(HDL(nps[cur->i1]),t,lcm),!ndl_equal(lcm,tl)) - && (ndl_lcm(HDL(nps[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) { - if ( !prev ) { - head = cur = NEXT(cur); - } else { - cur = NEXT(prev) = NEXT(cur); - } - } else { - prev = cur; - cur = NEXT(cur); - } - } - return head; -} - -ND_pairs crit_M( ND_pairs d1 ) -{ - ND_pairs e,d2,d3,dd,p; - unsigned int *id,*jd; - int itd,jtd; - - for ( dd = 0, e = d1; e; e = d3 ) { - if ( !(d2 = NEXT(e)) ) { - NEXT(e) = dd; - return e; - } - id = e->lcm; - itd = e->td; - for ( d3 = 0; d2; d2 = p ) { - p = NEXT(d2), - jd = d2->lcm; - jtd = d2->td; - if ( jtd == itd ) - if ( id == jd ); - else if ( ndl_reducible(jd,id) ) continue; - else if ( ndl_reducible(id,jd) ) goto delit; - else ; - else if ( jtd > itd ) - if ( ndl_reducible(jd,id) ) continue; - else ; - else if ( ndl_reducible(id,jd ) ) goto delit; - NEXT(d2) = d3; - d3 = d2; - } - NEXT(e) = dd; - dd = e; - continue; - /**/ - delit: NEXT(d2) = d3; - d3 = d2; - for ( ; p; p = d2 ) { - d2 = NEXT(p); - NEXT(p) = d3; - d3 = p; - } - } - return dd; -} - -ND_pairs crit_F( ND_pairs d1 ) -{ - ND_pairs rest, head; - ND_pairs last, p, r, w; - int s; - - for ( head = last = 0, p = d1; NEXT(p); ) { - r = w = equivalent_pairs(p,&rest); - s = r->sugar; - while ( w = NEXT(w) ) - if ( crit_2(w->i1,w->i2) ) { - r = w; - break; - } else if ( w->sugar < s ) { - r = w; - s = r->sugar; - } - if ( last ) NEXT(last) = r; - else head = r; - NEXT(last = r) = 0; - p = rest; - if ( !p ) return head; - } - if ( !last ) return p; - NEXT(last) = p; - return head; -} - -int crit_2( int dp1, int dp2 ) -{ - return ndl_disjoint(HDL(nps[dp1]),HDL(nps[dp2])); -} - -static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest ) -{ - ND_pairs w,p,r,s; - unsigned int *d; - int td; - - w = d1; - d = w->lcm; - td = w->td; - s = NEXT(w); - NEXT(w) = 0; - for ( r = 0; s; s = p ) { - p = NEXT(s); - if ( td == s->td && ndl_equal(d,s->lcm) ) { - NEXT(s) = w; - w = s; - } else { - NEXT(s) = r; - r = s; - } - } - *prest = r; - return w; -} - -NODE update_base(NODE nd,int ndp) -{ - unsigned int *dl, *dln; - NODE last, p, head; - int td,tdn; - - dl = HDL(nps[ndp]); - td = HTD(nps[ndp]); - for ( head = last = 0, p = nd; p; ) { - dln = HDL(nps[(int)BDY(p)]); - tdn = HTD(nps[(int)BDY(p)]); - if ( tdn >= td && ndl_reducible( dln, dl ) ) { - p = NEXT(p); - if ( last ) NEXT(last) = p; - } else { - if ( !last ) head = p; - p = NEXT(last = p); - } - } - head = append_one(head,ndp); - return head; -} - -ND_pairs nd_minp( ND_pairs d, ND_pairs *prest ) -{ - ND_pairs m,ml,p,l; - unsigned int *lcm; - int s,td,len,tlen,c; - - if ( !(p = NEXT(m = d)) ) { - *prest = p; - NEXT(m) = 0; - return m; - } - lcm = m->lcm; - s = m->sugar; - td = m->td; - len = nd_length(nps[m->i1])+nd_length(nps[m->i2]); - for ( ml = 0, l = m; p; p = NEXT(l = p) ) { - if (p->sugar < s) - goto find; - else if ( p->sugar == s ) { - if ( p->td < td ) - goto find; - else if ( p->td == td ) { - c = ndl_compare(p->lcm,lcm); - if ( c < 0 ) - goto find; - else if ( c == 0 ) { - tlen = nd_length(nps[p->i1])+nd_length(nps[p->i2]); - if ( tlen < len ) - goto find; - } - } - } - continue; -find: - ml = l; - m = p; - lcm = m->lcm; - s = m->sugar; - td = m->td; - len = tlen; - } - if ( !ml ) *prest = NEXT(m); - else { - NEXT(ml) = NEXT(m); - *prest = d; - } - NEXT(m) = 0; - return m; -} - -int nd_newps(ND a) -{ - if ( nd_psn == nd_pslen ) { - nd_pslen *= 2; - nps = (ND *)REALLOC((char *)nps,nd_pslen*sizeof(ND)); - } - nd_monic(a); - nps[nd_psn] = a; - return nd_psn++; -} - -NODE NODE_sortb(NODE f,int); -ND dptond(DP); -DP ndtodp(ND); - -NODE nd_setup(NODE f) -{ - int i; - NODE s,s0,f0; - -#if 0 - f0 = f = NODE_sortb(f,1); -#endif - nd_psn = length(f); nd_pslen = 2*nd_psn; - nps = (ND *)MALLOC(nd_pslen*sizeof(ND)); - nd_bpe = 4; - nd_epw = (sizeof(unsigned int)*8)/nd_bpe; - nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0); - if ( nd_bpe < 32 ) { - nd_mask0 = (1<id ) - error("nd_gr : unsupported order"); - switch ( ord->ord.simple ) { - case 0: - is_rlex = 1; - break; - case 1: - is_rlex = 0; - break; - default: - error("nd_gr : unsupported order"); - } - initd(ord); - nd_mod = m; - for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) { - ptod(CO,vv,(P)BDY(t),&b); - _dp_mod(b,m,0,&c); - if ( c ) { - NEXTNODE(fd0,fd); BDY(fd) = (pointer)c; - } - } - if ( fd0 ) NEXT(fd) = 0; - s = nd_setup(fd0); - x = nd_gb(s); -#if 0 - x = nd_reduceall(x,m); -#endif - for ( r0 = 0; x; x = NEXT(x) ) { - NEXTNODE(r0,r); - a = ndtodp(nps[(int)BDY(x)]); - _dtop_mod(CO,vv,a,(P *)&BDY(r)); - } - if ( r0 ) NEXT(r) = 0; - MKLIST(*rp,r0); -} - -void dltondl(int n,DL dl,unsigned int *r) -{ - unsigned int *d; - int i; - - d = dl->d; - bzero(r,nd_wpd*sizeof(unsigned int)); - if ( is_rlex ) - for ( i = 0; i < n; i++ ) - r[(n-1-i)/nd_epw] |= (d[i]<<((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe)); - else - for ( i = 0; i < n; i++ ) - r[i/nd_epw] |= d[i]<<((nd_epw-(i%nd_epw)-1)*nd_bpe); -} - -DL ndltodl(int n,int td,unsigned int *ndl) -{ - DL dl; - int *d; - int i; - - NEWDL(dl,n); - dl->td = td; - d = dl->d; - if ( is_rlex ) - for ( i = 0; i < n; i++ ) - d[i] = (ndl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe)) - &((1<>((nd_epw-(i%nd_epw)-1)*nd_bpe)) - &((1<c = ITOS(t->c); - m->td = t->dl->td; - dltondl(n,t->dl,m->dl); - } - NEXT(m) = 0; - MKND(n,m0,d); - d->nv = n; - d->sugar = p->sugar; - return d; -} - -DP ndtodp(ND p) -{ - DP d; - MP m0,m; - NM t; - int n; - - if ( !p ) - return 0; - n = NV(p); - m0 = 0; - for ( t = BDY(p); t; t = NEXT(t) ) { - NEXTMP(m0,m); - m->c = STOI(t->c); - m->dl = ndltodl(n,t->td,t->dl); - } - NEXT(m) = 0; - MKDP(n,m0,d); - d->sugar = p->sugar; - return d; -} - -void ndl_print(unsigned int *dl) -{ - int n; - int i; - - n = nd_nvar; - printf("<<"); - if ( is_rlex ) - for ( i = 0; i < n; i++ ) - printf(i==n-1?"%d":"%d,", - (dl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe)) - &((1<>((nd_epw-(i%nd_epw)-1)*nd_bpe)) - &((1<>"); -} - -void nd_print(ND p) -{ - NM m; - - if ( !p ) - printf("0\n"); - else { - for ( m = BDY(p); m; m = NEXT(m) ) { - printf("+%d*",m->c); - ndl_print(m->dl); - } - printf("\n"); - } -} - -void ndp_print(ND_pairs d) -{ - ND_pairs t; - - for ( t = d; t; t = NEXT(t) ) { - printf("%d,%d ",t->i1,t->i2); - } - printf("\n"); -} - -void nd_monic(ND p) -{ - if ( !p ) - return; - else - nd_mul_c(p,invm(HC(p),nd_mod)); -} - -void nd_mul_c(ND p,int mul) -{ - NM m; - int c,c1; - - if ( !p ) - return; - for ( m = BDY(p); m; m = NEXT(m) ) { - c1 = C(m); - DMAR(c1,mul,0,nd_mod,c); - C(m) = c; - } }