OpenXM_contrib2/asir2000/engine/nd.c - diff

Return to nd.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/nd.c between version 1.6 and 1.26

-version 1.6, 2003/07/25 04:19:54
+version 1.26, 2003/08/07 09:47:08
 Line 1
 Line 1
 Line 1
- /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.5 2003/07/24 03:45:41 noro Exp $ */
+ /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.25 2003/08/07 08:46:50 noro Exp $ */
  #include "ca.h"
  #include "inline.h"
- #define USE_NDV 1
  #if defined(__GNUC__)
  #define INLINE inline
  #elif defined(VISUAL)
-Line 35  typedef struct oNDV {
+Line 33  typedef struct oNDV {
 Line 35  typedef struct oNDV {
 Line 33  typedef struct oNDV {
  typedef struct oNM {
          struct oNM *next;
+         union {
+                 int m;
+                 Q z;
+         } c;
          int td;
-         int c;
          unsigned int dl[1];
  } *NM;
  typedef struct oNMV {
+         union {
+                 int m;
+                 Q z;
+         } c;
          int td;
-         int c;
          unsigned int dl[1];
  } *NMV;
+ typedef struct oRHist {
+         struct oRHist *next;
+         int index;
+         int td,sugar;
+         unsigned int dl[1];
+ } *RHist;
  typedef struct oND_pairs {
          struct oND_pairs *next;
          int i1,i2;
-Line 53  typedef struct oND_pairs {
+Line 64  typedef struct oND_pairs {
 Line 53  typedef struct oND_pairs {
 Line 64  typedef struct oND_pairs {
          unsigned int lcm[1];
  } *ND_pairs;
- static ND *nd_ps;
+ double nd_scale=2;
- static NDV *nd_psv;
  static unsigned int **nd_bound;
- int nd_mod,nd_nvar;
+ int 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;
+ static NDV *nd_ps;
+ static NDV *nd_psq;
+ int *nd_psl;
+ RHist *nd_psh;
  int nd_psn,nd_pslen;
+ RHist *nd_red;
+ int nd_red_len;
  int nd_found,nd_create,nd_notfirst;
- int *nd_psl;
+ int nm_adv;
+ int nmv_adv;
- void GC_gcollect();
+ extern int Top,Reverse;
- NODE append_one(NODE,int);
  #define HTD(d) ((d)->body->td)
  #define HDL(d) ((d)->body->dl)
- #define HC(d) ((d)->body->c)
+ #define HCM(d) ((d)->body->c.m)
+ #define HCQ(d) ((d)->body->c.z)
+ #define CM(a) ((a)->c.m)
+ #define CQ(a) ((a)->c.z)
+ #define DL(a) ((a)->dl)
+ #define TD(a) ((a)->td)
+ #define SG(a) ((a)->sugar)
+ #define LEN(a) ((a)->len)
+ #define NM_ADV(m) (m = (NM)(((char *)m)+nm_adv))
+ #define NEWRHist(r) \
+ ((r)=(RHist)MALLOC(sizeof(struct oRHist)+(nd_wpd-1)*sizeof(unsigned int)))
  #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 NEXTRHist(r,c) \
+ if(!(r)){NEWRHist(r);(c)=(r);}else{NEWRHist(NEXT(c));(c)=NEXT(c);}
  #define NEXTNM(r,c) \
  if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEXT(c);}
  #define NEXTNM2(r,c,s) \
-Line 94  if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
+Line 122  if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
 Line 94  if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
 Line 122  if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
  #define NEXTND_pairs(r,c) \
  if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);}
+ void nd_removecont(int mod,ND p);
+ void nd_removecont2(ND p1,ND p2);
+ void ndv_removecont(int mod,NDV p);
+ void ndv_dehomogenize(NDV p);
+ void ndv_mul_c_q(NDV p,Q mul);
+ void nd_mul_c_q(ND p,Q mul);
+ ND normalize_pbucket(int mod,PGeoBucket g);
+ int head_pbucket(int mod,PGeoBucket g);
+ int head_pbucket_q(PGeoBucket g);
+ void add_pbucket(int mod,PGeoBucket g,ND d,int l);
+ void free_pbucket(PGeoBucket b);
+ void mulq_pbucket(PGeoBucket g,Q c);
+ PGeoBucket create_pbucket();
+ ND nd_remove_head(ND p);
+ void GC_gcollect();
+ NODE append_one(NODE,int);
+ void removecont_array(Q *c,int n);
  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);
+ void nd_gr_trace(LIST f,LIST v,int m,int homo,struct order_spec *ord,LIST *rp);
- int nd_newps(ND a);
+ void nd_setup(int mod,NODE f);
+ void nd_setup_trace(int mod,NODE f);
+ int nd_newps(int mod,ND a);
+ int nd_newps_trace(int mod,ND nf,ND nfq);
  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 );
-Line 106  ND_pairs crit_F( ND_pairs d1 );
+Line 156  ND_pairs crit_F( ND_pairs d1 );
 Line 106  ND_pairs crit_F( ND_pairs d1 );
 Line 156  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);
+ NODE nd_gb(int m);
+ NODE nd_gb_trace(int m);
  void nd_free_private_storage();
  void _NM_alloc();
  void _ND_alloc();
  int ndl_td(unsigned int *d);
- ND nd_add(ND p1,ND p2);
+ int ndl_dehomogenize(unsigned int *p);
- ND nd_mul_nm(ND p,NM m0);
+ ND nd_add(int mod,ND p1,ND p2);
- ND nd_mul_ind_nm(int index,NM m0);
+ ND nd_add_q(ND p1,ND p2);
- ND nd_mul_term(ND p,int td,unsigned int *d);
+ ND nd_mul_nm(int mod,ND p,NM m0);
- int nd_sp(ND_pairs p,ND *nf);
+ ND nd_mul_ind_nm(int mod,int index,NM m0);
+ int nd_sp(int mod,ND_pairs p,ND *nf);
  int nd_find_reducer(ND g);
- int nd_nf(ND g,int full,ND *nf);
+ int nd_find_reducer_direct(ND g,NDV *ps,int len);
+ int nd_nf(int mod,ND g,int full,ND *nf);
  ND nd_reduce(ND p1,ND p2);
  ND nd_reduce_special(ND p1,ND p2);
+ NODE nd_reduceall(int m,NODE f);
  void nd_free(ND p);
+ void ndv_free(NDV p);
  void ndl_print(unsigned int *dl);
  void nd_print(ND p);
+ void nd_print_q(ND p);
+ void ndv_print(NDV p);
+ void ndv_print_q(NDV p);
  void ndp_print(ND_pairs d);
  int nd_length(ND p);
- void nd_monic(ND p);
+ void nd_mul_c(int mod,ND p,int mul);
- 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);
  unsigned int *dp_compute_bound(DP p);
- ND_pairs nd_reconstruct(ND_pairs);
+ ND_pairs nd_reconstruct(int mod,int trace,ND_pairs ndp);
+ void nd_reconstruct_direct(int mod,NDV *ps,int len);
  void nd_setup_parameters();
- ND nd_dup(ND p,int obpe);
+ void nd_realloc(ND p,int obpe);
  ND nd_copy(ND p);
  void ndl_dup(int obpe,unsigned int *d,unsigned int *r);
- #if USE_NDV
  #define NMV_ADV(m) (m = (NMV)(((char *)m)+nmv_adv))
  #define NEWNDV(d) ((d)=(NDV)MALLOC(sizeof(struct oNDV)))
- #define MKNDV(n,m,l,d) NEWNDV(d); (d)->nv=(n); BDY(d)=(m); (d)->len = l;
+ #define MKNDV(n,m,l,d) NEWNDV(d); NV(d)=(n); BDY(d)=(m); LEN(d) = l;
+ void ndv_mul_c(int mod,NDV p,int mul);
+ NDV ndtondv(int mod,ND p);
+ ND ndvtond(int mod,NDV p);
+ ND ndv_mul_nm(int mod,NDV p,NM m0);
+ void ndv_realloc(NDV p,int obpe,int oadv);
+ NDV dptondv(int,DP);
+ DP ndvtodp(int,NDV);
+ ND dptond(int,DP);
+ DP ndtodp(int,ND);
- int nmv_adv;
- int nmv_len;
- NDV ndv_red;
- ND ndv_add(ND p1,NDV p2);
- NDV ndtondv(ND p);
- void ndv_mul_nm(NDV pv,NM m,NDV r);
- #endif
  void nd_free_private_storage()
  {
          _nd_free_list = 0;
          _nm_free_list = 0;
          _ndp_free_list = 0;
-         nd_red = 0;
+         bzero(nd_red,sizeof(REDTAB_LEN*sizeof(RHist)));
          GC_gcollect();
  }
-Line 165  void _NM_alloc()
+Line 222  void _NM_alloc()
 Line 165  void _NM_alloc()
 Line 222  void _NM_alloc()
          NM p;
          int i;
-         for ( i = 0; i < 16; i++ ) {
+         for ( i = 0; i < 1024; i++ ) {
                  p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
                  p->next = _nm_free_list; _nm_free_list = p;
          }
-Line 187  void _NDP_alloc()
+Line 244  void _NDP_alloc()
 Line 187  void _NDP_alloc()
 Line 244  void _NDP_alloc()
          ND_pairs p;
          int i;
-         for ( i = 0; i < 10240; i++ ) {
+         for ( i = 0; i < 1024; 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;
-Line 271  INLINE int ndl_reducible(unsigned int *d1,unsigned int
+Line 328  INLINE int ndl_reducible(unsigned int *d1,unsigned int
 Line 271  INLINE int ndl_reducible(unsigned int *d1,unsigned int
 Line 328  INLINE int ndl_reducible(unsigned int *d1,unsigned int
          }
  }
+ /* returns the exponent of homo variable */
+ int ndl_dehomogenize(unsigned int *d)
+ {
+         unsigned int mask;
+         unsigned int h;
+         int i;
+         if ( is_rlex ) {
+                 if ( nd_bpe == 32 ) {
+                         h = d[0];
+                         for ( i = 1; i < nd_wpd; i++ )
+                                 d[i-1] = d[i];
+                         d[i-1] = 0;
+                         return h;
+                 } else {
+                         mask = (1<<(nd_epw*nd_bpe))-1;
+                         h = (d[0]>>((nd_epw-1)*nd_bpe))&nd_mask0;
+                         for ( i = 0; i < nd_wpd; i++ )
+                                 d[i] = ((d[i]<<nd_bpe)&mask)
+                                         |(i+1<nd_wpd?((d[i+1]>>((nd_epw-1)*nd_bpe))&nd_mask0):0);
+                         return h;
+                 }
+         } else {
+                 /* do nothing */
+                 return d[nd_wpd-1];
+         }
+ }
  void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d)
  {
          unsigned int t1,t2,u,u1,u2;
-Line 505  int ndl_disjoint(unsigned int *d1,unsigned int *d2)
+Line 591  int ndl_disjoint(unsigned int *d1,unsigned int *d2)
 Line 505  int ndl_disjoint(unsigned int *d1,unsigned int *d2)
 Line 591  int ndl_disjoint(unsigned int *d1,unsigned int *d2)
          }
  }
- 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;
-                 }
-         }
- }
- /* HDL(p1) = HDL(p2) */
- ND nd_reduce_special(ND p1,ND p2)
- {
-         int c,c1,c2,t,td,td2,mul;
-         NM m2,prev,head,cur,new;
-         if ( !p1 )
-                 return 0;
-         else {
-                 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;
-                         if ( !cur ) {
-                                 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;
-                                         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,m2->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:
-                                         bcopy(m2->dl,new->dl,nd_wpd*sizeof(unsigned int));
-                                         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;
-                 }
-         }
- }
  int ndl_check_bound2(int index,unsigned int *d2)
  {
          unsigned int u2;
-Line 773  int ndl_check_bound2(int index,unsigned int *d2)
+Line 660  int ndl_check_bound2(int index,unsigned int *d2)
 Line 773  int ndl_check_bound2(int index,unsigned int *d2)
 Line 660  int ndl_check_bound2(int index,unsigned int *d2)
          }
  }
