[BACK]Return to nd.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/nd.c between version 1.52 and 1.58

version 1.52, 2003/08/27 02:21:16 version 1.58, 2003/09/05 07:00:37
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.51 2003/08/27 01:53:29 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.57 2003/09/05 05:02:53 noro Exp $ */
   
 #include "ca.h"  #include "ca.h"
 #include "inline.h"  #include "inline.h"
Line 109  static ND _nd_free_list;
Line 109  static ND _nd_free_list;
 static ND_pairs _ndp_free_list;  static ND_pairs _ndp_free_list;
   
 static NDV *nd_ps;  static NDV *nd_ps;
 static NDV *nd_psq;  static NDV *nd_ps_trace;
 static RHist *nd_psh;  static RHist *nd_psh;
 static int nd_psn,nd_pslen;  static int nd_psn,nd_pslen;
   
Line 120  static int nm_adv;
Line 120  static int nm_adv;
 static int nmv_adv;  static int nmv_adv;
 static int nd_dcomp;  static int nd_dcomp;
   
 extern int Top,Reverse,dp_nelim;  extern int Top,Reverse,dp_nelim,do_weyl;
   extern int *current_weyl_weight_vector;
   
 /* fundamental macros */  /* fundamental macros */
 #define TD(d) (d[0])  #define TD(d) (d[0])
Line 185  if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT
Line 186  if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT
   
 /* macro for increasing pointer to NMV */  /* macro for increasing pointer to NMV */
 #define NMV_ADV(m) (m = (NMV)(((char *)m)+nmv_adv))  #define NMV_ADV(m) (m = (NMV)(((char *)m)+nmv_adv))
   #define NMV_PREV(m) (m = (NMV)(((char *)m)-nmv_adv))
   
 /* external functions */  /* external functions */
 void GC_gcollect();  void GC_gcollect();
Line 225  NODE nd_gb_trace(int m);
Line 227  NODE nd_gb_trace(int m);
 /* ndl functions */  /* ndl functions */
 int ndl_weight(unsigned int *d);  int ndl_weight(unsigned int *d);
 int ndl_weight_mask(unsigned int *d,int i);  int ndl_weight_mask(unsigned int *d,int i);
   void ndl_set_blockweight(unsigned int *d);
 void ndl_dehomogenize(unsigned int *p);  void ndl_dehomogenize(unsigned int *p);
 void ndl_reconstruct(int obpe,EPOS oepos,unsigned int *d,unsigned int *r);  void ndl_reconstruct(int obpe,EPOS oepos,unsigned int *d,unsigned int *r);
 INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2);  INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2);
Line 233  INLINE int ndl_block_compare(unsigned int *d1,unsigned
Line 236  INLINE int ndl_block_compare(unsigned int *d1,unsigned
 INLINE int ndl_equal(unsigned int *d1,unsigned int *d2);  INLINE int ndl_equal(unsigned int *d1,unsigned int *d2);
 INLINE void ndl_copy(unsigned int *d1,unsigned int *d2);  INLINE void ndl_copy(unsigned int *d1,unsigned int *d2);
 INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d);  INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d);
   INLINE void ndl_addto(unsigned int *d1,unsigned int *d2);
 INLINE void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d);  INLINE void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d);
 INLINE int ndl_hash_value(unsigned int *d);  INLINE int ndl_hash_value(unsigned int *d);
   
 /* normal forms */  /* normal forms */
 INLINE int nd_find_reducer(ND g);  INLINE int nd_find_reducer(ND g);
 INLINE int nd_find_reducer_direct(ND g,NDV *ps,int len);  INLINE int nd_find_reducer_direct(ND g,NDV *ps,int len);
 int nd_sp(int mod,ND_pairs p,ND *nf);  int nd_sp(int mod,int trace,ND_pairs p,ND *nf);
 int nd_nf(int mod,ND g,int full,ND *nf);  int nd_nf(int mod,ND g,NDV *ps,int full,ND *nf);
 int nd_nf_pbucket(int mod,ND g,int full,ND *nf);  int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *nf);
 int nd_nf_direct(int mod,ND g,BaseSet base,int full,ND *rp);  int nd_nf_direct(int mod,ND g,BaseSet base,int full,ND *rp);
   
 /* finalizers */  /* finalizers */
Line 271  void nd_reconstruct_direct(int mod,NDV *ps,int len);
Line 275  void nd_reconstruct_direct(int mod,NDV *ps,int len);
 void nd_setup(int mod,int trace,NODE f);  void nd_setup(int mod,int trace,NODE f);
 void nd_setup_parameters();  void nd_setup_parameters();
 BlockMask nd_create_blockmask(struct order_spec *ord);  BlockMask nd_create_blockmask(struct order_spec *ord);
   EPOS nd_create_epos(struct order_spec *ord);
 int nd_get_exporigin(struct order_spec *ord);  int nd_get_exporigin(struct order_spec *ord);
   
 /* ND functions */  /* ND functions */
Line 281  ND nd_remove_head(ND p);
Line 286  ND nd_remove_head(ND p);
 int nd_length(ND p);  int nd_length(ND p);
 void nd_append_red(unsigned int *d,int i);  void nd_append_red(unsigned int *d,int i);
 unsigned int *ndv_compute_bound(NDV p);  unsigned int *ndv_compute_bound(NDV p);
 unsigned int *dp_compute_bound(DP p);  
 ND nd_copy(ND p);  ND nd_copy(ND p);
 ND nd_add(int mod,ND p1,ND p2);  ND nd_add(int mod,ND p1,ND p2);
 ND nd_add_q(ND p1,ND p2);  ND nd_add_q(ND p1,ND p2);
 INLINE int nd_length(ND p);  INLINE int nd_length(ND p);
   
 /* NDV functions */  /* NDV functions */
   ND weyl_ndv_mul_nm(int mod,NM m0,NDV p);
   void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *tab,int tlen);
 void ndv_mul_c(int mod,NDV p,int mul);  void ndv_mul_c(int mod,NDV p,int mul);
 void ndv_mul_c_q(NDV p,Q mul);  void ndv_mul_c_q(NDV p,Q mul);
 void ndv_realloc(NDV p,int obpe,int oadv,EPOS oepos);  void ndv_realloc(NDV p,int obpe,int oadv,EPOS oepos);
 ND ndv_mul_nm(int mod,NDV p,NM m0);  ND ndv_mul_nm(int mod,NM m0,NDV p);
 void ndv_dehomogenize(NDV p,struct order_spec *spec);  void ndv_dehomogenize(NDV p,struct order_spec *spec);
 void ndv_removecont(int mod,NDV p);  void ndv_removecont(int mod,NDV p);
 void ndv_print(NDV p);  void ndv_print(NDV p);
