=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/distm.c,v retrieving revision 1.11 retrieving revision 1.16 diff -u -p -r1.11 -r1.16 --- OpenXM_contrib2/asir2000/engine/distm.c 2003/07/18 10:13:13 1.11 +++ OpenXM_contrib2/asir2000/engine/distm.c 2003/07/22 10:11:43 1.16 @@ -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.10 2002/01/28 00:54:43 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.15 2003/07/22 10:00:51 noro Exp $ */ #include "ca.h" #include "inline.h" @@ -896,6 +896,8 @@ void adddl_dup(int n,DL d1,DL d2,DL *dr) #define INLINE #endif +#define REDTAB_LEN 32003 + typedef struct oPGeoBucket { int m; struct oND *body[32]; @@ -921,16 +923,22 @@ typedef struct oND_pairs { unsigned int lcm[1]; } *ND_pairs; -static ND *nps; +static ND *nd_ps; +static unsigned int **nd_bound; int nd_mod,nd_nvar; int is_rlex; int nd_epw,nd_bpe,nd_wpd; +unsigned int nd_mask[32]; +unsigned int nd_mask0,nd_mask1; NM _nm_free_list; ND _nd_free_list; ND_pairs _ndp_free_list; +NM *nd_red; +int nd_red_len; extern int Top,Reverse; int nd_psn,nd_pslen; +int nd_found,nd_create,nd_notfirst; void GC_gcollect(); NODE append_one(NODE,int); @@ -948,6 +956,7 @@ if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEX #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 FREENDP(m) NEXT(m)=_ndp_free_list; _ndp_free_list=(m) #define FREEND(m) BDY(m)=(NM)_nd_free_list; _nd_free_list=(m) #define NEXTND_pairs(r,c) \ @@ -972,19 +981,32 @@ 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_sp(ND_pairs p); -ND nd_reducer(ND p1,ND p2); -ND nd_nf(NODE b,ND g,ND *ps,int full); +ND nd_mul_term(ND p,int td,unsigned int *d); +int nd_sp(ND_pairs p,ND *nf); +int nd_find_reducer(ND g,ND *red); +int nd_nf(ND g,int full,ND *nf); +ND nd_reduce(ND p1,ND p2); +ND nd_reduce_special(ND p1,ND p2); +void nd_free(ND p); 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_redlist(); +void nd_append_red(unsigned int *d,int td,int i); +unsigned int *nd_compute_bound(ND p); +ND_pairs nd_reconstruct(ND_pairs); +void nd_setup_parameters(); +ND nd_dup(ND p,int obpe); +void ndl_dup(int obpe,unsigned int *d,unsigned int *r); void nd_free_private_storage() { _nd_free_list = 0; _nm_free_list = 0; + nd_red = 0; GC_gcollect(); } @@ -993,7 +1015,7 @@ void _NM_alloc() NM p; int i; - for ( i = 0; i < 1024; i++ ) { + for ( i = 0; i < 16; i++ ) { p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int)); p->next = _nm_free_list; _nm_free_list = p; } @@ -1015,7 +1037,7 @@ void _NDP_alloc() ND_pairs p; int i; - for ( i = 0; i < 1024; 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; @@ -1038,7 +1060,7 @@ INLINE nd_length(ND p) INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2) { unsigned int u1,u2; - int i; + int i,j; switch ( nd_bpe ) { case 4: @@ -1055,6 +1077,17 @@ INLINE int ndl_reducible(unsigned int *d1,unsigned int } 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]; @@ -1079,14 +1112,19 @@ INLINE int ndl_reducible(unsigned int *d1,unsigned int return 1; break; default: - error("ndl_reducible : not implemented yet"); + 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; } } -INLINE void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d) +void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d) { unsigned int t1,t2,u,u1,u2; - int i; + int i,j; switch ( nd_bpe ) { case 4: @@ -1103,6 +1141,17 @@ INLINE void ndl_lcm(unsigned int *d1,unsigned *d2,unsi 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]; @@ -1128,43 +1177,26 @@ INLINE void ndl_lcm(unsigned int *d1,unsigned *d2,unsi } break; default: - error("ndl_lcm : not implemented yet"); + 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; } } -INLINE int ndl_td(unsigned int *d) +int ndl_td(unsigned int *d) { unsigned int t,u; - int i; + int i,j; - switch ( nd_bpe ) { - case 4: - for ( t = 0, i = 0; i < nd_wpd; i++ ) { - u = d[i]; - t += ((u&0xf0000000)>>28)+((u&0xf000000)>>24) - +((u&0xf00000)>>20)+((u&0xf0000)>>16) - +((u&0xf000)>>12)+((u&0xf00)>>8)+((u&0xf0)>>4)+(u&0xf); - } - break; - case 8: - for ( t = 0, i = 0; i < nd_wpd; i++ ) { - u = d[i]; - t += ((u&0xff000000)>>24)+((u&0xff0000)>>16) - +((u&0xff00)>>8)+(u&0xff); - } - break; - case 16: - for ( t = 0, i = 0; i < nd_wpd; i++ ) { - u = d[i]; - t += ((u&0xffff0000)>>16)+(u&0xffff); - } - break; - case 32: - for ( t = 0, i = 0; i < nd_wpd; i++ ) - t += d[i]; - break; - default: - error("ndl_td : not implemented yet"); + 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; } @@ -1173,10 +1205,10 @@ INLINE int ndl_compare(unsigned int *d1,unsigned int * { int i; - for ( i = 0; i < nd_wpd; i++ ) - if ( d1[i] > d2[i] ) + for ( i = 0; i < nd_wpd; i++, d1++, d2++ ) + if ( *d1 > *d2 ) return is_rlex ? -1 : 1; - else if ( d1[i] < d2[i] ) + else if ( *d1 < *d2 ) return is_rlex ? 1 : -1; return 0; } @@ -1195,11 +1227,12 @@ INLINE void ndl_add(unsigned int *d1,unsigned int *d2, { int i; - for ( i = 0; i < nd_wpd; i++ ) + for ( i = 0; i < nd_wpd; i++ ) { d[i] = d1[i]+d2[i]; + } } -INLINE void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d) +void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d) { int i; @@ -1207,10 +1240,10 @@ INLINE void ndl_sub(unsigned int *d1,unsigned int *d2, d[i] = d1[i]-d2[i]; } -INLINE int ndl_disjoint(unsigned int *d1,unsigned int *d2) +int ndl_disjoint(unsigned int *d1,unsigned int *d2) { unsigned int t1,t2,u,u1,u2; - int i; + int i,j; switch ( nd_bpe ) { case 4: @@ -1227,6 +1260,17 @@ INLINE int ndl_disjoint(unsigned int *d1,unsigned int } 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]; @@ -1251,117 +1295,140 @@ INLINE int ndl_disjoint(unsigned int *d1,unsigned int return 1; break; default: - error("ndl_disjoint : not implemented yet"); + 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_add(ND p1,ND p2) +ND nd_reduce(ND p1,ND p2) { - int n,c; - int t; - ND r; - NM m1,m2,mr0,mr,s; + int c,c1,c2,t,td,td2,mul; + NM m2,prev,head,cur,new; + unsigned int *d; if ( !p1 ) - return p2; - else if ( !p2 ) - return p1; + return 0; else { - for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) { - if ( m1->td > m2->td ) + 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 ( m1->td < m2->td ) + else if ( cur->td < td2 ) c = -1; else - c = ndl_compare(m1->dl,m2->dl); + c = ndl_compare(cur->dl,new->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); + 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 { - FREENM(s); + NEXT(prev) = NEXT(cur); + FREENM(cur); + cur = NEXT(prev); } - s = m2; m2 = NEXT(m2); FREENM(s); + m2 = NEXT(m2); break; case 1: - s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s); + prev = cur; + cur = NEXT(cur); break; case -1: - s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s); + 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; } } - 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; - } -} - -INLINE ND nd_mul_nm(ND p,NM m0) -{ - NM m,mr,mr0; - unsigned int *d,*dt,*dm; - int c,n,td; - 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); - C(mr) = (C(m)*c)%nd_mod; - mr->td = m->td+td; - ndl_add(m->dl,d,mr->dl); + FREENM(new); + if ( head ) { + BDY(p1) = head; + p1->sugar = MAX(p1->sugar,p2->sugar+td); + return p1; + } else { + FREEND(p1); + return 0; } - NEXT(mr) = 0; - MKND(NV(p),mr0,r); - r->sugar = p->sugar + td; - return r; + } } -ND nd_reduce(ND p1,ND p2) +/* HDL(p1) = HDL(p2) */ + +ND nd_reduce_special(ND p1,ND p2) { - int c,t,td,td2,mul; + int c,c1,c2,t,td,td2,mul; NM m2,prev,head,cur,new; - unsigned int *d; if ( !p1 ) return 0; else { - mul = ((nd_mod-HC(p1))*invm(HC(p2),nd_mod))%nd_mod; - td = HTD(p1)-HTD(p2); - d = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int)); - ndl_sub(HDL(p1),HDL(p2),d); + c2 = invm(HC(p2),nd_mod); + c1 = nd_mod-HC(p1); + DMAR(c1,c2,0,nd_mod,mul); 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); + td2 = new->td = m2->td; if ( !cur ) { - C(new) = (C(m2)*mul)%nd_mod; + c1 = C(m2); + DMAR(c1,mul,0,nd_mod,c2); + C(new) = c2; + bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int)); if ( !prev ) { prev = new; NEXT(prev) = 0; @@ -1380,10 +1447,12 @@ ND nd_reduce(ND p1,ND p2) else if ( cur->td < td2 ) c = -1; else - c = ndl_compare(cur->dl,new->dl); + c = ndl_compare(cur->dl,m2->dl); switch ( c ) { case 0: - t = (C(cur)+C(m2)*mul)%nd_mod; + c2 = C(m2); + c1 = C(cur); + DMAR(c2,mul,c1,nd_mod,t); if ( t ) C(cur) = t; else if ( !prev ) { @@ -1402,14 +1471,19 @@ ND nd_reduce(ND p1,ND p2) cur = NEXT(cur); break; case -1: + bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int)); if ( !prev ) { /* cur = head */ prev = new; - C(prev) = (C(m2)*mul)%nd_mod; + c2 = C(m2); + DMAR(c2,mul,0,nd_mod,c1); + C(prev) = c1; NEXT(prev) = head; head = prev; } else { - C(new) = (C(m2)*mul)%nd_mod; + c2 = C(m2); + DMAR(c2,mul,0,nd_mod,c1); + C(new) = c1; NEXT(prev) = new; NEXT(new) = cur; prev = new; @@ -1432,100 +1506,296 @@ ND nd_reduce(ND p1,ND p2) } } -ND nd_sp(ND_pairs p) +INLINE int ndl_check_bound(unsigned int *d) { + int i; + + for ( i = 0; i < nd_wpd; i++ ) + if ( d[i] & nd_mask1 ) + return 1; + return 0; +} + +int nd_sp(ND_pairs p,ND *rp) +{ NM m; ND p1,p2,t1,t2; - unsigned int *lcm; + unsigned int *lcm,*check; int td; - p1 = nps[p->i1]; - p2 = nps[p->i2]; + check = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int)); + p1 = nd_ps[p->i1]; + p2 = nd_ps[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; + ndl_add(nd_bound[p->i1],m->dl,check); + if ( ndl_check_bound(check) ) + return 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); + ndl_add(nd_bound[p->i2],m->dl,check); + if ( ndl_check_bound(check) ) { + nd_free(t1); + return 0; + } t2 = nd_mul_nm(p2,m); FREENM(m); - return nd_add(t1,t2); + *rp = nd_add(t1,t2); + return 1; } -ND nd_reducer(ND p1,ND p2) +int ndl_hash_value(int td,unsigned int *d) { + int i; + int r; + + r = td; + for ( i = 0; i < nd_wpd; i++ ) + r = ((r<<16)+d[i])%REDTAB_LEN; + return r; +} + +int nd_find_reducer(ND g, ND *rp) +{ NM m; - ND r; + ND r,p; + int i,c1,c2,c; + int d,k,append,index; + unsigned int *check; + NM t; + d = ndl_hash_value(HTD(g),HDL(g)); + for ( m = nd_red[d], k = 0; m; m = NEXT(m), k++ ) { + if ( HTD(g) == m->td && ndl_equal(HDL(g),m->dl) ) { + if ( k > 0 ) nd_notfirst++; + index = m->c; + append = 0; + nd_found++; + goto found; + } + } + + for ( i = 0; i < nd_psn; i++ ) { + p = nd_ps[i]; + if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) { + index = i; + append = 1; + nd_create++; + goto found; + } + } + return 0; + +found: NEWNM(m); - C(m) = ((nd_mod-HC(p1))*invm(HC(p2),nd_mod))%nd_mod; - m->td = HTD(p1)-HTD(p2); - ndl_sub(HDL(p1),HDL(p2),m->dl); + p = nd_ps[index]; + ndl_sub(HDL(g),HDL(p),m->dl); + + check = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int)); + ndl_add(nd_bound[index],m->dl,check); + if ( ndl_check_bound(check) ) { + FREENM(m); + return -1; + } + + 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); NEXT(m) = 0; - r = nd_mul_nm(p2,m); + *rp = r = nd_mul_nm(p,m); FREENM(m); - return r; + + if ( append ) nd_append_red(HDL(g),HTD(g),i); + return 1; } +ND nd_find_monic_reducer(ND g) +{ + int *d; + ND p,r; + int i; + + for ( i = 0; i < nd_psn; i++ ) { + p = nd_ps[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); + for ( mr0 = 0; 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(NODE b,ND g,ND *ps,int full) +/* ret=1 : success, ret=0 : overflow */ +int nd_nf(ND g,int full,ND *rp) { - ND u,p,d,red; - NODE l; - NM m,mrd; - int sugar,psugar,n,h_reducible; + ND p,d,red; + NM m,mrd,tail; + int n,sugar,psugar,stat; if ( !g ) { - return 0; + *rp = 0; + return 1; } sugar = g->sugar; - n = g->nv; + n = NV(g); for ( d = 0; g; ) { - for ( h_reducible = 0, l = b; l; l = NEXT(l) ) { - p = ps[(int)BDY(l)]; - if ( HTD(g)>=HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) { - h_reducible = 1; - psugar = HTD(g)-HTD(p) + p->sugar; -#if 0 - red = nd_reducer(g,p); - g = nd_add(g,red); + /* stat=1 : found, stat=0 : not found, stat=-1 : overflow */ + stat = nd_find_reducer(g,&red); + if ( stat == 1 ) { +#if 1 + g = nd_add(g,red); + sugar = MAX(sugar,red->sugar); #else - g = nd_reduce(g,p); + psugar = (HTD(g)-HTD(red))+red->sugar; + g = nd_reduce(g,red); + sugar = MAX(sugar,psugar); #endif - sugar = MAX(sugar,psugar); - if ( !g ) { - if ( d ) - d->sugar = sugar; - return d; - } - break; + } else if ( stat == -1 ) { + nd_free(g); + nd_free(d); + return 0; + } else if ( !full ) { + *rp = g; + return 1; + } else { + m = BDY(g); + if ( NEXT(m) ) { + BDY(g) = NEXT(m); NEXT(m) = 0; + } else { + FREEND(g); g = 0; } - } - if ( !h_reducible ) { - /* head term is not reducible */ - if ( !full ) { - if ( g ) - g->sugar = sugar; - return g; + if ( d ) { + NEXT(tail)=m; + tail=m; } else { - m = BDY(g); - if ( NEXT(m) ) { - BDY(g) = NEXT(m); NEXT(m) = 0; - } else { - FREEND(g); g = 0; - } - if ( d ) { - for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) ); - NEXT(mrd) = m; - } else { - MKND(n,m,d); - } + MKND(n,m,d); + tail = BDY(d); } } } if ( d ) d->sugar = sugar; - return d; + *rp = d; + return 1; } #else @@ -1632,7 +1902,7 @@ ND normalize_pbucket(PGeoBucket g) return r; } -ND nd_nf(NODE b,ND g,ND *ps,int full) +ND nd_nf(ND g,int full) { ND u,p,d,red; NODE l; @@ -1656,41 +1926,31 @@ ND nd_nf(NODE b,ND g,ND *ps,int full) return d; } g = bucket->body[h]; - for ( h_reducible = 0, l = b; l; l = NEXT(l) ) { - p = ps[(int)BDY(l)]; - if ( ndl_reducible(HDL(g),HDL(p)) ) { - h_reducible = 1; - psugar = HTD(g)-HTD(p) + p->sugar; - red = nd_reducer(g,p); - bucket->body[h] = nd_remove_head(g); - red = nd_remove_head(red); - add_pbucket(bucket,red); - sugar = MAX(sugar,psugar); - break; + 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; } - } - if ( !h_reducible ) { - /* head term is not reducible */ - if ( !full ) { - g = normalize_pbucket(bucket); - if ( g ) - g->sugar = sugar; - return g; + bucket->body[h] = g; + NEXT(m) = 0; + if ( d ) { + for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) ); + NEXT(mrd) = m; } 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); - } + MKND(n,m,d); } } } @@ -1699,7 +1959,7 @@ ND nd_nf(NODE b,ND g,ND *ps,int full) NODE nd_gb(NODE f) { - int i,nh; + int i,nh,sugar,stat; NODE r,g,gall; ND_pairs d; ND_pairs l; @@ -1711,24 +1971,35 @@ NODE nd_gb(NODE f) g = update_base(g,i); gall = append_one(gall,i); } + sugar = 0; while ( d ) { -#if 0 - ndp_print(d); -#endif +again: l = nd_minp(d,&d); - h = nd_sp(l); - nf = nd_nf(gall,h,nps,!Top); - if ( nf ) { + if ( l->sugar != sugar ) { + sugar = l->sugar; + fprintf(asir_out,"%d",sugar); + } + stat = nd_sp(l,&h); + if ( !stat ) { + NEXT(l) = d; d = l; + d = nd_reconstruct(d); + goto again; + } + stat = nd_nf(h,!Top,&nf); + if ( !stat ) { + NEXT(l) = d; d = l; + d = nd_reconstruct(d); + goto again; + } else 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); + FREENDP(l); } else { printf("."); fflush(stdout); + FREENDP(l); } } return g; @@ -1736,7 +2007,7 @@ NODE nd_gb(NODE f) ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t) { - ND_pairs d1,nd,cur,head,prev; + ND_pairs d1,nd,cur,head,prev,remove; if ( !g ) return d; d = crit_B(d,t); @@ -1746,11 +2017,13 @@ ND_pairs update_pairs( ND_pairs d, NODE /* of index */ prev = 0; cur = head = d1; while ( cur ) { if ( crit_2( cur->i1,cur->i2 ) ) { + remove = cur; if ( !prev ) { head = cur = NEXT(cur); } else { cur = NEXT(prev) = NEXT(cur); } + FREENDP(remove); } else { prev = cur; cur = NEXT(cur); @@ -1774,16 +2047,16 @@ ND_pairs nd_newpairs( NODE g, int t ) int td,ts,s; ND_pairs r,r0; - dl = HDL(nps[t]); - td = HTD(nps[t]); - ts = nps[t]->sugar - td; + dl = HDL(nd_ps[t]); + td = HTD(nd_ps[t]); + ts = nd_ps[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); + ndl_lcm(HDL(nd_ps[r->i1]),dl,r->lcm); r->td = ndl_td(r->lcm); - s = nps[r->i1]->sugar-HTD(nps[r->i1]); + s = nd_ps[r->i1]->sugar-HTD(nd_ps[r->i1]); r->sugar = MAX(s,ts) + r->td; } NEXT(r) = 0; @@ -1792,25 +2065,27 @@ ND_pairs nd_newpairs( NODE g, int t ) ND_pairs crit_B( ND_pairs d, int s ) { - ND_pairs cur,head,prev; + ND_pairs cur,head,prev,remove; unsigned int *t,*tl,*lcm; int td,tdl; if ( !d ) return 0; - t = HDL(nps[s]); + t = HDL(nd_ps[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)) ) { + && (ndl_lcm(HDL(nd_ps[cur->i1]),t,lcm),!ndl_equal(lcm,tl)) + && (ndl_lcm(HDL(nd_ps[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) { + remove = cur; if ( !prev ) { head = cur = NEXT(cur); } else { cur = NEXT(prev) = NEXT(cur); } + FREENDP(remove); } else { prev = cur; cur = NEXT(cur); @@ -1859,27 +2134,42 @@ ND_pairs crit_M( ND_pairs d1 ) NEXT(p) = d3; d3 = p; } + FREENDP(e); } return dd; } ND_pairs crit_F( ND_pairs d1 ) { - ND_pairs rest, head; + ND_pairs rest, head,remove; 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) ) + w = NEXT(w); + while ( w ) { if ( crit_2(w->i1,w->i2) ) { r = w; + w = NEXT(w); + while ( w ) { + remove = w; + w = NEXT(w); + FREENDP(remove); + } break; } else if ( w->sugar < s ) { + FREENDP(r); r = w; s = r->sugar; + w = NEXT(w); + } else { + remove = w; + w = NEXT(w); + FREENDP(remove); } + } if ( last ) NEXT(last) = r; else head = r; NEXT(last = r) = 0; @@ -1893,7 +2183,7 @@ ND_pairs crit_F( ND_pairs d1 ) int crit_2( int dp1, int dp2 ) { - return ndl_disjoint(HDL(nps[dp1]),HDL(nps[dp2])); + return ndl_disjoint(HDL(nd_ps[dp1]),HDL(nd_ps[dp2])); } static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest ) @@ -1927,11 +2217,11 @@ NODE update_base(NODE nd,int ndp) NODE last, p, head; int td,tdn; - dl = HDL(nps[ndp]); - td = HTD(nps[ndp]); + dl = HDL(nd_ps[ndp]); + td = HTD(nd_ps[ndp]); for ( head = last = 0, p = nd; p; ) { - dln = HDL(nps[(int)BDY(p)]); - tdn = HTD(nps[(int)BDY(p)]); + dln = HDL(nd_ps[(int)BDY(p)]); + tdn = HTD(nd_ps[(int)BDY(p)]); if ( tdn >= td && ndl_reducible( dln, dl ) ) { p = NEXT(p); if ( last ) NEXT(last) = p; @@ -1958,7 +2248,7 @@ ND_pairs nd_minp( ND_pairs d, ND_pairs *prest ) lcm = m->lcm; s = m->sugar; td = m->td; - len = nd_length(nps[m->i1])+nd_length(nps[m->i2]); + len = nd_length(nd_ps[m->i1])+nd_length(nd_ps[m->i2]); for ( ml = 0, l = m; p; p = NEXT(l = p) ) { if (p->sugar < s) goto find; @@ -1970,7 +2260,7 @@ ND_pairs nd_minp( ND_pairs d, ND_pairs *prest ) if ( c < 0 ) goto find; else if ( c == 0 ) { - tlen = nd_length(nps[p->i1])+nd_length(nps[p->i2]); + tlen = nd_length(nd_ps[p->i1])+nd_length(nd_ps[p->i2]); if ( tlen < len ) goto find; } @@ -1998,10 +2288,13 @@ int nd_newps(ND a) { if ( nd_psn == nd_pslen ) { nd_pslen *= 2; - nps = (ND *)REALLOC((char *)nps,nd_pslen*sizeof(ND)); + nd_ps = (ND *)REALLOC((char *)nd_ps,nd_pslen*sizeof(ND)); + nd_bound = (unsigned int **) + REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *)); } nd_monic(a); - nps[nd_psn] = a; + nd_ps[nd_psn] = a; + nd_bound[nd_psn] = nd_compute_bound(a); return nd_psn++; } @@ -2011,22 +2304,27 @@ DP ndtodp(ND); NODE nd_setup(NODE f) { - int i; + int i,td; NODE s,s0,f0; + nd_found = 0; + nd_notfirst = 0; + nd_create = 0; #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_ps = (ND *)MALLOC(nd_pslen*sizeof(ND)); + nd_bound = (unsigned int **)MALLOC(nd_pslen*sizeof(unsigned int *)); nd_bpe = 4; - nd_epw = (sizeof(unsigned int)*8)/nd_bpe; - nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0); + nd_setup_parameters(); nd_free_private_storage(); for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) { - nps[i] = dptond((DP)BDY(f)); - nd_monic(nps[i]); + nd_ps[i] = dptond((DP)BDY(f)); + nd_monic(nd_ps[i]); + nd_bound[i] = nd_compute_bound(nd_ps[i]); } + nd_red = (NM *)MALLOC(REDTAB_LEN*sizeof(NM)); for ( s0 = 0, i = 0; i < nd_psn; i++ ) { NEXTNODE(s0,s); BDY(s) = (pointer)i; } @@ -2072,11 +2370,13 @@ void nd_gr(LIST f,LIST v,int m,struct order_spec *ord, #endif for ( r0 = 0; x; x = NEXT(x) ) { NEXTNODE(r0,r); - a = ndtodp(nps[(int)BDY(x)]); + a = ndtodp(nd_ps[(int)BDY(x)]); _dtop_mod(CO,vv,a,(P *)&BDY(r)); } if ( r0 ) NEXT(r) = 0; MKLIST(*rp,r0); + fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n", + nd_found,nd_notfirst,nd_create); } void dltondl(int n,DL dl,unsigned int *r) @@ -2085,6 +2385,7 @@ void dltondl(int n,DL dl,unsigned int *r) 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)); @@ -2206,12 +2507,171 @@ void ndp_print(ND_pairs d) void nd_monic(ND p) { - int mul; + 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; - mul = invm(HC(p),nd_mod); - for ( m = BDY(p); m; m = NEXT(m) ) - C(m) = (C(m)*mul)%nd_mod; + for ( m = BDY(p); m; m = NEXT(m) ) { + c1 = C(m); + DMAR(c1,mul,0,nd_mod,c); + C(m) = c; + } +} + +void nd_free(ND p) +{ + NM t,s; + + if ( !p ) + return; + t = BDY(p); + while ( t ) { + s = NEXT(t); + FREENM(t); + t = s; + } + FREEND(p); +} + +void nd_append_red(unsigned int *d,int td,int i) +{ + NM m,m0; + int h; + + NEWNM(m); + h = ndl_hash_value(td,d); + m->c = i; + m->td = td; + bcopy(d,m->dl,nd_wpd*sizeof(unsigned int)); + NEXT(m) = nd_red[h]; + nd_red[h] = m; +} + +unsigned int *nd_compute_bound(ND p) +{ + unsigned int *d1,*d2,*t; + NM m; + + if ( !p ) + return 0; + d1 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int)); + d2 = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int)); + bcopy(HDL(p),d1,nd_wpd*sizeof(unsigned int)); + for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) { + ndl_lcm(m->dl,d1,d2); + t = d1; d1 = d2; d2 = t; + } + t = (unsigned int *)MALLOC_ATOMIC(nd_wpd*sizeof(unsigned int)); + bcopy(d1,t,nd_wpd*sizeof(unsigned int)); + return t; +} + +void nd_setup_parameters() { + int i; + + 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<i1 = t->i1; + s->i2 = t->i2; + s->td = t->td; + s->sugar = t->sugar; + ndl_dup(obpe,t->lcm,s->lcm); + } + if ( s0 ) NEXT(s) = 0; + prev_nm_free_list = 0; + prev_ndp_free_list = 0; + GC_gcollect(); + return s0; +} + +void ndl_dup(int obpe,unsigned int *d,unsigned int *r) +{ + int n,i,ei,oepw,cepw,cbpe; + + n = nd_nvar; + oepw = (sizeof(unsigned int)*8)/obpe; + cepw = nd_epw; + cbpe = nd_bpe; + if ( is_rlex ) + for ( i = 0; i < n; i++ ) { + ei = (d[(n-1-i)/oepw]>>((oepw-((n-1-i)%oepw)-1)*obpe)) + &((1<>((oepw-(i%oepw)-1)*obpe)) + &((1<td = m->td; + ndl_dup(obpe,m->dl,mr->dl); + } + NEXT(mr) = 0; + MKND(NV(p),mr0,r); + r->sugar = p->sugar; + return r; + } }