- int nd_sp(ND_pairs p,ND *rp)
+ int ndl_check_bound2_direct(unsigned int *d1,unsigned int *d2)
  {
-         NM m;
+         unsigned int u1,u2;
-         ND p1,p2,t1,t2;
+         int i,j,k;
-         unsigned int *lcm;
-         int td;
-         p1 = nd_ps[p->i1];
+         switch ( nd_bpe ) {
-         p2 = nd_ps[p->i2];
+                 case 4:
-         lcm = p->lcm;
+                         for ( i = 0; i < nd_wpd; i++ ) {
-         td = p->td;
+                                 u1 = d1[i]; u2 = d2[i];
-         NEWNM(m);
+                                 if ( ((u1>>28)&0xf)+((u2>>28)&0xf) >= 0x10 ) return 1;
-         C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl); NEXT(m) = 0;
+                                 if ( ((u1>>24)&0xf)+((u2>>24)&0xf) >= 0x10 ) return 1;
-         if ( ndl_check_bound2(p->i1,m->dl) )
+                                 if ( ((u1>>20)&0xf)+((u2>>20)&0xf) >= 0x10 ) return 1;
-                 return 0;
+                                 if ( ((u1>>16)&0xf)+((u2>>16)&0xf) >= 0x10 ) return 1;
-         t1 = nd_mul_nm(p1,m);
+                                 if ( ((u1>>12)&0xf)+((u2>>12)&0xf) >= 0x10 ) return 1;
-         C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl);
+                                 if ( ((u1>>8)&0xf)+((u2>>8)&0xf) >= 0x10 ) return 1;
-         if ( ndl_check_bound2(p->i2,m->dl) ) {
+                                 if ( ((u1>>4)&0xf)+((u2>>4)&0xf) >= 0x10 ) return 1;
-                 nd_free(t1);
+                                 if ( (u1&0xf)+(u2&0xf) >= 0x10 ) return 1;
-                 return 0;
+                         }
+                         return 0;
+                         break;
+                 case 6:
+                         for ( i = 0; i < nd_wpd; i++ ) {
+                                 u1 = d1[i]; u2 = d2[i];
+                                 if ( ((u1>>24)&0x3f)+((u2>>24)&0x3f) >= 0x40 ) return 1;
+                                 if ( ((u1>>18)&0x3f)+((u2>>18)&0x3f) >= 0x40 ) return 1;
+                                 if ( ((u1>>12)&0x3f)+((u2>>12)&0x3f) >= 0x40 ) return 1;
+                                 if ( ((u1>>6)&0x3f)+((u2>>6)&0x3f) >= 0x40 ) return 1;
+                                 if ( (u1&0x3f)+(u2&0x3f) >= 0x40 ) return 1;
+                         }
+                         return 0;
+                         break;
+                 case 8:
+                         for ( i = 0; i < nd_wpd; i++ ) {
+                                 u1 = d1[i]; u2 = d2[i];
+                                 if ( ((u1>>24)&0xff)+((u2>>24)&0xff) >= 0x100 ) return 1;
+                                 if ( ((u1>>16)&0xff)+((u2>>16)&0xff) >= 0x100 ) return 1;
+                                 if ( ((u1>>8)&0xff)+((u2>>8)&0xff) >= 0x100 ) return 1;
+                                 if ( (u1&0xff)+(u2&0xff) >= 0x100 ) return 1;
+                         }
+                         return 0;
+                         break;
+                 case 16:
+                         for ( i = 0; i < nd_wpd; i++ ) {
+                                 u1 = d1[i]; u2 = d2[i];
+                                         if ( ((u1>>16)&0xffff)+((u2>>16)&0xffff) > 0x10000 ) return 1;
+                                         if ( (u2&0xffff)+(u2&0xffff) > 0x10000 ) return 1;
+                         }
+                         return 0;
+                         break;
+                 case 32:
+                         for ( i = 0; i < nd_wpd; i++ )
+                                 if ( d1[i]+d2[i]<d1[i] ) return 1;
+                         return 0;
+                         break;
+                 default:
+                         for ( i = 0; i < nd_wpd; i++ ) {
+                                 u1 = d1[i]; u2 = d2[i];
+                                 k = (nd_epw-1)*nd_bpe;
+                                 for ( j = 0; j < nd_epw; j++, k -= nd_bpe )
+                                         if ( ((u1>>k)&nd_mask0)+((u2>>k)&nd_mask0) > nd_mask0 )
+                                                 return 1;
+                         }
+                         return 0;
+                         break;
          }
-         t2 = nd_mul_nm(p2,m);
-         FREENM(m);
-         *rp = nd_add(t1,t2);
-         return 1;
  }
  INLINE int ndl_hash_value(int td,unsigned int *d)