Line 536  void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned in
Line 542  void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned in
         }          }
 }  }
   
   void ndl_set_blockweight(unsigned int *d) {
           int l,j;
   
           if ( nd_blockmask ) {
                   l = nd_blockmask->n;
                   for ( j = 0; j < l; j++ )
                           d[j+1] = ndl_weight_mask(d,j);
           }
   }
   
 int ndl_weight(unsigned int *d)  int ndl_weight(unsigned int *d)
 {  {
         unsigned int t,u;          unsigned int t,u;
Line 606  int ndl_block_compare(unsigned int *d1,unsigned int *d
Line 622  int ndl_block_compare(unsigned int *d1,unsigned int *d
         return 0;          return 0;
 }  }
   
   /* TDH -> WW -> TD-> RL */
   
   int ndl_ww_lex_compare(unsigned int *d1,unsigned int *d2)
   {
           int i,m,e1,e2;
   
           if ( TD(d1) > TD(d2) ) return 1;
           else if ( TD(d1) < TD(d2) ) return -1;
           m = nd_nvar>>1;
           for ( i = 0, e1 = e2 = 0; i < m; i++ ) {
                   e1 += current_weyl_weight_vector[i]*(GET_EXP(d1,m+i)-GET_EXP(d1,i));
                   e2 += current_weyl_weight_vector[i]*(GET_EXP(d2,m+i)-GET_EXP(d2,i));
           }
           if ( e1 > e2 ) return 1;
           else if ( e1 < e2 ) return -1;
           return ndl_lex_compare(d1,d2);
   }
   
 INLINE int ndl_equal(unsigned int *d1,unsigned int *d2)  INLINE int ndl_equal(unsigned int *d1,unsigned int *d2)
 {  {
         int i;          int i;
Line 661  INLINE void ndl_add(unsigned int *d1,unsigned int *d2,
Line 695  INLINE void ndl_add(unsigned int *d1,unsigned int *d2,
 #endif  #endif
 }  }
   
   /* d1 += d2 */
   INLINE void ndl_addto(unsigned int *d1,unsigned int *d2)
   {
           int i;
   
   #if 1
           switch ( nd_wpd ) {
                   case 2:
                           TD(d1) += TD(d2);
                           d1[1] += d2[1];
                           break;
                   case 3:
                           TD(d1) += TD(d2);
                           d1[1] += d2[1];
                           d1[2] += d2[2];
                           break;
                   default:
                           for ( i = 0; i < nd_wpd; i++ ) d1[i] += d2[i];
                           break;
           }
   #else
           for ( i = 0; i < nd_wpd; i++ ) d1[i] += d2[i];
   #endif
   }
   
 INLINE void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)  INLINE void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)
 {  {
         int i;          int i;
Line 1038  ND nd_add_q(ND p1,ND p2)
Line 1097  ND nd_add_q(ND p1,ND p2)
 }  }
   
 /* ret=1 : success, ret=0 : overflow */  /* ret=1 : success, ret=0 : overflow */
 int nd_nf(int mod,ND g,int full,ND *rp)  int nd_nf(int mod,ND g,NDV *ps,int full,ND *rp)
 {  {
         ND d;          ND d;
         NM m,mrd,tail;          NM m,mrd,tail;
Line 1068  int nd_nf(int mod,ND g,int full,ND *rp)
Line 1127  int nd_nf(int mod,ND g,int full,ND *rp)
                                 nd_free(g); nd_free(d);                                  nd_free(g); nd_free(d);
                                 return 0;                                  return 0;
                         }                          }
                           p = ps[index];
                         if ( mod ) {                          if ( mod ) {
                                 p = nd_ps[index];  
                                 c1 = invm(HCM(p),mod); c2 = mod-HCM(g);                                  c1 = invm(HCM(p),mod); c2 = mod-HCM(g);
                                 DMAR(c1,c2,0,mod,c); CM(mul) = c;                                  DMAR(c1,c2,0,mod,c); CM(mul) = c;
                         } else {                          } else {
                                 p = nd_psq[index];  
                                 igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);                                  igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);
                                 chsgnq(cg,&CQ(mul));                                  chsgnq(cg,&CQ(mul));
                                 nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);                                  nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);
                         }                          }
                         g = nd_add(mod,g,ndv_mul_nm(mod,p,mul));                          g = nd_add(mod,g,ndv_mul_nm(mod,mul,p));
                         sugar = MAX(sugar,SG(p)+TD(DL(mul)));                          sugar = MAX(sugar,SG(p)+TD(DL(mul)));
                         if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {                          if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {
                                 nd_removecont2(d,g);                                  nd_removecont2(d,g);
Line 1106  int nd_nf(int mod,ND g,int full,ND *rp)
Line 1164  int nd_nf(int mod,ND g,int full,ND *rp)
         return 1;          return 1;
 }  }
   
 int nd_nf_pbucket(int mod,ND g,int full,ND *rp)  int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp)
 {  {
         int hindex,index;          int hindex,index;
         NDV p;          NDV p;
Line 1149  int nd_nf_pbucket(int mod,ND g,int full,ND *rp)
Line 1207  int nd_nf_pbucket(int mod,ND g,int full,ND *rp)
                                 *rp = 0;                                  *rp = 0;
                                 return 0;                                  return 0;
                         }                          }
                           p = ps[index];
                         if ( mod ) {                          if ( mod ) {
                                 p = nd_ps[index];  
                                 c1 = invm(HCM(p),mod); c2 = mod-HCM(g);                                  c1 = invm(HCM(p),mod); c2 = mod-HCM(g);
                                 DMAR(c1,c2,0,mod,c); CM(mul) = c;                                  DMAR(c1,c2,0,mod,c); CM(mul) = c;
                         } else {                          } else {
                                 p = nd_psq[index];  
                                 igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);                                  igcd_cofactor(HCQ(g),HCQ(p),&gcd,&cg,&cred);
                                 chsgnq(cg,&CQ(mul));                                  chsgnq(cg,&CQ(mul));
                                 nd_mul_c_q(d,cred);                                  nd_mul_c_q(d,cred);
Line 1162  int nd_nf_pbucket(int mod,ND g,int full,ND *rp)
Line 1219  int nd_nf_pbucket(int mod,ND g,int full,ND *rp)
                                 g = bucket->body[hindex];                                  g = bucket->body[hindex];
                                 gmag = (double)p_mag((P)HCQ(g));                                  gmag = (double)p_mag((P)HCQ(g));
                         }                          }
                         red = ndv_mul_nm(mod,p,mul);                          red = ndv_mul_nm(mod,mul,p);
                         bucket->body[hindex] = nd_remove_head(g);                          bucket->body[hindex] = nd_remove_head(g);
                         red = nd_remove_head(red);                          red = nd_remove_head(red);
                         add_pbucket(mod,bucket,red);                          add_pbucket(mod,bucket,red);
