[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.26 and 1.33

version 1.26, 2001/09/11 01:30:31 version 1.33, 2001/09/17 08:37:30
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.25 2001/09/10 05:55:14 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.32 2001/09/17 07:16:58 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 dptoca(DP,unsigned int **);
   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 193  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 228  int *b;
         }          }
 }  }
   
 /* create compressed poly */  /* [t,findex] -> tf -> compressed vector */
   
 void _dpmod_to_vect_compress(f,at,b)  void _tf_to_vect_compress(tf,at,b)
 DP f;  NODE tf;
 DL *at;  DL *at;
 CDP *b;  CDP *b;
 {  {
         int i,j,nv,len;          int i,j,k,nv,len;
           DL t,s,d1;
           DP f;
         MP m;          MP m;
         CDP r;          CDP r;
   
           t = (DL)BDY(tf);
           f = ps[(int)BDY(NEXT(tf))];
   
         nv = f->nv;          nv = f->nv;
         for ( m = BDY(f), len = 0; m; m = NEXT(m), len++ );          for ( m = BDY(f), len = 0; m; m = NEXT(m), len++ );
         r = (CDP)MALLOC(sizeof(struct oCDP));          r = (CDP)MALLOC(sizeof(struct oCDP));
         r->len = len;          r->len = len;
         r->body = (CM)MALLOC(sizeof(struct oCM)*len);          r->psindex = (int)BDY(NEXT(tf));
           r->body = (unsigned int *)MALLOC_ATOMIC(sizeof(unsigned int)*len);
   
           NEWDL(s,nv);
         for ( m = BDY(f), i = j = 0; m; m = NEXT(m), j++ ) {          for ( m = BDY(f), i = j = 0; m; m = NEXT(m), j++ ) {
                 for ( ; !eqdl(nv,m->dl,at[i]); i++ );                  d1 = m->dl;
                 r->body[j].index = i;                  s->td = t->td+d1->td;
                 r->body[j].c = ITOS(m->c);                  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] = i;
         }          }
         *b = r;          *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 296  DP f;
Line 294  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 486  LIST f,v;
Line 510  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 498  LIST *rp;
Line 522  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 519  int m;
Line 547  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 532  LIST *rp;
Line 560  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 543  LIST *rp;
Line 573  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 672  NODE f;
Line 704  NODE f;
   
 /* initial bases are monic */  /* initial bases are monic */
   
   unsigned int **psca;
   
 NODE gb_f4_mod(f,m)  NODE gb_f4_mod(f,m)
 NODE f;  NODE f;
 int m;  int m;
Line 682  int m;
Line 716  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;          int **spmat;
         CDP *redmat;          CDP *redmat;
Line 690  int m;
Line 724  int m;
         int rank,nred,nsp,nonzero,spcol;          int rank,nred,nsp,nonzero,spcol;
         int *indred,*isred;          int *indred,*isred;
         CDP ri;          CDP ri;
           int pscalen;
         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;
   
           /* initialize coeffcient array list of ps[] */
           pscalen = pslen;
           psca = (unsigned int **)MALLOC(pscalen*sizeof(unsigned int *));
   
         init_eg(&eg_symb); init_eg(&eg_elim1); init_eg(&eg_elim2);          init_eg(&eg_symb); init_eg(&eg_elim1); init_eg(&eg_elim2);
         for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {          for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
                 i = (int)BDY(r);                  i = (int)BDY(r);
                 d = updpairs(d,g,i);                  d = updpairs(d,g,i);
                 g = updbase(g,i);                  g = updbase(g,i);
                 gall = append_one(gall,i);                  gall = append_one(gall,i);
                   dptoca(ps[i],&psca[i]);
         }          }
         if ( gall )          if ( gall )
                 nv = ((DP)ps[(int)BDY(gall)])->nv;                  nv = ((DP)ps[(int)BDY(gall)])->nv;
