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

Diff for /OpenXM_contrib2/asir2000/builtin/gr.c between version 1.23 and 1.31

version 1.23, 2001/09/05 08:09:10 version 1.31, 2001/09/17 02:47:07
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/builtin/gr.c,v 1.22 2001/09/05 01:57:32 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.30 2001/09/17 01:18:34 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "parse.h"  #include "parse.h"
 #include "base.h"  #include "base.h"
 #include "ox.h"  #include "ox.h"
   
   #if defined(__GNUC__)
   #define INLINE inline
   #elif defined(VISUAL)
   #define INLINE __inline
   #else
   #define INLINE
   #endif
   
 #define ITOS(p) (((unsigned int)(p))&0x7fffffff)  #define ITOS(p) (((unsigned int)(p))&0x7fffffff)
 #define STOI(i) ((P)((unsigned int)(i)|0x80000000))  #define STOI(i) ((P)((unsigned int)(i)|0x80000000))
   
Line 85  extern int do_weyl;
Line 93  extern int do_weyl;
   
 extern DP_Print;  extern DP_Print;
   
   void _tf_to_vect_compress(NODE,DL *,CDP *);
   NODE mul_dllist(DL,DP);
 void dp_imul_d(DP,Q,DP *);  void dp_imul_d(DP,Q,DP *);
 void print_stat(void);  void print_stat(void);
 void init_stat(void);  void init_stat(void);
Line 182  int doing_f4;
Line 192  int doing_f4;
 NODE TraceList;  NODE TraceList;
 NODE AllTraceList;  NODE AllTraceList;
   
 int eqdl(nv,dl1,dl2)  INLINE int eqdl(nv,dl1,dl2)
 int nv;  int nv;
 DL dl1,dl2;  DL dl1,dl2;
 {  {
Line 217  int *b;
Line 227  int *b;
         }          }
 }  }
   
   /* create compressed poly */
   
   void _dpmod_to_vect_compress(f,at,b)
   DP f;
   DL *at;
   CDP *b;
   {
           int i,j,nv,len;
           MP m;
           CDP r;
   
           nv = f->nv;
           for ( m = BDY(f), len = 0; m; m = NEXT(m), len++ );
           r = (CDP)MALLOC(sizeof(struct oCDP));
           r->len = len;
           r->body = (CM)MALLOC_ATOMIC(sizeof(struct oCM)*len);
   
           for ( m = BDY(f), i = j = 0; m; m = NEXT(m), j++ ) {
                   for ( ; !eqdl(nv,m->dl,at[i]); i++ );
                   r->body[j].index = i;
                   r->body[j].c = ITOS(m->c);
           }
           *b = r;
   }
   
   /* [t,findex] -> tf -> compressed vector */
   
   void _tf_to_vect_compress(tf,at,b)
   NODE tf;
   DL *at;
   CDP *b;
   {
           int i,j,k,nv,len;
           DL t,s,d1;
           DP f;
           MP m;
           CDP r;
   
           t = (DL)BDY(tf);
           f = ps[(int)BDY(NEXT(tf))];
   
           nv = f->nv;
           for ( m = BDY(f), len = 0; m; m = NEXT(m), len++ );
           r = (CDP)MALLOC(sizeof(struct oCDP));
           r->len = len;
           r->body = (CM)MALLOC_ATOMIC(sizeof(struct oCM)*len);
   
           NEWDL(s,nv);
           for ( m = BDY(f), i = j = 0; m; m = NEXT(m), j++ ) {
                   d1 = m->dl;
                   s->td = t->td+d1->td;
                   for ( k = 0; k < nv; k++ )
                           s->d[k] = t->d[k]+d1->d[k];
                   for ( ; !eqdl(nv,s,at[i]); i++ );
                   r->body[j].index = i;
                   r->body[j].c = ITOS(m->c);
           }
           *b = r;
   }
   
   /* dense vector -> CDP  */
   void compress_vect(a,n,rp)
   int *a;
   int n;
   CDP *rp;
   {
           int i,j,nz;
           CDP r;
   
           for ( i = 0, nz = 0; i < n; i++ )
                   if ( a[i] ) nz++;
           *rp = r = (CDP)MALLOC(sizeof(struct oCDP));
           r->len = nz;
           r->body = (CM)MALLOC(sizeof(struct oCM)*nz);
           for ( i = 0, j = 0; i < n; i++ ) {
                   if ( a[i] ) {
                           r->body[j].index = i;
                           r->body[j].c = ITOS(a[i]);
                           j++;
                   }
           }
   }
   
 void dp_to_vect(f,at,b)  void dp_to_vect(f,at,b)
 DP f;  DP f;
 DL *at;  DL *at;