-Line 813  INLINE int ndl_hash_value(int td,unsigned int *d)
+Line 740  INLINE int ndl_hash_value(int td,unsigned int *d)
 Line 813  INLINE int ndl_hash_value(int td,unsigned int *d)
 Line 740  INLINE int ndl_hash_value(int td,unsigned int *d)
  INLINE int nd_find_reducer(ND g)
  {
-         NM m;
+         RHist r;
-         ND p;
          int d,k,i;
          d = ndl_hash_value(HTD(g),HDL(g));
-         for ( m = nd_red[d], k = 0; m; m = NEXT(m), k++ ) {
+         for ( r = nd_red[d], k = 0; r; r = NEXT(r), k++ ) {
-                 if ( HTD(g) == m->td && ndl_equal(HDL(g),m->dl) ) {
+                 if ( HTD(g) == TD(r) && ndl_equal(HDL(g),DL(r)) ) {
                          if ( k > 0 ) nd_notfirst++;
                          nd_found++;
-                         return m->c;
+                         return r->index;
                  }
          }
-         for ( i = 0; i < nd_psn; i++ ) {
+         if ( Reverse )
-                 p = nd_ps[i];
+                 for ( i = nd_psn-1; i >= 0; i-- ) {
-                 if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
+                         r = nd_psh[i];
-                         nd_create++;
+                         if ( HTD(g) >= TD(r) && ndl_reducible(HDL(g),DL(r)) ) {
-                         nd_append_red(HDL(g),HTD(g),i);
+                                 nd_create++;
-                         return i;
+                                 nd_append_red(HDL(g),HTD(g),i);
+                                 return i;
+                         }
                  }
-         }
+         else
+                 for ( i = 0; i < nd_psn; i++ ) {
+                         r = nd_psh[i];
+                         if ( HTD(g) >= TD(r) && ndl_reducible(HDL(g),DL(r)) ) {
+                                 nd_create++;
+                                 nd_append_red(HDL(g),HTD(g),i);
+                                 return i;
+                         }
+                 }
          return -1;
  }
- ND nd_find_monic_reducer(ND g)
+ INLINE int nd_find_reducer_direct(ND g,NDV *ps,int len)
  {
-         int *d;
-         ND p,r;
          int i;
+         NDV r;
-         for ( i = 0; i < nd_psn; i++ ) {
+         if ( Reverse )
-                 p = nd_ps[i];
+                 for ( i = len-1; i >= 0; i-- ) {
-                 if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
+                         r = ps[i];
-                         d = (int *)ALLOCA(nd_wpd*sizeof(int));
+                         if ( HTD(g) >= HTD(r) && ndl_reducible(HDL(g),HDL(r)) ) {
-                         ndl_sub(HDL(g),HDL(p),d);
+                                 nd_append_red(HDL(g),HTD(g),i);
-                         r = nd_mul_term(p,HTD(g)-HTD(p),d);
+                                 return i;
-                         return r;
+                         }
                  }
-         }
+         else
-         return 0;
+                 for ( i = 0; i < len; i++ ) {
+                         r = ps[i];
+                         if ( HTD(g) >= HTD(r) && ndl_reducible(HDL(g),HDL(r)) ) {
+                                 nd_append_red(HDL(g),HTD(g),i);
+                                 return i;
+                         }
+                 }
+         return -1;
  }
- ND nd_add(ND p1,ND p2)
+ ND nd_add(int mod,ND p1,ND p2)
  {
          int n,c;
          int t;
-Line 866  ND nd_add(ND p1,ND p2)
+Line 808  ND nd_add(ND p1,ND p2)
 Line 866  ND nd_add(ND p1,ND p2)
 Line 808  ND nd_add(ND p1,ND p2)
                  return p2;
          else if ( !p2 )
                  return p1;
+         else if ( !mod )
+                 return nd_add_q(p1,p2);
          else {
                  for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
-                         if ( m1->td > m2->td )
+                         if ( TD(m1) > TD(m2) )
                                  c = 1;
-                         else if ( m1->td < m2->td )
+                         else if ( TD(m1) < TD(m2) )
                                  c = -1;
                          else
-                                 c = ndl_compare(m1->dl,m2->dl);
+                                 c = ndl_compare(DL(m1),DL(m2));
                          switch ( c ) {
                                  case 0:
-                                         t = ((C(m1))+(C(m2))) - nd_mod;
+                                         t = ((CM(m1))+(CM(m2))) - mod;
                                          if ( t < 0 )
-                                                 t += nd_mod;
+                                                 t += mod;
                                          s = m1; m1 = NEXT(m1);
                                          if ( t ) {
-                                                 NEXTNM2(mr0,mr,s); C(mr) = (t);
+                                                 NEXTNM2(mr0,mr,s); CM(mr) = (t);
                                          } else {
                                                  FREENM(s);
                                          }
-Line 909  ND nd_add(ND p1,ND p2)
+Line 853  ND nd_add(ND p1,ND p2)
 Line 909  ND nd_add(ND p1,ND p2)
 Line 853  ND nd_add(ND p1,ND p2)
                  else
                          NEXT(mr) = 0;
                  BDY(p1) = mr0;
-                 p1->sugar = MAX(p1->sugar,p2->sugar);
+                 SG(p1) = MAX(SG(p1),SG(p2));
                  FREEND(p2);
                  return p1;
          }
  }
- ND nd_mul_nm(ND p,NM m0)
+ ND nd_add_q(ND p1,ND p2)
  {
-         NM m,mr,mr0;
+         int n,c;
-         unsigned int *d,*dt,*dm;
-         int c,n,td,i,c1,c2;
-         int *pt,*p1,*p2;
          ND r;
+         NM m1,m2,mr0,mr,s;
+         Q t;
-         if ( !p )
+         if ( !p1 )
-                 return 0;
+                 return p2;
+         else if ( !p2 )
+                 return p1;
          else {
-                 n = NV(p); m = BDY(p);
+                 for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
-                 d = m0->dl; td = m0->td; c = C(m0);
+                         if ( TD(m1) > TD(m2) )
-                 mr0 = 0;
+                                 c = 1;
-                 for ( ; m; m = NEXT(m) ) {
+                         else if ( TD(m1) < TD(m2) )
-                         NEXTNM(mr0,mr);
+                                 c = -1;
-                         c1 = C(m);
+                         else
-                         DMAR(c1,c,0,nd_mod,c2);
+                                 c = ndl_compare(DL(m1),DL(m2));
-                         C(mr) = c2;
+                         switch ( c ) {
-                         mr->td = m->td+td;
+                                 case 0:
-                         ndl_add(m->dl,d,mr->dl);
+                                         addq(CQ(m1),CQ(m2),&t);
+                                         s = m1; m1 = NEXT(m1);
+                                         if ( t ) {
+                                                 NEXTNM2(mr0,mr,s); CQ(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;
+                         }
                  }
-                 NEXT(mr) = 0;
+                 if ( !mr0 )
-                 MKND(NV(p),mr0,r);
+                         if ( m1 )
-                 r->sugar = p->sugar + td;
+                                 mr0 = m1;
-                 return r;
+                         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;
+                 SG(p1) = MAX(SG(p1),SG(p2));
+                 FREEND(p2);
+                 return p1;
          }
  }
- ND nd_mul_ind_nm(int index,NM m0)
+ #if 1
+ /* ret=1 : success, ret=0 : overflow */
+ int nd_nf(int mod,ND g,int full,ND *rp)
  {
-         register int c1,c2,c;
+         ND d;
-         register NM m,new,prev;
+         NM m,mrd,tail;
-         NM mr0;
+         NM mul;
-         unsigned int *d;
+         int n,sugar,psugar,sugar0,stat,index;
-         int n,td,i,len,d0,d1;
+         int c,c1,c2;
-         ND p,r;
+         RHist h;
+         NDV p,red;
+         Q cg,cred,gcd;
+         double hmag;
-         p = nd_ps[index];
+         if ( !g ) {
-         len = nd_psl[index];
+                 *rp = 0;
-         n = NV(p); m = BDY(p);
+                 return 1;
-         d = m0->dl; td = m0->td; c = C(m0);
+         }
+         if ( !mod )
+                 hmag = ((double)p_mag((P)HCQ(g)))*nd_scale;
-         NEWNM(mr0);
+         sugar0 = sugar = SG(g);
-         c1 = C(m); DMAR(c1,c,0,nd_mod,c2); C(mr0) = c2;
+         n = NV(g);
-         mr0->td = m->td+td; ndl_add(m->dl,d,mr0->dl);
+         mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
-         prev = mr0; m = NEXT(m);
+         for ( d = 0; g; ) {
-         len--;
+                 index = nd_find_reducer(g);
+                 if ( index >= 0 ) {
-         switch ( nd_wpd ) {
+                         h = nd_psh[index];
-                 case 1:
+                         ndl_sub(HDL(g),DL(h),DL(mul));
-                         d0 = d[0];
+                         TD(mul) = HTD(g)-TD(h);
-                         while ( len-- ) {
+ #if 0
-                                 c1 = C(m); DMAR(c1,c,0,nd_mod,c2);
+                         if ( d && (SG(p)+TD(mul)) > sugar ) {
-                                 NEWNM(new); C(new) = c2;
+                                 goto afo;
-                                 new->td = m->td+td; new->dl[0] = m->dl[0]+d0;
-                                 m = NEXT(m); NEXT(prev) = new; prev = new;
                          }
-                         break;
+ #endif
-                 case 2:
+                         if ( ndl_check_bound2(index,DL(mul)) ) {
-                         d0 = d[0]; d1 = d[1];
+                                 nd_free(g); nd_free(d);
-                         while ( len-- ) {
+                                 return 0;
-                                 c1 = C(m); DMAR(c1,c,0,nd_mod,c2);
-                                 NEWNM(new); C(new) = c2;
-                                 new->td = m->td+td;
-                                 new->dl[0] = m->dl[0]+d0;
-                                 new->dl[1] = m->dl[1]+d1;
-                                 m = NEXT(m); NEXT(prev) = new; prev = new;
                          }
-                         break;
+                         if ( mod ) {
-                 default:
+                                 p = nd_ps[index];
-                         while ( len-- ) {
+                                 c1 = invm(HCM(p),mod); c2 = mod-HCM(g);
-                                 c1 = C(m); DMAR(c1,c,0,nd_mod,c2);
+                                 DMAR(c1,c2,0,mod,c); CM(mul) = c;
-                                 NEWNM(new); C(new) = c2;
+                         } else {
-                                 new->td = m->td+td; ndl_add(m->dl,d,new->dl);
+                                 p = nd_psq[index];
-                                 m = NEXT(m); NEXT(prev) = new; prev = new;
+                                 igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);
+                                 chsgnq(cg,&CQ(mul));
+                                 nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);
                          }
-                         break;
+                         g = nd_add(mod,g,ndv_mul_nm(mod,p,mul));
+                         sugar = MAX(sugar,SG(p)+TD(mul));
+                         if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {
+                                 nd_removecont2(d,g);
+                                 hmag = ((double)p_mag((P)HCQ(g)))*nd_scale;
+                         }
+                 } else if ( !full ) {
+                         *rp = g;
+                         return 1;
+                 } else {
+ afo:
+                         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 )
-         NEXT(prev) = 0;
+                 SG(d) = sugar;
-         MKND(NV(p),mr0,r);
+         *rp = d;
-         r->sugar = p->sugar + td;
+         return 1;
-         return r;
  }
+ #else
- ND nd_mul_term(ND p,int td,unsigned int *d)
+ int nd_nf(int mod,ND g,int full,ND *rp)
  {
-         NM m,mr,mr0;
+         int hindex,index;
-         int c,n;
+         NDV p;
-         ND r;
+         ND u,d,red;
+         NODE l;
+         NM mul,m,mrd;
+         int sugar,psugar,n,h_reducible;
+         PGeoBucket bucket;
+         int c,c1,c2;
+         Q cg,cred,gcd,zzz;
+         RHist h;
-         if ( !p )
+         if ( !g ) {
-                 return 0;
+                 *rp = 0;
-         else {
+                 return 1;
-                 n = NV(p); m = BDY(p);
+         }
-                 for ( mr0 = 0; m; m = NEXT(m) ) {
+         sugar = SG(g);
-                         NEXTNM(mr0,mr);
+         n = NV(g);
-                         C(mr) = C(m);
+         bucket = create_pbucket();
-                         mr->td = m->td+td;
+         add_pbucket(mod,bucket,g,nd_length(g));
-                         ndl_add(m->dl,d,mr->dl);
+         d = 0;
+         mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
+         while ( 1 ) {
+                 hindex = mod?head_pbucket(mod,bucket):head_pbucket_q(bucket);
+                 if ( hindex < 0 ) {
+                         if ( d )
+                                 SG(d) = sugar;
+                         *rp = d;
+                         return 1;
                  }
-                 NEXT(mr) = 0;
+                 g = bucket->body[hindex];
-                 MKND(NV(p),mr0,r);
+                 index = nd_find_reducer(g);
-                 r->sugar = p->sugar + td;
+                 if ( index >= 0 ) {
-                 return r;
+                         h = nd_psh[index];
+                         ndl_sub(HDL(g),DL(h),DL(mul));
+                         TD(mul) = HTD(g)-TD(h);
+                         if ( ndl_check_bound2(index,DL(mul)) ) {
+                                 nd_free(d);
+                                 free_pbucket(bucket);
+                                 *rp = 0;
+                                 return 0;
+                         }
+                         if ( mod ) {
+                                 p = nd_ps[index];
+                                 c1 = invm(HCM(p),mod); c2 = mod-HCM(g);
+                                 DMAR(c1,c2,0,mod,c); CM(mul) = c;
+                         } else {
+                                 p = nd_psq[index];
+                                 igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);
+                                 chsgnq(cg,&CQ(mul));
+                                 nd_mul_c_q(d,cred);
+                                 mulq_pbucket(bucket,cred);
+                                 g = bucket->body[hindex];
+                         }
+                         red = ndv_mul_nm(mod,p,mul);
+                         bucket->body[hindex] = nd_remove_head(g);
+                         red = nd_remove_head(red);
+                         add_pbucket(mod,bucket,red,LEN(p)-1);
+                         sugar = MAX(sugar,SG(p)+TD(mul));
+                 } else if ( !full ) {
+                         g = normalize_pbucket(mod,bucket);
+                         if ( g )
+                                 SG(g) = sugar;
+                         *rp = g;
+                         return 1;
+                 } else {
+                         m = BDY(g);
+                         if ( NEXT(m) ) {
+                                 BDY(g) = NEXT(m); NEXT(m) = 0;
+                         } else {
+                                 FREEND(g); g = 0;
+                         }
+                         bucket->body[hindex] = g;
+                         NEXT(m) = 0;
+                         if ( d ) {
+                                 for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
+                                 NEXT(mrd) = m;
+                         } else {
+                                 MKND(n,m,d);
+                         }
+                 }
          }
  }
+ #endif
- #if 1
+ int nd_nf_direct(int mod,ND g,NDV *ps,int len,int full,ND *rp)
- /* ret=1 : success, ret=0 : overflow */
- int nd_nf(ND g,int full,ND *rp)
  {
-         ND p,d;
+         ND d;
          NM m,mrd,tail;
-         struct oNM mul;
+         NM mul;
-         int n,sugar,psugar,stat,index;
+         int n,sugar,psugar,sugar0,stat,index;
          int c,c1,c2;
- #if USE_NDV
+         RHist h;
-         NDV red;
+         NDV p,red;
- #else
+         Q cg,cred,gcd;
-         ND red;
+         double hmag;
- #endif
          if ( !g ) {
                  *rp = 0;
                  return 1;
          }
-         sugar = g->sugar;
+ #if 0
+         if ( !mod )
+                 hmag = ((double)p_mag((P)HCQ(g)))*nd_scale;
+ #else
+         /* XXX */
+         hmag = 0;
+ #endif
+         sugar0 = sugar = SG(g);
          n = NV(g);
+         mul = (NM)ALLOCA(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
          for ( d = 0; g; ) {
-                 index = nd_find_reducer(g);
+                 index = nd_find_reducer_direct(g,ps,len);
                  if ( index >= 0 ) {
-                         p = nd_ps[index];
+                         p = ps[index];
-                         ndl_sub(HDL(g),HDL(p),mul.dl);
+                         ndl_sub(HDL(g),HDL(p),DL(mul));
-                         mul.td = HTD(g)-HTD(p);
+                         TD(mul) = HTD(g)-HTD(p);
-                         if ( ndl_check_bound2(index,mul.dl) ) {
+                         if ( ndl_check_bound2_direct(HDL(p),DL(mul)) ) {
                                  nd_free(g); nd_free(d);
                                  return 0;
                          }
-                         c1 = invm(HC(p),nd_mod); c2 = nd_mod-HC(g);
+                         if ( mod ) {
-                         DMAR(c1,c2,0,nd_mod,c); mul.c = c;
+                                 c1 = invm(HCM(p),mod); c2 = mod-HCM(g);
- #if USE_NDV
+                                 DMAR(c1,c2,0,mod,c); CM(mul) = c;
-                         ndv_mul_nm(nd_psv[index],&mul,ndv_red);
+                         } else {
-                         g = ndv_add(g,ndv_red);
+                                 igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);
-                         sugar = MAX(sugar,ndv_red->sugar);
+                                 chsgnq(cg,&CQ(mul));
- #else
+                                 nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);
-                         red = nd_mul_ind_nm(index,&mul);
+                         }
-                         g = nd_add(g,red);
+                         g = nd_add(mod,g,ndv_mul_nm(mod,p,mul));
-                         sugar = MAX(sugar,red->sugar);
+                         sugar = MAX(sugar,SG(p)+TD(mul));
- #endif
+                         if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {
+                                 nd_removecont2(d,g);
+                                 hmag = ((double)p_mag((P)HCQ(g)))*nd_scale;
+                         }
                  } else if ( !full ) {
                          *rp = g;
                          return 1;
-Line 1086  int nd_nf(ND g,int full,ND *rp)
+Line 1155  int nd_nf(ND g,int full,ND *rp)
 Line 1086  int nd_nf(ND g,int full,ND *rp)
 Line 1155  int nd_nf(ND g,int full,ND *rp)
                  }
          }
          if ( d )
-                 d->sugar = sugar;
+                 SG(d) = sugar;
          *rp = d;
          return 1;
  }
- #else
  ND nd_remove_head(ND p)
  {
-Line 1115  PGeoBucket create_pbucket()
+Line 1183  PGeoBucket create_pbucket()
 Line 1115  PGeoBucket create_pbucket()
 Line 1183  PGeoBucket create_pbucket()
          return g;
  }
- void add_pbucket(PGeoBucket g,ND d)
+ void free_pbucket(PGeoBucket b) {
+         int i;
+         for ( i = 0; i <= b->m; i++ )
+                 if ( b->body[i] ) {
+                         nd_free(b->body[i]);
+                         b->body[i] = 0;
+                 }
+         GC_free(b);
+ }
+ void add_pbucket(int mod,PGeoBucket g,ND d,int l)
  {
-         int l,k,m;
+         int 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);
+         d = nd_add(mod,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);
+                 d = nd_add(mod,g->body[k+1],d);
          }
          g->body[k] = d;
          g->m = MAX(g->m,k);
  }
- int head_pbucket(PGeoBucket g)
+ void mulq_pbucket(PGeoBucket g,Q c)
  {
+         int k;
+         for ( k = 0; k <= g->m; k++ )
+                 nd_mul_c_q(g->body[k],c);
+ }
+ /* XXX not completed */
+ int head_pbucket(int mod,PGeoBucket g)
+ {
          int j,i,c,k,nv,sum;
          unsigned int *di,*dj;
          ND gi,gj;
-Line 1147  int head_pbucket(PGeoBucket g)
+Line 1234  int head_pbucket(PGeoBucket g)
 Line 1147  int head_pbucket(PGeoBucket g)
 Line 1234  int head_pbucket(PGeoBucket g)
                                  j = i;
                                  gj = g->body[j];
                                  dj = HDL(gj);
-                                 sum = HC(gj);
+                                 sum = HCM(gj);
                          } else {
                                  di = HDL(gi);
                                  nv = NV(gi);
-Line 1159  int head_pbucket(PGeoBucket g)
+Line 1246  int head_pbucket(PGeoBucket g)
 Line 1159  int head_pbucket(PGeoBucket g)
 Line 1246  int head_pbucket(PGeoBucket g)
                                          c = ndl_compare(di,dj);
                                  if ( c > 0 ) {
                                          if ( sum )
-                                                 HC(gj) = sum;
+                                                 HCM(gj) = sum;
                                          else
                                                  g->body[j] = nd_remove_head(gj);
                                          j = i;
                                          gj = g->body[j];
                                          dj = HDL(gj);
-                                         sum = HC(gj);
+                                         sum = HCM(gj);
                                  } else if ( c == 0 ) {
-                                         sum = sum+HC(gi)-nd_mod;
+                                         sum = sum+HCM(gi)-mod;
                                          if ( sum < 0 )
-                                                 sum += nd_mod;
+                                                 sum += mod;
                                          g->body[i] = nd_remove_head(gi);
                                  }
                          }
-Line 1177  int head_pbucket(PGeoBucket g)
+Line 1264  int head_pbucket(PGeoBucket g)
 Line 1177  int head_pbucket(PGeoBucket g)
 Line 1264  int head_pbucket(PGeoBucket g)
                  if ( j < 0 )
                          return -1;
                  else if ( sum ) {
-                         HC(gj) = sum;
+                         HCM(gj) = sum;
                          return j;
                  } else
                          g->body[j] = nd_remove_head(gj);
          }
  }
- ND normalize_pbucket(PGeoBucket g)
+ int head_pbucket_q(PGeoBucket g)
  {
+         int j,i,c,k,nv;
+         Q sum,t;
+         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 = HCQ(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 )
+                                                 HCQ(gj) = sum;
+                                         else
+                                                 g->body[j] = nd_remove_head(gj);
+                                         j = i;
+                                         gj = g->body[j];
+                                         dj = HDL(gj);
+                                         sum = HCQ(gj);
+                                 } else if ( c == 0 ) {
+                                         addq(sum,HCQ(gi),&t);
+                                         sum = t;
+                                         g->body[i] = nd_remove_head(gi);
+                                 }
+                         }
+                 }
+                 if ( j < 0 )
+                         return -1;
+                 else if ( sum ) {
+                         HCQ(gj) = sum;
+                         return j;
+                 } else
+                         g->body[j] = nd_remove_head(gj);
+         }
+ }
+ ND normalize_pbucket(int mod,PGeoBucket g)
+ {
          int i;
          ND r,t;
          r = 0;
          for ( i = 0; i <= g->m; i++ )
-                 r = nd_add(r,g->body[i]);
+                 r = nd_add(mod,r,g->body[i]);
          return r;
  }
- ND nd_nf(ND g,int full)
+ NODE nd_gb(int m)
  {
-         ND u,p,d,red;
+         int i,nh,sugar,stat;
-         NODE l;
+         NODE r,g,t;
-         NM m,mrd;
+         ND_pairs d;
-         int sugar,psugar,n,h_reducible,h;
+         ND_pairs l;
-         PGeoBucket bucket;
+         ND h,nf;
-         if ( !g ) {
+         g = 0; d = 0;
-                 return 0;
+         for ( i = 0; i < nd_psn; i++ ) {
+                 d = update_pairs(d,g,i);
+                 g = update_base(g,i);
          }
-         sugar = g->sugar;
+         sugar = 0;
-         n = g->nv;
+         while ( d ) {
-         bucket = create_pbucket();
+ again:
-         add_pbucket(bucket,g);
+                 l = nd_minp(d,&d);
-         d = 0;
+                 if ( SG(l) != sugar ) {
-         while ( 1 ) {
+                         sugar = SG(l);
-                 h = head_pbucket(bucket);
+                         fprintf(asir_out,"%d",sugar);
-                 if ( h < 0 ) {
-                         if ( d )
-                                 d->sugar = sugar;
-                         return d;
                  }
-                 g = bucket->body[h];
+                 stat = nd_sp(m,l,&h);
-                 red = nd_find_reducer(g);
+                 if ( !stat ) {
-                 if ( red ) {
+                         NEXT(l) = d; d = l;
-                         bucket->body[h] = nd_remove_head(g);
+                         d = nd_reconstruct(m,0,d);
-                         red = nd_remove_head(red);
+                         goto again;
-                         add_pbucket(bucket,red);
+                 }
-                         sugar = MAX(sugar,red->sugar);
+                 stat = nd_nf(m,h,!Top,&nf);
-                 } else if ( !full ) {
+                 if ( !stat ) {
-                         g = normalize_pbucket(bucket);
+                         NEXT(l) = d; d = l;
-                         if ( g )
+                         d = nd_reconstruct(m,0,d);
-                                 g->sugar = sugar;
+                         goto again;
-                         return g;
+                 } else if ( nf ) {
+                         printf("+"); fflush(stdout);
+                         nh = nd_newps(m,nf);
+                         d = update_pairs(d,g,nh);
+                         g = update_base(g,nh);
+                         FREENDP(l);
                  } else {
-                         m = BDY(g);
+                         printf("."); fflush(stdout);
-                         if ( NEXT(m) ) {
+                         FREENDP(l);
-                                 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);
-                         }
                  }
          }
+         if ( m )
+                 for ( t = g; t; t = NEXT(t) )
+                         BDY(t) = (pointer)nd_ps[(int)BDY(t)];
+         else
+                 for ( t = g; t; t = NEXT(t) )
+                         BDY(t) = (pointer)nd_psq[(int)BDY(t)];
+         return g;
  }
- #endif
- NODE nd_gb(NODE f)
+ NODE nd_gb_trace(int m)
  {
          int i,nh,sugar,stat;
-         NODE r,g,gall;
+         NODE r,g,t;
          ND_pairs d;
          ND_pairs l;
-         ND h,nf;
+         ND h,nf,nfq;
-         for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
+         g = 0; d = 0;
-                 i = (int)BDY(r);
+         for ( i = 0; i < nd_psn; i++ ) {
                  d = update_pairs(d,g,i);
                  g = update_base(g,i);
-                 gall = append_one(gall,i);
          }
          sugar = 0;
          while ( d ) {
  again:
                  l = nd_minp(d,&d);
-                 if ( l->sugar != sugar ) {
+                 if ( SG(l) != sugar ) {
-                         sugar = l->sugar;
+                         sugar = SG(l);
                          fprintf(asir_out,"%d",sugar);
                  }
-                 stat = nd_sp(l,&h);
+                 stat = nd_sp(m,l,&h);
                  if ( !stat ) {
                          NEXT(l) = d; d = l;
-                         d = nd_reconstruct(d);
+                         d = nd_reconstruct(m,1,d);
                          goto again;
                  }
-                 stat = nd_nf(h,!Top,&nf);
+                 stat = nd_nf(m,h,!Top,&nf);
                  if ( !stat ) {
                          NEXT(l) = d; d = l;
-                         d = nd_reconstruct(d);
+                         d = nd_reconstruct(m,1,d);
                          goto again;
                  } else if ( nf ) {
-                         printf("+"); fflush(stdout);
+                         /* overflow does not occur */
-                         nh = nd_newps(nf);
+                         nd_sp(0,l,&h);
-                         d = update_pairs(d,g,nh);
+                         nd_nf(0,h,!Top,&nfq);
-                         g = update_base(g,nh);
+                         if ( nfq ) {
-                         gall = append_one(gall,nh);
+                                 printf("+"); fflush(stdout);
-                         FREENDP(l);
+                                 nh = nd_newps_trace(m,nf,nfq);
+                                 d = update_pairs(d,g,nh);
+                                 g = update_base(g,nh);
+                         } else {
+                                 printf("*"); fflush(stdout);
+                         }
                  } else {
                          printf("."); fflush(stdout);
-                         FREENDP(l);
                  }
+                 FREENDP(l);
          }
+         for ( t = g; t; t = NEXT(t) )
+                 BDY(t) = (pointer)nd_psq[(int)BDY(t)];
          return g;
  }
+ int ndv_compare(NDV *p1,NDV *p2)
+ {
+         int td1,td2;
+         td1 = HTD(*p1); td2 = HTD(*p2);
+         if ( td1 > td2 ) return 1;
+         else if ( td1 < td2 ) return -1;
+         else return ndl_compare(HDL(*p1),HDL(*p2));
+ }
+ int ndv_compare_rev(NDV *p1,NDV *p2)
+ {
+         return -ndv_compare(p1,p2);
+ }
+ NODE nd_reduceall(int m,NODE f)
+ {
+         int i,j,n,stat;
+         NDV *w,*ps;
+         ND nf,g;
+         NODE t,a0,a;
+         for ( n = 0, t = f; t; t = NEXT(t), n++ );
+         ps = (NDV *)ALLOCA(n*sizeof(NDV));
+         for ( i = 0, t = f; i < n; i++, t = NEXT(t) )
+                 ps[i] = (NDV)BDY(t);
+         qsort(ps,n,sizeof(NDV),(int (*)(const void *,const void *))ndv_compare);
+         w = (NDV *)ALLOCA((n-1)*sizeof(NDV));
+         for ( i = 0; i < n; i++ ) {
+                 for ( j = 0; j < i; j++ )
+                         w[j] = (NDV)ps[j];
+                 for ( j = i+1; j < n; j++ )
+                         w[j-1] = ps[j];
+                 g = ndvtond(m,ps[i]);
+                 stat = nd_nf_direct(m,g,w,n-1,1,&nf);
+                 if ( !stat )
+                         nd_reconstruct_direct(m,ps,n);
+                 else if ( !nf ) {
+                         printf("."); fflush(stdout);
+                         ndv_free(ps[i]);
+                         for ( j = i+1; j < n; j++ )
+                                 ps[j-1] = ps[j];
+                         n--;
+                 } else {
+                         printf("."); fflush(stdout);
+                         ndv_free(ps[i]);
+                         nd_removecont(m,nf);
+                         ps[i] = ndtondv(m,nf);
+                         nd_free(nf);
+                 }
+         }
+         for ( a0 = 0, i = 0; i < n; i++ ) {
+                 NEXTNODE(a0,a);
+                 BDY(a) = (pointer)ps[i];
+         }
+         NEXT(a) = 0;
+         return a0;
+ }
  ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t)
  {
          ND_pairs d1,nd,cur,head,prev,remove;
-Line 1340  ND_pairs nd_newpairs( NODE g, int t )
+Line 1542  ND_pairs nd_newpairs( NODE g, int t )
 Line 1340  ND_pairs nd_newpairs( NODE g, int t )
 Line 1542  ND_pairs nd_newpairs( NODE g, int t )
          int td,ts,s;
          ND_pairs r,r0;
-         dl = HDL(nd_ps[t]);
+         dl = DL(nd_psh[t]);
-         td = HTD(nd_ps[t]);
+         td = TD(nd_psh[t]);
-         ts = nd_ps[t]->sugar - td;
+         ts = SG(nd_psh[t]) - 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(nd_ps[r->i1]),dl,r->lcm);
+                 ndl_lcm(DL(nd_psh[r->i1]),dl,r->lcm);
-                 r->td = ndl_td(r->lcm);
+                 TD(r) = ndl_td(r->lcm);
-                 s = nd_ps[r->i1]->sugar-HTD(nd_ps[r->i1]);
+                 s = SG(nd_psh[r->i1])-TD(nd_psh[r->i1]);
-                 r->sugar = MAX(s,ts) + r->td;
+                 SG(r) = MAX(s,ts) + TD(r);
          }
          NEXT(r) = 0;
          return r0;
-Line 1363  ND_pairs crit_B( ND_pairs d, int s )
+Line 1565  ND_pairs crit_B( ND_pairs d, int s )
 Line 1363  ND_pairs crit_B( ND_pairs d, int s )
 Line 1565  ND_pairs crit_B( ND_pairs d, int s )
          int td,tdl;
          if ( !d ) return 0;
-         t = HDL(nd_ps[s]);
+         t = DL(nd_psh[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(nd_ps[cur->i1]),t,lcm),!ndl_equal(lcm,tl))
+                         && (ndl_lcm(DL(nd_psh[cur->i1]),t,lcm),!ndl_equal(lcm,tl))
-                         && (ndl_lcm(HDL(nd_ps[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) {
+                         && (ndl_lcm(DL(nd_psh[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) {
                          remove = cur;
                          if ( !prev ) {
                                  head = cur = NEXT(cur);
-Line 1399  ND_pairs crit_M( ND_pairs d1 )
+Line 1601  ND_pairs crit_M( ND_pairs d1 )
 Line 1399  ND_pairs crit_M( ND_pairs d1 )
 Line 1601  ND_pairs crit_M( ND_pairs d1 )
                          return e;
                  }
                  id = e->lcm;
-                 itd = e->td;
+                 itd = TD(e);
                  for ( d3 = 0; d2; d2 = p ) {
                          p = NEXT(d2),
                          jd = d2->lcm;
-                         jtd = d2->td;
+                         jtd = TD(d2);
                          if ( jtd == itd  )
                                  if ( id == jd );
                                  else if ( ndl_reducible(jd,id) ) continue;
-Line 1440  ND_pairs crit_F( ND_pairs d1 )
+Line 1642  ND_pairs crit_F( ND_pairs d1 )
 Line 1440  ND_pairs crit_F( ND_pairs d1 )
 Line 1642  ND_pairs crit_F( ND_pairs d1 )
          for ( head = last = 0, p = d1; NEXT(p); ) {
                  r = w = equivalent_pairs(p,&rest);
-                 s = r->sugar;
+                 s = SG(r);
                  w = NEXT(w);
                  while ( w ) {
                          if ( crit_2(w->i1,w->i2) ) {
-Line 1452  ND_pairs crit_F( ND_pairs d1 )
+Line 1654  ND_pairs crit_F( ND_pairs d1 )
 Line 1452  ND_pairs crit_F( ND_pairs d1 )
 Line 1654  ND_pairs crit_F( ND_pairs d1 )
                                          FREENDP(remove);
                                  }
                                  break;
-                         } else if ( w->sugar < s ) {
+                         } else if ( SG(w) < s ) {
                                  FREENDP(r);
                                  r = w;
-                                 s = r->sugar;
+                                 s = SG(r);
                                  w = NEXT(w);
                          } else {
                                  remove = w;
-Line 1476  ND_pairs crit_F( ND_pairs d1 )
+Line 1678  ND_pairs crit_F( ND_pairs d1 )
 Line 1476  ND_pairs crit_F( ND_pairs d1 )
 Line 1678  ND_pairs crit_F( ND_pairs d1 )
  int crit_2( int dp1, int dp2 )
  {
-         return ndl_disjoint(HDL(nd_ps[dp1]),HDL(nd_ps[dp2]));
+         return ndl_disjoint(DL(nd_psh[dp1]),DL(nd_psh[dp2]));
  }
  static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest )
-Line 1487  static ND_pairs equivalent_pairs( ND_pairs d1, ND_pair
+Line 1689  static ND_pairs equivalent_pairs( ND_pairs d1, ND_pair
 Line 1487  static ND_pairs equivalent_pairs( ND_pairs d1, ND_pair
 Line 1689  static ND_pairs equivalent_pairs( ND_pairs d1, ND_pair
          w = d1;
          d = w->lcm;
-         td = w->td;
+         td = TD(w);
          s = NEXT(w);
          NEXT(w) = 0;
          for ( r = 0; s; s = p ) {
                  p = NEXT(s);
-                 if ( td == s->td && ndl_equal(d,s->lcm) ) {
+                 if ( td == TD(s) && ndl_equal(d,s->lcm) ) {
                          NEXT(s) = w;
                          w = s;
                  } else {
-Line 1510  NODE update_base(NODE nd,int ndp)
+Line 1712  NODE update_base(NODE nd,int ndp)
 Line 1510  NODE update_base(NODE nd,int ndp)
 Line 1712  NODE update_base(NODE nd,int ndp)
          NODE last, p, head;
          int td,tdn;
-         dl = HDL(nd_ps[ndp]);
+         dl = DL(nd_psh[ndp]);
-         td = HTD(nd_ps[ndp]);
+         td = TD(nd_psh[ndp]);
          for ( head = last = 0, p = nd; p; ) {
-                 dln = HDL(nd_ps[(int)BDY(p)]);
+                 dln = DL(nd_psh[(int)BDY(p)]);
-                 tdn = HTD(nd_ps[(int)BDY(p)]);
+                 tdn = TD(nd_psh[(int)BDY(p)]);
                  if ( tdn >= td && ndl_reducible( dln, dl ) ) {
                          p = NEXT(p);
                          if ( last ) NEXT(last) = p;
-Line 1539  ND_pairs nd_minp( ND_pairs d, ND_pairs *prest )
+Line 1741  ND_pairs nd_minp( ND_pairs d, ND_pairs *prest )
 Line 1539  ND_pairs nd_minp( ND_pairs d, ND_pairs *prest )
 Line 1741  ND_pairs nd_minp( ND_pairs d, ND_pairs *prest )
                  return m;
          }
          lcm = m->lcm;
-         s = m->sugar;
+         s = SG(m);
-         td = m->td;
+         td = TD(m);
          len = nd_psl[m->i1]+nd_psl[m->i2];
          for ( ml = 0, l = m; p; p = NEXT(l = p) ) {
-                 if (p->sugar < s)
+                 if (SG(p) < s)
                          goto find;
-                 else if ( p->sugar == s ) {
+                 else if ( SG(p) == s ) {
-                         if ( p->td < td )
+                         if ( TD(p) < td )
                                  goto find;
-                         else if ( p->td == td ) {
+                         else if ( TD(p) == td ) {
                                  c = ndl_compare(p->lcm,lcm);
                                  if ( c < 0 )
                                          goto find;
+ #if 0
                                  else if ( c == 0 ) {
                                          tlen = nd_psl[p->i1]+nd_psl[p->i2];
                                          if ( tlen < len )
                                                  goto find;
                                  }
+ #endif
                          }
                  }
                  continue;
-Line 1564  find:
+Line 1768  find:
 Line 1564  find:
 Line 1768  find:
                  ml = l;
                  m = p;
                  lcm = m->lcm;
-                 s = m->sugar;
+                 s = SG(m);
-                 td = m->td;
+                 td = TD(m);
                  len = tlen;
          }
          if ( !ml ) *prest = NEXT(m);
-Line 1577  find:
+Line 1781  find:
 Line 1577  find:
 Line 1781  find:
          return m;
  }
- int nd_newps(ND a)
+ int nd_newps(int mod,ND a)
  {
          int len;
+         RHist r;
+         NDV b;
          if ( nd_psn == nd_pslen ) {
                  nd_pslen *= 2;
-                 nd_ps = (ND *)REALLOC((char *)nd_ps,nd_pslen*sizeof(ND));
                  nd_psl = (int *)REALLOC((char *)nd_psl,nd_pslen*sizeof(int));
- #if USE_NDV
+                 nd_ps = (NDV *)REALLOC((char *)nd_ps,nd_pslen*sizeof(NDV));
-                 nd_psv = (NDV *)REALLOC((char *)nd_psv,nd_pslen*sizeof(NDV));
+                 nd_psq = (NDV *)REALLOC((char *)nd_psq,nd_pslen*sizeof(NDV));
- #endif
+                 nd_psh = (RHist *)REALLOC((char *)nd_psh,nd_pslen*sizeof(RHist));
                  nd_bound = (unsigned int **)
                          REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *));
          }
-         nd_monic(a);
+         nd_removecont(mod,a);
-         nd_ps[nd_psn] = a;
-         nd_psl[nd_psn] = nd_length(a);
          nd_bound[nd_psn] = nd_compute_bound(a);
- #if USE_NDV
+         NEWRHist(r); SG(r) = SG(a); TD(r) = HTD(a); ndl_copy(HDL(a),DL(r));
-         nd_psv[nd_psn] = ndtondv(a);
+         nd_psh[nd_psn] = r;
-         len = nd_psv[nd_psn]->len;
+         b = ndtondv(mod,a);
-         if ( len > nmv_len ) {
+         len = LEN(b);
-                 nmv_len = 2*len;
+         if ( mod )
-                 BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
+                 nd_ps[nd_psn] = b;
-         }
+         else
- #endif
+                 nd_psq[nd_psn] = b;
+         nd_psl[nd_psn] = len;
+         nd_free(a);
          return nd_psn++;
  }
- NODE NODE_sortb(NODE f,int);
+ int nd_newps_trace(int mod,ND nf,ND nfq)
- ND dptond(DP);
+ {
- DP ndtodp(ND);
+         int len;
+         RHist r;
+         NDV b;
- NODE nd_setup(NODE f)
+         if ( nd_psn == nd_pslen ) {
+                 nd_pslen *= 2;
+                 nd_psl = (int *)REALLOC((char *)nd_psl,nd_pslen*sizeof(int));
+                 nd_ps = (NDV *)REALLOC((char *)nd_ps,nd_pslen*sizeof(NDV));
+                 nd_psq = (NDV *)REALLOC((char *)nd_psq,nd_pslen*sizeof(NDV));
+                 nd_psh = (RHist *)REALLOC((char *)nd_psh,nd_pslen*sizeof(RHist));
+                 nd_bound = (unsigned int **)
+                         REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *));
+         }
+         nd_removecont(mod,nf);
+         nd_ps[nd_psn] = ndtondv(mod,nf);
+         nd_removecont(0,nfq);
+         nd_psq[nd_psn] = ndtondv(0,nfq);
+         nd_bound[nd_psn] = nd_compute_bound(nfq);
+         NEWRHist(r); SG(r) = SG(nf); TD(r) = HTD(nf); ndl_copy(HDL(nf),DL(r));
+         nd_psh[nd_psn] = r;
+         len = LEN(nd_psq[nd_psn]);
+         nd_psl[nd_psn] = len;
+         nd_free(nf); nd_free(nfq);
+         return nd_psn++;
+ }
+ void nd_setup(int mod,NODE f)
  {
          int i,j,td,len,max;
          NODE s,s0,f0;
          unsigned int *d;
+         RHist r;
+         NDV a;
-         nd_found = 0;
+         nd_found = 0; nd_notfirst = 0; nd_create = 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;
-         nd_ps = (ND *)MALLOC(nd_pslen*sizeof(ND));
          nd_psl = (int *)MALLOC(nd_pslen*sizeof(int));
- #if USE_NDV
+         nd_ps = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
-         nd_psv = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
+         nd_psq = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
- #endif
+         nd_psh = (RHist *)MALLOC(nd_pslen*sizeof(RHist));
          nd_bound = (unsigned int **)MALLOC(nd_pslen*sizeof(unsigned int *));
          for ( max = 0, i = 0, s = f; i < nd_psn; i++, s = NEXT(s) ) {
                  nd_bound[i] = d = dp_compute_bound((DP)BDY(s));
                  for ( j = 0; j < nd_nvar; j++ )
                          max = MAX(d[j],max);
          }
+         if ( !nd_red )
+                 nd_red = (RHist *)MALLOC(REDTAB_LEN*sizeof(RHist));
+         bzero(nd_red,REDTAB_LEN*sizeof(RHist));
          if ( max < 2 )
                  nd_bpe = 2;
-Line 1647  NODE nd_setup(NODE f)
+Line 1880  NODE nd_setup(NODE f)
 Line 1647  NODE nd_setup(NODE f)
 Line 1880  NODE nd_setup(NODE f)
                  nd_bpe = 16;
          else
                  nd_bpe = 32;
          nd_setup_parameters();
          nd_free_private_storage();
          len = 0;
          for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) {
-                 nd_ps[i] = dptond((DP)BDY(f));
+                 NEWRHist(r);
-                 nd_monic(nd_ps[i]);
+                 a = dptondv(mod,(DP)BDY(f));
-                 nd_psl[i] = nd_length(nd_ps[i]);
+                 ndv_removecont(mod,a);
- #if USE_NDV
+                 len = MAX(len,LEN(a));
-                 nd_psv[i] = ndtondv(nd_ps[i]);
+                 SG(r) = HTD(a); TD(r) = HTD(a); ndl_copy(HDL(a),DL(r));
-                 len = MAX(len,nd_psv[i]->len);
+                 if ( mod )
- #endif
+                         nd_ps[i] = a;
+                 else
+                         nd_psq[i] = a;
+                 nd_psh[i] = r;
          }
-         nd_red = (NM *)MALLOC(REDTAB_LEN*sizeof(NM));
+ }
- #if USE_NDV
-         nmv_len = 16*len;
+ void nd_setup_trace(int mod,NODE f)
-         NEWNDV(ndv_red);
+ {
-         BDY(ndv_red) = (NMV)MALLOC_ATOMIC(nmv_len*nmv_adv);
+         int i,j,td,len,max;
- #endif
+         NODE s,s0,f0;
-         for ( s0 = 0, i = 0; i < nd_psn; i++ ) {
+         unsigned int *d;
-                 NEXTNODE(s0,s); BDY(s) = (pointer)i;
+         RHist r;
+         NDV a;
+         nd_found = 0; nd_notfirst = 0; nd_create = 0;
+         nd_psn = length(f); nd_pslen = 2*nd_psn;
+         nd_psl = (int *)MALLOC(nd_pslen*sizeof(int));
+         nd_ps = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
+         nd_psq = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
+         nd_psh = (RHist *)MALLOC(nd_pslen*sizeof(RHist));
+         nd_bound = (unsigned int **)MALLOC(nd_pslen*sizeof(unsigned int *));
+         for ( max = 0, i = 0, s = f; i < nd_psn; i++, s = NEXT(s) ) {
+                 nd_bound[i] = d = dp_compute_bound((DP)BDY(s));
+                 for ( j = 0; j < nd_nvar; j++ )
+                         max = MAX(d[j],max);
          }
-         if ( s0 ) NEXT(s) = 0;
+         if ( !nd_red )
-         return s0;
+                 nd_red = (RHist *)MALLOC(REDTAB_LEN*sizeof(RHist));
+         bzero(nd_red,REDTAB_LEN*sizeof(RHist));
+         if ( max < 2 )
+                 nd_bpe = 2;
+         else if ( max < 4 )
+                 nd_bpe = 4;
+         else if ( max < 64 )
+                 nd_bpe = 6;
+         else if ( max < 256 )
+                 nd_bpe = 8;
+         else if ( max < 65536 )
+                 nd_bpe = 16;
+         else
+                 nd_bpe = 32;
+         nd_setup_parameters();
+         nd_free_private_storage();
+         len = 0;
+         for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) {
+                 a = dptondv(mod,(DP)BDY(f)); ndv_removecont(mod,a); nd_ps[i] = a;
+                 a = dptondv(0,(DP)BDY(f)); ndv_removecont(0,a); nd_psq[i] = a;
+                 NEWRHist(r);
+                 len = MAX(len,LEN(a));
+                 SG(r) = HTD(a); TD(r) = HTD(a); ndl_copy(HDL(a),DL(r));
+                 nd_psh[i] = r;
+         }
  }
  void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp)
-Line 1695  void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,
+Line 1971  void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,
 Line 1695  void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,
 Line 1971  void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,
                          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);
+                 NEXTNODE(fd0,fd); BDY(fd) = (pointer)b;
-                 if ( c ) {
-                         NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
-                 }
          }
          if ( fd0 ) NEXT(fd) = 0;
-         s = nd_setup(fd0);
+         nd_setup(m,fd0);
-         x = nd_gb(s);
+         x = nd_gb(m);
- #if 0
+         fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n",
-         x = nd_reduceall(x,m);
+                 nd_found,nd_notfirst,nd_create);
- #endif
+         x = nd_reduceall(m,x);
-         for ( r0 = 0; x; x = NEXT(x) ) {
+         for ( r0 = 0, t = x; t; t = NEXT(t) ) {
                  NEXTNODE(r0,r);
-                 a = ndtodp(nd_ps[(int)BDY(x)]);
+                 if ( m ) {
-                 _dtop_mod(CO,vv,a,(P *)&BDY(r));
+                         a = ndvtodp(m,BDY(t));
+                         _dtop_mod(CO,vv,a,(P *)&BDY(r));
+                 } else {
+                         a = ndvtodp(m,BDY(t));
+                         dtop(CO,vv,a,(P *)&BDY(r));
+                 }
          }
          if ( r0 ) NEXT(r) = 0;
          MKLIST(*rp,r0);
+ }
+ void nd_gr_trace(LIST f,LIST v,int m,int homo,struct order_spec *ord,LIST *rp)
+ {
+         struct order_spec ord1;
+         VL fv,vv,vc;
+         NODE fd,fd0,r,r0,t,x,s,xx,t0;
+         DP a,b,c,h;
+         NDV *w;
+         int len,i,j;
+         get_vars((Obj)f,&fv); pltovl(v,&vv);
+         nd_nvar = length(vv);
+         if ( ord->id )
+                 error("nd_gr : unsupported order");
+         initd(ord);
+         if ( homo ) {
+                 homogenize_order(ord,nd_nvar,&ord1);
+                 switch ( ord1.ord.simple ) {
+                         case 0: is_rlex = 1; break;
+                         case 1: is_rlex = 0; break;
+                         default: error("nd_gr : unsupported order");
+                 }
+                 for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
+                         ptod(CO,vv,(P)BDY(t),&c);
+                         if ( c ) {
+                                 dp_homo(c,&h); NEXTNODE(fd0,fd); BDY(fd) = (pointer)h;
+                         }
+                 }
+                 if ( fd0 ) NEXT(fd) = 0;
+                 initd(&ord1);
+                 nd_nvar++;
+                 nd_setup_trace(m,fd0);
+                 x = nd_gb_trace(m);
+                 /* dehomogenization */
+                 for ( t = x; t; t = NEXT(t) )
+                         ndv_dehomogenize((NDV)BDY(t));
+                 nd_nvar--;
+                 nd_setup_parameters();
+                 initd(ord);
+                 len = length(x);
+                 w = (NDV *)ALLOCA(len*sizeof(NDV));
+                 for ( i = 0, t = x; i < len; i++, t = NEXT(t) ) w[i] = BDY(t);
+                 for ( i = 0; i < len; i++ ) {
+                         for ( j = 0; j < i; j++ ) {
+                                 if ( w[i] && w[j] )
+                                         if ( ndl_reducible(HDL(w[i]),HDL(w[j])) ) w[i] = 0;
+                                         else if ( ndl_reducible(HDL(w[j]),HDL(w[i])) ) w[j] = 0;
+                         }
+                 }
+                 for ( i = len-1, t0 = 0; i >= 0; i-- ) {
+                         if ( w[i] ) { NEXTNODE(t0,t); BDY(t) = (pointer)w[i]; }
+                 }
+                 NEXT(t) = 0; x = t0;
+         } else {
+                 switch ( ord->ord.simple ) {
+                         case 0: is_rlex = 1; break;
+                         case 1: is_rlex = 0; break;
+                         default: error("nd_gr : unsupported order");
+                 }
+                 for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
+                         ptod(CO,vv,(P)BDY(t),&c);
+                         if ( c ) {
+                                 NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
+                         }
+                 }
+                 if ( fd0 ) NEXT(fd) = 0;
+                 /* setup over GF(m) */
+                 nd_setup_trace(m,fd0);
+                 x = nd_gb_trace(m);
+         }
          fprintf(asir_out,"found=%d,notfirst=%d,create=%d\n",
                  nd_found,nd_notfirst,nd_create);
+         x = nd_reduceall(0,x);
+         for ( r0 = 0; x; x = NEXT(x) ) {
+                 NEXTNODE(r0,r);
+                 a = ndvtodp(0,(NDV)BDY(x));
+                 dtop(CO,vv,a,(P *)&BDY(r));
+         }
+         if ( r0 ) NEXT(r) = 0;
+         MKLIST(*rp,r0);
  }
  void dltondl(int n,DL dl,unsigned int *r)
-Line 1742  DL ndltodl(int n,int td,unsigned int *ndl)
+Line 2099  DL ndltodl(int n,int td,unsigned int *ndl)
 Line 1742  DL ndltodl(int n,int td,unsigned int *ndl)
 Line 2099  DL ndltodl(int n,int td,unsigned int *ndl)
          int i;
          NEWDL(dl,n);
-         dl->td = td;
+         TD(dl) = td;
          d = dl->d;
          if ( is_rlex )
                  for ( i = 0; i < n; i++ )
-Line 1755  DL ndltodl(int n,int td,unsigned int *ndl)
+Line 2112  DL ndltodl(int n,int td,unsigned int *ndl)
 Line 1755  DL ndltodl(int n,int td,unsigned int *ndl)
 Line 2112  DL ndltodl(int n,int td,unsigned int *ndl)
          return dl;
  }
- ND dptond(DP p)
+ ND dptond(int mod,DP p)
  {
          ND d;
          NM m0,m;
-Line 1768  ND dptond(DP p)
+Line 2125  ND dptond(DP p)
 Line 1768  ND dptond(DP p)
 Line 2125  ND dptond(DP p)
          m0 = 0;
          for ( t = BDY(p); t; t = NEXT(t) ) {
                  NEXTNM(m0,m);
-                 m->c = ITOS(t->c);
+                 if ( mod )
-                 m->td = t->dl->td;
+                         CM(m) = ITOS(C(t));
-                 dltondl(n,t->dl,m->dl);
+                 else
+                         CQ(m) = (Q)C(t);
+                 TD(m) = TD(DL(t));
+                 dltondl(n,DL(t),DL(m));
          }
          NEXT(m) = 0;
          MKND(n,m0,d);
-         d->nv = n;
+         NV(d) = n;
-         d->sugar = p->sugar;
+         SG(d) = SG(p);
          return d;
  }
- DP ndtodp(ND p)
+ DP ndtodp(int mod,ND p)
  {
          DP d;
          MP m0,m;
-Line 1792  DP ndtodp(ND p)
+Line 2152  DP ndtodp(ND p)
 Line 1792  DP ndtodp(ND p)
 Line 2152  DP ndtodp(ND p)
          m0 = 0;
          for ( t = BDY(p); t; t = NEXT(t) ) {
                  NEXTMP(m0,m);
-                 m->c = STOI(t->c);
+                 if ( mod )
-                 m->dl = ndltodl(n,t->td,t->dl);
+                         C(m) = STOI(CM(t));
+                 else
+                         C(m) = (P)CQ(t);
+                 DL(m) = ndltodl(n,TD(t),DL(t));
          }
          NEXT(m) = 0;
          MKDP(n,m0,d);
-         d->sugar = p->sugar;
+         SG(d) = SG(p);
          return d;
  }
-Line 1829  void nd_print(ND p)
+Line 2192  void nd_print(ND p)
 Line 1829  void nd_print(ND p)
 Line 2192  void nd_print(ND p)
                  printf("0\n");
          else {
                  for ( m = BDY(p); m; m = NEXT(m) ) {
-                         printf("+%d*",m->c);
+                         printf("+%d*",CM(m));
-                         ndl_print(m->dl);
+                         ndl_print(DL(m));
                  }
                  printf("\n");
          }
  }
+ void nd_print_q(ND p)
+ {
+         NM m;
+         if ( !p )
+                 printf("0\n");
+         else {
+                 for ( m = BDY(p); m; m = NEXT(m) ) {
+                         printf("+");
+                         printexpr(CO,CQ(m));
+                         printf("*");
+                         ndl_print(DL(m));
+                 }
+                 printf("\n");
+         }
+ }
  void ndp_print(ND_pairs d)
  {
          ND_pairs t;
-Line 1846  void ndp_print(ND_pairs d)
+Line 2226  void ndp_print(ND_pairs d)
 Line 1846  void ndp_print(ND_pairs d)
 Line 2226  void ndp_print(ND_pairs d)
          printf("\n");
  }
- void nd_monic(ND p)
+ void nd_removecont(int mod,ND p)
  {
+         int i,n;
+         Q *w;
+         Q dvr,t;
+         NM m;
+         struct oVECT v;
+         N q,r;
+         if ( mod )
+                 nd_mul_c(mod,p,invm(HCM(p),mod));
+         else {
+                 for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
+                 w = (Q *)ALLOCA(n*sizeof(Q));
+                 v.len = n;
+                 v.body = (pointer *)w;
+                 for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ )
+                         w[i] = CQ(m);
+                 removecont_array(w,n);
+                 for ( m = BDY(p), i = 0; i < n; m = NEXT(m), i++ ) CQ(m) = w[i];
+         }
+ }
+ void nd_removecont2(ND p1,ND p2)
+ {
+         int i,n1,n2,n;
+         Q *w;
+         Q dvr,t;
+         NM m;
+         struct oVECT v;
+         N q,r;
+         if ( !p1 ) {
+                 nd_removecont(0,p2); return;
+         } else if ( !p2 ) {
+                 nd_removecont(0,p1); return;
+         }
+         n1 = nd_length(p1);
+         n2 = nd_length(p2);
+         n = n1+n2;
+         w = (Q *)ALLOCA(n*sizeof(Q));
+         v.len = n;
+         v.body = (pointer *)w;
+         for ( m = BDY(p1), i = 0; i < n1; m = NEXT(m), i++ )
+                 w[i] = CQ(m);
+         for ( m = BDY(p2); i < n; m = NEXT(m), i++ )
+                 w[i] = CQ(m);
+         removecont_array(w,n);
+         for ( m = BDY(p1), i = 0; i < n1; m = NEXT(m), i++ ) CQ(m) = w[i];
+         for ( m = BDY(p2); i < n; m = NEXT(m), i++ ) CQ(m) = w[i];
+ }
+ void ndv_removecont(int mod,NDV p)
+ {
+         int i,len;
+         Q *w;
+         Q dvr,t;
+         NMV m;
+         if ( mod )
+                 ndv_mul_c(mod,p,invm(HCM(p),mod));
+         else {
+                 len = p->len;
+                 w = (Q *)ALLOCA(len*sizeof(Q));
+                 for ( m = BDY(p), i = 0; i < len; NMV_ADV(m), i++ )
+                         w[i] = CQ(m);
+                 sortbynm(w,len);
+                 qltozl(w,len,&dvr);
+                 for ( m = BDY(p), i = 0; i < len; NMV_ADV(m), i++ ) {
+                         divq(CQ(m),dvr,&t); CQ(m) = t;
+                 }
+         }
+ }
+ void ndv_dehomogenize(NDV p)
+ {
+         int i,len,newnvar,newwpd,newadv;
+         Q *w;
+         Q dvr,t;
+         NMV m,r;
+         unsigned int *d;
+ #define NEWADV(m) (m = (NMV)(((char *)m)+newadv))
+         len = p->len;
+         newnvar = nd_nvar-1;
+         newwpd = newnvar/nd_epw+(newnvar%nd_epw?1:0);
+         for ( m = BDY(p), i = 0; i < len; NMV_ADV(m), i++ )
+                 TD(m) -= ndl_dehomogenize(DL(m));
+         if ( newwpd != nd_wpd ) {
+                 d = (unsigned int *)ALLOCA(newwpd*sizeof(unsigned int));
+                 newadv = sizeof(struct oNMV)+(newwpd-1)*sizeof(unsigned int);
+                 for ( m = r = BDY(p), i = 0; i < len; NMV_ADV(m), NEWADV(r), i++ ) {
+                         CQ(r) = CQ(m); TD(r) = TD(m); ndl_copy(DL(m),d); ndl_copy(d,DL(r));
+                 }
+         }
+         NV(p)--;
+ }
+ void removecont_array(Q *c,int n)
+ {
+         struct oVECT v;
+         Q d0,d1,a,u,u1,gcd;
+         int i;
+         N qn,rn,gn;
+         Q *q,*r;
+         q = (Q *)ALLOCA(n*sizeof(Q));
+         r = (Q *)ALLOCA(n*sizeof(Q));
+         v.id = O_VECT; v.len = n; v.body = (pointer *)c;
+         igcdv_estimate(&v,&d0);
+         for ( i = 0; i < n; i++ ) {
+                 divn(NM(c[i]),NM(d0),&qn,&rn);
+                 NTOQ(qn,SGN(c[i])*SGN(d0),q[i]);
+                 NTOQ(rn,SGN(c[i]),r[i]);
+         }
+         for ( i = 0; i < n; i++ )
+                 if ( r[i] )
+                         break;
+         if ( i < n ) {
+                 v.id = O_VECT; v.len = n; v.body = (pointer *)r;
+                 igcdv(&v,&d1);
+                 gcdn(NM(d0),NM(d1),&gn); NTOQ(gn,1,gcd);
+                 divsn(NM(d0),gn,&qn); NTOQ(qn,1,a);
+                 for ( i = 0; i < n; i++ ) {
+                         mulq(a,q[i],&u);
+                         if ( r[i] ) {
+                                 divsn(NM(r[i]),gn,&qn); NTOQ(qn,SGN(r[i]),u1);
+                                 addq(u,u1,&q[i]);
+                         } else
+                                 q[i] = u;
+                 }
+         }
+         for ( i = 0; i < n; i++ )
+                 c[i] = q[i];
+ }
+ void nd_mul_c(int mod,ND p,int mul)
+ {
+         NM m;
+         int c,c1;
          if ( !p )
                  return;
-         else
+         for ( m = BDY(p); m; m = NEXT(m) ) {
-                 nd_mul_c(p,invm(HC(p),nd_mod));
+                 c1 = CM(m);
+                 DMAR(c1,mul,0,mod,c);
+                 CM(m) = c;
+         }
  }
- void nd_mul_c(ND p,int mul)
+ void nd_mul_c_q(ND p,Q mul)
  {
          NM m;
-         int c,c1;
+         Q c;
          if ( !p )
                  return;
          for ( m = BDY(p); m; m = NEXT(m) ) {
-                 c1 = C(m);
+                 mulq(CQ(m),mul,&c); CQ(m) = c;
-                 DMAR(c1,mul,0,nd_mod,c);
-                 C(m) = c;
          }
  }
-Line 1883  void nd_free(ND p)
+Line 2403  void nd_free(ND p)
 Line 1883  void nd_free(ND p)
 Line 2403  void nd_free(ND p)
          FREEND(p);
  }
+ void ndv_free(NDV p)
+ {
+         GC_free(BDY(p));
+ }
  void nd_append_red(unsigned int *d,int td,int i)
  {
-         NM m,m0;
+         RHist m,m0;
          int h;
-         NEWNM(m);
+         NEWRHist(m);
          h = ndl_hash_value(td,d);
-         m->c = i;
+         m->index = i;
-         m->td = td;
+         TD(m) = td;
-         bcopy(d,m->dl,nd_wpd*sizeof(unsigned int));
+         ndl_copy(d,DL(m));
          NEXT(m) = nd_red[h];
          nd_red[h] = m;
  }
-Line 1901  unsigned int *dp_compute_bound(DP p)
+Line 2426  unsigned int *dp_compute_bound(DP p)
 Line 1901  unsigned int *dp_compute_bound(DP p)
 Line 2426  unsigned int *dp_compute_bound(DP p)
  {
          unsigned int *d,*d1,*d2,*t;
          MP m;
-         int i;
+         int i,l;
          if ( !p )
                  return 0;
          d1 = (unsigned int *)ALLOCA(nd_nvar*sizeof(unsigned int));
          d2 = (unsigned int *)ALLOCA(nd_nvar*sizeof(unsigned int));
          m = BDY(p);
-         bcopy(m->dl->d,d1,nd_nvar*sizeof(unsigned int));
+         bcopy(DL(m)->d,d1,nd_nvar*sizeof(unsigned int));
          for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) {
-                 d = m->dl->d;
+                 d = DL(m)->d;
                  for ( i = 0; i < nd_nvar; i++ )
                          d2[i] = d[i] > d1[i] ? d[i] : d1[i];
                  t = d1; d1 = d2; d2 = t;
          }
-         t = (unsigned int *)MALLOC_ATOMIC(nd_nvar*sizeof(unsigned int));
+         l = (nd_nvar+31);
+         t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));
+         bzero(t,l*sizeof(unsigned int));
          bcopy(d1,t,nd_nvar*sizeof(unsigned int));
          return t;
  }
-Line 1923  unsigned int *dp_compute_bound(DP p)
+Line 2450  unsigned int *dp_compute_bound(DP p)
 Line 1923  unsigned int *dp_compute_bound(DP p)
 Line 2450  unsigned int *dp_compute_bound(DP p)
  unsigned int *nd_compute_bound(ND p)
  {
          unsigned int *d1,*d2,*t;
-         int i;
+         int i,l;
          NM m;
          if ( !p )
-Line 1932  unsigned int *nd_compute_bound(ND p)
+Line 2459  unsigned int *nd_compute_bound(ND p)
 Line 1932  unsigned int *nd_compute_bound(ND p)
 Line 2459  unsigned int *nd_compute_bound(ND p)
          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);
+                 ndl_lcm(DL(m),d1,d2);
                  t = d1; d1 = d2; d2 = t;
          }
-         t = (unsigned int *)MALLOC_ATOMIC(nd_nvar*sizeof(unsigned int));
+         l = nd_nvar+31;
-         bzero(t,nd_nvar*sizeof(unsigned int));
+         t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));
+         bzero(t,l*sizeof(unsigned int));
          for ( i = 0; i < nd_nvar; i++ )
                  t[i] = (d1[i/nd_epw]>>((nd_epw-(i%nd_epw)-1)*nd_bpe))&nd_mask0;
          return t;
-Line 1958  void nd_setup_parameters() {
+Line 2486  void nd_setup_parameters() {
 Line 1958  void nd_setup_parameters() {
 Line 2486  void nd_setup_parameters() {
                  nd_mask[nd_epw-i-1] = (nd_mask0<<(i*nd_bpe));
                  nd_mask1 |= (1<<(nd_bpe-1))<<(i*nd_bpe);
          }
- #if USE_NDV
+         nm_adv = sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int);
          nmv_adv = sizeof(struct oNMV)+(nd_wpd-1)*sizeof(unsigned int);
- #endif
  }
- ND_pairs nd_reconstruct(ND_pairs d)
+ /* mod < 0 => realloc nd_ps and pd_psq */
+ ND_pairs nd_reconstruct(int mod,int trace,ND_pairs d)
  {
-         int i,obpe;
+         int i,obpe,oadv;
-         NM prev_nm_free_list,mr0,mr,m;
+         NM prev_nm_free_list;
+         RHist mr0,mr;
+         RHist r;
          ND_pairs s0,s,t,prev_ndp_free_list;
          obpe = nd_bpe;
+         oadv = nmv_adv;
          if ( obpe < 4 )
                  nd_bpe = 4;
          else if ( obpe < 6 )
-Line 1988  ND_pairs nd_reconstruct(ND_pairs d)
+Line 2520  ND_pairs nd_reconstruct(ND_pairs d)
 Line 1988  ND_pairs nd_reconstruct(ND_pairs d)
 Line 2520  ND_pairs nd_reconstruct(ND_pairs d)
          prev_ndp_free_list = _ndp_free_list;
          _nm_free_list = 0;
          _ndp_free_list = 0;
-         for ( i = 0; i < nd_psn; i++ ) {
+         if ( mod != 0 )
-                 nd_ps[i] = nd_dup(nd_ps[i],obpe);
+                 for ( i = nd_psn-1; i >= 0; i-- )
- #if USE_NDV
+                         ndv_realloc(nd_ps[i],obpe,oadv);
-                 nd_psv[i] = ndtondv(nd_ps[i]);
+         if ( !mod || trace )
- #endif
+                 for ( i = nd_psn-1; i >= 0; i-- )
-         }
+                         ndv_realloc(nd_psq[i],obpe,oadv);
          s0 = 0;
          for ( t = d; t; t = NEXT(t) ) {
                  NEXTND_pairs(s0,s);
                  s->i1 = t->i1;
                  s->i2 = t->i2;
-                 s->td = t->td;
+                 TD(s) = TD(t);
-                 s->sugar = t->sugar;
+                 SG(s) = SG(t);
                  ndl_dup(obpe,t->lcm,s->lcm);
          }
          for ( i = 0; i < REDTAB_LEN; i++ ) {
-                 for ( mr0 = 0, m = nd_red[i]; m; m = NEXT(m) ) {
+                 for ( mr0 = 0, r = nd_red[i]; r; r = NEXT(r) ) {
-                         NEXTNM(mr0,mr);
+                         NEXTRHist(mr0,mr);
-                         mr->c = m->c;
+                         mr->index = r->index;
-                         mr->td = m->td;
+                         SG(mr) = SG(r);
-                         ndl_dup(obpe,m->dl,mr->dl);
+                         TD(mr) = TD(r);
+                         ndl_dup(obpe,DL(r),DL(mr));
                  }
                  if ( mr0 ) NEXT(mr) = 0;
                  nd_red[i] = mr0;
          }
+         for ( i = 0; i < nd_psn; i++ ) {
+                 NEWRHist(r); SG(r) = SG(nd_psh[i]);
+                 TD(r) = TD(nd_psh[i]); ndl_dup(obpe,DL(nd_psh[i]),DL(r));
+                 nd_psh[i] = r;
+         }
          if ( s0 ) NEXT(s) = 0;
          prev_nm_free_list = 0;
          prev_ndp_free_list = 0;
- #if USE_NDV
-         BDY(ndv_red) = (NMV)REALLOC(BDY(ndv_red),nmv_len*nmv_adv);
- #endif
          GC_gcollect();
          return s0;
  }
+ void nd_reconstruct_direct(int mod,NDV *ps,int len)
+ {
+         int i,obpe,oadv;
+         NM prev_nm_free_list;
+         RHist mr0,mr;
+         RHist r;
+         ND_pairs s0,s,t,prev_ndp_free_list;
+         obpe = nd_bpe;
+         oadv = nmv_adv;
+         if ( obpe < 4 )
+                 nd_bpe = 4;
+         else if ( obpe < 6 )
+                 nd_bpe = 6;
+         else if ( obpe < 8 )
+                 nd_bpe = 8;
+         else if ( obpe < 16 )
+                 nd_bpe = 16;
+         else if ( obpe < 32 )
+                 nd_bpe = 32;
+         else
+                 error("nd_reconstruct_direct : exponent too large");
+         nd_setup_parameters();
+         prev_nm_free_list = _nm_free_list;
+         prev_ndp_free_list = _ndp_free_list;
+         _nm_free_list = 0;
+         _ndp_free_list = 0;
+         if ( mod != 0 )
+                 for ( i = len-1; i >= 0; i-- )
+                         ndv_realloc(ps[i],obpe,oadv);
+         prev_nm_free_list = 0;
+         prev_ndp_free_list = 0;
+         GC_gcollect();
+ }
  void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
  {
          int n,i,ei,oepw,cepw,cbpe;
-Line 2031  void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
+Line 2602  void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
 Line 2031  void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
 Line 2602  void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
          oepw = (sizeof(unsigned int)*8)/obpe;
          cepw = nd_epw;
          cbpe = nd_bpe;
+         for ( i = 0; i < nd_wpd; i++ )
+                 r[i] = 0;
          if ( is_rlex )
                  for ( i = 0; i < n; i++ ) {
                          ei = (d[(n-1-i)/oepw]>>((oepw-((n-1-i)%oepw)-1)*obpe))
-Line 2045  void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
+Line 2618  void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
 Line 2045  void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
 Line 2618  void ndl_dup(int obpe,unsigned int *d,unsigned int *r)
                  }
  }
- ND nd_dup(ND p,int obpe)
+ void nd_realloc(ND p,int obpe)
  {
          NM m,mr,mr0;
-         int c,n;
-         ND r;
-         if ( !p )
+         if ( p ) {
-                 return 0;
+                 m = BDY(p);
-         else {
-                 n = NV(p); m = BDY(p);
                  for ( mr0 = 0; m; m = NEXT(m) ) {
                          NEXTNM(mr0,mr);
-                         C(mr) = C(m);
+                         CM(mr) = CM(m);
-                         mr->td = m->td;
+                         TD(mr) = TD(m);
-                         ndl_dup(obpe,m->dl,mr->dl);
+                         ndl_dup(obpe,DL(m),DL(mr));
                  }
                  NEXT(mr) = 0;
-                 MKND(NV(p),mr0,r);
+                 BDY(p) = mr0;
-                 r->sugar = p->sugar;
-                 return r;
          }
  }
-Line 2080  ND nd_copy(ND p)
+Line 2647  ND nd_copy(ND p)
 Line 2080  ND nd_copy(ND p)
 Line 2647  ND nd_copy(ND p)
                  s = sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int);
                  for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                          NEXTNM(mr0,mr);
-                         C(mr) = C(m);
+                         CM(mr) = CM(m);
-                         mr->td = m->td;
+                         TD(mr) = TD(m);
-                         ndl_copy(m->dl,mr->dl);
+                         ndl_copy(DL(m),DL(mr));
                  }
                  NEXT(mr) = 0;
                  MKND(NV(p),mr0,r);
-                 r->sugar = p->sugar;
+                 SG(r) = SG(p);
                  return r;
          }
  }