Line 717  int m;
Line 757  int m;
                         }                          }
                 }                  }
                 /* s0 : all the terms appeared in symbolic reduction */                  /* s0 : all the terms appeared in symbolic reduction */
 #if 0  
                 for ( s = s0, nred = 0; s; s = NEXT(s) ) {                  for ( s = s0, nred = 0; s; s = NEXT(s) ) {
                         for ( j = psn-1; j >= 0; j-- )  
                                 if ( _dl_redble(BDY(ps[j])->dl,BDY(s),nv) )  
                                         break;  
                         if ( j >= 0 ) {  
                                 dltod(BDY(s),nv,&tdp);  
                                 dp_subd(tdp,ps[j],&sd);  
                                 for ( k = 0, i = 0; k < nv; k++ )  
                                         if ( BDY(sd)->dl->d[k] )  
                                                 i++;  
                                 fprintf(stderr,"%c ",i<=1 ? 'o' : 'x');  
                                 _dp_mod(sd,m,0,&sdm);  
                                 mulmd_dup(m,sdm,ps[j],&f2);  
                                 MKNODE(bt,f2,blist); blist = bt;  
                                 s = symb_merge(s,dp_dllist(f2),nv);  
                                 nred++;  
                         }  
                 }  
 #else  
                 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) )
                                         break;                                          break;
                         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++;
                         }                          }
                 }                  }
 #endif  
                 fprintf(stderr,"\n");  
   
                 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 779  int m;
Line 790  int m;
   
                 /* reducer matrix */                  /* reducer matrix */
                 /* indred : register the position of the head term */                  /* indred : register the position of the head term */
 #if 0  
                 reduce_reducers_mod_compress(blist,nred,at,col,m,&redmat,&indred);  
                 isred = (int *)MALLOC(col*sizeof(int));  
                 bzero(isred,col*sizeof(int));  
                 for ( i = 0; i < nred; i++ )  
                         isred[indred[i]] = 1;  
 #else  
                 redmat = (CDP *)MALLOC(nred*sizeof(CDP));                  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_compress(BDY(r),at,&redmat[i]);                          _tf_to_vect_compress(BDY(r),at,&redmat[i]);
                 /* XXX */  
 /*              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];
                         indred[i] = ri->body[0].index;                          indred[i] = ri->body[0];
                         isred[indred[i]] = 1;                          isred[indred[i]] = 1;
                 }                  }
 #endif  
   
                 spcol = col-nred;                  spcol = col-nred;
                 /* head terms not in ht */                  /* head terms not in ht */
Line 809  int m;
Line 811  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);                  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 = (int **)MALLOC(nsp*sizeof(int *));                  spmat = (int **)MALLOC(nsp*sizeof(int *));
                 w = (int *)MALLOC(col*sizeof(int));                  w = (int *)MALLOC_ATOMIC(col*sizeof(int));
   
                 /* skip reducers in blist */                  /* skip reducers in blist */
                 for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ );                  for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ );
Line 836  int m;
Line 840  int m;
                 /* update nsp */                  /* update nsp */
                 nsp = i;                  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);
   
Line 878  int m;
Line 888  int m;
                                 NEXT(mp) = 0;                                  NEXT(mp) = 0;
                                 MKDP(nv,mp0,nf); nf->sugar = dm->sugar;                                  MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
                                 nh = newps_mod(nf,m);                                  nh = newps_mod(nf,m);
                                   if ( nh == pscalen ) {
                                           psca = (unsigned int **)
                                                   REALLOC(psca,2*pscalen*sizeof(unsigned int *));
                                           pscalen *= 2;
                                   }
                                   dptoca(ps[nh],&psca[nh]);
                                 d = updpairs(d,g,nh);                                  d = updpairs(d,g,nh);
                                 g = updbase(g,nh);                                  g = updbase(g,nh);
                                 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);
Line 2458  VECT *rp;
Line 2479  VECT *rp;
         MKVECT(r,n); *rp = r;          MKVECT(r,n); *rp = r;
         for ( i = 0; i < n; i++ )          for ( i = 0; i < n; i++ )
                 mulq((Q)BDY(w)[i],(Q)c,(Q *)&BDY(r)[i]);                  mulq((Q)BDY(w)[i],(Q)c,(Q *)&BDY(r)[i]);
   }
   
   void dptoca(p,rp)
   DP p;
   unsigned int **rp;
   {
           int i;
           MP m;
           unsigned int *r;
   
           if ( !p )
                   *rp = 0;
           else {
                   for ( m = BDY(p), i = 0; m; m = NEXT(m), i++ );
                   *rp = r = (unsigned int *)MALLOC_ATOMIC(i*sizeof(unsigned int));
                   for ( m = BDY(p), i = 0; m; m = NEXT(m), i++ )
                           r[i] = ITOS(C(m));
           }
 }  }
   

Legend:
Removed from v.1.26  
changed lines
  Added in v.1.33

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