Line 248  DP f;
Line 341  DP f;
         return mp0;          return mp0;
 }  }
   
   NODE mul_dllist(d,f)
   DL d;
   DP f;
   {
           MP m;
           NODE mp,mp0;
           DL t,d1;
           int i,nv;
   
           if ( !f )
                   return 0;
           nv = NV(f);
           mp0 = 0;
           for ( m = BDY(f); m; m = NEXT(m) ) {
                   NEXTNODE(mp0,mp);
                   NEWDL(t,nv);
                   d1 = m->dl;
                   t->td = d->td+d1->td;
                   for ( i = 0; i < nv; i++ )
                           t->d[i] = d->d[i]+d1->d[i];
                   BDY(mp) = (pointer)t;
           }
           NEXT(mp) = 0;
           return mp0;
   }
   
 void pdl(f)  void pdl(f)
 NODE f;  NODE f;
 {  {
Line 438  LIST f,v;
Line 557  LIST f,v;
 struct order_spec *ord;  struct order_spec *ord;
 LIST *rp;  LIST *rp;
 {  {
         int i,mindex,m,nochk;          int i,mindex,m,nochk,homogen;
         struct order_spec ord1;          struct order_spec ord1;
         VL fv,vv,vc;          VL fv,vv,vc;
         NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;          NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;
Line 450  LIST *rp;
Line 569  LIST *rp;
         if ( ord->id && NVars != ord->nv )          if ( ord->id && NVars != ord->nv )
                 error("dp_f4_main : invalid order specification");                  error("dp_f4_main : invalid order specification");
         initd(ord);          initd(ord);
         for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {          for ( fd0 = 0, t = BDY(f), homogen = 1; t; t = NEXT(t) ) {
                 NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));                  NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
                   if ( homogen )
                           homogen = dp_homogeneous(BDY(fd));
         }          }
         if ( fd0 ) NEXT(fd) = 0;          if ( fd0 ) NEXT(fd) = 0;
         setup_arrays(fd0,0,&s);          setup_arrays(fd0,0,&s);
         x = gb_f4(s);          x = gb_f4(s);
         reduceall(x,&xx); x = xx;          if ( !homogen ) {
                   reduceall(x,&xx); x = xx;
           }
         for ( r0 = 0; x; x = NEXT(x) ) {          for ( r0 = 0; x; x = NEXT(x) ) {
                 NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);                  NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);
                 dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));                  dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
Line 471  int m;
Line 594  int m;
 struct order_spec *ord;  struct order_spec *ord;
 LIST *rp;  LIST *rp;
 {  {
         int i;          int i,homogen;
         struct order_spec ord1;          struct order_spec ord1;
         VL fv,vv,vc;          VL fv,vv,vc;
         DP b,c,c1;          DP b,c,c1;
Line 484  LIST *rp;
Line 607  LIST *rp;
         if ( ord->id && NVars != ord->nv )          if ( ord->id && NVars != ord->nv )
                 error("dp_f4_mod_main : invalid order specification");                  error("dp_f4_mod_main : invalid order specification");
         initd(ord);          initd(ord);
         for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {          for ( fd0 = 0, t = BDY(f), homogen = 1; t; t = NEXT(t) ) {
                 ptod(CO,vv,(P)BDY(t),&b);                  ptod(CO,vv,(P)BDY(t),&b);
                   if ( homogen )
                           homogen = dp_homogeneous(b);
                 _dp_mod(b,m,0,&c);                  _dp_mod(b,m,0,&c);
                 _dp_monic(c,m,&c1);                  _dp_monic(c,m,&c1);
                 if ( c ) {                  if ( c ) {
Line 495  LIST *rp;
Line 620  LIST *rp;
         if ( fd0 ) NEXT(fd) = 0;          if ( fd0 ) NEXT(fd) = 0;
         setup_arrays(fd0,m,&s);          setup_arrays(fd0,m,&s);
         x = gb_f4_mod(s,m);          x = gb_f4_mod(s,m);
         reduceall_mod(x,m,&xx); x = xx;          if ( !homogen ) {
                   reduceall_mod(x,m,&xx); x = xx;
           }
         for ( r0 = 0; x; x = NEXT(x) ) {          for ( r0 = 0; x; x = NEXT(x) ) {
                 NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));                  NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
         }          }
Line 634  int m;
Line 761  int m;
         DP_pairs d,dm,dr,t;          DP_pairs d,dm,dr,t;
         DP h,nf,f1,f2,f21,f21r,sp,sp1,sd,sdm,tdp;          DP h,nf,f1,f2,f21,f21r,sp,sp1,sd,sdm,tdp;
         MP mp,mp0;          MP mp,mp0;
         NODE blist,bt,nt;          NODE blist,bt,nt,bt1,dt,rhtlist;
         DL *ht,*at,*st;          DL *ht,*at,*st;
         int **spmat,**redmat;          int **spmat;
         int *colstat,*w;          CDP *redmat;
           int *colstat,*w,*w1;
         int rank,nred,nsp,nonzero,spcol;          int rank,nred,nsp,nonzero,spcol;
         int *indred,*isred,*ri;          int *indred,*isred;
           CDP ri;
         struct oEGT tmp0,tmp1,tmp2,eg_split_symb,eg_split_elim1,eg_split_elim2;          struct oEGT tmp0,tmp1,tmp2,eg_split_symb,eg_split_elim1,eg_split_elim2;
         extern struct oEGT eg_symb,eg_elim1,eg_elim2;          extern struct oEGT eg_symb,eg_elim1,eg_elim2;
   
Line 666  int m;
Line 795  int m;
                                 s0 = symb_merge(s0,dp_dllist(sp),nv);                                  s0 = symb_merge(s0,dp_dllist(sp),nv);
                         }                          }
                 }                  }
                 /* s0 : all the terms appeared in symbolic redunction */                  /* s0 : all the terms appeared in symbolic reduction */
                 for ( s = s0, nred = 0; s; s = NEXT(s) ) {                  for ( s = s0, nred = 0; s; s = NEXT(s) ) {
                         for ( r = gall; r; r = NEXT(r) )                          for ( r = gall; r; r = NEXT(r) )
                                 if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )                                  if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )
Line 674  int m;
Line 803  int m;
                         if ( r ) {                          if ( r ) {
                                 dltod(BDY(s),nv,&tdp);                                  dltod(BDY(s),nv,&tdp);
                                 dp_subd(tdp,ps[(int)BDY(r)],&sd);                                  dp_subd(tdp,ps[(int)BDY(r)],&sd);
                                 _dp_mod(sd,m,0,&sdm);                                  dt = mul_dllist(BDY(sd)->dl,ps[(int)BDY(r)]);
                                 mulmd_dup(m,sdm,ps[(int)BDY(r)],&f2);                                  /* list of [t,f] */
                                 MKNODE(bt,f2,blist); blist = bt;                                  bt1 = mknode(2,BDY(sd)->dl,BDY(r));
                                 s = symb_merge(s,dp_dllist(f2),nv);                                  MKNODE(bt,bt1,blist); blist = bt;
                                   symb_merge(s,dt,nv);
                                 nred++;                                  nred++;
                         }                          }
                 }                  }
   
                 get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);  
                 init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);  
   
                 /* the first nred polys in blist are reducers */                  /* the first nred polys in blist are reducers */
                 /* row = the number of all the polys */                  /* row = the number of all the polys */
                 for ( r = blist, row = 0; r; r = NEXT(r), row++ );                  for ( r = blist, row = 0; r; r = NEXT(r), row++ );
   
                 /* head terms of reducers */  
                 ht = (DL *)MALLOC(nred*sizeof(DL));  
                 for ( r = blist, i = 0; i < nred; r = NEXT(r), i++ )  
                         ht[i] = BDY((DP)BDY(r))->dl;  
   
                 /* col = number of all terms */                  /* col = number of all terms */
                 for ( s = s0, col = 0; s; s = NEXT(s), col++ );                  for ( s = s0, col = 0; s; s = NEXT(s), col++ );
   
Line 706  int m;
Line 828  int m;
                 nsp = row-nred;                  nsp = row-nred;
   
                 /* reducer matrix */                  /* reducer matrix */
                 redmat = (int **)almat(nred,col);                  /* indred : register the position of the head term */
                   redmat = (CDP *)MALLOC(nred*sizeof(CDP));
                 for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ )                  for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ )
                         _dpmod_to_vect(BDY(r),at,redmat[i]);                          _tf_to_vect_compress(BDY(r),at,&redmat[i]);
                 /* XXX */                  /* XXX */
 /*              reduce_reducers_mod(redmat,nred,col,m); */  /*              reduce_reducers_mod(redmat,nred,col,m); */
                 /* register the position of the head term */                  /* register the position of the head term */
                 indred = (int *)MALLOC(nred*sizeof(int));                  indred = (int *)MALLOC_ATOMIC(nred*sizeof(int));
                 bzero(indred,nred*sizeof(int));                  bzero(indred,nred*sizeof(int));
                 isred = (int *)MALLOC(col*sizeof(int));                  isred = (int *)MALLOC_ATOMIC(col*sizeof(int));
                 bzero(isred,col*sizeof(int));                  bzero(isred,col*sizeof(int));
                 for ( i = 0; i < nred; i++ ) {                  for ( i = 0; i < nred; i++ ) {
                         ri = redmat[i];                          ri = redmat[i];
                         for ( j = 0; j < col && !ri[j]; j++ );                          indred[i] = ri->body[0].index;
                         indred[i] = j;                          isred[indred[i]] = 1;
                         isred[j] = 1;  
                 }                  }
   
                 spcol = col-nred;                  spcol = col-nred;