Line 1251  int nd_nf_direct(int mod,ND g,BaseSet base,int full,ND
Line 1308  int nd_nf_direct(int mod,ND g,BaseSet base,int full,ND
                                 chsgnq(cg,&CQ(mul));                                  chsgnq(cg,&CQ(mul));
                                 nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);                                  nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);
                         }                          }
                         g = nd_add(mod,g,ndv_mul_nm(mod,p,mul));                          g = nd_add(mod,g,ndv_mul_nm(mod,mul,p));
                         sugar = MAX(sugar,SG(p)+TD(DL(mul)));                          sugar = MAX(sugar,SG(p)+TD(DL(mul)));
                         if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {                          if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {
                                 nd_removecont2(d,g);                                  nd_removecont2(d,g);
Line 1301  int nd_check_candidate(NODE input,NODE cand)
Line 1358  int nd_check_candidate(NODE input,NODE cand)
         for ( t = input; t; t = NEXT(t) ) {          for ( t = input; t; t = NEXT(t) ) {
 again:  again:
                 d = dptond(0,(DP)BDY(t));                  d = dptond(0,(DP)BDY(t));
                 stat = nd_nf(0,d,0,&nf);                  stat = nd_nf(0,d,nd_ps,0,&nf);
                 if ( !stat ) {                  if ( !stat ) {
                         nd_reconstruct(0,0,0);                          nd_reconstruct(0,0,0);
                         goto again;                          goto again;
Line 1494  again:
Line 1551  again:
                         sugar = SG(l);                          sugar = SG(l);
                         fprintf(asir_out,"%d",sugar);                          fprintf(asir_out,"%d",sugar);
                 }                  }
                 stat = nd_sp(m,l,&h);                  stat = nd_sp(m,0,l,&h);
                 if ( !stat ) {                  if ( !stat ) {
                         NEXT(l) = d; d = l;                          NEXT(l) = d; d = l;
                         d = nd_reconstruct(m,0,d);                          d = nd_reconstruct(m,0,d);
                         goto again;                          goto again;
                 }                  }
 #if USE_GEOBUCKET  #if USE_GEOBUCKET
                 stat = m?nd_nf_pbucket(m,h,!Top,&nf):nd_nf(m,h,!Top,&nf);                  stat = m?nd_nf_pbucket(m,h,nd_ps,!Top,&nf):nd_nf(m,h,nd_ps,!Top,&nf);
 #else  #else
                 stat = nd_nf(m,h,!Top,&nf);                  stat = nd_nf(m,h,nd_ps,!Top,&nf);
 #endif  #endif
                 if ( !stat ) {                  if ( !stat ) {
                         NEXT(l) = d; d = l;                          NEXT(l) = d; d = l;
Line 1512  again:
Line 1569  again:
                 } else if ( nf ) {                  } else if ( nf ) {
                         if ( checkonly ) return 0;                          if ( checkonly ) return 0;
                         printf("+"); fflush(stdout);                          printf("+"); fflush(stdout);
                         nh = m?nd_newps(m,nf,0):nd_newps(m,0,nf);                          nh = nd_newps(m,nf,0);
                         d = update_pairs(d,g,nh);                          d = update_pairs(d,g,nh);
                         g = update_base(g,nh);                          g = update_base(g,nh);
                         FREENDP(l);                          FREENDP(l);
Line 1521  again:
Line 1578  again:
                         FREENDP(l);                          FREENDP(l);
                 }                  }
         }          }
         if ( m )          for ( t = g; t; t = NEXT(t) ) BDY(t) = (pointer)nd_ps[(int)BDY(t)];
                 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;          return g;
 }  }
   
Line 1549  again:
Line 1603  again:
                         sugar = SG(l);                          sugar = SG(l);
                         fprintf(asir_out,"%d",sugar);                          fprintf(asir_out,"%d",sugar);
                 }                  }
                 stat = nd_sp(m,l,&h);                  stat = nd_sp(m,0,l,&h);
                 if ( !stat ) {                  if ( !stat ) {
                         NEXT(l) = d; d = l;                          NEXT(l) = d; d = l;
                         d = nd_reconstruct(m,1,d);                          d = nd_reconstruct(m,1,d);
                         goto again;                          goto again;
                 }                  }
 #if USE_GEOBUCKET  #if USE_GEOBUCKET
                 stat = nd_nf_pbucket(m,h,!Top,&nf);                  stat = nd_nf_pbucket(m,h,nd_ps,!Top,&nf);
 #else  #else
                 stat = nd_nf(m,h,!Top,&nf);                  stat = nd_nf(m,h,nd_ps,!Top,&nf);
 #endif  #endif
                 if ( !stat ) {                  if ( !stat ) {
                         NEXT(l) = d; d = l;                          NEXT(l) = d; d = l;
Line 1566  again:
Line 1620  again:
                         goto again;                          goto again;
                 } else if ( nf ) {                  } else if ( nf ) {
                         /* overflow does not occur */                          /* overflow does not occur */
                         nd_sp(0,l,&h);                          nd_sp(0,1,l,&h);
                         nd_nf(0,h,!Top,&nfq);                          nd_nf(0,h,nd_ps_trace,!Top,&nfq);
                         if ( nfq ) {                          if ( nfq ) {
                                 printf("+"); fflush(stdout);                                  printf("+"); fflush(stdout);
                                 nh = nd_newps(m,nf,nfq);                                  nh = nd_newps(m,nf,nfq);
Line 1584  again:
Line 1638  again:
                 FREENDP(l);                  FREENDP(l);
         }          }
         for ( t = g; t; t = NEXT(t) )          for ( t = g; t; t = NEXT(t) )
                 BDY(t) = (pointer)nd_psq[(int)BDY(t)];                  BDY(t) = (pointer)nd_ps_trace[(int)BDY(t)];
         return g;          return g;
 }  }
   
Line 1610  NODE nd_reduceall(int m,NODE f)
Line 1664  NODE nd_reduceall(int m,NODE f)
         for ( n = 0, t = f; t; t = NEXT(t), n++ );          for ( n = 0, t = f; t; t = NEXT(t), n++ );
         ps = (NDV *)ALLOCA(n*sizeof(NDV));          ps = (NDV *)ALLOCA(n*sizeof(NDV));
         bound = (unsigned int **)ALLOCA(n*sizeof(unsigned int *));          bound = (unsigned int **)ALLOCA(n*sizeof(unsigned int *));
         for ( i = 0, t = f; i < n; i++, t = NEXT(t) ) {          for ( i = 0, t = f; i < n; i++, t = NEXT(t) ) ps[i] = (NDV)BDY(t);
                 ps[i] = (NDV)BDY(t);  
                 bound[i] = ndv_compute_bound(ps[i]);  
         }  
         qsort(ps,n,sizeof(NDV),(int (*)(const void *,const void *))ndv_compare);          qsort(ps,n,sizeof(NDV),(int (*)(const void *,const void *))ndv_compare);
           for ( i = 0; i < n; i++ ) bound[i] = ndv_compute_bound(ps[i]);
         base.ps = (NDV *)ALLOCA((n-1)*sizeof(NDV));          base.ps = (NDV *)ALLOCA((n-1)*sizeof(NDV));
         base.bound = (unsigned int **)ALLOCA((n-1)*sizeof(unsigned int *));          base.bound = (unsigned int **)ALLOCA((n-1)*sizeof(unsigned int *));
         base.len = n-1;          base.len = n-1;
