[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.61 and 1.70

version 1.61, 2012/12/17 07:20:44 version 1.70, 2017/01/08 03:05:39
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.60 2010/11/09 16:23:45 ohara Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.69 2015/09/03 23:05:35 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
Line 75  void Pnewbytearray(),Pmemoryplot_to_coord();
Line 75  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 98  void Pmatc();
Line 100  void Pmatc();
 void Pnd_det();  void Pnd_det();
 void Plu_mat();  void Plu_mat();
 void Pmat_col();  void Pmat_col();
   void Plusolve_prep();
   void Plusolve_main();
   
 struct ftab array_tab[] = {  struct ftab array_tab[] = {
         {"lu_mat",Plu_mat,1},          {"lu_mat",Plu_mat,1},
Line 106  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 144  struct ftab array_tab[] = {
Line 149  struct ftab array_tab[] = {
         {"mat_swap_row_destructive",Pmat_swap_row_destructive,3},          {"mat_swap_row_destructive",Pmat_swap_row_destructive,3},
         {"mat_swap_col_destructive",Pmat_swap_col_destructive,3},          {"mat_swap_col_destructive",Pmat_swap_col_destructive,3},
         {"mat_col",Pmat_col,2},          {"mat_col",Pmat_col,2},
           {"lusolve_prep",Plusolve_prep,1},
           {"lusolve_main",Plusolve_main,1},
         {0,0,0},          {0,0,0},
 };  };
   
   typedef struct _ent { int j; unsigned int e; } ent;
   
   ent *get_row(FILE *,int *l);
   void put_row(FILE *out,int l,ent *a);
   int lu_elim(int *l,ent **a,int k,int i,int mul,int mod);
   
   static int *ul,*ll;
   static ent **u,**l;
   static int modulus;
   
   void Plusolve_prep(NODE arg,Q *rp)
   {
           char *fname;
           FILE *in;
           int len,i,rank;
           int *rhs;
   
           fname = BDY((STRING)ARG0(arg));
           in = fopen(fname,"r");
           modulus = getw(in);
           len = getw(in);
           ul = (int *)MALLOC_ATOMIC(len*sizeof(int));
           u = (ent **)MALLOC(len*sizeof(ent *));
           ll = (int *)MALLOC_ATOMIC(len*sizeof(int));
           l = (ent **)MALLOC(len*sizeof(ent *));
           for ( i = 0; i < len; i++ ) {
                   u[i] = get_row(in,&ul[i]);
           }
           for ( i = 0; i < len; i++ ) {
                   l[i] = get_row(in,&ll[i]);
           }
           fclose(in);
           *rp = ONE;
   }
   
   void Plusolve_main(NODE arg,VECT *rp)
   {
           Q *d,*p;
           VECT v,r;
           int len,i;
           int *rhs;
   
           v = (VECT)ARG0(arg); len = v->len;
           d = (Q *)BDY(v);
           rhs = (int *)MALLOC_ATOMIC(len*sizeof(int));
           for ( i = 0; i < len; i++ ) rhs[i] = QTOS(d[i]);
           solve_l(ll,l,len,rhs,modulus);
           solve_u(ul,u,len,rhs,modulus);
           NEWVECT(r); r->len = len;
           r->body = (pointer *)MALLOC(len*sizeof(pointer));
           p = (Q *)r->body;
           for ( i = 0; i < len; i++ )
                   STOQ(rhs[i],p[i]);
           *rp = r;
   }
   
   ent *get_row(FILE *in,int *l)
   {
           int len,i;
           ent *a;
   
           *l = len = getw(in);
           a = (ent *)MALLOC_ATOMIC(len*sizeof(ent));
           for ( i = 0; i < len; i++ ) {
                   a[i].j = getw(in);
                   a[i].e = getw(in);
           }
           return a;
   }
   
   int lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod)
   {
           int i,j,k,s,mul;
           unsigned int inv;
           int *ll2;
   
           ll2 = (int *)MALLOC_ATOMIC(n*sizeof(int));
           for ( i = 0; i < n; i++ ) ll2[i] = 0;
           for ( i = 0; i < n; i++ ) {
                   fprintf(stderr,"i=%d\n",i);
                   inv = invm(u[i][0].e,mod);
                   for ( k = i+1; k < n; k++ )
                           if ( u[k][0].j == n-i ) {
                                   s = u[k][0].e;
                                   DMAR(s,inv,0,mod,mul);
                                   lu_elim(ul,u,k,i,mul,mod);
                                   lu_append(ll,l,ll2,k,i,mul);
                           }
           }
   }
   
   #define INITLEN 10
   
   lu_append(int *l,ent **a,int *l2,int k,int i,int mul)
   {
           int len;
           ent *p;
   
           len = l[k];
           if ( !len ) {
                   a[k] = p = (ent *)MALLOC_ATOMIC(INITLEN*sizeof(ent));
                   p[0].j = i; p[0].e = mul;
                   l[k] = 1; l2[k] = INITLEN;
           } else {
                   if ( l2[k] == l[k] ) {
                           l2[k] *= 2;
                           a[k] = REALLOC(a[k],l2[k]*sizeof(ent));
                   }
                   p =a[k];
                   p[l[k]].j = i; p[l[k]].e = mul;
                   l[k]++;
           }
   }
   
   /* a[k] = a[k]-mul*a[i] */
   
   int lu_elim(int *l,ent **a,int k,int i,int mul,int mod)
   {
           ent *ak,*ai,*w;
           int lk,li,j,m,p,q,r,s,t,j0;
   
           ak = a[k]; ai = a[i]; lk = l[k]; li = l[i];
           w = (ent *)alloca((lk+li)*sizeof(ent));
           p = 0; q = 0; j = 0;
           mul = mod-mul;
           while ( p < lk && q < li ) {
                   if ( ak[p].j > ai[q].j ) {
                           w[j] = ak[p]; j++; p++;
                   } else if ( ak[p].j < ai[q].j ) {
                           w[j].j = ai[q].j;
                           t = ai[q].e;
                           DMAR(t,mul,0,mod,r);
                           w[j].e = r;
                           j++; q++;
                   } else {
                           t = ai[q].e; s = ak[p].e;
                           DMAR(t,mul,s,mod,r);
                           if ( r ) {
                                   w[j].j = ai[q].j; w[j].e = r; j++;
                           }
                           p++; q++;
                   }
           }
           if ( q == li )
                   while ( p < lk ) {
                           w[j] = ak[p]; j++; p++;
                   }
           else if ( p == lk )
                   while ( q < li ) {
                           w[j].j = ai[q].j;
                           t = ai[q].e;
                           DMAR(t,mul,0,mod,r);
                           w[j].e = r;
                           j++; q++;
                   }
           if ( j <= lk ) {
                   for ( m = 0; m < j; m++ ) ak[m] = w[m];
           } else {
                   a[k] = ak = (ent *)MALLOC_ATOMIC(j*sizeof(ent));
                   for ( m = 0; m < j; m++ ) ak[m] = w[m];
           }
           l[k] = j;
   }
   
   int solve_l(int *ll,ent **l,int n,int *rhs,int mod)
   {
           int j,k,s,len;
           ent *p;
   
           for ( j = 0; j < n; j++ ) {
                   len = ll[j]; p = l[j];
                   for ( k = 0, s = 0; k < len; k++ )
                           s = dmar(p[k].e,rhs[p[k].j],s,mod);
                   rhs[j] -=  s;
                   if ( rhs[j] < 0 ) rhs[j] += mod;
           }
   }
   
   int solve_u(int *ul,ent **u,int n,int *rhs,int mod)
   {
           int j,k,s,len,inv;
           ent *p;
   
           for ( j = n-1; j >= 0; j-- ) {
                   len = ul[j]; p = u[j];
                   for ( k = 1, s = 0; k < len; k++ )
                           s = dmar(p[k].e,rhs[p[k].j],s,mod);
                   rhs[j] -=  s;
                   if ( rhs[j] < 0 ) rhs[j] += mod;
                   inv = invm((unsigned int)p[0].e,mod);
                   rhs[j] = dmar(rhs[j],inv,0,mod);
           }
   }
   
 int comp_obj(Obj *a,Obj *b)  int comp_obj(Obj *a,Obj *b)
 {  {
         return arf_comp(CO,*a,*b);          return arf_comp(CO,*a,*b);
Line 411  void Pnewvect(NODE arg,VECT *rp)
Line 612  void Pnewvect(NODE arg,VECT *rp)
 }  }
   
 void Pvect(NODE arg,VECT *rp) {  void Pvect(NODE arg,VECT *rp) {
         int len,i,r;          int len,i;
         VECT vect;          VECT vect;
         pointer *vb;          pointer *vb;
         NODE tn;          NODE tn;
Line 697  void Pvtol(NODE arg,LIST *rp)
Line 898  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 708  void Pvtol(NODE arg,LIST *rp)
Line 913  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 923  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
Line 1137  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
         int *ri,*ci;          int *ri,*ci;
         VECT rind,cind;          VECT rind,cind;
         Q dn,q;          Q dn,q;
         int i,j,k,l,row,col,t,rank;          int i,row,col,t,rank;
         int is_hensel = 0;          int is_hensel = 0;
         char *key;          char *key;
         Obj value;          Obj value;
Line 961  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
Line 1175  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 1266  RESET:
Line 1520  RESET:
         }          }
 }  }
   
   void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm);
   
 /* XXX broken */  /* XXX broken */
 int lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)  void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
 {  {
         Q **a0,**b;          Q **a0,**b;
         Q *aiq;          Q *aiq;
Line 1369  int lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
Line 1625  int lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
         }          }
 }  }
   
 int nmat(N **m,int n)  void nmat(N **m,int n)
 {  {
         int i,j;          int i,j;
   
Line 2103  void red_by_vect_sf(int m,unsigned int *p,unsigned int
Line 2359  void red_by_vect_sf(int m,unsigned int *p,unsigned int
                         *p = _addsf(_mulsf(*r,hc),*p);                          *p = _addsf(_mulsf(*r,hc),*p);
 }  }
   
   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);
   }
   
   
 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 2191  int generic_gauss_elim_mod(int **mat0,int row,int col,
Line 2455  int generic_gauss_elim_mod(int **mat0,int row,int col,
         return rank;          return rank;
 }  }
   
   int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,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 ) {
                           colstat[j] = 0;
                           continue;
                   } else
                           colstat[j] = 1;
                   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++;
           }
           for ( j = col-1, l = rank-1; j >= 0; j-- )
                   if ( colstat[j] ) {
                           pivot = mat[l];
                           for ( i = 0; i < l; i++ ) {
                                   t = mat[i];
                                   t[j] %= md;
                                   if ( a = t[j] )
                                           red_by_vect(md,t+j,pivot+j,md-a,col-j);
                           }
                           l--;
                   }
           for ( j = 0, l = 0; l < rank; j++ )
                   if ( colstat[j] ) {
                           t = mat[l];
                           for ( k = j; k < col; k++ )
                                   if ( t[k] >= (unsigned int)md )
                                           t[k] %= md;
                           l++;
                   }
           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 3458  void Pnd_det(NODE arg,P *rp)
Line 3819  void Pnd_det(NODE arg,P *rp)
                 nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);                  nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);
 }  }
   
 void Pmat_col(NODE arg,P *rp)  void Pmat_col(NODE arg,VECT *rp)
 {  {
         int i,j,n;          int i,j,n;
         pointer t;  
         MAT mat;          MAT mat;
         VECT vect;          VECT vect;
   
Legend:
Removed from v.1.61  
changed lines
  Added in v.1.70
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>