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

Diff for /OpenXM_contrib2/asir2000/engine/_distm.c between version 1.3 and 1.14

version 1.3, 2000/08/22 05:04:05 version 1.14, 2009/03/16 16:43:02
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/engine/_distm.c,v 1.2 2000/08/21 08:31:27 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/_distm.c,v 1.13 2002/07/24 07:38:55 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "inline.h"  #include "inline.h"
   
 #define NV(p) ((p)->nv)  
 #define C(p) ((p)->c)  
 #if 0  
 #define ITOS(p) (((unsigned int)(p))>>1)  
 #define STOI(i) ((P)((((unsigned int)(i))<<1)|1))  
 #else  
 #define ITOS(p) (((unsigned int)(p))&0x7fffffff)  
 #define STOI(i) ((P)((unsigned int)(i)|0x80000000))  
 #endif  
   
 extern int (*cmpdl)();  extern int (*cmpdl)();
 extern int do_weyl;  extern int do_weyl;
   
 void _dptodp();  
 void adddl_dup();  
 void _mulmdm_dup();  
 void _free_dp();  
 void _comm_mulmd_dup();  
 void _weyl_mulmd_dup();  
 void _weyl_mulmmm_dup();  
 void _weyl_mulmdm_dup();  
   
 MP _mp_free_list;  MP _mp_free_list;
 DP _dp_free_list;  DP _dp_free_list;
 DL _dl_free_list;  DL _dl_free_list;
 int current_dl_length;  int current_dl_length;
   
 #define _NEWDL_NOINIT(d,n) if ((n)!= current_dl_length){_dl_free_list=0; current_dl_length=(n);} if(!_dl_free_list)_DL_alloc(); (d)=_dl_free_list; _dl_free_list = *((DL *)_dl_free_list)  void GC_gcollect();
 #define _NEWDL(d,n) if ((n)!= current_dl_length){_dl_free_list=0; current_dl_length=(n);} if(!_dl_free_list)_DL_alloc(); (d)=_dl_free_list; _dl_free_list = *((DL *)_dl_free_list); bzero((d),(((n)+1)*sizeof(int)))  
 #define _NEWMP(m) if(!_mp_free_list)_MP_alloc(); (m)=_mp_free_list; _mp_free_list = NEXT(_mp_free_list)  
 #define _MKDP(n,m,d) if(!_dp_free_list)_DP_alloc(); (d)=_dp_free_list; _dp_free_list = (DP)BDY(_dp_free_list); (d)->nv=(n); BDY(d)=(m)  
   
 #define _NEXTMP(r,c) \  void _free_private_storage()
 if(!(r)){_NEWMP(r);(c)=(r);}else{_NEWMP(NEXT(c));(c)=NEXT(c);}  {
           _mp_free_list = 0;
           _dp_free_list = 0;
           _dl_free_list = 0;
           GC_gcollect();
   }
   
 #define _NEXTMP2(r,c,s) \  
 if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}  
   
 #define FREEDL(m) *((DL *)m)=_dl_free_list; _dl_free_list=(m)  
 #define FREEMP(m) NEXT(m)=_mp_free_list; _mp_free_list=(m)  
 #define FREEDP(m) BDY(m)=(MP)_dp_free_list; _dp_free_list=(m)  
   
   
 void _DL_alloc()  void _DL_alloc()
 {  {
         int *p;          int *p;
Line 101  void _DL_alloc()
Line 76  void _DL_alloc()
   
 /*      fprintf(stderr,"DL_alloc : %d \n",++DL_alloc_count); */  /*      fprintf(stderr,"DL_alloc : %d \n",++DL_alloc_count); */
         dl_len = (current_dl_length+1);          dl_len = (current_dl_length+1);
         p = (int *)GC_malloc(128*dl_len*sizeof(int));  #if SIZEOF_LONG == 8
           if ( dl_len & 1 )
                   dl_len += 1;
   #endif
         for ( i = 0; i < 128; i++, p += dl_len ) {          for ( i = 0; i < 128; i++, p += dl_len ) {
                   p = (int *)GC_malloc(dl_len*sizeof(int));
                 *(DL *)p = _dl_free_list;                  *(DL *)p = _dl_free_list;
                 _dl_free_list = (DL)p;                  _dl_free_list = (DL)p;
         }          }
Line 115  void _MP_alloc()
Line 94  void _MP_alloc()
         static int MP_alloc_count;          static int MP_alloc_count;
   
 /*      fprintf(stderr,"MP_alloc : %d \n",++MP_alloc_count); */  /*      fprintf(stderr,"MP_alloc : %d \n",++MP_alloc_count); */
         p = (MP)GC_malloc(1024*sizeof(struct oMP));  
         for ( i = 0; i < 1024; i++ ) {          for ( i = 0; i < 1024; i++ ) {
                 p[i].next = _mp_free_list; _mp_free_list = &p[i];                  p = (MP)GC_malloc(sizeof(struct oMP));
                   p->next = _mp_free_list; _mp_free_list = p;
         }          }
 }  }
   
