[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.239 and 1.240

version 1.239, 2017/09/14 01:34:53 version 1.240, 2017/09/15 01:52:51
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.238 2017/08/31 02:36:21 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.239 2017/09/14 01:34:53 noro Exp $ */
   
 #include "nd.h"  #include "nd.h"
   
Line 100  void nd_mul_c_lf(ND p,GZ mul);
Line 100  void nd_mul_c_lf(ND p,GZ mul);
 void ndv_mul_c_lf(NDV p,GZ mul);  void ndv_mul_c_lf(NDV p,GZ mul);
 NODE nd_f4_red_main(int m,ND_pairs sp0,int nsp,UINT *s0vect,int col,  NODE nd_f4_red_main(int m,ND_pairs sp0,int nsp,UINT *s0vect,int col,
         NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred,ND_pairs *nz);          NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred,ND_pairs *nz);
 NODE nd_f4_red_mod128_main(int m,ND_pairs sp0,int nsp,UINT *s0vect,int col,  NODE nd_f4_red_mod64_main(int m,ND_pairs sp0,int nsp,UINT *s0vect,int col,
         NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred,ND_pairs *nz);          NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred,ND_pairs *nz);
 NODE nd_f4_red_lf_main(int m,ND_pairs sp0,int nsp,int trace,UINT *s0vect,int col,  NODE nd_f4_red_lf_main(int m,ND_pairs sp0,int nsp,int trace,UINT *s0vect,int col,
         NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred);          NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred);
Line 5823  int nd_to_vect(int mod,UINT *s0,int n,ND d,UINT *r)
Line 5823  int nd_to_vect(int mod,UINT *s0,int n,ND d,UINT *r)
 }  }
   
 #if defined(__GNUC__)  #if defined(__GNUC__)
 int nd_to_vect128(int mod,UINT *s0,int n,ND d,U128 *r)  
   #define MOD128(a,c,m) ((a)=(((c)!=0||((a)>=(m)))?(((((U128)(c))<<64)+(a))%(m)):(a)))
   
   int nd_to_vect64(int mod,UINT *s0,int n,ND d,U64 *r)
 {  {
     NM m;      NM m;
     UINT *t,*s;      UINT *t,*s;
Line 5833  int nd_to_vect128(int mod,UINT *s0,int n,ND d,U128 *r)
Line 5836  int nd_to_vect128(int mod,UINT *s0,int n,ND d,U128 *r)
     for ( i = 0, s = s0, m = BDY(d); m; m = NEXT(m) ) {      for ( i = 0, s = s0, m = BDY(d); m; m = NEXT(m) ) {
         t = DL(m);          t = DL(m);
         for ( ; !ndl_equal(t,s); s += nd_wpd, i++ );          for ( ; !ndl_equal(t,s); s += nd_wpd, i++ );
         r[i] = (U128)CM(m);          r[i] = (U64)CM(m);
     }      }
     for ( i = 0; !r[i]; i++ );      for ( i = 0; !r[i]; i++ );
     return i;      return i;
Line 6235  int ndv_reduce_vect(int m,UINT *svect,int col,IndArray
Line 6238  int ndv_reduce_vect(int m,UINT *svect,int col,IndArray
 }  }
   
 #if defined(__GNUC__)  #if defined(__GNUC__)
 int ndv_reduce_vect128(int m,U128 *svect,int col,IndArray *imat,NM_ind_pair *rp0,int nred)  
   int ndv_reduce_vect64(int m,U64 *svect,U64 *cvect,int col,IndArray *imat,NM_ind_pair *rp0,int nred)
 {  {
     int i,j,k,len,pos,prev;      int i,j,k,len,pos,prev;
     U64 c,c1;      U64 a,c,c1,c2;
     U128 a;  
     IndArray ivect;      IndArray ivect;
     unsigned char *ivc;      unsigned char *ivc;
     unsigned short *ivs;      unsigned short *ivs;
Line 6249  int ndv_reduce_vect128(int m,U128 *svect,int col,IndAr
Line 6252  int ndv_reduce_vect128(int m,U128 *svect,int col,IndAr
     NODE rp;      NODE rp;
     int maxrs;      int maxrs;
   
       for ( i = 0; i < col; i++ ) cvect[i] = 0;
     maxrs = 0;      maxrs = 0;
     for ( i = 0; i < nred; i++ ) {      for ( i = 0; i < nred; i++ ) {
         ivect = imat[i];          ivect = imat[i];
         k = ivect->head; svect[k] %= m;          k = ivect->head;
           a = svect[k]; c = cvect[k];
           MOD128(a,c,m);
           svect[k] = a; cvect[k] = 0;
         if ( c = svect[k] ) {          if ( c = svect[k] ) {
             maxrs = MAX(maxrs,rp0[i]->sugar);              maxrs = MAX(maxrs,rp0[i]->sugar);
             c = m-c; redv = nd_ps[rp0[i]->index];              c = m-c; redv = nd_ps[rp0[i]->index];
Line 6263  int ndv_reduce_vect128(int m,U128 *svect,int col,IndAr
Line 6270  int ndv_reduce_vect128(int m,U128 *svect,int col,IndAr
                     ivc = ivect->index.c;                      ivc = ivect->index.c;
                     for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) {                      for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) {
                         pos = prev+ivc[j]; c1 = CM(mr); prev = pos;                          pos = prev+ivc[j]; c1 = CM(mr); prev = pos;
                                                 if ( c1 ) svect[pos] += c1*c;                          if ( c1 ) {
                             c2 = svect[pos]+c1*c;
                             if ( c2 < svect[pos] ) cvect[pos]++;
                             svect[pos] = c2;
                           }
                     }                      }
                     break;                      break;
                 case 2:                  case 2:
                     ivs = ivect->index.s;                      ivs = ivect->index.s;
                     for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) {                      for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) {
                         pos = prev+ivs[j]; c1 = CM(mr); prev = pos;                          pos = prev+ivs[j]; c1 = CM(mr); prev = pos;
                                                 if ( c1 ) svect[pos] += c1*c;                          if ( c1 ) {
                             c2 = svect[pos]+c1*c;
                             if ( c2 < svect[pos] ) cvect[pos]++;
                             svect[pos] = c2;
                           }
                     }                      }
                     break;                      break;
                 case 4:                  case 4:
                     ivi = ivect->index.i;                      ivi = ivect->index.i;
                     for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) {                      for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) {
                         pos = prev+ivi[j]; c1 = CM(mr); prev = pos;                          pos = prev+ivi[j]; c1 = CM(mr); prev = pos;
                                                 if ( c1 ) svect[pos] += c1*c;                          if ( c1 ) {
                             c2 = svect[pos]+c1*c;
                             if ( c2 < svect[pos] ) cvect[pos]++;
                             svect[pos] = c2;
                           }
                     }                      }
                     break;                      break;
             }              }
         }          }
     }      }
     for ( i = 0; i < col; i++ ) svect[i] %= m;      for ( i = 0; i < col; i++ ) {
         a = svect[i]; c = cvect[i]; MOD128(a,c,m); svect[i] = a;
       }
     return maxrs;      return maxrs;
 }  }
 #endif  #endif