Line 1666  ND_pairs update_pairs( ND_pairs d, NODE /* of index */
Line 1718  ND_pairs update_pairs( ND_pairs d, NODE /* of index */
         d1 = nd_newpairs(g,t);          d1 = nd_newpairs(g,t);
         d1 = crit_M(d1);          d1 = crit_M(d1);
         d1 = crit_F(d1);          d1 = crit_F(d1);
         prev = 0; cur = head = d1;          if ( do_weyl )
         while ( cur ) {                  head = d1;
                 if ( crit_2( cur->i1,cur->i2 ) ) {          else {
                         remove = cur;                  prev = 0; cur = head = d1;
                         if ( !prev ) head = cur = NEXT(cur);                  while ( cur ) {
                         else cur = NEXT(prev) = NEXT(cur);                          if ( crit_2( cur->i1,cur->i2 ) ) {
                         FREENDP(remove);                                  remove = cur;
                 } else {                                  if ( !prev ) head = cur = NEXT(cur);
                         prev = cur; cur = NEXT(cur);                                  else cur = NEXT(prev) = NEXT(cur);
                                   FREENDP(remove);
                           } else {
                                   prev = cur; cur = NEXT(cur);
                           }
                 }                  }
         }          }
         if ( !d )          if ( !d )
Line 1738  ND_pairs crit_B( ND_pairs d, int s )
Line 1794  ND_pairs crit_B( ND_pairs d, int s )
         return head;          return head;
 }  }
   
 /* XXX : check is necessary */  
   
 ND_pairs crit_M( ND_pairs d1 )  ND_pairs crit_M( ND_pairs d1 )
 {  {
         ND_pairs e,d2,d3,dd,p;          ND_pairs e,d2,d3,dd,p;
Line 1899  int nd_newps(int mod,ND a,ND aq)
Line 1953  int nd_newps(int mod,ND a,ND aq)
         RHist r;          RHist r;
         NDV b;          NDV b;
   
           if ( aq ) {
                   /* trace lifting */
                   if ( !rem(NM(HCQ(aq)),mod) )
                           return -1;
           }
         if ( nd_psn == nd_pslen ) {          if ( nd_psn == nd_pslen ) {
                 nd_pslen *= 2;                  nd_pslen *= 2;
                 nd_ps = (NDV *)REALLOC((char *)nd_ps,nd_pslen*sizeof(NDV));                  nd_ps = (NDV *)REALLOC((char *)nd_ps,nd_pslen*sizeof(NDV));
                 nd_psq = (NDV *)REALLOC((char *)nd_psq,nd_pslen*sizeof(NDV));                  nd_ps_trace = (NDV *)REALLOC((char *)nd_ps_trace,nd_pslen*sizeof(NDV));
                 nd_psh = (RHist *)REALLOC((char *)nd_psh,nd_pslen*sizeof(RHist));                  nd_psh = (RHist *)REALLOC((char *)nd_psh,nd_pslen*sizeof(RHist));
                 nd_bound = (unsigned int **)                  nd_bound = (unsigned int **)
                         REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *));                          REALLOC((char *)nd_bound,nd_pslen*sizeof(unsigned int *));
         }          }
         if ( a && aq ) {  
                 /* trace lifting */  
                 if ( !rem(NM(HCQ(aq)),mod) )  
                         return -1;  
         }  
         NEWRHist(r); nd_psh[nd_psn] = r;          NEWRHist(r); nd_psh[nd_psn] = r;
           nd_removecont(mod,a); nd_ps[nd_psn] = ndtondv(mod,a);
         if ( aq ) {          if ( aq ) {
                 nd_removecont(0,aq);                  nd_removecont(0,aq); nd_ps_trace[nd_psn] = ndtondv(0,aq);
                 nd_psq[nd_psn] = ndtondv(0,aq);                  nd_bound[nd_psn] = ndv_compute_bound(nd_ps_trace[nd_psn]);
                 nd_bound[nd_psn] = ndv_compute_bound(nd_psq[nd_psn]);  
                 SG(r) = SG(aq); ndl_copy(HDL(aq),DL(r));                  SG(r) = SG(aq); ndl_copy(HDL(aq),DL(r));
                   nd_free(a); nd_free(aq);
           } else {
                   nd_bound[nd_psn] = ndv_compute_bound(nd_ps[nd_psn]);
                   SG(r) = SG(a); ndl_copy(HDL(a),DL(r));
                   nd_free(a);
         }          }
         if ( a ) {  
                 nd_removecont(mod,a);  
                 nd_ps[nd_psn] = ndtondv(mod,a);  
                 if ( !aq ) {  
                         nd_bound[nd_psn] = ndv_compute_bound(nd_ps[nd_psn]);  
                         SG(r) = SG(a); ndl_copy(HDL(a),DL(r));  
                 }  
         }  
         nd_free(a); nd_free(aq);  
         return nd_psn++;          return nd_psn++;
 }  }
   
Line 1938  void nd_setup(int mod,int trace,NODE f)
Line 1988  void nd_setup(int mod,int trace,NODE f)
         unsigned int *d;          unsigned int *d;
         RHist r;          RHist r;
         NDV a;          NDV a;
           MP t;
   
         nd_found = 0; nd_notfirst = 0; nd_create = 0;          nd_found = 0; nd_notfirst = 0; nd_create = 0;
   
         nd_psn = length(f); nd_pslen = 2*nd_psn;          nd_psn = length(f); nd_pslen = 2*nd_psn;
         nd_ps = (NDV *)MALLOC(nd_pslen*sizeof(NDV));          nd_ps = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
         nd_psq = (NDV *)MALLOC(nd_pslen*sizeof(NDV));          nd_ps_trace = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
         nd_psh = (RHist *)MALLOC(nd_pslen*sizeof(RHist));          nd_psh = (RHist *)MALLOC(nd_pslen*sizeof(RHist));
         nd_bound = (unsigned int **)MALLOC(nd_pslen*sizeof(unsigned int *));          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 )          if ( !nd_red )
                 nd_red = (RHist *)MALLOC(REDTAB_LEN*sizeof(RHist));                  nd_red = (RHist *)MALLOC(REDTAB_LEN*sizeof(RHist));
         bzero(nd_red,REDTAB_LEN*sizeof(RHist));          bzero(nd_red,REDTAB_LEN*sizeof(RHist));
   
           for ( max = 0, s = f; s; s = NEXT(s) )
                   for ( t = BDY((DP)BDY(s)); t; t = NEXT(t) ) {
                           d = t->dl->d;
                           for ( j = 0; j < nd_nvar; j++ ) max = MAX(d[j],max);
                   }
   
         if ( max < 2 ) nd_bpe = 2;          if ( max < 2 ) nd_bpe = 2;
         else if ( max < 4 ) nd_bpe = 4;          else if ( max < 4 ) nd_bpe = 4;
         else if ( max < 64 ) nd_bpe = 6;          else if ( max < 64 ) nd_bpe = 6;