- #if USE_NDV
+ int nd_sp(int mod,ND_pairs p,ND *rp)
- void ndv_mul_nm(NDV p,NM m0,NDV r)
  {
-         NMV m,mr,mr0;
+         NM m;
-         unsigned int *d,*dt,*dm;
+         NDV p1,p2;
-         int c,n,td,i,c1,c2,len;
+         ND t1,t2;
+         unsigned int *lcm;
+         int td;
+         if ( mod ) {
+                 p1 = nd_ps[p->i1]; p2 = nd_ps[p->i2];
+         } else {
+                 p1 = nd_psq[p->i1]; p2 = nd_psq[p->i2];
+         }
+         lcm = p->lcm;
+         td = TD(p);
+         NEWNM(m);
+         CQ(m) = HCQ(p2);
+         TD(m) = td-HTD(p1); ndl_sub(lcm,HDL(p1),DL(m));
+         if ( ndl_check_bound2(p->i1,DL(m)) )
+                 return 0;
+         t1 = ndv_mul_nm(mod,p1,m);
+         if ( mod )
+                 CM(m) = mod-HCM(p1);
+         else
+                 chsgnq(HCQ(p1),&CQ(m));
+         TD(m) = td-HTD(p2); ndl_sub(lcm,HDL(p2),DL(m));
+         if ( ndl_check_bound2(p->i2,DL(m)) ) {
+                 nd_free(t1);
+                 return 0;
+         }
+         t2 = ndv_mul_nm(mod,p2,m);
+         *rp = nd_add(mod,t1,t2);
+         FREENM(m);
+         return 1;
+ }
+ void ndv_mul_c(int mod,NDV p,int mul)
+ {
+         NMV m;
+         int c,c1,len,i;
          if ( !p )
-                 /* XXX */
+                 return;
-                 r->len = 0;
+         len = LEN(p);
-         else {
+         for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
-                 n = NV(p); m = BDY(p); len = p->len;
+                 c1 = CM(m);
-                 d = m0->dl; td = m0->td; c = C(m0);
+                 DMAR(c1,mul,0,mod,c);
-                 mr = BDY(r);
+                 CM(m) = c;
-                 for ( i = 0; i < len; i++, NMV_ADV(m), NMV_ADV(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);
-                 }
-                 NV(r) = NV(p);
-                 r->len = p->len;
-                 r->sugar = p->sugar + td;
          }
  }
