OpenXM_contrib2/asir2000/engine/dist.c - diff

Return to dist.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/dist.c between version 1.21 and 1.26

-version 1.21, 2002/01/30 01:09:07
+version 1.26, 2003/08/22 08:14:45
 Line 45
 Line 45
 Line 45
   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
   *
-  * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.20 2002/01/28 00:54:43 noro Exp $
+  * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.25 2003/07/09 07:11:08 noro Exp $
  */
  #include "ca.h"
 Line 75  int dp_nelim,dp_fcoeffs;
 Line 75  int dp_nelim,dp_fcoeffs;
 Line 75  int dp_nelim,dp_fcoeffs;
  struct order_spec dp_current_spec;
  int *dp_dl_work;
- int has_fcoef(DP f)
+ void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr);
+ void comm_quod(VL vl,DP p1,DP p2,DP *pr);
+ void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr);
+ void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr);
+ int has_sfcoef(DP f)
  {
          MP t;
          if ( !f )
                  return 0;
          for ( t = BDY(f); t; t = NEXT(t) )
-                 if ( has_fcoef_p(t->c) )
+                 if ( has_sfcoef_p(t->c) )
                          break;
          return t ? 1 : 0;
  }
- int has_fcoef_p(P f)
+ int has_sfcoef_p(P f)
  {
          DCP dc;
          if ( !f )
                  return 0;
          else if ( NUM(f) )
-                 return (NID((Num)f) == N_LM
+                 return (NID((Num)f) == N_GFS) ? 1 : 0;
-                         || NID((Num)f) == N_GF2N
-                         || NID((Num)f) == N_GFPN
-                         || NID((Num)f) == N_GFS) ? 1 : 0;
          else {
                  for ( dc = DC(f); dc; dc = NEXT(dc) )
-                         if ( has_fcoef_p(COEF(dc)) )
+                         if ( has_sfcoef_p(COEF(dc)) )
                                  return 1;
                  return 0;
          }
-Line 204  void ptod(VL vl,VL dvl,P p,DP *pr)
+Line 206  void ptod(VL vl,VL dvl,P p,DP *pr)
 Line 204  void ptod(VL vl,VL dvl,P p,DP *pr)
 Line 206  void ptod(VL vl,VL dvl,P p,DP *pr)
                  }
          }
  #if 0
-         if ( !dp_fcoeffs && has_fcoef(*pr) )
+         if ( !dp_fcoeffs && has_sfcoef(*pr) )
-                 dp_fcoeffs = 1;
+                 dp_fcoeffs = N_GFS;
  #endif
  }