Line 1968  void nd_setup(int mod,int trace,NODE f)
Line 2021  void nd_setup(int mod,int trace,NODE f)
                 NEWRHist(r);                  NEWRHist(r);
                 a = dptondv(mod,(DP)BDY(f)); ndv_removecont(mod,a);                  a = dptondv(mod,(DP)BDY(f)); ndv_removecont(mod,a);
                 SG(r) = HTD(a); ndl_copy(HDL(a),DL(r));                  SG(r) = HTD(a); ndl_copy(HDL(a),DL(r));
                   nd_ps[i] = a;
                 if ( trace ) {                  if ( trace ) {
                         nd_ps[i] = a;  
                         a = dptondv(0,(DP)BDY(f)); ndv_removecont(0,a);                          a = dptondv(0,(DP)BDY(f)); ndv_removecont(0,a);
                         nd_psq[i] = a;                          nd_ps_trace[i] = a;
                 } else {  
                         if ( mod ) nd_ps[i] = a;  
                         else nd_psq[i] = a;  
                 }                  }
                   nd_bound[i] = ndv_compute_bound(a);
                 nd_psh[i] = r;                  nd_psh[i] = r;
         }          }
 }  }
Line 2122  void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru
Line 2172  void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru
 void dltondl(int n,DL dl,unsigned int *r)  void dltondl(int n,DL dl,unsigned int *r)
 {  {
         unsigned int *d;          unsigned int *d;
         int i,j,l,s,ord_l,ord_o;          int i,j,l,s,ord_l;
         struct order_pair *op;          struct order_pair *op;
   
         d = dl->d;          d = dl->d;
Line 2131  void dltondl(int n,DL dl,unsigned int *r)
Line 2181  void dltondl(int n,DL dl,unsigned int *r)
                 l = nd_blockmask->n;                  l = nd_blockmask->n;
                 op = nd_blockmask->order_pair;                  op = nd_blockmask->order_pair;
                 for ( j = 0, s = 0; j < l; j++ ) {                  for ( j = 0, s = 0; j < l; j++ ) {
                         ord_o = op[j].order;  
                         ord_l = op[j].length;                          ord_l = op[j].length;
                         if ( !ord_o )                          for ( i = 0; i < ord_l; i++, s++ ) PUT_EXP(r,s,d[s]);
                                 for ( i = 0; i < ord_l; i++ )  
                                         PUT_EXP(r,s+ord_l-i-1,d[s+i]);  
                         else  
                                 for ( i = 0; i < ord_l; i++ )  
                                         PUT_EXP(r,s+i,d[s+i]);  
                         s += ord_l;  
                 }                  }
                 TD(r) = ndl_weight(r);                  TD(r) = ndl_weight(r);
                 for ( j = 0; j < l; j++ )                  for ( j = 0; j < l; j++ )
                         r[j+1] = ndl_weight_mask(r,j);                          r[j+1] = ndl_weight_mask(r,j);
         } else {          } else {
                 if ( nd_isrlex )                  for ( i = 0; i < n; i++ ) PUT_EXP(r,i,d[i]);
                         for ( i = 0; i < n; i++ ) PUT_EXP(r,n-1-i,d[i]);  
                 else  
                         for ( i = 0; i < n; i++ ) PUT_EXP(r,i,d[i]);  
                 TD(r) = ndl_weight(r);                  TD(r) = ndl_weight(r);
         }          }
 }  }
Line 2157  DL ndltodl(int n,unsigned int *ndl)
Line 2197  DL ndltodl(int n,unsigned int *ndl)
 {  {
         DL dl;          DL dl;
         int *d;          int *d;
         int i,j,l,s,ord_l,ord_o;          int i,j,l,s,ord_l;
         struct order_pair *op;          struct order_pair *op;
   
         NEWDL(dl,n);          NEWDL(dl,n);
Line 2167  DL ndltodl(int n,unsigned int *ndl)
Line 2207  DL ndltodl(int n,unsigned int *ndl)
                 l = nd_blockmask->n;                  l = nd_blockmask->n;
                 op = nd_blockmask->order_pair;                  op = nd_blockmask->order_pair;
                 for ( j = 0, s = 0; j < l; j++ ) {                  for ( j = 0, s = 0; j < l; j++ ) {
                         ord_o = op[j].order;  
                         ord_l = op[j].length;                          ord_l = op[j].length;
                         if ( !ord_o )                          for ( i = 0; i < ord_l; i++, s++ ) d[s] = GET_EXP(ndl,s);
                                 for ( i = 0; i < ord_l; i++ )  
                                         d[s+i] = GET_EXP(ndl,s+ord_l-i-1);  
                         else  
                                 for ( i = 0; i < ord_l; i++ )  
                                         d[s+i] = GET_EXP(ndl,s+i);  
                         s += ord_l;  
                 }                  }
         } else {          } else {
                 if ( nd_isrlex )                  for ( i = 0; i < n; i++ ) d[i] = GET_EXP(ndl,i);
                         for ( i = 0; i < n; i++ )  
                                 d[i] = GET_EXP(ndl,n-1-i);  
                 else  
                         for ( i = 0; i < n; i++ )  
                                 d[i] = GET_EXP(ndl,i);  
         }          }
         return dl;          return dl;
 }  }
Line 2237  DP ndtodp(int mod,ND p)
Line 2265  DP ndtodp(int mod,ND p)
 void ndl_print(unsigned int *dl)  void ndl_print(unsigned int *dl)
 {  {
         int n;          int n;
         int i,j,l,ord_o,ord_l,s,s0;          int i,j,l,ord_l,s,s0;
         struct order_pair *op;          struct order_pair *op;
   
         n = nd_nvar;          n = nd_nvar;
Line 2246  void ndl_print(unsigned int *dl)
Line 2274  void ndl_print(unsigned int *dl)
                 l = nd_blockmask->n;                  l = nd_blockmask->n;
                 op = nd_blockmask->order_pair;                  op = nd_blockmask->order_pair;
                 for ( j = 0, s = s0 = 0; j < l; j++ ) {                  for ( j = 0, s = s0 = 0; j < l; j++ ) {
                         ord_o = op[j].order;  
                         ord_l = op[j].length;                          ord_l = op[j].length;
                         if ( !ord_o )                          for ( i = 0; i < ord_l; i++, s++ )
                                 for ( i = 0, s0 += ord_l; i < ord_l; i++, s++ )                                  printf(s==n-1?"%d":"%d,",GET_EXP(dl,s));
                                         printf(s==n-1?"%d":"%d,",GET_EXP(dl,s0-i-1));  
                         else  
                                 for ( i = 0; i < ord_l; i++, s++ )  
                                         printf(s==n-1?"%d":"%d,",GET_EXP(dl,s));  
                 }                  }
         } else {          } else {
                 if ( nd_isrlex )                  for ( i = 0; i < n; i++ ) printf(i==n-1?"%d":"%d,",GET_EXP(dl,i));
                         for ( i = 0; i < n; i++ ) printf(i==n-1?"%d":"%d,",GET_EXP(dl,n-1-i));  
                 else  
                         for ( i = 0; i < n; i++ ) printf(i==n-1?"%d":"%d,",GET_EXP(dl,i));  
         }          }
         printf(">>");          printf(">>");
 }  }
