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

Diff for /OpenXM_contrib2/asir2000/builtin/array.c between version 1.65 and 1.75

version 1.65, 2013/12/20 02:02:23 version 1.75, 2017/09/17 02:34:02
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/array.c,v 1.64 2013/11/05 02:55:02 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.74 2017/09/15 01:52:51 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
Line 68 
Line 68 
 extern int DP_Print; /* XXX */  extern int DP_Print; /* XXX */
   
   
 void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm();  void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(), Ptriangleq();
 void Pinvmat();  void Pinvmat();
 void Pnewbytearray(),Pmemoryplot_to_coord();  void Pnewbytearray(),Pmemoryplot_to_coord();
   
 void Pgeneric_gauss_elim();  void Pgeneric_gauss_elim();
 void Pgeneric_gauss_elim_mod();  void Pgeneric_gauss_elim_mod();
   
   void Pindep_rows_mod();
   
 void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();  void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();
 void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov();  void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov();
 void Pgeninv_sf_swap();  void Pgeninv_sf_swap();
Line 108  struct ftab array_tab[] = {
Line 110  struct ftab array_tab[] = {
         {"mat_to_gfmmat",Pmat_to_gfmmat,2},          {"mat_to_gfmmat",Pmat_to_gfmmat,2},
         {"generic_gauss_elim",Pgeneric_gauss_elim,1},          {"generic_gauss_elim",Pgeneric_gauss_elim,1},
         {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},          {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
           {"indep_rows_mod",Pindep_rows_mod,2},
         {"newvect",Pnewvect,-2},          {"newvect",Pnewvect,-2},
         {"vect",Pvect,-99999999},          {"vect",Pvect,-99999999},
         {"vector",Pnewvect,-2},          {"vector",Pnewvect,-2},
Line 148  struct ftab array_tab[] = {
Line 151  struct ftab array_tab[] = {
         {"mat_col",Pmat_col,2},          {"mat_col",Pmat_col,2},
         {"lusolve_prep",Plusolve_prep,1},          {"lusolve_prep",Plusolve_prep,1},
         {"lusolve_main",Plusolve_main,1},          {"lusolve_main",Plusolve_main,1},
           {"triangleq",Ptriangleq,1},
         {0,0,0},          {0,0,0},
 };  };
   
Line 155  typedef struct _ent { int j; unsigned int e; } ent;
Line 159  typedef struct _ent { int j; unsigned int e; } ent;
   
 ent *get_row(FILE *,int *l);  ent *get_row(FILE *,int *l);
 void put_row(FILE *out,int l,ent *a);  void put_row(FILE *out,int l,ent *a);
 int lu_elim(int *l,ent **a,int k,int i,int mul,int mod);  void lu_elim(int *l,ent **a,int k,int i,int mul,int mod);
   void lu_append(int *,ent **,int *,int,int,int);
   void solve_l(int *,ent **,int,int *,int);
   void solve_u(int *,ent **,int,int *,int);
   
   
 static int *ul,*ll;  static int *ul,*ll;
 static ent **u,**l;  static ent **u,**l;
 static int modulus;  static int modulus;
Line 221  ent *get_row(FILE *in,int *l)
Line 229  ent *get_row(FILE *in,int *l)
         return a;          return a;
 }  }
   
 int lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod)  void lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod)
 {  {
         int i,j,k,s,mul;          int i,j,k,s,mul;
         unsigned int inv;          unsigned int inv;
Line 244  int lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int
Line 252  int lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int
   
 #define INITLEN 10  #define INITLEN 10
   
 lu_append(int *l,ent **a,int *l2,int k,int i,int mul)  void lu_append(int *l,ent **a,int *l2,int k,int i,int mul)
 {  {
         int len;          int len;
         ent *p;          ent *p;
Line 267  lu_append(int *l,ent **a,int *l2,int k,int i,int mul)
Line 275  lu_append(int *l,ent **a,int *l2,int k,int i,int mul)
   
 /* a[k] = a[k]-mul*a[i] */  /* a[k] = a[k]-mul*a[i] */
   
 int lu_elim(int *l,ent **a,int k,int i,int mul,int mod)  void lu_elim(int *l,ent **a,int k,int i,int mul,int mod)
 {  {
         ent *ak,*ai,*w;          ent *ak,*ai,*w;
         int lk,li,j,m,p,q,r,s,t,j0;          int lk,li,j,m,p,q,r,s,t,j0;
Line 315  int lu_elim(int *l,ent **a,int k,int i,int mul,int mod
Line 323  int lu_elim(int *l,ent **a,int k,int i,int mul,int mod
         l[k] = j;          l[k] = j;
 }  }
   
 int solve_l(int *ll,ent **l,int n,int *rhs,int mod)  void solve_l(int *ll,ent **l,int n,int *rhs,int mod)
 {  {
         int j,k,s,len;          int j,k,s,len;
         ent *p;          ent *p;
Line 329  int solve_l(int *ll,ent **l,int n,int *rhs,int mod)
Line 337  int solve_l(int *ll,ent **l,int n,int *rhs,int mod)
         }          }
 }  }
   
 int solve_u(int *ul,ent **u,int n,int *rhs,int mod)  void solve_u(int *ul,ent **u,int n,int *rhs,int mod)
 {  {
         int j,k,s,len,inv;          int j,k,s,len,inv;
         ent *p;          ent *p;
Line 895  void Pvtol(NODE arg,LIST *rp)
Line 903  void Pvtol(NODE arg,LIST *rp)
         pointer *a;          pointer *a;
         int len,i;          int len,i;
   
           if ( OID(ARG0(arg)) == O_LIST ) {
                   *rp = ARG0(arg);
                   return;
           }
         asir_assert(ARG0(arg),O_VECT,"vtol");          asir_assert(ARG0(arg),O_VECT,"vtol");
         v = (VECT)ARG0(arg); len = v->len; a = BDY(v);          v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
         for ( i = len - 1, n = 0; i >= 0; i-- ) {          for ( i = len - 1, n = 0; i >= 0; i-- ) {
Line 906  void Pvtol(NODE arg,LIST *rp)
Line 918  void Pvtol(NODE arg,LIST *rp)
 void Pltov(NODE arg,VECT *rp)  void Pltov(NODE arg,VECT *rp)
 {  {
         NODE n;          NODE n;
         VECT v;          VECT v,v0;
         int len,i;          int len,i;
   
           if ( OID(ARG0(arg)) == O_VECT ) {
                   v0 = (VECT)ARG0(arg); len = v0->len;
                   MKVECT(v,len);
                   for ( i = 0; i < len; i++ ) {
                           BDY(v)[i] = BDY(v0)[i];
                   }
                   *rp = v;
                   return;
           }
         asir_assert(ARG0(arg),O_LIST,"ltov");          asir_assert(ARG0(arg),O_LIST,"ltov");
         n = (NODE)BDY((LIST)ARG0(arg));          n = (NODE)BDY((LIST)ARG0(arg));
         len = length(n);          len = length(n);
Line 1159  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
Line 1180  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
         MKLIST(*rp,n0);          MKLIST(*rp,n0);
 }  }
   
   void Pindep_rows_mod(NODE arg,VECT *rp)
   {
           MAT m,mat;
           VECT rind;
           Q **tmat;
           int **wmat,**row0;
           Q *rib;
           int *rowstat,*p;
           Q q;
           int md,i,j,k,l,row,col,t,rank;
   
           asir_assert(ARG0(arg),O_MAT,"indep_rows_mod");
           asir_assert(ARG1(arg),O_N,"indep_rows_mod");
           m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
           row = m->row; col = m->col; tmat = (Q **)m->body;
           wmat = (int **)almat(row,col);
   
           row0 = (int **)ALLOCA(row*sizeof(int *));
           for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
   
           rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
           for ( i = 0; i < row; i++ )
                   for ( j = 0; j < col; j++ )
                           if ( q = (Q)tmat[i][j] ) {
                                   t = rem(NM(q),md);
                                   if ( t && SGN(q) < 0 )
                                           t = (md - t) % md;
                                   wmat[i][j] = t;
                           } else
                                   wmat[i][j] = 0;
           rank = indep_rows_mod(wmat,row,col,md,rowstat);
   
           MKVECT(rind,rank);
           rib = (Q *)rind->body;
           for ( j = 0; j < rank; j++ ) {
                   STOQ(rowstat[j],rib[j]);
           }
       *rp = rind;
   }
   
 /*  /*
         input : a row x col matrix A          input : a row x col matrix A
                 A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...                  A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
Line 2295  void red_by_vect(int m,unsigned int *p,unsigned int *r
Line 2356  void red_by_vect(int m,unsigned int *p,unsigned int *r
                 }                  }
 }  }
   
   #if defined(__GNUC__) && SIZEOF_LONG==8
   /* 64bit vector += UNIT vector(normalized) */
   
   void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *r,unsigned int hc,int len)
   {
     U64 t;
   
     /* (p[0],c[0]) is normalized */
     *p++ = 0; *c++ = 0; r++; len--;
     for ( ; len; len--, r++, p++, c++ )
       if ( *r ) {
         t = (*p)+(*r)*hc;
         if ( t < *p ) (*c)++;
         *p = t;
       }
   }
   #endif
   
 void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)  void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
 {  {
         *p++ = 0; r++; len--;          *p++ = 0; r++; len--;
Line 2303  void red_by_vect_sf(int m,unsigned int *p,unsigned int
Line 2382  void red_by_vect_sf(int m,unsigned int *p,unsigned int
                         *p = _addsf(_mulsf(*r,hc),*p);                          *p = _addsf(_mulsf(*r,hc),*p);
 }  }
   
   extern GZ current_mod_lf;
   extern int current_mod_lf_size;
   
   void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len)
   {
           mpz_set_ui(*p++,0); r++; len--;
           for ( ; len; len--, r++, p++ ) {
          mpz_addmul(*p,*r,hc);
   #if 0
          if ( mpz_size(*p) > current_mod_lf_size )
            mpz_mod(*p,*p,BDY(current_mod_lf));
   #endif
       }
   }
   
   
 extern unsigned int **psca;  extern unsigned int **psca;
   
 void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,  void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
Line 2451  int generic_gauss_elim_mod2(int **mat0,int row,int col
Line 2546  int generic_gauss_elim_mod2(int **mat0,int row,int col
         return rank;          return rank;
 }  }
   
   int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat)
   {
           int i,j,k,l,inv,a,rank;
           unsigned int *t,*pivot,*pk;
           unsigned int **mat;
   
           for ( i = 0; i < row; i++ ) rowstat[i] = i;
           mat = (unsigned int **)mat0;
           for ( rank = 0, j = 0; j < col; j++ ) {
                   for ( i = rank; i < row; i++ )
                           mat[i][j] %= md;
                   for ( i = rank; i < row; i++ )
                           if ( mat[i][j] )
                                   break;
                   if ( i == row ) continue;
                   if ( i != rank ) {
                           t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                           k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
                   }
                   pivot = mat[rank];
                   inv = invm(pivot[j],md);
                   for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                           if ( *pk ) {
                                   if ( *pk >= (unsigned int)md )
                                           *pk %= md;
                                   DMAR(*pk,inv,0,md,*pk)
                           }
                   for ( i = rank+1; i < row; i++ ) {
                           t = mat[i];
                           if ( a = t[j] )
                                   red_by_vect(md,t+j,pivot+j,md-a,col-j);
                   }
                   rank++;
           }
           return rank;
   }
   
 int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)  int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)
 {  {
         int i,j,k,l,inv,a,rank;          int i,j,k,l,inv,a,rank;
Line 2977  void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
Line 3109  void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
         TOGFMMAT(row,col,wmat,*rp);          TOGFMMAT(row,col,wmat,*rp);
 }  }
   
 void Pgeninvm_swap(arg,rp)  void Pgeninvm_swap(NODE arg,LIST *rp)
 NODE arg;  
 LIST *rp;  
 {  {
         MAT m;          MAT m;
         pointer **mat;          pointer **mat;
Line 3025  LIST *rp;
Line 3155  LIST *rp;
         }          }
 }  }
   
 gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp)  int gauss_elim_geninv_mod_swap(unsigned int **mat,int row,int col,unsigned int md,
 unsigned int **mat;      unsigned int ***invmatp,int **indexp)
 int row,col;  
 unsigned int md;  
 unsigned int ***invmatp;  
 int **indexp;  
 {  {
         int i,j,k,inv,a,n,m;          int i,j,k,inv,a,n,m;
         unsigned int *t,*pivot,*s;          unsigned int *t,*pivot,*s;
Line 3737  void Pmat_col(NODE arg,VECT *rp)
Line 3863  void Pmat_col(NODE arg,VECT *rp)
                 BDY(vect)[i] = BDY(mat)[i][j];                  BDY(vect)[i] = BDY(mat)[i][j];
         }          }
         *rp = vect;          *rp = vect;
   }
   
   NODE triangleq(NODE e)
   {
     int n,i,k;
     V v;
     VL vl;
     P *p;
     NODE r,r1;
   
     n = length(e);
     p = (P *)MALLOC(n*sizeof(P));
     for ( i = 0; i < n; i++, e = NEXT(e) ) p[i] = (P)BDY(e);
     i = 0;
     while ( 1 ) {
       for ( ; i < n && !p[i]; i++ );
       if ( i == n ) break;
       if ( OID(p[i]) == O_N ) return 0;
       v = p[i]->v;
       for ( k = i+1; k < n; k++ )
         if ( p[k] ) {
           if ( OID(p[k]) == O_N ) return 0;
           if ( p[k]->v == v ) p[k] = 0;
         }
       i++;
     }
     for ( r = 0, i = 0; i < n; i++ ) {
       if ( p[i] ) {
         MKNODE(r1,p[i],r); r = r1;
       }
     }
     return r;
   }
   
   void Ptriangleq(NODE arg,LIST *rp)
   {
     NODE ret;
   
     asir_assert(ARG0(arg),O_LIST,"sparseleq");
     ret = triangleq(BDY((LIST)ARG0(arg)));
     MKLIST(*rp,ret);
 }  }
Legend:
Removed from v.1.65  
changed lines
  Added in v.1.75
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>