-Line 383  void symb_addd(DP p1,DP p2,DP *pr)
+Line 385  void symb_addd(DP p1,DP p2,DP *pr)
 Line 383  void symb_addd(DP p1,DP p2,DP *pr)
 Line 385  void symb_addd(DP p1,DP p2,DP *pr)
  NODE symb_merge(NODE m1,NODE m2,int n)
  {
          NODE top,prev,cur,m,t;
+         int c,i;
+         DL d1,d2;
          if ( !m1 )
                  return m2;
-Line 403  NODE symb_merge(NODE m1,NODE m2,int n)
+Line 407  NODE symb_merge(NODE m1,NODE m2,int n)
 Line 403  NODE symb_merge(NODE m1,NODE m2,int n)
 Line 407  NODE symb_merge(NODE m1,NODE m2,int n)
                  prev = top; cur = NEXT(top);
                  /* BDY(prev) > BDY(m) always holds */
                  while ( cur && m ) {
+                         d1 = (DL)BDY(cur);
+                         d2 = (DL)BDY(m);
+ #if 1
                          switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
+ #else
+                         /* XXX only valid for DRL */
+                         if ( d1->td > d2->td )
+                                 c = 1;
+                         else if ( d1->td < d2->td )
+                                 c = -1;
+                         else {
+                                 for ( i = n-1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
+                                 if ( i < 0 )
+                                         c = 0;
+                                 else if ( d1->d[i] < d2->d[i] )
+                                         c = 1;
+                                 else
+                                         c = -1;
+                         }
+                         switch ( c ) {
+ #endif
                                  case 0:
                                          m = NEXT(m);
                                          prev = cur; cur = NEXT(cur);
-Line 424  NODE symb_merge(NODE m1,NODE m2,int n)
+Line 448  NODE symb_merge(NODE m1,NODE m2,int n)
 Line 424  NODE symb_merge(NODE m1,NODE m2,int n)
 Line 448  NODE symb_merge(NODE m1,NODE m2,int n)
          }
  }
+ void _adddl(int n,DL d1,DL d2,DL d3)
+ {
+         int i;
+         d3->td = d1->td+d2->td;
+         for ( i = 0; i < n; i++ )
+                 d3->d[i] = d1->d[i]+d2->d[i];
+ }
+ /* m1 <- m1 U dl*f, destructive */
+ NODE mul_dllist(DL dl,DP f);
+ NODE symb_mul_merge(NODE m1,DL dl,DP f,int n)
+ {
+         NODE top,prev,cur,n1;
+         DP g;
+         DL t,s;
+         MP m;
+         if ( !m1 )
+                 return mul_dllist(dl,f);
+         else if ( !f )
+                 return m1;
+         else {
+                 m = BDY(f);
+                 NEWDL_NOINIT(t,n);
+                 _adddl(n,m->dl,dl,t);
+                 top = m1; prev = 0; cur = m1;
+                 while ( m ) {
+                         switch ( (*cmpdl)(n,(DL)BDY(cur),t) ) {
+                                 case 0:
+                                         prev = cur; cur = NEXT(cur);
+                                         if ( !cur ) {
+                                                 MKDP(n,m,g);
+                                                 NEXT(prev) = mul_dllist(dl,g);
+                                                 return;
+                                         }
+                                         m = NEXT(m);
+                                         if ( m ) _adddl(n,m->dl,dl,t);
+                                         break;
+                                 case 1:
+                                         prev = cur; cur = NEXT(cur);
+                                         if ( !cur ) {
+                                                 MKDP(n,m,g);
+                                                 NEXT(prev) = mul_dllist(dl,g);
+                                                 return;
+                                         }
+                                         break;
+                                 case -1:
+                                         NEWDL_NOINIT(s,n);
+                                         s->td = t->td;
+                                         bcopy(t->d,s->d,n*sizeof(int));
+                                         NEWNODE(n1);
+                                         n1->body = (pointer)s;
+                                         NEXT(n1) = cur;
+                                         if ( !prev ) {
+                                                 top = n1; cur = n1;
+                                         } else {
+                                                 NEXT(prev) = n1; prev = n1;
+                                         }
+                                         m = NEXT(m);
+                                         if ( m ) _adddl(n,m->dl,dl,t);
+                                         break;
+                         }
+                 }
+                 return top;
+         }
+ }
  DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
  {
          DLBUCKET top,prev,cur,m,t;
-Line 534  void comm_muld(VL vl,DP p1,DP p2,DP *pr)
+Line 628  void comm_muld(VL vl,DP p1,DP p2,DP *pr)
 Line 534  void comm_muld(VL vl,DP p1,DP p2,DP *pr)
 Line 628  void comm_muld(VL vl,DP p1,DP p2,DP *pr)
          }
  }
+ /* discard terms which is not a multiple of dl */
+ void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr)
+ {
+         MP m;
+         DP s,t,u;
+         int i,l,l1;
+         static MP *w;
+         static int wlen;
+         if ( !p1 || !p2 )
+                 *pr = 0;
+         else if ( OID(p1) <= O_P )
+                 muldc_trunc(vl,p2,(P)p1,dl,pr);
+         else if ( OID(p2) <= O_P )
+                 muldc_trunc(vl,p1,(P)p2,dl,pr);
+         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_trunc(vl,p1,w[i],dl,&t); addd(vl,s,t,&u); s = u;
+                 }
+                 bzero(w,l*sizeof(MP));
+                 *pr = s;
+         }
+ }
+ void comm_quod(VL vl,DP p1,DP p2,DP *pr)
+ {
+         MP m,m0;
+         DP s,t;
+         int i,n,sugar;
+         DL d1,d2,d;
+         Q a,b;
+         if ( !p2 )
+                 error("comm_quod : invalid input");
+         if ( !p1 )
+                 *pr = 0;
+         else {
+                 n = NV(p1);
+                 d2 = BDY(p2)->dl;
+                 m0 = 0;
+                 sugar = p1->sugar;
+                 while ( p1 ) {
+                         d1 = BDY(p1)->dl;
+                         NEWDL(d,n);
+                         d->td = d1->td - d2->td;
+                         for ( i = 0; i < n; i++ )
+                                 d->d[i] = d1->d[i]-d2->d[i];
+                         NEXTMP(m0,m);
+                         m->dl = d;
+                         divq((Q)BDY(p1)->c,(Q)BDY(p2)->c,&a); chsgnq(a,&b);
+                         C(m) = (P)b;
+                         muldm_trunc(vl,p2,m,d2,&t);
+                         addd(vl,p1,t,&s); p1 = s;
+                         C(m) = (P)a;
+                 }
+                 if ( m0 ) {
+                         NEXT(m) = 0; MKDP(n,m0,*pr);
+                 } else
+                         *pr = 0;
+                 /* XXX */
+                 if ( *pr )
+                         (*pr)->sugar = sugar - d2->td;
+         }
+ }
  void muldm(VL vl,DP p,MP m0,DP *pr)
  {
          MP m,mr,mr0;
-Line 559  void muldm(VL vl,DP p,MP m0,DP *pr)
+Line 732  void muldm(VL vl,DP p,MP m0,DP *pr)
 Line 559  void muldm(VL vl,DP p,MP m0,DP *pr)
 Line 732  void muldm(VL vl,DP p,MP m0,DP *pr)
          }
  }