Line 2490  void nd_append_red(unsigned int *d,int i)
Line 2510  void nd_append_red(unsigned int *d,int i)
         nd_red[h] = m;          nd_red[h] = m;
 }  }
   
 unsigned int *dp_compute_bound(DP p)  
 {  
         unsigned int *d,*d1,*d2,*t;  
         MP m;  
         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);  
         d = DL(m)->d;  
         for ( i = 0; i < nd_nvar; i++ ) d1[i] = d[i];  
         for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) {  
                 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;  
         }  
         l = (nd_nvar+31);  
         t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));  
         for ( i = 0; i < nd_nvar; i++ ) t[i] = d1[i];  
         for ( ; i < l; i++ ) t[i] = 0;  
         return t;  
 }  
   
 unsigned int *ndv_compute_bound(NDV p)  unsigned int *ndv_compute_bound(NDV p)
 {  {
         unsigned int *d1,*d2,*t;          unsigned int *d1,*d2,*t;
         int i,l,len;          unsigned int u;
           int i,j,k,l,len,ind;
         NMV m;          NMV m;
   
         if ( !p )          if ( !p )
Line 2534  unsigned int *ndv_compute_bound(NDV p)
Line 2529  unsigned int *ndv_compute_bound(NDV p)
         }          }
         l = nd_nvar+31;          l = nd_nvar+31;
         t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));          t = (unsigned int *)MALLOC_ATOMIC(l*sizeof(unsigned int));
         for ( i = 0; i < nd_nvar; i++ ) t[i] = GET_EXP(d1,i);          for ( i = nd_exporigin, ind = 0; i < nd_wpd; i++ ) {
         for ( ; i < l; i++ ) t[i] = 0;                  u = d1[i];
                   k = (nd_epw-1)*nd_bpe;
                   for ( j = 0; j < nd_epw; j++, k -= nd_bpe, ind++ )
                           t[ind] = (u>>k)&nd_mask0;
           }
           for ( ; ind < l; ind++ ) t[ind] = 0;
         return t;          return t;
 }  }
   
Line 2554  int nd_get_exporigin(struct order_spec *ord)
Line 2554  int nd_get_exporigin(struct order_spec *ord)
 }  }
   
 void nd_setup_parameters() {  void nd_setup_parameters() {
         int i,n,elen;          int i,j,n,elen,ord_o,ord_l,l,s;
           struct order_pair *op;
   
         nd_epw = (sizeof(unsigned int)*8)/nd_bpe;          nd_epw = (sizeof(unsigned int)*8)/nd_bpe;
         elen = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);          elen = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
   
         nd_exporigin = nd_get_exporigin(nd_ord);          nd_exporigin = nd_get_exporigin(nd_ord);
         nd_wpd = nd_exporigin+elen;          nd_wpd = nd_exporigin+elen;
         nd_epos = (EPOS)MALLOC_ATOMIC(nd_nvar*sizeof(struct oEPOS));  
         for ( i = 0; i < nd_nvar; i++ ) {  
                 nd_epos[i].i = nd_exporigin + i/nd_epw;  
                 nd_epos[i].s = (nd_epw-(i%nd_epw)-1)*nd_bpe;  
         }  
         if ( nd_bpe < 32 ) {          if ( nd_bpe < 32 ) {
                 nd_mask0 = (1<<nd_bpe)-1;                  nd_mask0 = (1<<nd_bpe)-1;
         } else {          } else {
Line 2579  void nd_setup_parameters() {
Line 2576  void nd_setup_parameters() {
         }          }
         nm_adv = sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int);          nm_adv = sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int);
         nmv_adv = sizeof(struct oNMV)+(nd_wpd-1)*sizeof(unsigned int);          nmv_adv = sizeof(struct oNMV)+(nd_wpd-1)*sizeof(unsigned int);
           nd_epos = nd_create_epos(nd_ord);
         nd_blockmask = nd_create_blockmask(nd_ord);          nd_blockmask = nd_create_blockmask(nd_ord);
 }  }
   
