[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.55 and 1.57

version 1.55, 2003/09/03 07:33:35 version 1.57, 2003/09/05 05:02:53
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.54 2003/08/31 07:42:23 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.56 2003/09/04 08:35:09 noro Exp $ */
   
 #include "ca.h"  #include "ca.h"
 #include "inline.h"  #include "inline.h"
Line 185  if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT
Line 185  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 226  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 272  void nd_reconstruct_direct(int mod,NDV *ps,int len);
Line 274  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 282  ND nd_remove_head(ND p);
Line 285  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);
Line 291  INLINE int nd_length(ND p);
Line 293  INLINE int nd_length(ND p);
 /* NDV functions */  /* NDV functions */
 ND weyl_ndv_mul_nm(int mod,NM m0,NDV p);  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 weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *tab,int tlen);
 void weyl_mul_nm_nmv_q(int n,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);
Line 540  void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned in
Line 541  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 1764  ND_pairs crit_B( ND_pairs d, int s )
Line 1775  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 1960  void nd_setup(int mod,int trace,NODE f)
Line 1969  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;
   
Line 1968  void nd_setup(int mod,int trace,NODE f)
Line 1978  void nd_setup(int mod,int trace,NODE f)
         nd_ps_trace = (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 1990  void nd_setup(int mod,int trace,NODE f)
Line 2002  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;                  nd_ps[i] = a;
                 if ( trace ) {                  if ( trace ) {
                         a = dptondv(0,(DP)BDY(f)); ndv_removecont(0,a);                          a = dptondv(0,(DP)BDY(f)); ndv_removecont(0,a);
                         nd_ps_trace[i] = a;                          nd_ps_trace[i] = a;
                 }                  }
                   nd_bound[i] = ndv_compute_bound(a);
                 nd_psh[i] = r;                  nd_psh[i] = r;
         }          }
 }  }
Line 2141  void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru
Line 2153  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 2150  void dltondl(int n,DL dl,unsigned int *r)
Line 2162  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 2176  DL ndltodl(int n,unsigned int *ndl)
Line 2178  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 2186  DL ndltodl(int n,unsigned int *ndl)
Line 2188  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 2256  DP ndtodp(int mod,ND p)
Line 2246  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 2265  void ndl_print(unsigned int *dl)
Line 2255  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 2509  void nd_append_red(unsigned int *d,int i)
Line 2491  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 2553  unsigned int *ndv_compute_bound(NDV p)
Line 2510  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 2573  int nd_get_exporigin(struct order_spec *ord)
Line 2535  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 2598  void nd_setup_parameters() {
Line 2557  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 2701  void nd_reconstruct_direct(int mod,NDV *ps,int len)
Line 2661  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 2717  void ndl_reconstruct(int obpe,EPOS oepos,unsigned int 
Line 2677  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 2783  int nd_sp(int mod,int trace,ND_pairs p,ND *rp)
Line 2729  int nd_sp(int mod,int trace,ND_pairs p,ND *rp)
         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)) )
                   return 0;
         t1 = ndv_mul_nm(mod,m,p1);          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));
