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

Diff for /OpenXM_contrib2/asir2000/engine/nd.c between version 1.215 and 1.229

version 1.215, 2013/12/20 02:02:24 version 1.229, 2016/11/17 05:33:20
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.214 2013/09/27 07:00:45 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.228 2016/08/08 07:18:10 noro Exp $ */
   
 #include "nd.h"  #include "nd.h"
   
Line 8  int diag_period = 6;
Line 8  int diag_period = 6;
 int weight_check = 1;  int weight_check = 1;
 int (*ndl_compare_function)(UINT *a1,UINT *a2);  int (*ndl_compare_function)(UINT *a1,UINT *a2);
 int nd_dcomp;  int nd_dcomp;
   int nd_rref2;
 NM _nm_free_list;  NM _nm_free_list;
 ND _nd_free_list;  ND _nd_free_list;
 ND_pairs _ndp_free_list;  ND_pairs _ndp_free_list;
 NODE nd_hcf;  NODE nd_hcf;
   
   Obj nd_top_weight;
   
 static NODE nd_subst;  static NODE nd_subst;
 static VL nd_vc;  static VL nd_vc;
 static int nd_ntrans;  static int nd_ntrans;
Line 53  static int nd_found,nd_create,nd_notfirst;
Line 56  static int nd_found,nd_create,nd_notfirst;
 static int nmv_adv;  static int nmv_adv;
 static int nd_demand;  static int nd_demand;
 static int nd_module,nd_ispot,nd_mpos,nd_pot_nelim;  static int nd_module,nd_ispot,nd_mpos,nd_pot_nelim;
   static int nd_module_rank,nd_poly_weight_len;
   static int *nd_poly_weight,*nd_module_weight;
 static NODE nd_tracelist;  static NODE nd_tracelist;
 static NODE nd_alltracelist;  static NODE nd_alltracelist;
 static int nd_gentrace,nd_gensyz,nd_nora,nd_newelim,nd_intersect;  static int nd_gentrace,nd_gensyz,nd_nora,nd_newelim,nd_intersect;