Line 2607  ND_pairs nd_reconstruct(int mod,int trace,ND_pairs d)
Line 2605  ND_pairs nd_reconstruct(int mod,int trace,ND_pairs d)
         prev_ndp_free_list = _ndp_free_list;          prev_ndp_free_list = _ndp_free_list;
         _nm_free_list = 0;          _nm_free_list = 0;
         _ndp_free_list = 0;          _ndp_free_list = 0;
         if ( mod != 0 )          for ( i = nd_psn-1; i >= 0; i-- ) ndv_realloc(nd_ps[i],obpe,oadv,oepos);
                 for ( i = nd_psn-1; i >= 0; i-- ) ndv_realloc(nd_ps[i],obpe,oadv,oepos);          if ( trace )
         if ( !mod || trace )                  for ( i = nd_psn-1; i >= 0; i-- )
                 for ( i = nd_psn-1; i >= 0; i-- ) ndv_realloc(nd_psq[i],obpe,oadv,oepos);                          ndv_realloc(nd_ps_trace[i],obpe,oadv,oepos);
         s0 = 0;          s0 = 0;
         for ( t = d; t; t = NEXT(t) ) {          for ( t = d; t; t = NEXT(t) ) {
                 NEXTND_pairs(s0,s);                  NEXTND_pairs(s0,s);
Line 2682  void nd_reconstruct_direct(int mod,NDV *ps,int len)
Line 2680  void nd_reconstruct_direct(int mod,NDV *ps,int len)
   
 void ndl_reconstruct(int obpe,EPOS oepos,unsigned int *d,unsigned int *r)  void ndl_reconstruct(int obpe,EPOS oepos,unsigned int *d,unsigned int *r)
 {  {
         int n,i,ei,oepw,omask0,j,s,ord_l,ord_o,l;          int n,i,ei,oepw,omask0,j,s,ord_l,l;
         struct order_pair *op;          struct order_pair *op;
 #define GET_EXP_OLD(d,a) (((d)[oepos[a].i]>>oepos[a].s)&omask0)  #define GET_EXP_OLD(d,a) (((d)[oepos[a].i]>>oepos[a].s)&omask0)
 #define PUT_EXP_OLD(r,a,e) ((r)[oepos[a].i] |= ((e)<<oepos[a].s))  #define PUT_EXP_OLD(r,a,e) ((r)[oepos[a].i] |= ((e)<<oepos[a].s))
Line 2698  void ndl_reconstruct(int obpe,EPOS oepos,unsigned int 
Line 2696  void ndl_reconstruct(int obpe,EPOS oepos,unsigned int 
                 for ( i = 1; i < nd_exporigin; i++ )                  for ( i = 1; i < nd_exporigin; i++ )
                         r[i] = d[i];                          r[i] = d[i];
                 for ( j = 0, s = 0; j < l; j++ ) {                  for ( j = 0, s = 0; j < l; j++ ) {
                         ord_o = op[j].order;  
                         ord_l = op[j].length;                          ord_l = op[j].length;
                         if ( !ord_o )                          for ( i = 0; i < ord_l; i++, s++ ) {
                                 for ( i = 0; i < ord_l; i++ ) {                                  ei =  GET_EXP_OLD(d,s);
                                         ei =  GET_EXP_OLD(d,s+ord_l-i-1);                                  PUT_EXP(r,s,ei);
                                         PUT_EXP(r,s+ord_l-i-1,ei);                          }
                                 }  
                         else  
                                 for ( i = 0; i < ord_l; i++ ) {  
                                         ei =  GET_EXP_OLD(d,s+i);  
                                         PUT_EXP(r,s+i,ei);  
                                 }  
                         s += ord_l;  
                 }                  }
         } else {          } else {
                 if ( nd_isrlex )                  for ( i = 0; i < n; i++ ) {
                         for ( i = 0; i < n; i++ ) {                          ei = GET_EXP_OLD(d,i);
                                 ei = GET_EXP_OLD(d,n-1-i);                          PUT_EXP(r,i,ei);
                                 PUT_EXP(r,n-1-i,ei);                  }
                         }  
                 else  
                         for ( i = 0; i < n; i++ ) {  
                                 ei = GET_EXP_OLD(d,i);  
                                 PUT_EXP(r,i,ei);  
                         }  
         }          }
 }  }
   
Line 2747  ND nd_copy(ND p)
Line 2731  ND nd_copy(ND p)
         }          }
 }  }
   
 int nd_sp(int mod,ND_pairs p,ND *rp)  int nd_sp(int mod,int trace,ND_pairs p,ND *rp)
 {  {
         NM m;          NM m;
         NDV p1,p2;          NDV p1,p2;
Line 2755  int nd_sp(int mod,ND_pairs p,ND *rp)
Line 2739  int nd_sp(int mod,ND_pairs p,ND *rp)
         unsigned int *lcm;          unsigned int *lcm;
         int td;          int td;
   
         if ( mod ) {          if ( trace ) {
                 p1 = nd_ps[p->i1]; p2 = nd_ps[p->i2];                  p1 = nd_ps_trace[p->i1]; p2 = nd_ps_trace[p->i2];
         } else {          } else {
                 p1 = nd_psq[p->i1]; p2 = nd_psq[p->i2];                  p1 = nd_ps[p->i1]; p2 = nd_ps[p->i2];
         }          }
         lcm = LCM(p);          lcm = LCM(p);
         NEWNM(m);          NEWNM(m);
         CQ(m) = HCQ(p2);          CQ(m) = HCQ(p2);
         ndl_sub(lcm,HDL(p1),DL(m));          ndl_sub(lcm,HDL(p1),DL(m));
         if ( ndl_check_bound2(p->i1,DL(m)) ) return 0;          if ( ndl_check_bound2(p->i1,DL(m)) )
         t1 = ndv_mul_nm(mod,p1,m);                  return 0;
           t1 = ndv_mul_nm(mod,m,p1);
         if ( mod ) CM(m) = mod-HCM(p1);          if ( mod ) CM(m) = mod-HCM(p1);
         else chsgnq(HCQ(p1),&CQ(m));          else chsgnq(HCQ(p1),&CQ(m));
         ndl_sub(lcm,HDL(p2),DL(m));          ndl_sub(lcm,HDL(p2),DL(m));
Line 2773  int nd_sp(int mod,ND_pairs p,ND *rp)
Line 2758  int nd_sp(int mod,ND_pairs p,ND *rp)
                 nd_free(t1);                  nd_free(t1);
                 return 0;                  return 0;
         }          }
         t2 = ndv_mul_nm(mod,p2,m);          t2 = ndv_mul_nm(mod,m,p2);
         *rp = nd_add(mod,t1,t2);          *rp = nd_add(mod,t1,t2);
         FREENM(m);          FREENM(m);
         return 1;          return 1;