Line 128  void _DP_alloc()
Line 107  void _DP_alloc()
         static int DP_alloc_count;          static int DP_alloc_count;
   
 /*      fprintf(stderr,"DP_alloc : %d \n",++DP_alloc_count); */  /*      fprintf(stderr,"DP_alloc : %d \n",++DP_alloc_count); */
         p = (DP)GC_malloc(1024*sizeof(struct oDP));  
         for ( i = 0; i < 1024; i++ ) {          for ( i = 0; i < 1024; i++ ) {
                 p[i].body = (MP)_dp_free_list; _dp_free_list = &p[i];                  p = (DP)GC_malloc(sizeof(struct oDP));
                   p->body = (MP)_dp_free_list; _dp_free_list = p;
         }          }
 }  }
   
 /* merge p1 and p2 into pr */  /* merge p1 and p2 into pr */
   
 void _addmd_destructive(mod,p1,p2,pr)  void _addmd_destructive(int mod,DP p1,DP p2,DP *pr)
 int mod;  
 DP p1,p2,*pr;  
 {  {
         int n;          int n;
         MP m1,m2,mr,mr0,s;          MP m1,m2,mr,mr0,s;
Line 159  DP p1,p2,*pr;
Line 136  DP p1,p2,*pr;
                                         if ( t ) {                                          if ( t ) {
                                                 _NEXTMP2(mr0,mr,s); C(mr) = STOI(t);                                                  _NEXTMP2(mr0,mr,s); C(mr) = STOI(t);
                                         } else {                                          } else {
                                                 FREEDL(s->dl); FREEMP(s);                                                  _FREEDL(s->dl); _FREEMP(s);
                                         }                                          }
                                         s = m2; m2 = NEXT(m2); FREEDL(s->dl); FREEMP(s);                                          s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
                                         break;                                          break;
                                 case 1:                                  case 1:
                                         s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);                                          s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);
Line 188  DP p1,p2,*pr;
Line 165  DP p1,p2,*pr;
                 _MKDP(NV(p1),mr0,*pr);                  _MKDP(NV(p1),mr0,*pr);
                 if ( *pr )                  if ( *pr )
                         (*pr)->sugar = MAX(p1->sugar,p2->sugar);                          (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                 FREEDP(p1); FREEDP(p2);                  _FREEDP(p1); _FREEDP(p2);
         }          }
 }  }
   
 void _mulmd_dup(mod,p1,p2,pr)  void _mulmd_dup(int mod,DP p1,DP p2,DP *pr)
 int mod;  
 DP p1,p2,*pr;  
 {  {
         if ( !do_weyl )          if ( !do_weyl )
                 _comm_mulmd_dup(mod,p1,p2,pr);                  _comm_mulmd_dup(mod,p1,p2,pr);
Line 202  DP p1,p2,*pr;
Line 177  DP p1,p2,*pr;
                 _weyl_mulmd_dup(mod,p1,p2,pr);                  _weyl_mulmd_dup(mod,p1,p2,pr);
 }  }
   
 void _comm_mulmd_dup(mod,p1,p2,pr)  void _comm_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
 int mod;  
 DP p1,p2,*pr;  
 {  {
         MP m;          MP m;
         DP s,t,u;          DP s,t,u;
Line 236  DP p1,p2,*pr;
Line 209  DP p1,p2,*pr;
         }          }
 }  }
   
 void _weyl_mulmd_dup(mod,p1,p2,pr)  void _weyl_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
 int mod;  
 DP p1,p2,*pr;  
 {  {
         MP m;          MP m;
         DP s,t,u;          DP s,t,u;
         int i,l,l1;          int i,l;
         static MP *w;          static MP *w;
         static int wlen;          static int wlen;
   
         if ( !p1 || !p2 )          if ( !p1 || !p2 )
                 *pr = 0;                  *pr = 0;
         else {          else {
                 for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );                  for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
                 if ( l > wlen ) {                  if ( l > wlen ) {
                         if ( w ) GC_free(w);                          if ( w ) GC_free(w);
                         w = (MP *)MALLOC(l*sizeof(MP));                          w = (MP *)MALLOC(l*sizeof(MP));
                         wlen = l;                          wlen = l;
                 }                  }
                 for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )                  for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
                         w[i] = m;                          w[i] = m;
                 for ( s = 0, i = l-1; i >= 0; i-- ) {                  for ( s = 0, i = l-1; i >= 0; i-- ) {
                         _weyl_mulmdm_dup(mod,p1,w[i],&t); _addmd_destructive(mod,s,t,&u); s = u;                          _weyl_mulmdm_dup(mod,w[i],p2,&t); _addmd_destructive(mod,s,t,&u); s = u;
                 }                  }
                 bzero(w,l*sizeof(MP));                  bzero(w,l*sizeof(MP));
                 *pr = s;                  *pr = s;
         }          }
 }  }
   
 void _mulmdm_dup(mod,p,m0,pr)  void _mulmdm_dup(int mod,DP p,MP m0,DP *pr)
 int mod;  
 DP p;  
 MP m0;  
 DP *pr;  
 {  {
         MP m,mr,mr0;          MP m,mr,mr0;
         DL d,dt,dm;          DL d,dt,dm;
         int c,n,r,i;          int c,n,i,c1,c2;
         int *pt,*p1,*p2;          int *pt,*p1,*p2;
   
         if ( !p )          if ( !p )
Line 282  DP *pr;
Line 249  DP *pr;
                 for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);                  for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
                         m; m = NEXT(m) ) {                          m; m = NEXT(m) ) {
                         _NEXTMP(mr0,mr);                          _NEXTMP(mr0,mr);
                         C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod));                          c1 = ITOS(C(m));
                           DMAR(c1,c,0,mod,c2);
                           C(mr) = (P)STOI(c2);
                         _NEWDL_NOINIT(dt,n); mr->dl = dt;                          _NEWDL_NOINIT(dt,n); mr->dl = dt;
                         dm = m->dl;                          dm = m->dl;
                         dt->td = d->td + dm->td;                          dt->td = d->td + dm->td;