Line 85  NDV ndvgztondv(NDV p);
Line 90  NDV ndvgztondv(NDV p);
 ND ndtondgz(ND p);  ND ndtondgz(ND p);
 ND ndgztond(ND p);  ND ndgztond(ND p);
   
 extern int Denominator,DP_Multiple;  void Pdp_set_weight(NODE,VECT *);
   void Pox_cmo_rpc(NODE,Obj *);
   
   extern int Denominator,DP_Multiple,MaxDeg;
   
   #define BLEN (8*sizeof(unsigned long))
   
   typedef struct matrix {
     int row,col;
     unsigned long **a;
   } *matrix;
   
   
 void nd_free_private_storage()  void nd_free_private_storage()
 {  {
     _nm_free_list = 0;      _nm_free_list = 0;
Line 107  void _NM_alloc()
Line 123  void _NM_alloc()
     }      }
 }  }
   
   matrix alloc_matrix(int row,int col)
   {
     unsigned long **a;
     int i,len,blen;
     matrix mat;
   
     mat = (matrix)MALLOC(sizeof(struct matrix));
     mat->row = row;
     mat->col = col;
     mat->a = a = (unsigned long **)MALLOC(row*sizeof(unsigned long *));
     return mat;
   }
   
   
 void _ND_alloc()  void _ND_alloc()
 {  {
     ND p;      ND p;
Line 497  int ndl_block_compare(UINT *d1,UINT *d2)
Line 527  int ndl_block_compare(UINT *d1,UINT *d2)
   
 int ndl_matrix_compare(UINT *d1,UINT *d2)  int ndl_matrix_compare(UINT *d1,UINT *d2)
 {  {
     int i,j,s;      int i,j,s,row;
     int *v;      int *v;
           Q **mat;
       Q *w;
       Q t,t1,t2;
   
     for ( j = 0; j < nd_nvar; j++ )      for ( j = 0; j < nd_nvar; j++ )
         nd_work_vector[j] = GET_EXP(d1,j)-GET_EXP(d2,j);          nd_work_vector[j] = GET_EXP(d1,j)-GET_EXP(d2,j);
           if ( nd_top_weight ) {
                   if ( OID(nd_top_weight) == O_VECT ) {
                       mat = (Q **)&BDY((VECT)nd_top_weight);
                       row = 1;
                   } else {
                           mat = (Q **)BDY((MAT)nd_top_weight);
                       row = ((MAT)nd_top_weight)->row;
                   }
                   for ( i = 0; i < row; i++ ) {
                       w = (Q *)mat[i];
                           for ( j = 0, t = 0; j < nd_nvar; j++ ) {
                                   STOQ(nd_work_vector[j],t1);
                               mulq(w[j],t1,&t2);
                                   addq(t,t2,&t1);
                                   t = t1;
                           }
                       if ( t ) {
                           s = SGN(t);
                           if ( s > 0 ) return 1;
                           else if ( s < 0 ) return -1;
                       }
                   }
           }
     for ( i = 0; i < nd_matrix_len; i++ ) {      for ( i = 0; i < nd_matrix_len; i++ ) {
         v = nd_matrix[i];          v = nd_matrix[i];
         for ( j = 0, s = 0; j < nd_nvar; j++ )          for ( j = 0, s = 0; j < nd_nvar; j++ )
Line 509  int ndl_matrix_compare(UINT *d1,UINT *d2)
Line 565  int ndl_matrix_compare(UINT *d1,UINT *d2)
         if ( s > 0 ) return 1;          if ( s > 0 ) return 1;
         else if ( s < 0 ) return -1;          else if ( s < 0 ) return -1;
     }      }
           if ( !ndl_equal(d1,d2) )
                   error("afo");
     return 0;      return 0;
 }  }
   
Line 591  int ndl_ww_lex_compare(UINT *d1,UINT *d2)
Line 649  int ndl_ww_lex_compare(UINT *d1,UINT *d2)
     return ndl_lex_compare(d1,d2);      return ndl_lex_compare(d1,d2);
 }  }
   
   int ndl_module_weight_compare(UINT *d1,UINT *d2)
   {
     int s,j;
   
     if ( nd_nvar != nd_poly_weight_len )
       error("invalid module weight : the length of polynomial weight != the number of variables");
     s = 0;
     for ( j = 0; j < nd_nvar; j++ )
        s += (GET_EXP(d1,j)-GET_EXP(d2,j))*nd_poly_weight[j];
     if ( MPOS(d1) >= 1 && MPOS(d2) >= 1 ) {
       s += nd_module_weight[MPOS(d1)-1]-nd_module_weight[MPOS(d2)-1];
     }
     if ( s > 0 ) return 1;
     else if ( s < 0 ) return -1;
     else return 0;
   }
   
 int ndl_module_grlex_compare(UINT *d1,UINT *d2)  int ndl_module_grlex_compare(UINT *d1,UINT *d2)
 {  {
     int i,c;      int i,c;
   
       if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {      if ( nd_ispot ) {
                 if ( nd_pot_nelim && MPOS(d1)>=nd_pot_nelim+1 && MPOS(d2) >= nd_pot_nelim+1 ) {                  if ( nd_pot_nelim && MPOS(d1)>=nd_pot_nelim+1 && MPOS(d2) >= nd_pot_nelim+1 ) {
             if ( TD(d1) > TD(d2) ) return 1;              if ( TD(d1) > TD(d2) ) return 1;
Line 621  int ndl_module_glex_compare(UINT *d1,UINT *d2)
Line 697  int ndl_module_glex_compare(UINT *d1,UINT *d2)
 {  {
     int i,c;      int i,c;
   
       if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {      if ( nd_ispot ) {
         if ( MPOS(d1) < MPOS(d2) ) return 1;          if ( MPOS(d1) < MPOS(d2) ) return 1;
         else if ( MPOS(d1) > MPOS(d2) ) return -1;          else if ( MPOS(d1) > MPOS(d2) ) return -1;
Line 639  int ndl_module_lex_compare(UINT *d1,UINT *d2)
Line 716  int ndl_module_lex_compare(UINT *d1,UINT *d2)
 {  {
     int i,c;      int i,c;
   
       if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {      if ( nd_ispot ) {
         if ( MPOS(d1) < MPOS(d2) ) return 1;          if ( MPOS(d1) < MPOS(d2) ) return 1;
         else if ( MPOS(d1) > MPOS(d2) ) return -1;          else if ( MPOS(d1) > MPOS(d2) ) return -1;
Line 655  int ndl_module_block_compare(UINT *d1,UINT *d2)
Line 733  int ndl_module_block_compare(UINT *d1,UINT *d2)
 {  {
     int i,c;      int i,c;
   
       if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {      if ( nd_ispot ) {
         if ( MPOS(d1) < MPOS(d2) ) return 1;          if ( MPOS(d1) < MPOS(d2) ) return 1;
         else if ( MPOS(d1) > MPOS(d2) ) return -1;          else if ( MPOS(d1) > MPOS(d2) ) return -1;
Line 671  int ndl_module_matrix_compare(UINT *d1,UINT *d2)
Line 750  int ndl_module_matrix_compare(UINT *d1,UINT *d2)
 {  {
     int i,c;      int i,c;
   
       if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {      if ( nd_ispot ) {
         if ( MPOS(d1) < MPOS(d2) ) return 1;          if ( MPOS(d1) < MPOS(d2) ) return 1;
         else if ( MPOS(d1) > MPOS(d2) ) return -1;          else if ( MPOS(d1) > MPOS(d2) ) return -1;
Line 687  int ndl_module_composite_compare(UINT *d1,UINT *d2)
Line 767  int ndl_module_composite_compare(UINT *d1,UINT *d2)
 {  {
     int i,c;      int i,c;
   
       if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {      if ( nd_ispot ) {
         if ( MPOS(d1) > MPOS(d2) ) return 1;          if ( MPOS(d1) > MPOS(d2) ) return 1;
         else if ( MPOS(d1) < MPOS(d2) ) return -1;          else if ( MPOS(d1) < MPOS(d2) ) return -1;
Line 1924  NODE nd_gb(int m,int ishomo,int checkonly,int gensyz,i
Line 2005  NODE nd_gb(int m,int ishomo,int checkonly,int gensyz,i
     while ( d ) {      while ( d ) {
 again:  again:
         l = nd_minp(d,&d);          l = nd_minp(d,&d);
           if ( MaxDeg > 0 && SG(l) > MaxDeg ) break;
         if ( SG(l) != sugar ) {          if ( SG(l) != sugar ) {
             if ( ishomo ) {              if ( ishomo ) {
                 diag_count = 0;                  diag_count = 0;
Line 2182  NODE nd_gb_trace(int m,int ishomo,int **indp)
Line 2264  NODE nd_gb_trace(int m,int ishomo,int **indp)
     while ( d ) {      while ( d ) {
 again:  again:
         l = nd_minp(d,&d);          l = nd_minp(d,&d);
           if ( MaxDeg > 0 && SG(l) > MaxDeg ) break;
         if ( SG(l) != sugar ) {          if ( SG(l) != sugar ) {
 #if 1  #if 1
             if ( ishomo ) {              if ( ishomo ) {
Line 2279  again:
Line 2362  again:
         FREENDP(l);          FREENDP(l);
     }      }
     if ( nd_nalg ) {      if ( nd_nalg ) {
         print_eg("monic",&eg_monic);          if ( DP_Print ) {
         print_eg("invdalg",&eg_invdalg);            print_eg("monic",&eg_monic);
         print_eg("le",&eg_le);            print_eg("invdalg",&eg_invdalg);
             print_eg("le",&eg_le);
           }
     }      }
         conv_ilist(nd_demand,1,g,indp);          conv_ilist(nd_demand,1,g,indp);
     if ( DP_Print ) { printf("nd_gb_trace done.\n"); fflush(stdout); }      if ( DP_Print ) { printf("nd_gb_trace done.\n"); fflush(stdout); }
Line 2703  int ndv_newps(int m,NDV a,NDV aq,int f4)
Line 2788  int ndv_newps(int m,NDV a,NDV aq,int f4)
     nd_ps[nd_psn] = a;      nd_ps[nd_psn] = a;
     if ( aq ) {      if ( aq ) {
         nd_ps_trace[nd_psn] = aq;          nd_ps_trace[nd_psn] = aq;
                 nd_ps_gz[nd_psn] = ndvtondvgz(aq);                  if ( !nd_vc ) nd_ps_gz[nd_psn] = ndvtondvgz(aq);
         register_hcf(aq);          register_hcf(aq);
         nd_bound[nd_psn] = ndv_compute_bound(aq);          nd_bound[nd_psn] = ndv_compute_bound(aq);
         SG(r) = SG(aq); ndl_copy(HDL(aq),DL(r));          SG(r) = SG(aq); ndl_copy(HDL(aq),DL(r));
Line 2711  int ndv_newps(int m,NDV a,NDV aq,int f4)
Line 2796  int ndv_newps(int m,NDV a,NDV aq,int f4)
         if ( !m ) register_hcf(a);          if ( !m ) register_hcf(a);
         nd_bound[nd_psn] = ndv_compute_bound(a);          nd_bound[nd_psn] = ndv_compute_bound(a);
         SG(r) = SG(a); ndl_copy(HDL(a),DL(r));          SG(r) = SG(a); ndl_copy(HDL(a),DL(r));
                 if ( !m ) nd_ps_gz[nd_psn] = ndvtondvgz(a);                  if ( !m && !nd_vc ) nd_ps_gz[nd_psn] = ndvtondvgz(a);
     }      }
     if ( nd_demand ) {      if ( nd_demand ) {
         if ( aq ) {          if ( aq ) {
Line 2789  int ndv_setup(int mod,int trace,NODE f,int dont_sort,i
Line 2874  int ndv_setup(int mod,int trace,NODE f,int dont_sort,i
         hc = HCU(w[i].p);          hc = HCU(w[i].p);
         if ( trace ) {          if ( trace ) {
             a = nd_ps_trace[i] = ndv_dup(0,w[i].p);              a = nd_ps_trace[i] = ndv_dup(0,w[i].p);
                         nd_ps_gz[i] = ndvtondvgz(a);                          if ( !nd_vc ) nd_ps_gz[i] = ndvtondvgz(a);
             if ( !dont_removecont) ndv_removecont(0,a);              if ( !dont_removecont) ndv_removecont(0,a);
             register_hcf(a);              register_hcf(a);
             am = nd_ps[i] = ndv_dup(mod,a);              am = nd_ps[i] = ndv_dup(mod,a);
Line 2799  int ndv_setup(int mod,int trace,NODE f,int dont_sort,i
Line 2884  int ndv_setup(int mod,int trace,NODE f,int dont_sort,i
             ndv_removecont(mod,am);              ndv_removecont(mod,am);
         } else {          } else {
             a = nd_ps[i] = ndv_dup(mod,w[i].p);              a = nd_ps[i] = ndv_dup(mod,w[i].p);
                         if ( !mod ) nd_ps_gz[i] = ndvtondvgz(a);                          if ( !mod && !nd_vc ) nd_ps_gz[i] = ndvtondvgz(a);
             if ( mod || !dont_removecont ) ndv_removecont(mod,a);              if ( mod || !dont_removecont ) ndv_removecont(mod,a);
             if ( !mod ) register_hcf(a);              if ( !mod ) register_hcf(a);
         }          }
Line 2967  void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int 
Line 3052  void nd_gr(LIST f,LIST v,int m,int homo,int retdp,int 
     ndv_alloc = 0;      ndv_alloc = 0;
 #endif  #endif
     get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc);      get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc);
       if ( m && nd_vc )
          error("nd_{gr,f4} : computation over Fp(X) is unsupported. Use dp_gr_mod_main().");
     for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );      for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );
     switch ( ord->id ) {      switch ( ord->id ) {
         case 1:          case 1:
Line 3481  void nd_gr_trace(LIST f,LIST v,int trace,int homo,int 
Line 3568  void nd_gr_trace(LIST f,LIST v,int trace,int homo,int 
         for ( t = fd0; t; t = NEXT(t) )          for ( t = fd0; t; t = NEXT(t) )
             ndv_homogenize((NDV)BDY(t),obpe,oadv,oepos,ompos);              ndv_homogenize((NDV)BDY(t),obpe,oadv,oepos,ompos);
     }      }
       if ( MaxDeg > 0 ) nocheck = 1;
     while ( 1 ) {      while ( 1 ) {
                 tl1 = tl2 = tl3 = tl4 = 0;                  tl1 = tl2 = tl3 = tl4 = 0;
         if ( Demand )          if ( Demand )
Line 5193  NODE ndv_reducebase(NODE x,int *perm)
Line 5281  NODE ndv_reducebase(NODE x,int *perm)
   
 void nd_init_ord(struct order_spec *ord)  void nd_init_ord(struct order_spec *ord)
 {  {
     nd_module = (ord->id >= 256);    nd_module = (ord->id >= 256);
     if ( nd_module ) {
       nd_dcomp = -1;
       nd_ispot = ord->ispot;
       nd_pot_nelim = ord->pot_nelim;
       nd_poly_weight_len = ord->nv;
       nd_poly_weight = ord->top_weight;
       nd_module_rank = ord->module_rank;
       nd_module_weight = ord->module_top_weight;
     }
         nd_matrix = 0;          nd_matrix = 0;
         nd_matrix_len = 0;          nd_matrix_len = 0;
     switch ( ord->id ) {      switch ( ord->id ) {
Line 5249  void nd_init_ord(struct order_spec *ord)
Line 5346  void nd_init_ord(struct order_spec *ord)
   
         /* module order */          /* module order */
         case 256:          case 256:
             nd_ispot = ord->ispot;  
             nd_pot_nelim = ord->pot_nelim;  
             nd_dcomp = -1;  
             switch ( ord->ord.simple ) {              switch ( ord->ord.simple ) {
                 case 0:                  case 0:
                     nd_isrlex = 1;                      nd_isrlex = 1;
Line 5271  void nd_init_ord(struct order_spec *ord)
Line 5365  void nd_init_ord(struct order_spec *ord)
             break;              break;
         case 257:          case 257:
             /* block order */              /* block order */
             nd_ispot = ord->ispot;  
             nd_pot_nelim = ord->pot_nelim;  
             nd_dcomp = -1;  
             nd_isrlex = 0;              nd_isrlex = 0;
             ndl_compare_function = ndl_module_block_compare;              ndl_compare_function = ndl_module_block_compare;
             break;              break;
         case 258:          case 258:
             /* matrix order */              /* matrix order */
             nd_ispot = ord->ispot;  
             nd_pot_nelim = ord->pot_nelim;  
             nd_dcomp = -1;  
             nd_isrlex = 0;              nd_isrlex = 0;
             nd_matrix_len = ord->ord.matrix.row;              nd_matrix_len = ord->ord.matrix.row;
             nd_matrix = ord->ord.matrix.matrix;              nd_matrix = ord->ord.matrix.matrix;
Line 5289  void nd_init_ord(struct order_spec *ord)
Line 5377  void nd_init_ord(struct order_spec *ord)
             break;              break;
         case 259:          case 259:
             /* composite order */              /* composite order */
             nd_ispot = ord->ispot;  
             nd_pot_nelim = ord->pot_nelim;  
             nd_dcomp = -1;  
             nd_isrlex = 0;              nd_isrlex = 0;
             nd_worb_len = ord->ord.composite.length;              nd_worb_len = ord->ord.composite.length;
             nd_worb = ord->ord.composite.w_or_b;              nd_worb = ord->ord.composite.w_or_b;
Line 5476  int nd_to_vect_q(UINT *s0,int n,ND d,Q *r)
Line 5561  int nd_to_vect_q(UINT *s0,int n,ND d,Q *r)
     return i;      return i;
 }  }
   
   unsigned long *nd_to_vect_2(UINT *s0,int n,int *s0hash,ND p)
   {
       NM m;
       unsigned long *v;
       int i,j,h,size;
           UINT *s,*t;
   
           size = sizeof(unsigned long)*(n+BLEN-1)/BLEN;
       v = (unsigned long *)MALLOC_ATOMIC_IGNORE_OFF_PAGE(size);
       bzero(v,size);
       for ( i = j = 0, s = s0, m = BDY(p); m; j++, m = NEXT(m) ) {
                   t = DL(m);
                   h = ndl_hash_value(t);
           for ( ; h != s0hash[i] || !ndl_equal(t,s); s += nd_wpd, i++ );
               v[i/BLEN] |= 1L <<(i%BLEN);
       }
       return v;
   }
   
   int nd_nm_to_vect_2(UINT *s0,int n,int *s0hash,NDV p,NM m,unsigned long *v)
   {
       NMV mr;
       UINT *d,*t,*s;
       int i,j,len,h,head;
   
       d = DL(m);
       len = LEN(p);
       t = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));
       for ( i = j = 0, s = s0, mr = BDY(p); j < len; j++, NMV_ADV(mr) ) {
           ndl_add(d,DL(mr),t);
                   h = ndl_hash_value(t);
           for ( ; h != s0hash[i] || !ndl_equal(t,s); s += nd_wpd, i++ );
                   if ( j == 0 ) head = i;
               v[i/BLEN] |= 1L <<(i%BLEN);
       }
       return head;
   }
   
 Q *nm_ind_pair_to_vect(int mod,UINT *s0,int n,NM_ind_pair pair)  Q *nm_ind_pair_to_vect(int mod,UINT *s0,int n,NM_ind_pair pair)
 {  {
     NM m;      NM m;
Line 5875  NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea
Line 5998  NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea
     }      }
 }  }
   
   NDV vect_to_ndv_2(unsigned long *vect,int col,UINT *s0vect)
   {
       int j,k,len;
       UINT *p;
       NDV r;
       NMV mr0,mr;
   
       for ( j = 0, len = 0; j < col; j++ ) if ( vect[j/BLEN] & (1L<<(j%BLEN)) ) len++;
       if ( !len ) return 0;
       else {
           mr0 = (NMV)MALLOC_ATOMIC_IGNORE_OFF_PAGE(nmv_adv*len);
           mr = mr0;
           p = s0vect;
           for ( j = 0; j < col; j++, p += nd_wpd )
                     if ( vect[j/BLEN] & (1L<<(j%BLEN)) ) {
               ndl_copy(p,DL(mr)); CM(mr) = 1; NMV_ADV(mr);
             }
           MKNDV(nd_nvar,mr0,len,r);
           return r;
       }
   }
   
 /* for preprocessed vector */  /* for preprocessed vector */
   
 NDV vect_to_ndv_q(Q *vect,int spcol,int col,int *rhead,UINT *s0vect)  NDV vect_to_ndv_q(Q *vect,int spcol,int col,int *rhead,UINT *s0vect)
Line 6063  NODE nd_f4(int m,int **indp)
Line 6208  NODE nd_f4(int m,int **indp)
         get_eg(&eg0);          get_eg(&eg0);
         l = nd_minsugarp(d,&d);          l = nd_minsugarp(d,&d);
         sugar = SG(l);          sugar = SG(l);
           if ( MaxDeg > 0 && sugar > MaxDeg ) break;
         if ( nd_nzlist ) {          if ( nd_nzlist ) {
             for ( tn = nd_nzlist; tn; tn = NEXT(tn) ) {              for ( tn = nd_nzlist; tn; tn = NEXT(tn) ) {
                 node = BDY((LIST)BDY(tn));                  node = BDY((LIST)BDY(tn));
                             if ( QTOS((Q)ARG0(node)) == sugar ) break;                              if ( QTOS((Q)ARG0(node)) == sugar ) break;
             }              }
             if ( !tn ) error("nd_f4 : inconsistent non-zero list");                          if ( tn ) {
                         for ( t = l, ll0 = 0; t; t = NEXT(t) ) {                                  for ( t = l, ll0 = 0; t; t = NEXT(t) ) {
                 for ( tn = BDY((LIST)ARG1(node)); tn; tn = NEXT(tn) ) {                          for ( tn = BDY((LIST)ARG1(node)); tn; tn = NEXT(tn) ) {
                                   i1s = QTOS((Q)ARG0(BDY((LIST)BDY(tn))));                                          i1s = QTOS((Q)ARG0(BDY((LIST)BDY(tn))));
                                   i2s = QTOS((Q)ARG1(BDY((LIST)BDY(tn))));                                          i2s = QTOS((Q)ARG1(BDY((LIST)BDY(tn))));
                                   if ( t->i1 == i1s && t->i2 == i2s ) break;                                          if ( t->i1 == i1s && t->i2 == i2s ) break;
                                 }                                          }
                             if ( tn ) {                                  if ( tn ) {
                                     if ( !ll0 ) ll0 = t;                                          if ( !ll0 ) ll0 = t;
                                         else NEXT(ll) = t;                                                  else NEXT(ll) = t;
                                         ll = t;                                                  ll = t;
                                 }                                          }
             }                  }
                         if ( ll0 ) NEXT(ll) = 0;                                  if ( ll0 ) NEXT(ll) = 0;
                     l = ll0;                          l = ll0;
                           } else l = 0;
         }          }
         bucket = create_pbucket();          bucket = create_pbucket();
         stat = nd_sp_f4(m,0,l,bucket);          stat = nd_sp_f4(m,0,l,bucket);
Line 6174  NODE nd_f4_trace(int m,int **indp)
Line 6321  NODE nd_f4_trace(int m,int **indp)
         get_eg(&eg0);          get_eg(&eg0);
         l = nd_minsugarp(d,&d);          l = nd_minsugarp(d,&d);
         sugar = SG(l);          sugar = SG(l);
           if ( MaxDeg > 0 && sugar > MaxDeg ) break;
         bucket = create_pbucket();          bucket = create_pbucket();
         stat = nd_sp_f4(m,0,l,bucket);          stat = nd_sp_f4(m,0,l,bucket);
         if ( !stat ) {          if ( !stat ) {
Line 6347  NODE nd_f4_pseudo_trace(int m,int **indp)
Line 6495  NODE nd_f4_pseudo_trace(int m,int **indp)
     return g;      return g;
 }  }
   
   int rref(matrix mat,int *sugar)
   {
     int row,col,i,j,k,l,s,wcol,wj;
     unsigned long bj;
     unsigned long **a;
     unsigned long *ai,*ak,*as,*t;
     int *pivot;
   
     row = mat->row;
     col = mat->col;
     a = mat->a;
     wcol = (col+BLEN-1)/BLEN;
     pivot = (int *)MALLOC_ATOMIC(row*sizeof(int));
     i = 0;
     for ( j = 0; j < col; j++ ) {
           wj = j/BLEN; bj = 1L<<(j%BLEN);
       for ( k = i; k < row; k++ )
             if ( a[k][wj] & bj ) break;
       if ( k == row ) continue;
           pivot[i] = j;
       if ( k != i ) {
            t = a[i]; a[i] = a[k]; a[k] = t;
            s = sugar[i]; sugar[i] = sugar[k]; sugar[k] = s;
           }
           ai = a[i];
       for ( k = i+1; k < row; k++ ) {
             ak = a[k];
             if ( ak[wj] & bj ) {
               for ( l = wj; l < wcol; l++ )
                     ak[l] ^= ai[l];
               sugar[k] = MAX(sugar[k],sugar[i]);
             }
           }
           i++;
     }
     for ( k = i-1; k >= 0; k-- ) {
       j = pivot[k]; wj = j/BLEN; bj = 1L<<(j%BLEN);
           ak = a[k];
       for ( s = 0; s < k; s++ ) {
             as = a[s];
         if ( as[wj] & bj ) {
           for ( l = wj; l < wcol; l++ )
                     as[l] ^= ak[l];
               sugar[s] = MAX(sugar[s],sugar[k]);
             }
           }
     }
     return i;
   }
   
   void print_matrix(matrix mat)
   {
     int row,col,i,j;
     unsigned long *ai;
   
     row = mat->row;
     col = mat->col;
     printf("%d x %d\n",row,col);
     for ( i = 0; i < row; i++ ) {
           ai = mat->a[i];
       for ( j = 0; j < col; j++ ) {
             if ( ai[j/BLEN] & (1L<<(j%BLEN)) ) putchar('1');
             else putchar('0');
           }
           putchar('\n');
     }
   }
   
   NDV vect_to_ndv_2(unsigned long *vect,int col,UINT *s0vect);
   
   void red_by_vect_2(matrix mat,int *sugar,unsigned long *v,int rhead,int rsugar)
   {
     int row,col,wcol,wj,i,j;
     unsigned long bj;
     unsigned long *ai;
     unsigned long **a;
     int len;
     int *pos;
   
     row = mat->row;
     col = mat->col;
     wcol = (col+BLEN-1)/BLEN;
     pos = (int *)ALLOCA(wcol*sizeof(int));
     bzero(pos,wcol*sizeof(int));
     for ( i = j = 0; i < wcol; i++ )
       if ( v[i] ) pos[j++] = i;;
     len = j;
     wj = rhead/BLEN;
     bj = 1L<<rhead%BLEN;
     a = mat->a;
     for ( i = 0; i < row; i++ ) {
           ai = a[i];
       if ( ai[wj]&bj ) {
             for ( j = 0; j < len; j++ )
               ai[pos[j]] ^= v[pos[j]];
             sugar[i] = MAX(sugar[i],rsugar);
           }
     }
   }
   
   NODE nd_f4_red_2(ND_pairs sp0,UINT *s0vect,int col,NODE rp0,ND_pairs *nz)
   {
       int nsp,nred,i,i0,k,rank,row;
       NODE r0,rp;
       ND_pairs sp;
           ND spol;
           NM_ind_pair rt;
       int *s0hash;
           UINT *s;
           int *pivot,*sugar,*head;
           matrix mat;
       NM m;
       NODE r;
           struct oEGT eg0,eg1,eg2,eg_elim1,eg_elim2;
           int rhead,rsugar,size;
       unsigned long *v;
   
       get_eg(&eg0);
   init_eg(&eg_search);
       for ( sp = sp0, nsp = 0; sp; sp = NEXT(sp), nsp++ );
       nred = length(rp0);
       mat = alloc_matrix(nsp,col);
       s0hash = (int *)ALLOCA(col*sizeof(int));
       for ( i = 0, s = s0vect; i < col; i++, s += nd_wpd )
           s0hash[i] = ndl_hash_value(s);
   
           sugar = (int *)ALLOCA(nsp*sizeof(int));
           for ( i = 0, sp = sp0; sp; sp = NEXT(sp) ) {
                   nd_sp(2,0,sp,&spol);
                   if ( spol ) {
                 mat->a[i] = nd_to_vect_2(s0vect,col,s0hash,spol);
                     sugar[i] = SG(spol);
                     i++;
                   }
           }
           mat->row = i;
       if ( DP_Print ) {
         fprintf(asir_out,"%dx%d,",mat->row,mat->col); fflush(asir_out);
       }
           size = ((col+BLEN-1)/BLEN)*sizeof(unsigned long);
           v = CALLOC((col+BLEN-1)/BLEN,sizeof(unsigned long));
       for ( rp = rp0, i = 0; rp; rp = NEXT(rp), i++ ) {
                   rt = (NM_ind_pair)BDY(rp);
                   bzero(v,size);
           rhead = nd_nm_to_vect_2(s0vect,col,s0hash,nd_ps[rt->index],rt->mul,v);
                   rsugar = SG(nd_ps[rt->index])+TD(DL(rt->mul));
               red_by_vect_2(mat,sugar,v,rhead,rsugar);
           }
   
       get_eg(&eg1);
       init_eg(&eg_elim1); add_eg(&eg_elim1,&eg0,&eg1);
           rank = rref(mat,sugar);
   
       for ( i = 0, r0 = 0; i < rank; i++ ) {
         NEXTNODE(r0,r);
             BDY(r) = (pointer)vect_to_ndv_2(mat->a[i],col,s0vect);
         SG((NDV)BDY(r)) = sugar[i];
       }
       if ( r0 ) NEXT(r) = 0;
       get_eg(&eg2);
       init_eg(&eg_elim2); add_eg(&eg_elim2,&eg1,&eg2);
       if ( DP_Print ) {
           fprintf(asir_out,"elim1=%fsec,elim2=%fsec\n",
                     eg_elim1.exectime+eg_elim1.gctime,eg_elim2.exectime+eg_elim2.gctime);
           fflush(asir_out);
           }
       return r0;
   }
   
   
 NODE nd_f4_red(int m,ND_pairs sp0,int trace,UINT *s0vect,int col,NODE rp0,ND_pairs *nz)  NODE nd_f4_red(int m,ND_pairs sp0,int trace,UINT *s0vect,int col,NODE rp0,ND_pairs *nz)
 {  {
     IndArray *imat;      IndArray *imat;
Line 6358  NODE nd_f4_red(int m,ND_pairs sp0,int trace,UINT *s0ve
Line 6676  NODE nd_f4_red(int m,ND_pairs sp0,int trace,UINT *s0ve
     UINT *s;      UINT *s;
     int *s0hash;      int *s0hash;
   
       if ( m == 2 && nd_rref2 )
              return nd_f4_red_2(sp0,s0vect,col,rp0,nz);
   
 init_eg(&eg_search);  init_eg(&eg_search);
     for ( sp = sp0, nsp = 0; sp; sp = NEXT(sp), nsp++ );      for ( sp = sp0, nsp = 0; sp; sp = NEXT(sp), nsp++ );
     nred = length(rp0);      nred = length(rp0);
Line 6366  init_eg(&eg_search);
Line 6687  init_eg(&eg_search);
     for ( i = 0; i < col; i++ ) rhead[i] = 0;      for ( i = 0; i < col; i++ ) rhead[i] = 0;
   
     /* construction of index arrays */      /* construction of index arrays */
         fprintf(stderr,"%dx%d,",nsp+nred,col);      if ( DP_Print ) {
             fprintf(stderr,"%dx%d,",nsp+nred,col);
       }
     rvect = (NM_ind_pair *)ALLOCA(nred*sizeof(NM_ind_pair));      rvect = (NM_ind_pair *)ALLOCA(nred*sizeof(NM_ind_pair));
     s0hash = (int *)ALLOCA(col*sizeof(int));      s0hash = (int *)ALLOCA(col*sizeof(int));
     for ( i = 0, s = s0vect; i < col; i++, s += nd_wpd )      for ( i = 0, s = s0vect; i < col; i++, s += nd_wpd )
Line 6380  init_eg(&eg_search);
Line 6703  init_eg(&eg_search);
         r0 = nd_f4_red_main(m,sp0,nsp,s0vect,col,rvect,rhead,imat,nred,nz);          r0 = nd_f4_red_main(m,sp0,nsp,s0vect,col,rvect,rhead,imat,nred,nz);
     else      else
         r0 = nd_f4_red_gz_main(sp0,nsp,trace,s0vect,col,rvect,rhead,imat,nred);          r0 = nd_f4_red_gz_main(sp0,nsp,trace,s0vect,col,rvect,rhead,imat,nred);
 print_eg("search",&eg_search);      if ( DP_Print ) print_eg("search",&eg_search);
     return r0;      return r0;
 }  }
   
Line 6978  int nd_gauss_elim_gz(GZ **mat0,int *sugar,int row,int 
Line 7301  int nd_gauss_elim_gz(GZ **mat0,int *sugar,int row,int 
 {  {
     int i,j,t,c,rank,inv;      int i,j,t,c,rank,inv;
     int *ci,*ri;      int *ci,*ri;
     GZ *dn;      GZ dn;
     MAT m,nm;      MAT m,nm;
   
     NEWMAT(m); m->row = row; m->col = col; m->body = (pointer **)mat0;      NEWMAT(m); m->row = row; m->col = col; m->body = (pointer **)mat0;
         dn = (GZ *)MALLOC(row*sizeof(GZ));      rank = gz_generic_gauss_elim(m,&nm,&dn,&ri,&ci);
     rank = gz_generic_gauss_elim2(m,&nm,dn,&ri,&ci);  
     for ( i = 0; i < row; i++ )      for ( i = 0; i < row; i++ )
         for ( j = 0; j < col; j++ )          for ( j = 0; j < col; j++ )
             mat0[i][j] = 0;              mat0[i][j] = 0;
     c = col-rank;      c = col-rank;
     for ( i = 0; i < rank; i++ ) {      for ( i = 0; i < rank; i++ ) {
         mat0[i][ri[i]] = dn[i];          mat0[i][ri[i]] = dn;
         for ( j = 0; j < c; j++ )          for ( j = 0; j < c; j++ )
             mat0[i][ci[j]] = (GZ)BDY(nm)[i][j];              mat0[i][ci[j]] = (GZ)BDY(nm)[i][j];
     }      }
Line 7354  void nd_det(int mod,MAT f,P *rp)
Line 7676  void nd_det(int mod,MAT f,P *rp)
     if ( mod ) ndv_mod(mod,d);      if ( mod ) ndv_mod(mod,d);
     chsgnq(ONE,&mone);      chsgnq(ONE,&mone);
     for ( j = 0, sgn = 1; j < n; j++ ) {      for ( j = 0, sgn = 1; j < n; j++ ) {
         if ( DP_Print ) fprintf(stderr,".",j);          if ( DP_Print ) {
             fprintf(stderr,".",j);
           }
         for ( i = j; i < n && !dm[i][j]; i++ );          for ( i = j; i < n && !dm[i][j]; i++ );
         if ( i == n ) {          if ( i == n ) {
             *rp = 0;              *rp = 0;
Line 7413  void nd_det(int mod,MAT f,P *rp)
Line 7737  void nd_det(int mod,MAT f,P *rp)
         }          }
         d = mjj;          d = mjj;
     }      }
     if ( DP_Print ) fprintf(stderr,"\n",k);      if ( DP_Print ) {
         fprintf(stderr,"\n",k);
       }
     if ( sgn < 0 )      if ( sgn < 0 )
         if ( mod )          if ( mod )
             ndv_mul_c(mod,d,mod-1);              ndv_mul_c(mod,d,mod-1);
Line 7948  MAT nd_btog(LIST f,LIST v,int mod,struct order_spec *o
Line 8274  MAT nd_btog(LIST f,LIST v,int mod,struct order_spec *o
         pi0 = QTOS((Q)ARG0(pi)); pi1 = QTOS((Q)ARG1(pi));          pi0 = QTOS((Q)ARG0(pi)); pi1 = QTOS((Q)ARG1(pi));
         p[pi0] = c = (ND *)MALLOC(nb*sizeof(ND));          p[pi0] = c = (ND *)MALLOC(nb*sizeof(ND));
         ptomp(mod,(P)ARG2(pi),&inv);          ptomp(mod,(P)ARG2(pi),&inv);
           ((MQ)inv)->cont = invm(((MQ)inv)->cont,mod);
         u = ptond(CO,vv,(P)ONE);          u = ptond(CO,vv,(P)ONE);
         HCM(u) = ((MQ)inv)->cont;          HCM(u) = ((MQ)inv)->cont;
         c[pi1] = u;          c[pi1] = u;
Line 8019  VECT nd_btog_one(LIST f,LIST v,int mod,struct order_sp
Line 8346  VECT nd_btog_one(LIST f,LIST v,int mod,struct order_sp
     pi = BDY((LIST)BDY(t));      pi = BDY((LIST)BDY(t));
         pi0 = QTOS((Q)ARG0(pi)); pi1 = QTOS((Q)ARG1(pi));          pi0 = QTOS((Q)ARG0(pi)); pi1 = QTOS((Q)ARG1(pi));
         if ( pi1 == pos ) {          if ( pi1 == pos ) {
           ptomp(mod,(P)ARG2(pi),&inv);            ptomp(mod,(P)ARG2(pi),&inv);
             ((MQ)inv)->cont = invm(((MQ)inv)->cont,mod);
           u = ptond(CO,vv,(P)ONE);            u = ptond(CO,vv,(P)ONE);
           HCM(u) = ((MQ)inv)->cont;            HCM(u) = ((MQ)inv)->cont;
           p[pi0] = u;            p[pi0] = u;

Legend:
Removed from v.1.215  
changed lines
  Added in v.1.229

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