- ND ndv_add(ND p1,NDV p2)
+ void ndv_mul_c_q(NDV p,Q mul)
  {
-         int c,c1,c2,t,td,td2,mul,len,i;
+         NMV m;
-         NM prev,head,cur,new;
+         Q c;
-         unsigned int *d;
+         int len,i;
-         NMV m2;
-         if ( !p1 )
+         if ( !p )
+                 return;
+         len = LEN(p);
+         for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
+                 mulq(CQ(m),mul,&c); CQ(m) = c;
+         }
+ }
+ ND ndv_mul_nm(int mod,NDV p,NM m0)
+ {
+         NM mr,mr0;
+         NMV m;
+         unsigned int *d,*dt,*dm;
+         int c,n,td,i,c1,c2,len;
+         Q q;
+         ND r;
+         if ( !p )
                  return 0;
          else {
-                 prev = 0; head = cur = BDY(p1);
+                 n = NV(p); m = BDY(p);
-                 NEWNM(new);
+                 d = DL(m0); td = TD(m0);
-                 len = p2->len;
+                 len = LEN(p);
-                 for ( m2 = BDY(p2), i = 0; i < len; ) {
+                 mr0 = 0;
-                         td2 = new->td = m2->td;
+                 if ( mod ) {
-                         if ( !cur ) {
+                         c = CM(m0);
-                                 C(new) = C(m2);
+                         for ( i = 0; i < len; i++, NMV_ADV(m) ) {
-                                 ndl_copy(m2->dl,new->dl);
+                                 NEXTNM(mr0,mr);
-                                 if ( !prev ) {
+                                 c1 = CM(m);
-                                         prev = new;
+                                 DMAR(c1,c,0,mod,c2);
-                                         NEXT(prev) = 0;
+                                 CM(mr) = c2;
-                                         head = prev;
+                                 TD(mr) = TD(m)+td;
-                                 } else {
+                                 ndl_add(DL(m),d,DL(mr));
-                                         NEXT(prev) = new;
-                                         NEXT(new) = 0;
-                                         prev = new;
-                                 }
-                                 NMV_ADV(m2); i++;
-                                 NEWNM(new);
-                                 continue;
                          }
-                         if ( cur->td > td2 )
+                 } else {
-                                 c = 1;
+                         q = CQ(m0);
-                         else if ( cur->td < td2 )
+                         for ( i = 0; i < len; i++, NMV_ADV(m) ) {
-                                 c = -1;
+                                 NEXTNM(mr0,mr);
-                         else
+                                 mulq(CQ(m),q,&CQ(mr));
-                                 c = ndl_compare(cur->dl,m2->dl);
+                                 TD(mr) = TD(m)+td;
-                         switch ( c ) {
+                                 ndl_add(DL(m),d,DL(mr));
-                                 case 0:
-                                         t = C(m2)+C(cur)-nd_mod;
-                                         if ( t < 0 )
-                                                 t += nd_mod;
-                                         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);
-                                         }
-                                         NMV_ADV(m2); i++;
-                                         break;
-                                 case 1:
-                                         prev = cur;
-                                         cur = NEXT(cur);
-                                         break;
-                                 case -1:
-                                         ndl_copy(m2->dl,new->dl);
-                                         if ( !prev ) {
-                                                 /* cur = head */
-                                                 prev = new;
-                                                 C(prev) = C(m2);
-                                                 NEXT(prev) = head;
-                                                 head = prev;
-                                         } else {
-                                                 C(new) = C(m2);
-                                                 NEXT(prev) = new;
-                                                 NEXT(new) = cur;
-                                                 prev = new;
-                                         }
-                                         NEWNM(new);
-                                         NMV_ADV(m2); i++;
-                                         break;
                          }
                  }
-                 FREENM(new);
+                 NEXT(mr) = 0;
-                 if ( head ) {
+                 MKND(NV(p),mr0,r);
-                         BDY(p1) = head;
+                 SG(r) = SG(p) + td;
-                         p1->sugar = MAX(p1->sugar,p2->sugar+td);
+                 return r;
-                         return p1;
+         }
-                 } else {
+ }
-                         FREEND(p1);
-                         return 0;
+ void ndv_realloc(NDV p,int obpe,int oadv)
+ {
+         NMV m,mr,mr0,t;
+         int len,i,k;
+ #define NMV_OPREV(m) (m = (NMV)(((char *)m)-oadv))
+ #define NMV_PREV(m) (m = (NMV)(((char *)m)-nmv_adv))
+         if ( p ) {
+                 m = BDY(p); len = LEN(p);
+                 if ( nmv_adv > oadv )
+                         mr0 = (NMV)REALLOC(BDY(p),len*nmv_adv);
+                 else
+                         mr0 = BDY(p);
+                 m = (NMV)((char *)mr0+(len-1)*oadv);
+                 mr = (NMV)((char *)mr0+(len-1)*nmv_adv);
+                 t = (NMV)ALLOCA(nmv_adv);
+                 for ( i = 0; i < len; i++, NMV_OPREV(m), NMV_PREV(mr) ) {
+                         CQ(t) = CQ(m);
+                         TD(t) = TD(m);
+                         for ( k = 0; k < nd_wpd; k++ ) DL(t)[k] = 0;
+                         ndl_dup(obpe,DL(m),DL(t));
+                         CQ(mr) = CQ(t);
+                         TD(mr) = TD(t);
+                         ndl_copy(DL(t),DL(mr));
                  }
+                 BDY(p) = mr0;
          }
  }
- NDV ndtondv(ND p)
+ NDV ndtondv(int mod,ND p)
  {
          NDV d;
          NMV m,m0;
-Line 2219  NDV ndtondv(ND p)
+Line 2805  NDV ndtondv(ND p)
 Line 2219  NDV ndtondv(ND p)
 Line 2805  NDV ndtondv(ND p)
          if ( !p )
                  return 0;
          len = nd_length(p);
-         m0 = m = (NMV)MALLOC_ATOMIC(len*nmv_adv);
+         if ( mod )
+                 m0 = m = (NMV)MALLOC_ATOMIC(len*nmv_adv);
+         else
+                 m0 = m = (NMV)MALLOC(len*nmv_adv);
          for ( t = BDY(p), i = 0; t; t = NEXT(t), i++, NMV_ADV(m) ) {
-                 m->td = t->td;
+                 TD(m) = TD(t);
-                 bcopy(t->dl,m->dl,nd_wpd*sizeof(unsigned int));
+                 ndl_copy(DL(t),DL(m));
-                 m->c = t->c;
+                 CQ(m) = CQ(t);
          }
          MKNDV(NV(p),m0,len,d);
-         d->sugar = p->sugar;
+         SG(d) = SG(p);
          return d;
  }
+ ND ndvtond(int mod,NDV p)
+ {
+         ND d;
+         NM m,m0;
+         NMV t;
+         int i,len;
+         if ( !p )
+                 return 0;
+         m0 = 0;
+         len = p->len;
+         for ( t = BDY(p), i = 0; i < len; NMV_ADV(t), i++ ) {
+                 NEXTNM(m0,m);
+                 TD(m) = TD(t);
+                 ndl_copy(DL(t),DL(m));
+                 CQ(m) = CQ(t);
+         }
+         NEXT(m) = 0;
+         MKND(NV(p),m0,d);
+         SG(d) = SG(p);
+         return d;
+ }
+ NDV dptondv(int mod,DP p)
+ {
+         NDV d;
+         NMV m0,m;
+         MP t;
+         DP q;
+         int l,i,n;
+         if ( !p )
+                 return 0;
+         for ( t = BDY(p), l = 0; t; t = NEXT(t), l++ );
+         if ( mod ) {
+                 _dp_mod(p,mod,0,&q); p = q;
+                 m0 = m = (NMV)MALLOC_ATOMIC(l*nmv_adv);
+         } else
+                 m0 = m = (NMV)MALLOC(l*nmv_adv);
+         n = NV(p);
+         for ( t = BDY(p), i = 0; i < l; i++, t = NEXT(t), NMV_ADV(m) ) {
+                 if ( mod )
+                         CM(m) = ITOS(C(t));
+                 else
+                         CQ(m) = (Q)C(t);
+                 TD(m) = TD(DL(t));
+                 dltondl(n,DL(t),DL(m));
+         }
+         MKNDV(n,m0,l,d);
+         SG(d) = SG(p);
+         return d;
+ }
+ DP ndvtodp(int mod,NDV p)
+ {
+         DP d;
+         MP m0,m;
+         NMV t;
+         int len,i,n;
+         if ( !p )
+                 return 0;
+         m0 = 0;
+         len = LEN(p);
+         n = NV(p);
+         for ( t = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) {
+                 NEXTMP(m0,m);
+                 if ( mod )
+                         C(m) = STOI(CM(t));
+                 else
+                         C(m) = (P)CQ(t);
+                 DL(m) = ndltodl(n,TD(t),DL(t));
+         }
+         NEXT(m) = 0;
+         MKDP(NV(p),m0,d);
+         SG(d) = SG(p);
+         return d;
+ }
  void ndv_print(NDV p)
  {
          NMV m;
-Line 2238  void ndv_print(NDV p)
+Line 2906  void ndv_print(NDV p)
 Line 2238  void ndv_print(NDV p)
 Line 2906  void ndv_print(NDV p)
          if ( !p )
                  printf("0\n");
          else {
-                 len = p->len;
+                 len = LEN(p);
                  for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
-                         printf("+%d*",m->c);
+                         printf("+%d*",CM(m));
-                         ndl_print(m->dl);
+                         ndl_print(DL(m));
                  }
                  printf("\n");
          }
  }
- #endif
+ void ndv_print_q(NDV p)
+ {
+         NMV m;
+         int i,len;
+         if ( !p )
+                 printf("0\n");
+         else {
+                 len = LEN(p);
+                 for ( m = BDY(p), i = 0; i < len; i++, NMV_ADV(m) ) {
+                         printf("+");
+                         printexpr(CO,CQ(m));
+                         printf("*");
+                         ndl_print(DL(m));
+                 }
+                 printf("\n");
+         }
+ }
+ int nd_equal(ND a,ND b)
+ {
+         NM ma,mb;
+         if ( (!a && b)||(a&& !b)||(NV(a) != NV(b))||(SG(a) != SG(b)) ) return 0;
+         for ( ma = BDY(a), mb = BDY(b); ma && mb; ma = NEXT(ma), mb = NEXT(mb) ) {
+                 if ( (TD(ma) != TD(mb)) || !ndl_equal(DL(ma),DL(mb)) ) return 0;
+         }
+         if ( !ma && !mb ) return 1;
+         else return 0;
+ }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>