Line 729  int m;
Line 851  int m;
                 for ( j = 0, k = 0; j < col; j++ )                  for ( j = 0, k = 0; j < col; j++ )
                         if ( !isred[j] )                          if ( !isred[j] )
                                 st[k++] = at[j];                                  st[k++] = at[j];
                   get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
                   init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
   
                   get_eg(&tmp1);
                 /* spoly matrix; stored in reduced form; terms in ht[] are omitted */                  /* spoly matrix; stored in reduced form; terms in ht[] are omitted */
                 spmat = almat(nsp,spcol);                  spmat = (int **)MALLOC(nsp*sizeof(int *));
                 w = (int *)MALLOC(col*sizeof(int));                  w = (int *)MALLOC_ATOMIC(col*sizeof(int));
                 for ( ; i < row; r = NEXT(r), i++ ) {  
                   /* skip reducers in blist */
                   for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ );
                   for ( i = 0; r; r = NEXT(r) ) {
                         bzero(w,col*sizeof(int));                          bzero(w,col*sizeof(int));
                         _dpmod_to_vect(BDY(r),at,w);                          _dpmod_to_vect(BDY(r),at,w);
                         reduce_sp_by_red_mod(w,redmat,indred,nred,col,m);                          reduce_sp_by_red_mod_compress(w,redmat,indred,nred,col,m);
                         for ( j = 0, k = 0; j < col; j++ )                          for ( j = 0; j < col; j++ )
                                 if ( !isred[j] )                                  if ( w[j] )
                                         spmat[i-nred][k++] = w[j];                                          break;
                           if ( j < col ) {
                                   w1 = (int *)MALLOC_ATOMIC(spcol*sizeof(int));
                                   for ( j = 0, k = 0; j < col; j++ )
                                           if ( !isred[j] )
                                                   w1[k++] = w[j];
                                   spmat[i] = w1;
                                   i++;
                           }
                 }                  }
                   /* update nsp */
                   nsp = i;
   
                   /* XXX free redmat explicitly */
                   for ( k = 0; k < nred; k++ ) {
                           GC_free(BDY(redmat[k]));
                           GC_free(redmat[k]);
                   }
   
                 get_eg(&tmp0); add_eg(&eg_elim1,&tmp1,&tmp0);                  get_eg(&tmp0); add_eg(&eg_elim1,&tmp1,&tmp0);
                 init_eg(&eg_split_elim1); add_eg(&eg_split_elim1,&tmp1,&tmp0);                  init_eg(&eg_split_elim1); add_eg(&eg_split_elim1,&tmp1,&tmp0);
   
                 colstat = (int *)MALLOC_ATOMIC(spcol*sizeof(int));                  colstat = (int *)MALLOC_ATOMIC(spcol*sizeof(int));
                   bzero(colstat,spcol*sizeof(int));
                 for ( i = 0, nonzero=0; i < nsp; i++ )                  for ( i = 0, nonzero=0; i < nsp; i++ )
                         for ( j = 0; j < spcol; j++ )                          for ( j = 0; j < spcol; j++ )
                                 if ( spmat[i][j] )                                  if ( spmat[i][j] )
Line 767  int m;
Line 912  int m;
                         print_eg("Elim2",&eg_split_elim2);                          print_eg("Elim2",&eg_split_elim2);
                         fprintf(asir_out,"\n");                          fprintf(asir_out,"\n");
                 }                  }
   
                   if ( !rank )
                           continue;
   
                 for ( j = 0, i = 0; j < spcol; j++ )                  for ( j = 0, i = 0; j < spcol; j++ )
                         if ( colstat[j] ) {                          if ( colstat[j] ) {
                                 mp0 = 0;                                  mp0 = 0;
Line 784  int m;
Line 933  int m;
                                 gall = append_one(gall,nh);                                  gall = append_one(gall,nh);
                                 i++;                                  i++;
                         }                          }
   
                   /* XXX free spmat[] explicitly */
                   for ( j = 0; j < nsp; j++ ) {
                           GC_free(spmat[j]);
                   }
         }          }
         if ( DP_Print ) {          if ( DP_Print ) {
                 print_eg("Symb",&eg_symb);                  print_eg("Symb",&eg_symb);
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.31
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>