Line 6546  NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea
Line 6567  NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea
 }  }
   
 #if defined(__GNUC__)  #if defined(__GNUC__)
 NDV vect128o_ndv(U128 *vect,int spcol,int col,int *rhead,UINT *s0vect)  NDV vect64_to_ndv(U64 *vect,int spcol,int col,int *rhead,UINT *s0vect)
 {  {
     int j,k,len;      int j,k,len;
     UINT *p;      UINT *p;
Line 7206  init_eg(&eg_search);
Line 7227  init_eg(&eg_search);
         imat[i] = nm_ind_pair_to_vect_compress(trace,s0vect,col,s0hash,rvect[i]);          imat[i] = nm_ind_pair_to_vect_compress(trace,s0vect,col,s0hash,rvect[i]);
         rhead[imat[i]->head] = 1;          rhead[imat[i]->head] = 1;
     }      }
 #if defined(__GNUC__)  
     if ( m >= (1<<15) )  
         r0 = nd_f4_red_mod128_main(m,sp0,nsp,s0vect,col,rvect,rhead,imat,nred,nz);  
     else  
 #endif  
     if ( m > 0 )      if ( m > 0 )
   #if defined(__GNUC__)
           r0 = nd_f4_red_mod64_main(m,sp0,nsp,s0vect,col,rvect,rhead,imat,nred,nz);
   #else
         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);
   #endif
     else if ( m == -1 )      else if ( m == -1 )
         r0 = nd_f4_red_sf_main(m,sp0,nsp,s0vect,col,rvect,rhead,imat,nred,nz);          r0 = nd_f4_red_sf_main(m,sp0,nsp,s0vect,col,rvect,rhead,imat,nred,nz);
     else if ( m == -2 )      else if ( m == -2 )