Line 2839  ND weyl_ndv_mul_nm(int mod,NM m0,NDV p) {
Line 2786  ND weyl_ndv_mul_nm(int mod,NM m0,NDV p) {
         tab = (NM *)ALLOCA(tlen*sizeof(NM));          tab = (NM *)ALLOCA(tlen*sizeof(NM));
         psum = (NM *)ALLOCA(tlen*sizeof(NM));          psum = (NM *)ALLOCA(tlen*sizeof(NM));
         for ( i = 0; i < tlen; i++ ) psum[i] = 0;          for ( i = 0; i < tlen; i++ ) psum[i] = 0;
         for ( i = l-1, m1 = BDY(p)+nmv_adv*(l-1); i >= 0; i--, m1 -= nmv_adv ) {          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) */                  /* m0(NM) * m1(NMV) => tab(NM) */
                 if ( mod )                  weyl_mul_nm_nmv(n,mod,m0,m1,tab,tlen);
                         weyl_mul_nm_nmv(mod,n,m0,m1,tab,tlen);  
                 else  
                         weyl_mul_nm_nmv_q(n,m0,m1,tab,tlen);  
                 for ( j = 0; j < tlen; j++ ) {                  for ( j = 0; j < tlen; j++ ) {
                         if ( tab[j] ) {                          if ( tab[j] ) {
                                 NEXT(tab[j]) = psum[j]; psum[j] = tab[j];                                  NEXT(tab[j]) = psum[j]; psum[j] = tab[j];
Line 2857  ND weyl_ndv_mul_nm(int mod,NM m0,NDV p) {
Line 2802  ND weyl_ndv_mul_nm(int mod,NM m0,NDV p) {
                         MKND(n,psum[i],j,s);                          MKND(n,psum[i],j,s);
                         r = nd_add(mod,r,s);                          r = nd_add(mod,r,s);
                 }                  }
         if ( s ) SG(s) = SG(p)+TD(d0);          if ( r ) SG(r) = SG(p)+TD(d0);
         return s;          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)  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,min;          int i,n2,j,s,curlen,homo,h,a,b,k,l,u,min;
         unsigned int *d0,*d1,*d,*ctab;          unsigned int *d0,*d1,*d,*dt,*ctab;
           Q *ctab_q;
           Q q,q1;
         unsigned int c0,c1,c;          unsigned int c0,c1,c;
         NM *p;          NM *p;
         NM m,t;          NM m,t;
Line 2873  void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *ta
Line 2823  void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *ta
         if ( !m0 || !m1 ) return;          if ( !m0 || !m1 ) return;
         d0 = DL(m0); d1 = DL(m1); n2 = n>>1;          d0 = DL(m0); d1 = DL(m1); n2 = n>>1;
         NEWNM(m); d = DL(m);          NEWNM(m); d = DL(m);
         c0 = CM(m0); c1 = CM(m1); DMAR(c1,c,0,mod,c); CM(m) = c;          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;          for ( i = 0; i < nd_wpd; i++ ) d[i] = 0;
         homo = n&1 ? 1 : 0;          homo = n&1 ? 1 : 0;
         if ( homo ) {          if ( homo ) {
Line 2881  void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *ta
Line 2834  void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *ta
                 h = GET_EXP(d0,n-1)+GET_EXP(d1,n-1);                  h = GET_EXP(d0,n-1)+GET_EXP(d1,n-1);
                 PUT_EXP(DL(m),n-1,h);                  PUT_EXP(DL(m),n-1,h);
                 TD(DL(m)) = h;                  TD(DL(m)) = h;
                 /* XXX other weights */                  if ( nd_blockmask ) ndl_set_blockweight(DL(m));
         }          }
         tab[0] = m;          tab[0] = m;
         curlen = 1;  
         s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);  
         NEWNM(m); d = DL(m);          NEWNM(m); d = DL(m);
         for ( i = 0; i < n2; i++ ) {          for ( i = 0, curlen = 1; i < n2; i++ ) {
                 a = GET_EXP(d0,i); b = GET_EXP(d1,n2+i);                  a = GET_EXP(d0,i); b = GET_EXP(d1,n2+i);
                 k = GET_EXP(d0,n2+i); l = GET_EXP(d1,i);                  k = GET_EXP(d0,n2+i); l = GET_EXP(d1,i);
                 /* xi^a*(Di^k*xi^l)*Di^b */                  /* xi^a*(Di^k*xi^l)*Di^b */
                 a += l; b += k;                  a += l; b += k;
                   s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
                 if ( !k || !l ) {                  if ( !k || !l ) {
                         for ( j = 0; j < curlen; j++ )                          for ( j = 0; j < curlen; j++ )
                                 if ( m = tab[j] ) {                                  if ( t = tab[j] ) {
                                         d = DL(m);                                          dt = DL(t);
                                         PUT_EXP(d,i,a); PUT_EXP(d,n2+i,b); TD(d) += s;                                          PUT_EXP(dt,i,a); PUT_EXP(dt,n2+i,b); TD(dt) += s;
                                         /* XXX other weights */                                          if ( nd_blockmask ) ndl_set_blockweight(dt);
                                 }                                  }
                         curlen *= k+1;                          curlen *= k+1;
                         continue;                          continue;
                 }                  }
                 min = MIN(k,l);                  min = MIN(k,l);
                 ctab = (unsigned int *)ALLOCA((min+1)*sizeof(unsigned int));                  if ( mod ) {
                 mkwcm(k,l,mod,ctab);                          ctab = (unsigned int *)ALLOCA((min+1)*sizeof(unsigned int));
                 for ( j = 0; j < nd_wpd; i++ ) d[j] = 0;                          mkwcm(k,l,mod,ctab);
                 p = tab+curlen;                  } else {
                 for ( j = 1; j <= min; j++ ) {                          ctab_q = (Q *)ALLOCA((min+1)*sizeof(Q));
                         PUT_EXP(d,i,a-j); PUT_EXP(d,n2+i,b-j);                          mkwc(k,l,ctab_q);
                         h = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));  
                         if ( homo ) {  
                                 TD(d) = s;  
                                 PUT_EXP(d,n-1,h);  
                         } else TD(d) = h;  
                         /* XXX other weights */  
                         c = ctab[j];  
                         for ( k = 0; k < curlen; k++, p++ ) {  
                                 NEWNM(t);  
                                 ndl_add(DL(tab[k]),d,DL(t));  
                                 c0 = CM(tab[k]); DMAR(c0,c,0,mod,c1); CM(t) = c1;  
                                 *p = t;  
                         }  
                 }                  }
                 /* destructive for j = 0 */                  for ( j = min; j >= 0; j-- ) {
                 PUT_EXP(d,i,a); PUT_EXP(d,n2+i,b);                          for ( u = 0; u < nd_wpd; u++ ) d[u] = 0;
                 h = s-(MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i));  
                 if ( homo ) {  
                         TD(d) = s;  
                         PUT_EXP(d,n-1,h);  
                 } else TD(d) = h;  
                 /* XXX other weights */  
                 c = ctab[0];  
                 p = tab;  
                 for ( k = 0; k < curlen; k++, p++ ) {  
                         ndl_addto(DL(tab[k]),d);  
                         c0 = CM(tab[k]); DMAR(c0,c,0,mod,c1); CM(tab[k]) = c1;  
                 }  
                 curlen *= k+1;  
         }  
         FREENM(m);  
 }  
   
 void weyl_mul_nm_nmv_q(int n,NM m0,NMV m1,NM *tab,int tlen)  
 {  
         int i,n2,j,s,curlen,homo,h,a,b,k,l,min;  
         unsigned int *d0,*d1,*d;  
         Q *ctab;  
         Q 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);  
         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;  
                 /* XXX other weights */  
         }  
         tab[0] = m;  
         curlen = 1;  
         s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);  
         NEWNM(m); d = DL(m);  
         for ( i = 0; 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;  
                 if ( !k || !l ) {  
                         for ( j = 0; j < curlen; j++ )  
                                 if ( m = tab[j] ) {  
                                         d = DL(m);  
                                         PUT_EXP(d,i,a); PUT_EXP(d,n2+i,b); TD(d) += s;  
                                         /* XXX other weights */  
                                 }  
                         curlen *= k+1;  
                         continue;  
                 }  
                 min = MIN(k,l);  
                 ctab = (Q *)ALLOCA((min+1)*sizeof(Q));  
                 mkwc(k,l,ctab);  
                 for ( j = 0; j < nd_wpd; i++ ) d[j] = 0;  
                 p = tab+curlen;  
                 for ( j = 1; j <= min; j++ ) {  
                         PUT_EXP(d,i,a-j); PUT_EXP(d,n2+i,b-j);                          PUT_EXP(d,i,a-j); PUT_EXP(d,n2+i,b-j);
                         h = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));                          h = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i);
                         if ( homo ) {                          if ( homo ) {
                                 TD(d) = s;                                  TD(d) = s;
                                 PUT_EXP(d,n-1,h);                                  PUT_EXP(d,n-1,s-h);
                         } else TD(d) = h;                          } else TD(d) = h;
                         /* XXX other weights */                          if ( nd_blockmask ) ndl_set_blockweight(d);
                         c = ctab[j];                          if ( mod ) c = ctab[j];
                         for ( k = 0; k < curlen; k++, p++ ) {                          else q = ctab_q[j];
                                 NEWNM(t);                          p = tab+curlen*j;
                                 ndl_add(DL(tab[k]),d,DL(t));                          if ( j == 0 ) {
                                 mulq(CQ(tab[k]),c,&CQ(t));                                  for ( u = 0; u < curlen; u++, p++ ) {
                                 *p = t;                                          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;
                                           }
                                   }
                         }                          }
                 }                  }
                 /* destructive for j = 0 */  
                 PUT_EXP(d,i,a); PUT_EXP(d,n2+i,b);  
                 h = s-(MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i));  
                 if ( homo ) {  
                         TD(d) = s;  
                         PUT_EXP(d,n-1,h);  
                 } else TD(d) = h;  
                 /* XXX other weights */  
                 c = ctab[0];  
                 p = tab;  
                 for ( k = 0; k < curlen; k++, p++ ) {  
                         ndl_addto(DL(tab[k]),d);  
                         mulq(CQ(tab[k]),c,&CQ(tab[k]));  
                 }  
                 curlen *= k+1;                  curlen *= k+1;
         }          }
         FREENM(m);          FREENM(m);
Line 3072  void ndv_realloc(NDV p,int obpe,int oadv,EPOS oepos)
Line 2952  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 3297  BlockMask nd_create_blockmask(struct order_spec *ord)
Line 3176  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.55  
changed lines
  Added in v.1.57

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