Line 295  DP *pr;
Line 264  DP *pr;
         }          }
 }  }
   
 void _weyl_mulmdm_dup(mod,p,m0,pr)  void _weyl_mulmdm_dup(int mod,MP m0,DP p,DP *pr)
 int mod;  
 DP p;  
 MP m0;  
 DP *pr;  
 {  {
         DP r,t,t1;          DP r,t,t1;
         MP m;          MP m;
         int n,l,i;          DL d0;
         static MP *w;          int n,n2,l,i,j,tlen;
           static MP *w,*psum;
           static struct cdlm *tab;
         static int wlen;          static int wlen;
           static int rtlen;
   
         if ( !p )          if ( !p )
                 *pr = 0;                  *pr = 0;
Line 318  DP *pr;
Line 286  DP *pr;
                 }                  }
                 for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )                  for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
                         w[i] = m;                          w[i] = m;
                 for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {                  n = NV(p); n2 = n>>1;
                         _weyl_mulmmm_dup(mod,w[i],m0,n,&t);                  d0 = m0->dl;
                         _addmd_destructive(mod,r,t,&t1); r = t1;  
                   for ( i = 0, tlen = 1; i < n2; i++ )
                           tlen *= d0->d[n2+i]+1;
                   if ( tlen > rtlen ) {
                           if ( tab ) GC_free(tab);
                           if ( psum ) GC_free(psum);
                           rtlen = tlen;
                           tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm));
                           psum = (MP *)MALLOC(rtlen*sizeof(MP));
                 }                  }
                 bzero(w,l*sizeof(MP));                  bzero(psum,tlen*sizeof(MP));
                   for ( i = l-1; i >= 0; i-- ) {
                           bzero(tab,tlen*sizeof(struct cdlm));
                           _weyl_mulmmm_dup(mod,m0,w[i],n,tab,tlen);
                           for ( j = 0; j < tlen; j++ ) {
                                   if ( tab[j].c ) {
                                           _NEWMP(m); m->dl = tab[j].d;
                                           C(m) = STOI(tab[j].c); NEXT(m) = psum[j];
                                           psum[j] = m;
                                   }
                           }
                   }
                   for ( j = tlen-1, r = 0; j >= 0; j-- )
                           if ( psum[j] ) {
                                   _MKDP(n,psum[j],t); _addmd_destructive(mod,r,t,&t1); r = t1;
                           }
                 if ( r )                  if ( r )
                         r->sugar = p->sugar + m0->dl->td;                          r->sugar = p->sugar + m0->dl->td;
                 *pr = r;                  *pr = r;
         }          }
 }  }
   
 /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */  /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
   
 void _weyl_mulmmm_dup(mod,m0,m1,n,pr)  void _weyl_mulmmm_dup(int mod,MP m0,MP m1,int n,struct cdlm *rtab,int rtablen)
 int mod;  
 MP m0,m1;  
 int n;  
 DP *pr;  
 {  {
         MP m,mr,mr0;          int c,c0,c1;
         DP r,t,t1;          DL d,d0,d1,dt;
         int c,c0,c1,cc;          int i,j,a,b,k,l,n2,s,min,curlen;
         DL d,d0,d1;          struct cdlm *p;
         int i,j,a,b,k,l,n2,s,min,h;          static int *ctab;
         static int *tab;          static struct cdlm *tab;
         static int tablen;          static int tablen;
           static struct cdlm *tmptab;
           static int tmptablen;
   
         if ( !m0 || !m1 )          if ( !m0 || !m1 ) {
                 *pr = 0;                  rtab[0].c = 0;
         else {                  rtab[0].d = 0;
                 c0 = ITOS(C(m0)); c1 = ITOS(C(m1));                  return;
                 c = dmar(c0,c1,0,mod);          }
                 d0 = m0->dl; d1 = m1->dl;          c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
                 n2 = n>>1;          c = dmar(c0,c1,0,mod);
           d0 = m0->dl; d1 = m1->dl;
           n2 = n>>1;
           curlen = 1;
   
                 _NEWDL(d,n);          _NEWDL(d,n);
           if ( n & 1 )
                   /* offset of h-degree */
                   d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
           else
                   d->td = 0;
           rtab[0].c = c;
           rtab[0].d = d;
   
           if ( rtablen > tmptablen ) {
                   if ( tmptab ) GC_free(tmptab);
                   tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm));
                   tmptablen = rtablen;
           }
   
           for ( i = 0; i < n2; i++ ) {
                   a = d0->d[i]; b = d1->d[n2+i];
                   k = d0->d[n2+i]; l = d1->d[i];
   
                   /* degree of 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, p = rtab; j < curlen; j++, p++ ) {
                                   if ( p->c ) {
                                           dt = p->d;
                                           dt->d[i] = a;
                                           dt->d[n2+i] = b;
                                           dt->td += s;
                                   }
                           }
                           curlen *= k+1;
                           continue;
                   }
                   if ( k+1 > tablen ) {
                           if ( tab ) GC_free(tab);
                           if ( ctab ) GC_free(ctab);
                           tablen = k+1;
                           tab = (struct cdlm *)MALLOC(tablen*sizeof(struct cdlm));
                           ctab = (int *)MALLOC(tablen*sizeof(int));
                   }
                   /* compute xi^a*(Di^k*xi^l)*Di^b */
                   min = MIN(k,l);
                   mkwcm(k,l,mod,ctab);
                   bzero(tab,(k+1)*sizeof(struct cdlm));
                   /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
                 if ( n & 1 )                  if ( n & 1 )
                         /* offset of h-degree */                          for ( j = 0; j <= min; j++ ) {
                         d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];                                  _NEWDL(d,n);
                                   d->d[i] = a-j; d->d[n2+i] = b-j;
                                   d->td = s;
                                   d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
                                   tab[j].d = d;
                                   tab[j].c = ctab[j];
                           }
                 else                  else
                         d->td = 0;                          for ( j = 0; j <= min; j++ ) {
                 _NEWMP(mr); mr->c = STOI(c); mr->dl = d; NEXT(mr) = 0;                                  _NEWDL(d,n);
                 _MKDP(n,mr,r); r->sugar = d->td;                                  d->d[i] = a-j; d->d[n2+i] = b-j;
                                   d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i);  /* XXX */
                                   tab[j].d = d;
                                   tab[j].c = ctab[j];
                           }
   #if 0
                   _comm_mulmd_tab(mod,n,rtab,curlen,tab,k+1,tmptab);
                   for ( j = 0; j < curlen; j++ )
                           if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
                   for ( j = 0; j <= min; j++ )
                           if ( tab[j].d ) { _FREEDL(tab[j].d); }
                   curlen *= k+1;
                   bcopy(tmptab,rtab,curlen*sizeof(struct cdlm));
   #else
                   _comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
                   for ( j = 0; j <= min; j++ )
                           if ( tab[j].d ) { _FREEDL(tab[j].d); }
                   curlen *= k+1;
   #endif
           }
   }
   
                 /* homogenized computation; dx-xd=h^2 */  /* direct product of two cdlm tables
                 for ( i = 0; i < n2; i++ ) {    rt[] = [
                         a = d0->d[i]; b = d1->d[n2+i];      t[0]*t1[0],...,t[n-1]*t1[0],
                         k = d0->d[n2+i]; l = d1->d[i];      t[0]*t1[1],...,t[n-1]*t1[1],
                         /* degree of xi^a*(Di^k*xi^l)*Di^b */      ...
                         s = a+k+l+b;      t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
                         /* compute xi^a*(Di^k*xi^l)*Di^b */    ]
                         min = MIN(k,l);  */
   
                         if ( min+1 > tablen ) {  void _comm_mulmd_tab(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1,struct cdlm *rt)
                                 if ( tab ) GC_free(tab);  {
                                 tab = (int *)MALLOC((min+1)*sizeof(int));          int i,j;
                                 tablen = min+1;          struct cdlm *p;
           int c;
           DL d;
   
           bzero(rt,n*n1*sizeof(struct cdlm));
           for ( j = 0, p = rt; j < n1; j++ ) {
                   c = t1[j].c;
                   d = t1[j].d;
                   if ( !c )
                           break;
                   for ( i = 0; i < n; i++, p++ ) {
                           if ( t[i].c ) {
                                   p->c = dmar(t[i].c,c,0,mod);
                                   _adddl_dup(nv,t[i].d,d,&p->d);
                         }                          }
                         mkwcm(k,l,mod,tab);  
                         if ( n & 1 )  
                                 for ( mr0 = 0, j = 0; j <= min; j++ ) {  
                                         _NEXTMP(mr0,mr); _NEWDL(d,n);  
                                         d->d[i] = l-j+a; d->d[n2+i] = k-j+b;  
                                         d->td = s;  
                                         d->d[n-1] = s-(d->d[i]+d->d[n2+i]);  
                                         mr->c = STOI(tab[j]); mr->dl = d;  
                                 }  
                         else  
                                 for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {  
                                         _NEXTMP(mr0,mr); _NEWDL(d,n);  
                                         d->d[i] = l-j+a; d->d[n2+i] = k-j+b;  
                                         d->td = d->d[i]+d->d[n2+i]; /* XXX */  
                                         s = MAX(s,d->td); /* XXX */  
                                         mr->c = STOI(tab[j]); mr->dl = d;  
                                 }  
                         bzero(tab,(min+1)*sizeof(int));  
                         if ( mr0 )  
                                 NEXT(mr) = 0;  
                         _MKDP(n,mr0,t);  
                         if ( t )  
                                 t->sugar = s;  
                         _comm_mulmd_dup(mod,r,t,&t1);  
                         _free_dp(r); _free_dp(t);  
                         r = t1;  
                 }                  }
                 *pr = r;  
         }          }
 }  }
   
 void _dp_red_mod_destructive(p1,p2,mod,rp)  void _comm_mulmd_tab_destructive(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1)
 DP p1,p2;  
 int mod;  
 DP *rp;  
 {  {
         int i,n;          int i,j;
         DL d1,d2,d;          struct cdlm *p;
         MP m;          int c;
         DP t,s;          DL d;
         int c,c1;  
   
         n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;          for ( j = 1, p = t+n; j < n1; j++ ) {
         _NEWDL(d,n); d->td = d1->td - d2->td;                  c = t1[j].c;
         for ( i = 0; i < n; i++ )                  d = t1[j].d;
                 d->d[i] = d1->d[i]-d2->d[i];                  if ( !c )
         c = invm(ITOS(BDY(p2)->c),mod); c1 = dmar(c,ITOS(BDY(p1)->c),0,mod);                          break;
         _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;                  for ( i = 0; i < n; i++, p++ ) {
         _MKDP(n,m,s); s->sugar = d->td;                          if ( t[i].c ) {
         _mulmd_dup(mod,s,p2,&t); _free_dp(s);                                  p->c = dmar(t[i].c,c,0,mod);
         _addmd_destructive(mod,p1,t,rp);                                  _adddl_dup(nv,t[i].d,d,&p->d);
                           }
                   }
           }
           c = t1[0].c;
           d = t1[0].d;
           for ( i = 0, p = t; i < n; i++, p++ )
                   if ( t[i].c ) {
                           p->c = dmar(t[i].c,c,0,mod);
                           /* t[i].d += d */
                           adddl_destructive(nv,t[i].d,d);
                   }
 }  }
   
 void _dp_sp_mod_dup(p1,p2,mod,rp)  void dlto_dl(DL d,DL *dr)
 DP p1,p2;  
 int mod;  
 DP *rp;  
 {  {
         int i,n,td;          int i,n;
         int *w;          DL t;
         DL d1,d2,d;  
         MP m;  
         DP t,s,u;  
   
         n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;          n = current_dl_length;
         w = (int *)ALLOCA(n*sizeof(int));          _NEWDL(t,n); *dr = t;
         for ( i = 0, td = 0; i < n; i++ ) {          t->td = d->td;
                 w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];  
         }  
         _NEWDL(d,n); d->td = td - d1->td;  
         for ( i = 0; i < n; i++ )          for ( i = 0; i < n; i++ )
                 d->d[i] = w[i] - d1->d[i];                  t->d[i] = d->d[i];
         _NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;  
         _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p1,&t); _free_dp(s);  
         _NEWDL(d,n); d->td = td - d2->td;  
         for ( i = 0; i < n; i++ )  
                 d->d[i] = w[i] - d2->d[i];  
         _NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;  
         _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p2,&u); _free_dp(s);  
         _addmd_destructive(mod,t,u,rp);  
 }  }
   
 void _dltodl(d,dr)  void _dltodl(DL d,DL *dr)
 DL d;  
 DL *dr;  
 {  {
         int i,n;          int i,n;
         DL t;          DL t;
Line 471  DL *dr;
Line 511  DL *dr;
                 t->d[i] = d->d[i];                  t->d[i] = d->d[i];
 }  }
   
 void adddl_dup(n,d1,d2,dr)  void _adddl_dup(int n,DL d1,DL d2,DL *dr)
 int n;  
 DL d1,d2;  
 DL *dr;  
 {  {
         DL dt;          DL dt;
         int i;          int i;
   
         _NEWDL(dt,n); *dr = dt;          _NEWDL(dt,n);
           *dr = dt;
         dt->td = d1->td + d2->td;          dt->td = d1->td + d2->td;
         for ( i = 0; i < n; i++ )          for ( i = 0; i < n; i++ )
                 dt->d[i] = d1->d[i]+d2->d[i];                  dt->d[i] = d1->d[i]+d2->d[i];
 }  }
   
 void _free_dp(f)  void _free_dlarray(DL *a,int n)
 DP f;  
 {  {
           int i;
   
           for ( i = 0; i < n; i++ ) { _FREEDL(a[i]); }
   }
   
   void _free_dp(DP f)
   {
         MP m,s;          MP m,s;
   
         if ( !f )          if ( !f )
                 return;                  return;
         m = f->body;          m = f->body;
         while ( m ) {          while ( m ) {
                 s = m; m = NEXT(m); FREEDL(s->dl); FREEMP(s);                  s = m; m = NEXT(m); _FREEDL(s->dl); _FREEMP(s);
         }          }
         FREEDP(f);          _FREEDP(f);
 }  }
   
 void _dptodp(p,r)  void dpto_dp(DP p,DP *r)
 DP p;  
 DP *r;  
 {  {
         MP m,mr0,mr;          MP m,mr0,mr;
           DL t;
   
         if ( !p )          if ( !p )
                 *r = 0;                  *r = 0;
         else {          else {
                   /* XXX : dummy call to set current_dl_length */
                   _NEWDL_NOINIT(t,NV(p));
   
                 for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {                  for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
                           _NEXTMP(mr0,mr);
                           dlto_dl(m->dl,&mr->dl);
                           mr->c = m->c;
                   }
                   NEXT(mr) = 0;
                   _MKDP(p->nv,mr0,*r);
                   (*r)->sugar = p->sugar;
           }
   }
   
   void _dptodp(DP p,DP *r)
   {
           MP m,mr0,mr;
   
           if ( !p )
                   *r = 0;
           else {
                   for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
                         NEXTMP(mr0,mr);                          NEXTMP(mr0,mr);
                         _dltodl(m->dl,&mr->dl);                          _dltodl(m->dl,&mr->dl);
                         mr->c = m->c;                          mr->c = m->c;
Line 519  DP *r;
Line 583  DP *r;
         }          }
 }  }
   
 void _dp_nf_mod_destructive(b,g,ps,mod,full,rp)  /*
 NODE b;   * destructive merge of two list
 DP g;   *
 DP *ps;   * p1, p2 : list of DL
 int mod,full;   * return : a merged list
 DP *rp;   */
   
   NODE _symb_merge(NODE m1,NODE m2,int n)
 {  {
         DP u,p,d,s,t;          NODE top,prev,cur,m,t;
         NODE l;  
         MP m,mr,mrd;  
         int sugar,psugar,n,h_reducible,i;  
   
         if ( !g ) {          if ( !m1 )
                 *rp = 0; return;                  return m2;
           else if ( !m2 )
                   return m1;
           else {
                   switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
                           case 0:
                                   top = m1; _FREEDL((DL)BDY(m2)); m = NEXT(m2);
                                   break;
                           case 1:
                                   top = m1; m = m2;
                                   break;
                           case -1:
                                   top = m2; m = m1;
                                   break;
                   }
                   prev = top; cur = NEXT(top);
                   /* BDY(prev) > BDY(m) always holds */
                   while ( cur && m ) {
                           switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
                                   case 0:
                                           _FREEDL(BDY(m)); m = NEXT(m);
                                           prev = cur; cur = NEXT(cur);
                                           break;
                                   case 1:
                                           t = NEXT(cur); NEXT(cur) = m; m = t;
                                           prev = cur; cur = NEXT(cur);
                                           break;
                                   case -1:
                                           NEXT(prev) = m; m = cur;
                                           prev = NEXT(prev); cur = NEXT(prev);
                                           break;
                           }
                   }
                   if ( !cur )
                           NEXT(prev) = m;
                   return top;
         }          }
         sugar = g->sugar;  }
         n = g->nv;  
         for ( d = 0; g; ) {  /* merge p1 and p2 into pr */
                 for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {  
                         if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {  void _addd_destructive(VL vl,DP p1,DP p2,DP *pr)
                                 h_reducible = 1;  {
                                 psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;          int n;
                                 _dp_red_mod_destructive(g,p,mod,&u); g = u;          MP m1,m2,mr,mr0,s;
                                 sugar = MAX(sugar,psugar);          P t;
                                 if ( !g ) {  
                                         if ( d )          if ( !p1 )
                                                 d->sugar = sugar;                  *pr = p2;
                                         _dptodp(d,rp); _free_dp(d); return;          else if ( !p2 )
                   *pr = p1;
           else {
                   for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
                           switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                                   case 0:
                                           addp(vl,C(m1),C(m2),&t);
                                           s = m1; m1 = NEXT(m1);
                                           if ( t ) {
                                                   _NEXTMP2(mr0,mr,s); C(mr) = t;
                                           } else {
                                                   _FREEDL(s->dl); _FREEMP(s);
                                           }
                                           s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
                                           break;
                                   case 1:
                                           s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);
                                           break;
                                   case -1:
                                           s = m2; m2 = NEXT(m2); _NEXTMP2(mr0,mr,s);
                                           break;
                           }
                   if ( !mr0 )
                           if ( m1 )
                                   mr0 = m1;
                           else if ( m2 )
                                   mr0 = m2;
                           else {
                                   *pr = 0;
                                   return;
                           }
                   else if ( m1 )
                           NEXT(mr) = m1;
                   else if ( m2 )
                           NEXT(mr) = m2;
                   else
                           NEXT(mr) = 0;
                   _MKDP(NV(p1),mr0,*pr);
                   if ( *pr )
                           (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                   _FREEDP(p1); _FREEDP(p2);
           }
   }
   
   void _muld_dup(VL vl,DP p1,DP p2,DP *pr)
   {
           if ( !do_weyl )
                   _comm_muld_dup(vl,p1,p2,pr);
           else
                   _weyl_muld_dup(vl,p1,p2,pr);
   }
   
   void _comm_muld_dup(VL vl,DP p1,DP p2,DP *pr)
   {
           MP m;
           DP s,t,u;
           int i,l,l1;
           static MP *w;
           static int wlen;
   
           if ( !p1 || !p2 )
                   *pr = 0;
           else {
                   for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
                   for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                   if ( l1 < l ) {
                           t = p1; p1 = p2; p2 = t;
                           l = l1;
                   }
                   if ( l > wlen ) {
                           if ( w ) GC_free(w);
                           w = (MP *)MALLOC(l*sizeof(MP));
                           wlen = l;
                   }
                   for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                           w[i] = m;
                   for ( s = 0, i = l-1; i >= 0; i-- ) {
                           _muldm_dup(vl,p1,w[i],&t); _addd_destructive(vl,s,t,&u); s = u;
                   }
                   bzero(w,l*sizeof(MP));
                   *pr = s;
           }
   }
   
   void _weyl_muld_dup(VL vl,DP p1,DP p2,DP *pr)
   {
           MP m;
           DP s,t,u;
           int i,l;
           static MP *w;
           static int wlen;
   
           if ( !p1 || !p2 )
                   *pr = 0;
           else {
                   for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
                   if ( l > wlen ) {
                           if ( w ) GC_free(w);
                           w = (MP *)MALLOC(l*sizeof(MP));
                           wlen = l;
                   }
                   for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
                           w[i] = m;
                   for ( s = 0, i = l-1; i >= 0; i-- ) {
                           _weyl_muldm_dup(vl,w[i],p2,&t); _addd_destructive(vl,s,t,&u); s = u;
                   }
                   bzero(w,l*sizeof(MP));
                   *pr = s;
           }
   }
   
   void _muldm_dup(VL vl,DP p,MP m0,DP *pr)
   {
           MP m,mr,mr0;
           DL d,dt,dm;
           P c;
           int n,i;
           int *pt,*p1,*p2;
   
           if ( !p )
                   *pr = 0;
           else {
                   for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
                           m; m = NEXT(m) ) {
                           _NEXTMP(mr0,mr);
                           mulp(vl,C(m),c,&C(mr));
                           _NEWDL_NOINIT(dt,n); mr->dl = dt;
                           dm = m->dl;
                           dt->td = d->td + dm->td;
                           for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
                                   *pt++ = *p1++ + *p2++;
                   }
                   NEXT(mr) = 0; _MKDP(NV(p),mr0,*pr);
                   if ( *pr )
                           (*pr)->sugar = p->sugar + m0->dl->td;
           }
   }
   
   void _weyl_muldm_dup(VL vl,MP m0,DP p,DP *pr)
   {
           DP r,t,t1;
           MP m;
           DL d0;
           int n,n2,l,i,j,tlen;
           static MP *w,*psum;
           static struct cdl *tab;
           static int wlen;
           static int rtlen;
   
           if ( !p )
                   *pr = 0;
           else {
                   for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
                   if ( l > wlen ) {
                           if ( w ) GC_free(w);
                           w = (MP *)MALLOC(l*sizeof(MP));
                           wlen = l;
                   }
                   for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
                           w[i] = m;
                   n = NV(p); n2 = n>>1;
                   d0 = m0->dl;
   
                   for ( i = 0, tlen = 1; i < n2; i++ )
                           tlen *= d0->d[n2+i]+1;
                   if ( tlen > rtlen ) {
                           if ( tab ) GC_free(tab);
                           if ( psum ) GC_free(psum);
                           rtlen = tlen;
                           tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
                           psum = (MP *)MALLOC(rtlen*sizeof(MP));
                   }
                   bzero(psum,tlen*sizeof(MP));
                   for ( i = l-1; i >= 0; i-- ) {
                           bzero(tab,tlen*sizeof(struct cdl));
                           _weyl_mulmm_dup(vl,m0,w[i],n,tab,tlen);
                           for ( j = 0; j < tlen; j++ ) {
                                   if ( tab[j].c ) {
                                           _NEWMP(m); m->dl = tab[j].d;
                                           C(m) = tab[j].c; NEXT(m) = psum[j];
                                           psum[j] = m;
                                 }                                  }
                                 break;  
                         }                          }
                 }                  }
                 if ( !h_reducible ) {                  for ( j = tlen-1, r = 0; j >= 0; j-- )
                         /* head term is not reducible */                          if ( psum[j] ) {
                         if ( !full ) {                                  _MKDP(n,psum[j],t); _addd_destructive(vl,r,t,&t1); r = t1;
                                 if ( g )                          }
                                         g->sugar = sugar;                  if ( r )
                                 _dptodp(g,rp); _free_dp(g); return;                          r->sugar = p->sugar + m0->dl->td;
                         } else {                  *pr = r;
                                 m = BDY(g);          }
                                 if ( NEXT(m) ) {  }
                                         BDY(g) = NEXT(m); NEXT(m) = 0;  
                                 } else {  /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
                                         FREEDP(g); g = 0;  
   void _weyl_mulmm_dup(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
   {
           P c;
           DL d,d0,d1,dt;
           int i,j,a,b,k,l,n2,s,min,curlen;
           struct cdl *p;
           static Q *ctab;
           static struct cdl *tab;
           static int tablen;
           static struct cdl *tmptab;
           static int tmptablen;
   
           if ( !m0 || !m1 ) {
                   rtab[0].c = 0;
                   rtab[0].d = 0;
                   return;
           }
           mulp(vl,C(m0),C(m1),&c);
           d0 = m0->dl; d1 = m1->dl;
           n2 = n>>1;
           curlen = 1;
   
           _NEWDL(d,n);
           if ( n & 1 )
                   /* offset of h-degree */
                   d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
           else
                   d->td = 0;
           rtab[0].c = c;
           rtab[0].d = d;
   
           if ( rtablen > tmptablen ) {
                   if ( tmptab ) GC_free(tmptab);
                   tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
                   tmptablen = rtablen;
           }
   
           for ( i = 0; i < n2; i++ ) {
                   a = d0->d[i]; b = d1->d[n2+i];
                   k = d0->d[n2+i]; l = d1->d[i];
   
                   /* degree of 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, p = rtab; j < curlen; j++, p++ ) {
                                   if ( p->c ) {
                                           dt = p->d;
                                           dt->d[i] = a;
                                           dt->d[n2+i] = b;
                                           dt->td += s;
                                 }                                  }
                                 if ( d ) {  
                                         for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );  
                                         NEXT(mrd) = m;  
                                 } else {  
                                         _MKDP(n,m,d);  
                                 }  
                         }                          }
                           curlen *= k+1;
                           continue;
                 }                  }
                   if ( k+1 > tablen ) {
                           if ( tab ) GC_free(tab);
                           if ( ctab ) GC_free(ctab);
                           tablen = k+1;
                           tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
                           ctab = (Q *)MALLOC(tablen*sizeof(P));
                   }
                   /* compute xi^a*(Di^k*xi^l)*Di^b */
                   min = MIN(k,l);
                   mkwc(k,l,ctab);
                   bzero(tab,(k+1)*sizeof(struct cdl));
                   /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
                   if ( n & 1 )
                           for ( j = 0; j <= min; j++ ) {
                                   _NEWDL(d,n);
                                   d->d[i] = a-j; d->d[n2+i] = b-j;
                                   d->td = s;
                                   d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
                                   tab[j].d = d;
                                   tab[j].c = (P)ctab[j];
                           }
                   else
                           for ( j = 0; j <= min; j++ ) {
                                   _NEWDL(d,n);
                                   d->d[i] = a-j; d->d[n2+i] = b-j;
                                   d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
                                   tab[j].d = d;
                                   tab[j].c = (P)ctab[j];
                           }
   #if 0
                   _comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
                   for ( j = 0; j < curlen; j++ )
                           if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
                   for ( j = 0; j <= min; j++ )
                           if ( tab[j].d ) { _FREEDL(tab[j].d); }
                   curlen *= k+1;
                   bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
   #else
                   _comm_muld_tab_destructive(vl,n,rtab,curlen,tab,k+1);
                   for ( j = 0; j <= min; j++ )
                           if ( tab[j].d ) { _FREEDL(tab[j].d); }
                   curlen *= k+1;
   #endif
         }          }
         if ( d )  }
                 d->sugar = sugar;  
         _dptodp(d,rp); _free_dp(d);  /* direct product of two cdl tables
     rt[] = [
       t[0]*t1[0],...,t[n-1]*t1[0],
       t[0]*t1[1],...,t[n-1]*t1[1],
       ...
       t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
     ]
   */
   
   void _comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
   {
           int i,j;
           struct cdl *p;
           P c;
           DL d;
   
           bzero(rt,n*n1*sizeof(struct cdl));
           for ( j = 0, p = rt; j < n1; j++ ) {
                   c = t1[j].c;
                   d = t1[j].d;
                   if ( !c )
                           break;
                   for ( i = 0; i < n; i++, p++ ) {
                           if ( t[i].c ) {
                                   mulp(vl,t[i].c,c,&p->c);
                                   _adddl_dup(nv,t[i].d,d,&p->d);
                           }
                   }
           }
   }
   
   void _comm_muld_tab_destructive(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1)
   {
           int i,j;
           struct cdl *p;
           P c;
           DL d;
   
           for ( j = 1, p = t+n; j < n1; j++ ) {
                   c = t1[j].c;
                   d = t1[j].d;
                   if ( !c )
                           break;
                   for ( i = 0; i < n; i++, p++ ) {
                           if ( t[i].c ) {
                                   mulp(vl,t[i].c,c,&p->c);
                                   _adddl_dup(nv,t[i].d,d,&p->d);
                           }
                   }
           }
           c = t1[0].c;
           d = t1[0].d;
           for ( i = 0, p = t; i < n; i++, p++ )
                   if ( t[i].c ) {
                           mulp(vl,t[i].c,c,&p->c);
                           /* t[i].d += d */
                           adddl_destructive(nv,t[i].d,d);
                   }
 }  }
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.14
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>