[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.56

version 1.55, 2003/09/03 07:33:35 version 1.56, 2003/09/04 08:35:09
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.55 2003/09/03 07:33:35 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 291  INLINE int nd_length(ND p);
Line 292  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 1764  ND_pairs crit_B( ND_pairs d, int s )
Line 1764  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 2164  void dltondl(int n,DL dl,unsigned int *r)
Line 2162  void dltondl(int n,DL dl,unsigned int *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 2197  DL ndltodl(int n,unsigned int *ndl)
Line 2192  DL ndltodl(int n,unsigned int *ndl)
                         s += ord_l;                          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 2275  void ndl_print(unsigned int *dl)
Line 2265  void ndl_print(unsigned int *dl)
                                         printf(s==n-1?"%d":"%d,",GET_EXP(dl,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 2581  void nd_setup_parameters() {
Line 2568  void nd_setup_parameters() {
         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));          nd_epos = (EPOS)MALLOC_ATOMIC(nd_nvar*sizeof(struct oEPOS));
         for ( i = 0; i < nd_nvar; i++ ) {          if ( nd_isrlex ) {
                 nd_epos[i].i = nd_exporigin + i/nd_epw;                  for ( i = 0; i < nd_nvar; i++ ) {
                 nd_epos[i].s = (nd_epw-(i%nd_epw)-1)*nd_bpe;                          nd_epos[i].i = nd_exporigin + (nd_nvar-1-i)/nd_epw;
                           nd_epos[i].s = (nd_epw-((nd_nvar-1-i)%nd_epw)-1)*nd_bpe;
                   }
           } else {
                   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;
Line 2732  void ndl_reconstruct(int obpe,EPOS oepos,unsigned int 
Line 2726  void ndl_reconstruct(int obpe,EPOS oepos,unsigned int 
                         s += ord_l;                          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 2771  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 2828  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 2844  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 2865  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 2885  void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *ta
Line 2880  void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *ta
         }          }
         tab[0] = m;          tab[0] = m;
         curlen = 1;          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; 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 */                                          /* XXX other weights */
                                 }                                  }
                         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);
                   } else {
                           ctab_q = (Q *)ALLOCA((min+1)*sizeof(Q));
                           mkwc(k,l,ctab_q);
                   }
                 p = tab+curlen;                  p = tab+curlen;
                 for ( j = 1; j <= min; j++ ) {                  for ( j = 1; j <= min; j++ ) {
                           for ( u = 0; u < nd_wpd; u++ ) d[u] = 0;
                         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 */                          /* XXX other weights */
                         c = ctab[j];                          if ( mod ) c = ctab[j];
                         for ( k = 0; k < curlen; k++, p++ ) {                          else q = ctab_q[j];
                                 NEWNM(t);                          for ( u = 0; u < curlen; u++, p++ ) {
                                 ndl_add(DL(tab[k]),d,DL(t));                                  if ( tab[u] ) {
                                 c0 = CM(tab[k]); DMAR(c0,c,0,mod,c1); CM(t) = c1;                                          NEWNM(t);
                                 *p = 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 */                  /* destructive for j = 0 */
                   for ( u = 0; u < nd_wpd; u++ ) d[u] = 0;
                 PUT_EXP(d,i,a); PUT_EXP(d,n2+i,b);                  PUT_EXP(d,i,a); PUT_EXP(d,n2+i,b);
                 h = s-(MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i));                  h = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,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 */                  /* XXX other weights */
                 c = ctab[0];                  if ( mod ) c = ctab[0];
                   else q = ctab_q[0];
                 p = tab;                  p = tab;
                 for ( k = 0; k < curlen; k++, p++ ) {                  for ( u = 0; u < curlen; u++, p++ ) {
                         ndl_addto(DL(tab[k]),d);                          if ( tab[u] ) {
                         c0 = CM(tab[k]); DMAR(c0,c,0,mod,c1); CM(tab[k]) = c1;                                  ndl_addto(DL(tab[u]),d);
                 }                                  if ( mod ) {
                 curlen *= k+1;                                          c0 = CM(tab[u]); DMAR(c0,c,0,mod,c1); CM(tab[u]) = c1;
         }                                  } else {
         FREENM(m);                                          mulq(CQ(tab[u]),q,&q1); CQ(tab[u]) = q1;
 }  
   
 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);  
                         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));  
                                 mulq(CQ(tab[k]),c,&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 3004  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);
Legend:
Removed from v.1.55  
changed lines
  Added in v.1.56
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>