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

Diff for /OpenXM_contrib2/asir2000/engine/dist.c between version 1.33 and 1.36

version 1.33, 2005/11/16 23:42:53 version 1.36, 2005/11/25 07:18:32
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/dist.c,v 1.32 2004/06/15 16:14:50 ohara Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.35 2005/11/25 02:43:39 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
   
Line 1784  int ni_next(int *a,int n)
Line 1784  int ni_next(int *a,int n)
         int i,j,k,kj;          int i,j,k,kj;
   
         /* find the first nonzero a[j] */          /* find the first nonzero a[j] */
         for ( j = 0; a[j] == 0; j++ );          for ( j = 0; j < n && a[j] == 0; j++ );
         /* find the first zero a[k] after a[j] */          /* find the first zero a[k] after a[j] */
         for ( k = j; k < n && a[k] == 1; k++ );          for ( k = j; k < n && a[k] == 1; k++ );
         if ( k == n ) return 0;          if ( k == n ) return 0;
Line 1893  NBP shuffle_mul_nbm(NBM a,NBM b)
Line 1893  NBP shuffle_mul_nbm(NBM a,NBM b)
         return u;          return u;
 }  }
   
   int nbmtoxky(NBM a,int *b)
   {
           int d,i,j,k;
           int *p;
   
           d = a->d; p = a->b;
           for ( i = j = 0, k = 1; i < d; i++ ) {
                   if ( !NBM_GET(p,i) ) {
                           b[j++] = k;
                           k = 1;
                   } else k++;
           }
           return j;
   }
   
   NBP harmonic_mul_nbm(NBM a,NBM b)
   {
           int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
           int i,j,k,ia,ib,s;
           int *wa,*wb,*w,*wab,*wa1,*wmb;
           Q c,c1;
           NBM wm,tm;
           NODE r,t1,t,p;
           NBP u;
   
           da = a->d; db = b->d; d = da+db;
           wa = (int *)ALLOCA(da*sizeof(int));
           wb = (int *)ALLOCA(db*sizeof(int));
           la = nbmtoxky(a,wa);
           lb = nbmtoxky(b,wb);
           mulq(a->c,b->c,&c);
           /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
           /* lmax : total length */
           lmax = la+lb;
           lmin = la>lb?la:lb;
           w = (int *)ALLOCA(lmax*sizeof(int));
           /* position of a+b */
           wab = (int *)ALLOCA(lmax*sizeof(int));
           /* position of a */
           wa1 = (int *)ALLOCA(lmax*sizeof(int));
           NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
           for ( l = lmin, r = 0; l <= lmax; l++ ) {
                   lab = lmax - l;
                   la1 = la - lab;
                   lb1 = lb - lab;
                   lab1 = l-lab;
                   /* partion l into three parts: a, b, a+b */
                   /* initialize wab */
                   for ( i = 0; i < lab; i++ ) wab[i] = 1;
                   for ( ; i < l; i++ ) wab[i] = 0;
                   do {
                           /* initialize wa1 */
                           for ( i = 0; i < la1; i++ ) wa1[i] = 1;
                           for ( ; i < lab1; i++ ) wa1[i] = 0;
                           do {
                                   ia = 0; ib = 0;
                                   for ( i = j = 0; i < l; i++ )
                                           if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
                                           else if ( wa1[j++] ) w[i] = wa[ia++];
                                           else w[i] = wb[ib++];
                                   for ( i = j = 0; i < l; i++ ) {
                                           for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
                                           NBM_CLR(wmb,j); j++;
                                   }
                                   wm->d = j; wm->c = c;
                                   for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
                                           tm = (NBM)BDY(t);
                                           s = comp_nbm(tm,wm);
                                           if ( s < 0 ) {
                                                   /* insert */
                                                   MKNODE(t1,wm,t);
                                                   if ( !p ) r = t1;
                                                   else NEXT(p) = t1;
                                                   NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                                                   break;
                                           } else if ( s == 0 ) {
                                                   /* add coefs */
                                                   addq(tm->c,c,&c1);
                                                   if ( c1 ) tm->c = c1;
                                                   else NEXT(p) = NEXT(t);
                                                   break;
                                           }
                                   }
                                   if ( !t ) {
                                           /* append */
                                           MKNODE(t1,wm,t);
                                           if ( !p ) r = t1;
                                           else NEXT(p) = t1;
                                           NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
                                   }
                           } while ( ni_next(wa1,lab1) );
                   } while ( ni_next(wab,l) );
           }
           MKNBP(u,r);
           return u;
   }
   
 void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)  void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
 {  {
         NODE b1,b2,br,br0;          NODE b1,b2,br,br0;
Line 1922  void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
Line 2019  void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
                                         NEXTNODE(br0,br); BDY(br) = BDY(b2);                                          NEXTNODE(br0,br); BDY(br) = BDY(b2);
                                         b2 = NEXT(b2); break;                                          b2 = NEXT(b2); break;
                         }                          }
                         if ( !br0 )                  }
                                 if ( b1 )                  if ( !br0 )
                                         br0 = b1;                          if ( b1 )
                                 else if ( b2 )                                  br0 = b1;
                                         br0 = b2;  
                                 else {  
                                         *rp = 0;  
                                         return;  
                                 }  
                         else if ( b1 )  
                                 NEXT(br) = b1;  
                         else if ( b2 )                          else if ( b2 )
                                   br0 = b2;
                           else {
                                   *rp = 0;
                                   return;
                           }
                   else if ( b1 )
                           NEXT(br) = b1;
                   else if ( b2 )
                                 NEXT(br) = b2;                                  NEXT(br) = b2;
                         else                  else
                                 NEXT(br) = 0;                          NEXT(br) = 0;
                         MKNBP(*rp,br0);                  MKNBP(*rp,br0);
                 }  
         }          }
 }  }
   