+ void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr)
+ {
+         MP m,mr,mr0;
+         P c;
+         DL d,tdl;
+         int n,i;
+         if ( !p )
+                 *pr = 0;
+         else {
+                 n = NV(p);
+                 NEWDL(tdl,n);
+                 for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl;
+                         m; m = NEXT(m) ) {
+                         _adddl(n,m->dl,d,tdl);
+                         for ( i = 0; i < n; i++ )
+                                 if ( tdl->d[i] < dl->d[i] )
+                                         break;
+                         if ( i < n )
+                                 continue;
+                         NEXTMP(mr0,mr);
+                         mr->dl = tdl;
+                         NEWDL(tdl,n);
+                         if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
+                                 mulq((Q)C(m),(Q)c,(Q *)&C(mr));
+                         else
+                                 mulp(vl,C(m),c,&C(mr));
+                 }
+                 if ( mr0 ) {
+                         NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
+                 } else
+                         *pr = 0;
+                 if ( *pr )
+                         (*pr)->sugar = p->sugar + m0->dl->td;
+         }
+ }
  void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
  {
          MP m;
-Line 799  void muldc(VL vl,DP p,P c,DP *pr)
+Line 1009  void muldc(VL vl,DP p,P c,DP *pr)
 Line 799  void muldc(VL vl,DP p,P c,DP *pr)
 Line 1009  void muldc(VL vl,DP p,P c,DP *pr)
          }
  }