Line 2804  void ndv_mul_c_q(NDV p,Q mul)
Line 2789  void ndv_mul_c_q(NDV p,Q mul)
         }          }
 }  }
   
 ND ndv_mul_nm(int mod,NDV p,NM m0)  ND weyl_ndv_mul_nm(int mod,NM m0,NDV p) {
           int n2,i,j,l,n,tlen;
           unsigned int *d0;
           NM *tab,*psum;
           ND s,r;
           NM t;
           NMV m1;
   
           if ( !p ) return 0;
           n = NV(p); n2 = n>>1;
           d0 = DL(m0);
           l = LEN(p);
           for ( i = 0, tlen = 1; i < n2; i++ ) tlen *= (GET_EXP(d0,n2+i)+1);
           tab = (NM *)ALLOCA(tlen*sizeof(NM));
           psum = (NM *)ALLOCA(tlen*sizeof(NM));
           for ( i = 0; i < tlen; i++ ) psum[i] = 0;
           m1 = (NMV)(((char *)BDY(p))+nmv_adv*(l-1));
           for ( i = l-1; i >= 0; i--, NMV_PREV(m1) ) {
                   /* m0(NM) * m1(NMV) => tab(NM) */
                   weyl_mul_nm_nmv(n,mod,m0,m1,tab,tlen);
                   for ( j = 0; j < tlen; j++ ) {
                           if ( tab[j] ) {
                                   NEXT(tab[j]) = psum[j]; psum[j] = tab[j];
                           }
                   }
           }
           for ( i = tlen-1, r = 0; i >= 0; i-- )
                   if ( psum[i] ) {
                           for ( j = 0, t = psum[i]; t; t = NEXT(t), j++ );
                           MKND(n,psum[i],j,s);
                           r = nd_add(mod,r,s);
                   }
           if ( r ) SG(r) = SG(p)+TD(d0);
           return r;
   }
   
   /* product of monomials */
   /* XXX block order is not handled correctly */
   
   void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *tab,int tlen)
 {  {
           int i,n2,j,s,curlen,homo,h,a,b,k,l,u,min;
           unsigned int *d0,*d1,*d,*dt,*ctab;
           Q *ctab_q;
           Q q,q1;
           unsigned int c0,c1,c;
           NM *p;
           NM m,t;
   
           for ( i = 0; i < tlen; i++ ) tab[i] = 0;
           if ( !m0 || !m1 ) return;
           d0 = DL(m0); d1 = DL(m1); n2 = n>>1;
           NEWNM(m); d = DL(m);
           if ( mod ) {
                   c0 = CM(m0); c1 = CM(m1); DMAR(c0,c1,0,mod,c); CM(m) = c;
           } else
                   mulq(CQ(m0),CQ(m1),&CQ(m));
           for ( i = 0; i < nd_wpd; i++ ) d[i] = 0;
           homo = n&1 ? 1 : 0;
           if ( homo ) {
                   /* offset of h-degree */
                   h = GET_EXP(d0,n-1)+GET_EXP(d1,n-1);
                   PUT_EXP(DL(m),n-1,h);
                   TD(DL(m)) = h;
                   if ( nd_blockmask ) ndl_set_blockweight(DL(m));
           }
           tab[0] = m;
           NEWNM(m); d = DL(m);
           for ( i = 0, curlen = 1; i < n2; i++ ) {
                   a = GET_EXP(d0,i); b = GET_EXP(d1,n2+i);
                   k = GET_EXP(d0,n2+i); l = GET_EXP(d1,i);
                   /* xi^a*(Di^k*xi^l)*Di^b */
                   a += l; b += k;
                   s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
                   if ( !k || !l ) {
                           for ( j = 0; j < curlen; j++ )
                                   if ( t = tab[j] ) {
                                           dt = DL(t);
                                           PUT_EXP(dt,i,a); PUT_EXP(dt,n2+i,b); TD(dt) += s;
                                           if ( nd_blockmask ) ndl_set_blockweight(dt);
                                   }
                           curlen *= k+1;
                           continue;
                   }
                   min = MIN(k,l);
                   if ( mod ) {
                           ctab = (unsigned int *)ALLOCA((min+1)*sizeof(unsigned int));
                           mkwcm(k,l,mod,ctab);
                   } else {
                           ctab_q = (Q *)ALLOCA((min+1)*sizeof(Q));
                           mkwc(k,l,ctab_q);
                   }
                   for ( j = min; j >= 0; j-- ) {
                           for ( u = 0; u < nd_wpd; u++ ) d[u] = 0;
                           PUT_EXP(d,i,a-j); PUT_EXP(d,n2+i,b-j);
                           h = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i);
                           if ( homo ) {
                                   TD(d) = s;
                                   PUT_EXP(d,n-1,s-h);
                           } else TD(d) = h;
                           if ( nd_blockmask ) ndl_set_blockweight(d);
                           if ( mod ) c = ctab[j];
                           else q = ctab_q[j];
                           p = tab+curlen*j;
                           if ( j == 0 ) {
                                   for ( u = 0; u < curlen; u++, p++ ) {
                                           if ( tab[u] ) {
                                                   ndl_addto(DL(tab[u]),d);
                                                   if ( mod ) {
                                                           c0 = CM(tab[u]); DMAR(c0,c,0,mod,c1); CM(tab[u]) = c1;
                                                   } else {
                                                           mulq(CQ(tab[u]),q,&q1); CQ(tab[u]) = q1;
                                                   }
                                           }
                                   }
                           } else {
                                   for ( u = 0; u < curlen; u++, p++ ) {
                                           if ( tab[u] ) {
                                                   NEWNM(t);
                                                   ndl_add(DL(tab[u]),d,DL(t));
                                                   if ( mod ) {
                                                           c0 = CM(tab[u]); DMAR(c0,c,0,mod,c1); CM(t) = c1;
                                                   } else
                                                           mulq(CQ(tab[u]),q,&CQ(t));
                                                   *p = t;
                                           }
                                   }
                           }
                   }
                   curlen *= k+1;
           }
           FREENM(m);
   }
   
   ND ndv_mul_nm(int mod,NM m0,NDV p)
   {
         NM mr,mr0;          NM mr,mr0;
         NMV m;          NMV m;
         unsigned int *d,*dt,*dm;          unsigned int *d,*dt,*dm;
Line 2814  ND ndv_mul_nm(int mod,NDV p,NM m0)
Line 2933  ND ndv_mul_nm(int mod,NDV p,NM m0)
         ND r;          ND r;
   
         if ( !p ) return 0;          if ( !p ) return 0;
           else if ( do_weyl )
                   return weyl_ndv_mul_nm(mod,m0,p);
         else {          else {
                 n = NV(p); m = BDY(p);                  n = NV(p); m = BDY(p);
                 d = DL(m0);                  d = DL(m0);
Line 2850  void ndv_realloc(NDV p,int obpe,int oadv,EPOS oepos)
Line 2971  void ndv_realloc(NDV p,int obpe,int oadv,EPOS oepos)
         int len,i,k;          int len,i,k;
   
 #define NMV_OPREV(m) (m = (NMV)(((char *)m)-oadv))  #define NMV_OPREV(m) (m = (NMV)(((char *)m)-oadv))
 #define NMV_PREV(m) (m = (NMV)(((char *)m)-nmv_adv))  
   
         if ( p ) {          if ( p ) {
                 m = BDY(p); len = LEN(p);                  m = BDY(p); len = LEN(p);
Line 3037  void nd_init_ord(struct order_spec *ord)
Line 3157  void nd_init_ord(struct order_spec *ord)
                                         nd_isrlex = 0;                                          nd_isrlex = 0;
                                         ndl_compare_function = ndl_lex_compare;                                          ndl_compare_function = ndl_lex_compare;
                                         break;                                          break;
                                   case 11:
                                           /* XXX */
                                           nd_dcomp = 0;
                                           nd_isrlex = 1;
                                           ndl_compare_function = ndl_ww_lex_compare;
                                           break;
                                 default:                                  default:
                                         error("nd_gr : unsupported order");                                          error("nd_gr : unsupported order");
                         }                          }
Line 3075  BlockMask nd_create_blockmask(struct order_spec *ord)
Line 3201  BlockMask nd_create_blockmask(struct order_spec *ord)
                 for ( j = 0; j < l; j++, s++ ) PUT_EXP(t,s,nd_mask0);                  for ( j = 0; j < l; j++, s++ ) PUT_EXP(t,s,nd_mask0);
         }          }
         return bm;          return bm;
   }
   
   EPOS nd_create_epos(struct order_spec *ord)
   {
           int i,j,l,s,ord_l,ord_o;
           EPOS epos;
           struct order_pair *op;
   
           epos = (EPOS)MALLOC_ATOMIC(nd_nvar*sizeof(struct oEPOS));
           switch ( ord->id ) {
                   case 0:
                           if ( nd_isrlex ) {
                                   for ( i = 0; i < nd_nvar; i++ ) {
                                           epos[i].i = nd_exporigin + (nd_nvar-1-i)/nd_epw;
                                           epos[i].s = (nd_epw-((nd_nvar-1-i)%nd_epw)-1)*nd_bpe;
                                   }
                           } else {
                                   for ( i = 0; i < nd_nvar; i++ ) {
                                           epos[i].i = nd_exporigin + i/nd_epw;
                                           epos[i].s = (nd_epw-(i%nd_epw)-1)*nd_bpe;
                                   }
                           }
                           break;
                   case 1:
                           /* block order */
                           l = ord->ord.block.length;
                           op = ord->ord.block.order_pair;
                           for ( j = 0, s = 0; j < l; j++ ) {
                                   ord_o = op[j].order;
                                   ord_l = op[j].length;
                                   if ( !ord_o )
                                           for ( i = 0; i < ord_l; i++ ) {
                                                   epos[s+i].i = nd_exporigin + (s+ord_l-i-1)/nd_epw;
                                                   epos[s+i].s = (nd_epw-((s+ord_l-i-1)%nd_epw)-1)*nd_bpe;
                                           }
                                   else
                                           for ( i = 0; i < ord_l; i++ ) {
                                                   epos[s+i].i = nd_exporigin + (s+i)/nd_epw;
                                                   epos[s+i].s = (nd_epw-((s+i)%nd_epw)-1)*nd_bpe;
                                           }
                                   s += ord_l;
                           }
                           break;
                   case 2:
                           error("nd_create_epos : matrix order is not supported yet.");
           }
           return epos;
 }  }

Legend:
Removed from v.1.52  
changed lines
  Added in v.1.58

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