Line 1959  void chsgnnbp(NBP p,NBP *rp)
Line 2056  void chsgnnbp(NBP p,NBP *rp)
                 NEXTNODE(r0,r);                  NEXTNODE(r0,r);
                 m = (NBM)BDY(b);                  m = (NBM)BDY(b);
                 NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnq(m->c,&m1->c);                  NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnq(m->c,&m1->c);
                 BDY(r) = m;                  BDY(r) = m1;
         }          }
         if ( r0 ) NEXT(r) = 0;          if ( r0 ) NEXT(r) = 0;
         MKNBP(*rp,r0);          MKNBP(*rp,r0);
Line 1970  void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
Line 2067  void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
   
 void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)  void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
 {  {
         NODE b;          NODE b,n;
         NBP r,t,s;          NBP r,t,s;
           NBM m;
   
         if ( !p1 || !p2 ) *rp = 0;          if ( !p1 || !p2 ) {
         else if ( length(BDY(p1)) < length(BDY(p2)) ) {                  *rp = 0; return;
           }
           if ( OID(p1) != O_NBP ) {
                   if ( !NUM(p1) || !RATN(p1) ) error("mulnbp : invalid argument");
                   NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p1;
                   MKNODE(n,m,0); MKNBP(p1,n);
           }
           if ( OID(p2) != O_NBP ) {
                   if ( !NUM(p2) || !RATN(p2) ) error("mulnbp : invalid argument");
                   NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p2;
                   MKNODE(n,m,0); MKNBP(p2,n);
           }
           if ( length(BDY(p1)) < length(BDY(p2)) ) {
                 for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {                  for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
                         mulnbmnbp(vl,(NBM)BDY(b),p2,&t);                          mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
                         addnbp(vl,r,t,&s); r = s;                          addnbp(vl,r,t,&s); r = s;
Line 2046  void pwrnbp(VL vl,NBP a,Q q,NBP *c)
Line 2156  void pwrnbp(VL vl,NBP a,Q q,NBP *c)
         }          }
 }  }
   
 void shuffle_mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);  
 void shuffle_mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);  
 void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);  
   
 void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)  void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
 {  {
         NODE b;          NODE b1,b2,n;
         NBP r,t,s;          NBP r,t,s;
           NBM m;
   
         if ( !p1 || !p2 ) *rp = 0;          if ( !p1 || !p2 ) {
         else if ( length(BDY(p1)) < length(BDY(p2)) ) {                  *rp = 0; return;
                 for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {  
                         shuffle_mulnbmnbp(vl,(NBM)BDY(b),p2,&t);  
                         addnbp(vl,r,t,&s); r = s;  
                 }  
                 *rp = r;  
         } else {  
                 for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {  
                         shuffle_mulnbpnbm(vl,p1,(NBM)BDY(b),&t);  
                         addnbp(vl,r,t,&s); r = s;  
                 }  
                 *rp = r;  
         }          }
 }          if ( OID(p1) != O_NBP ) {
                   if ( !NUM(p1) || !RATN(p1) ) error("mulnbp : invalid argument");
 void shuffle_mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)                  NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p1;
 {                  MKNODE(n,m,0); MKNBP(p1,n);
         NODE b;          }
         NBP t,s,r;          if ( OID(p2) != O_NBP ) {
                   if ( !NUM(p2) || !RATN(p2) ) error("mulnbp : invalid argument");
         if ( !p ) *rp = 0;                  NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p2;
         else {                  MKNODE(n,m,0); MKNBP(p2,n);
                 r = 0;          }
                 for ( b = BDY(p); b; b = NEXT(b) ) {          for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
                         t = shuffle_mul_nbm(m,(NBM)BDY(b));                  for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
                           t = shuffle_mul_nbm(m,(NBM)BDY(b2));
                         addnbp(vl,r,t,&s); r = s;                          addnbp(vl,r,t,&s); r = s;
                 }                  }
                 *rp = r;          *rp = r;
         }  
 }  }
   
 void shuffle_mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)  void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
 {  {
         NODE b;          NODE b1,b2,n;
         NBP t,s,r;          NBP r,t,s;
           NBM m;
   
         if ( !p ) *rp = 0;          if ( !p1 || !p2 ) {
         else {                  *rp = 0; return;
                 r = 0;          }
                 for ( b = BDY(p); b; b = NEXT(b) ) {          if ( OID(p1) != O_NBP ) {
                         t = shuffle_mul_nbm((NBM)BDY(b),m);                  if ( !NUM(p1) || !RATN(p1) ) error("mulnbp : invalid argument");
                   NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p1;
                   MKNODE(n,m,0); MKNBP(p1,n);
           }
           if ( OID(p2) != O_NBP ) {
                   if ( !NUM(p2) || !RATN(p2) ) error("mulnbp : invalid argument");
                   NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p2;
                   MKNODE(n,m,0); MKNBP(p2,n);
           }
           for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
                   for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
                           t = harmonic_mul_nbm(m,(NBM)BDY(b2));
                         addnbp(vl,r,t,&s); r = s;                          addnbp(vl,r,t,&s); r = s;
                 }                  }
                 *rp = r;          *rp = r;
         }  
 }  }
Legend:
Removed from v.1.33  
changed lines
  Added in v.1.36
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>