Line 7318  NODE nd_f4_red_main(int m,ND_pairs sp0,int nsp,UINT *s
Line 7338  NODE nd_f4_red_main(int m,ND_pairs sp0,int nsp,UINT *s
 #if defined(__GNUC__)  #if defined(__GNUC__)
 /* for Fp, 2^15=<p<2^29 */  /* for Fp, 2^15=<p<2^29 */
   
 NODE nd_f4_red_mod128_main(int m,ND_pairs sp0,int nsp,UINT *s0vect,int col,  NODE nd_f4_red_mod64_main(int m,ND_pairs sp0,int nsp,UINT *s0vect,int col,
         NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred,ND_pairs *nz)          NM_ind_pair *rvect,int *rhead,IndArray *imat,int nred,ND_pairs *nz)
 {  {
     int spcol,sprow,a;      int spcol,sprow,a;
Line 7326  NODE nd_f4_red_mod128_main(int m,ND_pairs sp0,int nsp,
Line 7346  NODE nd_f4_red_mod128_main(int m,ND_pairs sp0,int nsp,
     NODE r0,r;      NODE r0,r;
     ND_pairs sp;      ND_pairs sp;
     ND spol;      ND spol;
     U128 **spmat;      U64 **spmat;
     U128 *svect,*v;      U64 *svect,*cvect;
       U64 *v;
     int *colstat;      int *colstat;
     struct oEGT eg0,eg1,eg2,eg_f4,eg_f4_1,eg_f4_2;      struct oEGT eg0,eg1,eg2,eg_f4,eg_f4_1,eg_f4_2;
     int maxrs;      int maxrs;
     int *spsugar;      int *spsugar;
     ND_pairs *spactive;      ND_pairs *spactive;
     int **spmat32;  
   
     spcol = col-nred;      spcol = col-nred;
     get_eg(&eg0);      get_eg(&eg0);
     /* elimination (1st step) */      /* elimination (1st step) */
     spmat = (U128 **)MALLOC(nsp*sizeof(U128 *));      spmat = (U64 **)MALLOC(nsp*sizeof(U64 *));
     svect = (U128 *)MALLOC(col*sizeof(U128));      svect = (U64 *)MALLOC(col*sizeof(U64));
       cvect = (U64 *)MALLOC(col*sizeof(U64));
     spsugar = (int *)MALLOC(nsp*sizeof(int));      spsugar = (int *)MALLOC(nsp*sizeof(int));
     spactive = !nz?0:(ND_pairs *)MALLOC(nsp*sizeof(ND_pairs));      spactive = !nz?0:(ND_pairs *)MALLOC(nsp*sizeof(ND_pairs));
     for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) {      for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) {
         nd_sp(m,0,sp,&spol);          nd_sp(m,0,sp,&spol);
         if ( !spol ) continue;          if ( !spol ) continue;
         nd_to_vect128(m,s0vect,col,spol,svect);          nd_to_vect64(m,s0vect,col,spol,svect);
         maxrs = ndv_reduce_vect128(m,svect,col,imat,rvect,nred);          maxrs = ndv_reduce_vect64(m,svect,cvect,col,imat,rvect,nred);
         for ( i = 0; i < col; i++ ) if ( svect[i] ) break;          for ( i = 0; i < col; i++ ) if ( svect[i] ) break;
         if ( i < col ) {          if ( i < col ) {
             spmat[sprow] = v = (U128 *)MALLOC_ATOMIC(spcol*sizeof(U128));              spmat[sprow] = v = (U64 *)MALLOC_ATOMIC(spcol*sizeof(U64));
             for ( j = k = 0; j < col; j++ )              for ( j = k = 0; j < col; j++ )
                 if ( !rhead[j] ) v[k++] = (UINT)svect[j];                  if ( !rhead[j] ) v[k++] = (UINT)svect[j];
             spsugar[sprow] = MAX(maxrs,SG(spol));              spsugar[sprow] = MAX(maxrs,SG(spol));
Line 7369  NODE nd_f4_red_mod128_main(int m,ND_pairs sp0,int nsp,
Line 7390  NODE nd_f4_red_mod128_main(int m,ND_pairs sp0,int nsp,
   
     /* elimination (2nd step) */      /* elimination (2nd step) */
     colstat = (int *)MALLOC(spcol*sizeof(int));      colstat = (int *)MALLOC(spcol*sizeof(int));
     rank = nd_gauss_elim_mod128(spmat,spsugar,spactive,sprow,spcol,m,colstat);      rank = nd_gauss_elim_mod64(spmat,spsugar,spactive,sprow,spcol,m,colstat);
     r0 = 0;      r0 = 0;
     for ( i = 0; i < rank; i++ ) {      for ( i = 0; i < rank; i++ ) {
         NEXTNODE(r0,r); BDY(r) =          NEXTNODE(r0,r); BDY(r) =
           (pointer)vect128o_ndv(spmat[i],spcol,col,rhead,s0vect);            (pointer)vect64_to_ndv(spmat[i],spcol,col,rhead,s0vect);
         SG((NDV)BDY(r)) = spsugar[i];          SG((NDV)BDY(r)) = spsugar[i];
         GCFREE(spmat[i]);          GCFREE(spmat[i]);
     }      }
Line 7977  int nd_gauss_elim_mod(int **mat0,int *sugar,ND_pairs *
Line 7998  int nd_gauss_elim_mod(int **mat0,int *sugar,ND_pairs *
 }  }
   
 #if defined(__GNUC__)  #if defined(__GNUC__)
 int nd_gauss_elim_mod128(U128 **mat,int *sugar,ND_pairs *spactive,int row,int col,int md,int *colstat)  
   int nd_gauss_elim_mod64(U64 **mat,int *sugar,ND_pairs *spactive,int row,int col,int md,int *colstat)
 {  {
     int i,j,k,l,rank,s;    int i,j,k,l,rank,s;
     U64 inv;    U64 inv;
     U128 a;    U64 a;
     U128 *t,*pivot,*pk;    UINT c;
     ND_pairs pair;    U64 *t,*pivot,*pk;
     UINT *ck;
     UINT **cmat;
     UINT *ct;
     ND_pairs pair;
   
     for ( rank = 0, j = 0; j < col; j++ ) {    cmat = (UINT **)MALLOC(row*sizeof(UINT *));
         for ( i = rank; i < row; i++ )    for ( i = 0; i < row; i++ ) {
             mat[i][j] %= md;      cmat[i] = MALLOC_ATOMIC(col*sizeof(UINT));
         for ( i = rank; i < row; i++ )      bzero(cmat[i],col*sizeof(UINT));
             if ( mat[i][j] )    }
                 break;  
         if ( i == row ) {    for ( rank = 0, j = 0; j < col; j++ ) {
             colstat[j] = 0;      for ( i = rank; i < row; i++ ) {
             continue;        a = mat[i][j]; c = cmat[i][j];
         } else        MOD128(a,c,md);
             colstat[j] = 1;        mat[i][j] = a; cmat[i][j] = 0;
         if ( i != rank ) {  
             t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;  
             s = sugar[i]; sugar[i] = sugar[rank]; sugar[rank] = s;  
             if ( spactive ) {  
                 pair = spactive[i]; spactive[i] = spactive[rank];  
                 spactive[rank] = pair;  
             }  
         }  
         pivot = mat[rank];  
         s = sugar[rank];  
         inv = invm((UINT)pivot[j],md);  
         for ( k = j, pk = pivot+k; k < col; k++, pk++ )  
             if ( a = *pk ) {  
                 if ( a >= md ) a %= md;  
                 *pk = ((U64)a*inv)%md;  
             }  
         for ( i = rank+1; i < row; i++ ) {  
             t = mat[i];  
             if ( a = (t[j]%md) ) {  
                 sugar[i] = MAX(sugar[i],s);  
                 red_by_vect128(md,t+j,pivot+j,(int)(md-a),col-j);  
             }  
         }  
         rank++;  
     }      }
     for ( j = col-1, l = rank-1; j >= 0; j-- )      for ( i = rank; i < row; i++ )
         if ( colstat[j] ) {        if ( mat[i][j] )
             pivot = mat[l];          break;
             for ( k = j; k < col; k++ )      if ( i == row ) {
               if ( pivot[k] >= md ) pivot[k] %= md;        colstat[j] = 0;
             s = sugar[l];        continue;
             for ( i = 0; i < l; i++ ) {      } else
                 t = mat[i];        colstat[j] = 1;
                 t[j] %= md;      if ( i != rank ) {
                 if ( a = (t[j]%md) ) {        t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                     sugar[i] = MAX(sugar[i],s);        ct = cmat[i]; cmat[i] = cmat[rank]; cmat[rank] = ct;
                     red_by_vect128(md,t+j,pivot+j,(int)(md-a),col-j);        s = sugar[i]; sugar[i] = sugar[rank]; sugar[rank] = s;
                 }        if ( spactive ) {
             }          pair = spactive[i]; spactive[i] = spactive[rank];
             l--;          spactive[rank] = pair;
         }
       }
       /* column j is normalized */
       s = sugar[rank];
       inv = invm((UINT)mat[rank][j],md);
       /* normalize pivot row */
       for ( k = j, pk = mat[rank]+j, ck = cmat[rank]+j; k < col; k++, pk++, ck++ ) {
         a = *pk; c = *ck; MOD128(a,c,md); *pk = (a*inv)%md; *ck = 0;
       }
       for ( i = rank+1; i < row; i++ ) {
         if ( a = mat[i][j] ) {
           sugar[i] = MAX(sugar[i],s);
           red_by_vect64(md,mat[i]+j,cmat[i]+j,mat[rank]+j,(int)(md-a),col-j);
         }
       }
       rank++;
     }
     for ( j = col-1, l = rank-1; j >= 0; j-- )
       if ( colstat[j] ) {
         for ( k = j, pk = mat[l]+j, ck = cmat[l]+j; k < col; k++, pk++, ck++ ) {
           a = *pk; c = *ck; MOD128(a,c,md); *pk = a; *ck = 0;
         }
         s = sugar[l];
         for ( i = 0; i < l; i++ ) {
           a = mat[i][j]; c = cmat[i][j]; MOD128(a,c,md); mat[i][j] = a; cmat[i][j] = 0;
           if ( a ) {
             sugar[i] = MAX(sugar[i],s);
             red_by_vect64(md,mat[i]+j,cmat[i]+j,mat[l]+j,(int)(md-a),col-j);
         }          }
     return rank;        }
         l--;
       }
     for ( i = 0; i < row; i++ ) GCFREE(cmat[i]);
     GCFREE(cmat);
     return rank;
 }  }
 #endif  #endif
   
Legend:
Removed from v.1.239  
changed lines
  Added in v.1.240
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>