+ void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr)
+ {
+         MP m,mr,mr0;
+         DL mdl;
+         int i,n;
+         if ( !p || !c ) {
+                 *pr = 0; return;
+         }
+         n = NV(p);
+         for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
+                 mdl = m->dl;
+                 for ( i = 0; i < n; i++ )
+                         if ( mdl->d[i] < dl->d[i] )
+                                 break;
+                 if ( i < n )
+                         break;
+                 NEXTMP(mr0,mr);
+                 if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
+                         mulq((Q)C(m),(Q)c,(Q *)&C(mr));
+                 else
+                         mulp(vl,C(m),c,&C(mr));
+                 mr->dl = m->dl;
+         }
+         NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
+         if ( *pr )
+                 (*pr)->sugar = p->sugar;
+ }
  void divsdc(VL vl,DP p,P c,DP *pr)
  {
          MP m,mr,mr0;
-Line 897  int cmpdl_gradlex(int n,DL d1,DL d2)
+Line 1136  int cmpdl_gradlex(int n,DL d1,DL d2)
 Line 897  int cmpdl_gradlex(int n,DL d1,DL d2)
 Line 1136  int cmpdl_gradlex(int n,DL d1,DL d2)
  int cmpdl_revgradlex(int n,DL d1,DL d2)
  {
-         register int i;
+         register int i,c;
          register int *p1,*p2;
          if ( d1->td > d2->td )
-Line 905  int cmpdl_revgradlex(int n,DL d1,DL d2)
+Line 1144  int cmpdl_revgradlex(int n,DL d1,DL d2)
 Line 905  int cmpdl_revgradlex(int n,DL d1,DL d2)
 Line 1144  int cmpdl_revgradlex(int n,DL d1,DL d2)
          else if ( d1->td < d2->td )
                  return -1;
          else {
-                 for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
+                 i = n-1;
-                         i >= 0 && *p1 == *p2; i--, p1--, p2-- );
+                 p1 = d1->d+n-1;
-                 return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
+                 p2 = d2->d+n-1;
+                 while ( i >= 7 ) {
+                         c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                         c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                         c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                         c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                         c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                         c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                         c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                         c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                         i -= 8;
+                 }
+                 switch ( i ) {
+                         case 6:
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 return 0;
+                         case 5:
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 return 0;
+                         case 4:
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 return 0;
+                         case 3:
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 return 0;
+                         case 2:
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 return 0;
+                         case 1:
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 return 0;
+                         case 0:
+                                 c = (*p1--) - (*p2--); if ( c ) goto LAST;
+                                 return 0;
+                         default:
+                                 return 0;
+                 }
+ LAST:
+                 if ( c > 0 ) return -1;
+                 else return 1;
          }
  }
-Line 1238  int cmpdl_matrix(int n,DL d1,DL d2)
+Line 1537  int cmpdl_matrix(int n,DL d1,DL d2)
 Line 1238  int cmpdl_matrix(int n,DL d1,DL d2)
 Line 1537  int cmpdl_matrix(int n,DL d1,DL d2)
                          return -1;
          }
          return 0;
+ }
+ GeoBucket create_bucket()
+ {
+         GeoBucket g;
+         g = CALLOC(1,sizeof(struct oGeoBucket));
+         g->m = 32;
+         return g;
+ }
+ void add_bucket(GeoBucket g,NODE d,int nv)
+ {
+         int l,k,m;
+         l = length(d);
+         for ( k = 0, m = 1; l > m; k++, m <<= 1 );
+         /* 2^(k-1) < l <= 2^k */
+         d = symb_merge(g->body[k],d,nv);
+         for ( ; length(d) > (1<<(k)); k++ ) {
+                 g->body[k] = 0;
+                 d = symb_merge(g->body[k+1],d,nv);
+         }
+         g->body[k] = d;
+         g->m = MAX(g->m,k);
+ }
+ DL remove_head_bucket(GeoBucket g,int nv)
+ {
+         int j,i,c,m;
+         DL d;
+         j = -1;
+         m = g->m;
+         for ( i = 0; i <= m; i++ ) {
+                 if ( !g->body[i] )
+                         continue;
+                 if ( j < 0 ) j = i;
+                 else {
+                         c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
+                         if ( c > 0 )
+                                 j = i;
+                         else if ( c == 0 )
+                                 g->body[i] = NEXT(g->body[i]);
+                 }
+         }
+         if ( j < 0 )
+                 return 0;
+         else {
+                 d = g->body[j]->body;
+                 g->body[j] = NEXT(g->body[j]);
+                 return d;
+         }
  }

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