[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.49 and 1.62

version 1.49, 2005/12/21 23:18:15 version 1.62, 2013/06/14 05:55:24
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.48 2005/11/27 05:37:53 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.61 2012/12/17 07:20:44 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
 #include "parse.h"  #include "parse.h"
 #include "inline.h"  #include "inline.h"
   
   #include <sys/types.h>
   #include <sys/stat.h>
   #if !defined(_MSC_VER)
   #include <unistd.h>
   #endif
   
 #define F4_INTRAT_PERIOD 8  #define F4_INTRAT_PERIOD 8
   
 #if 0  #if 0
Line 90  void Pvect();
Line 96  void Pvect();
 void Pmat();  void Pmat();
 void Pmatc();  void Pmatc();
 void Pnd_det();  void Pnd_det();
   void Plu_mat();
   void Pmat_col();
   
 struct ftab array_tab[] = {  struct ftab array_tab[] = {
           {"lu_mat",Plu_mat,1},
         {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},          {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},
         {"lu_gfmmat",Plu_gfmmat,2},          {"lu_gfmmat",Plu_gfmmat,2},
         {"mat_to_gfmmat",Pmat_to_gfmmat,2},          {"mat_to_gfmmat",Pmat_to_gfmmat,2},
Line 134  struct ftab array_tab[] = {
Line 143  struct ftab array_tab[] = {
         {"nbpoly_up2",Pnbpoly_up2,2},          {"nbpoly_up2",Pnbpoly_up2,2},
         {"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},
         {0,0,0},          {0,0,0},
 };  };
   
Line 144  int comp_obj(Obj *a,Obj *b)
Line 154  int comp_obj(Obj *a,Obj *b)
   
 static FUNC generic_comp_obj_func;  static FUNC generic_comp_obj_func;
 static NODE generic_comp_obj_arg;  static NODE generic_comp_obj_arg;
   static NODE generic_comp_obj_option;
   
 int generic_comp_obj(Obj *a,Obj *b)  int generic_comp_obj(Obj *a,Obj *b)
 {  {
Line 151  int generic_comp_obj(Obj *a,Obj *b)
Line 162  int generic_comp_obj(Obj *a,Obj *b)
   
         BDY(generic_comp_obj_arg)=(pointer)(*a);          BDY(generic_comp_obj_arg)=(pointer)(*a);
         BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);          BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);
         r = (Q)bevalf(generic_comp_obj_func,generic_comp_obj_arg);          r = (Q)bevalf_with_opts(generic_comp_obj_func,generic_comp_obj_arg,generic_comp_obj_option);
         if ( !r )          if ( !r )
                 return 0;                  return 0;
         else          else
Line 198  void Pqsort(NODE arg,LIST *rp)
Line 209  void Pqsort(NODE arg,LIST *rp)
                         func = (FUNC)v->priv;                          func = (FUNC)v->priv;
                 }                  }
                 generic_comp_obj_func = func;                  generic_comp_obj_func = func;
                 MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);                  MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);
                   generic_comp_obj_option = current_option;
                 qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);                  qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);
         }          }
     if (OID(t) == O_LIST) {      if (OID(t) == O_LIST) {
Line 332  void Psepmat_destructive(NODE arg,LIST *rp)
Line 344  void Psepmat_destructive(NODE arg,LIST *rp)
                         sgn = SGN(ent);                          sgn = SGN(ent);
                         divn(nm,mod,&quo,&rem);                          divn(nm,mod,&quo,&rem);
 /*                      if ( quo != nm && rem != nm ) */  /*                      if ( quo != nm && rem != nm ) */
 /*                              GC_free(nm); */  /*                              GCFREE(nm); */
 /*                      GC_free(ent); */  /*                      GCFREE(ent); */
                         NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]);                          NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]);
                 }                  }
         MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);          MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);
Line 385  void Pnewvect(NODE arg,VECT *rp)
Line 397  void Pnewvect(NODE arg,VECT *rp)
         if ( argc(arg) == 2 ) {          if ( argc(arg) == 2 ) {
                 list = (LIST)ARG1(arg);                  list = (LIST)ARG1(arg);
                 asir_assert(list,O_LIST,"newvect");                  asir_assert(list,O_LIST,"newvect");
   #if 0
                 for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );                  for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
                 if ( r > len ) {                  if ( r > len ) {
                         *rp = vect;                          *rp = vect;
                         return;                          return;
                 }                  }
   #endif
                 for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )                  for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )
                         vb[i] = (pointer)BDY(tn);                          vb[i] = (pointer)BDY(tn);
         }          }
Line 397  void Pnewvect(NODE arg,VECT *rp)
Line 411  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 444  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
Line 458  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
         char *str;          char *str;
         LIST list;          LIST list;
         NODE tn;          NODE tn;
           int ac;
           struct stat sbuf;
           char *fname;
           FILE *fp;
   
         asir_assert(ARG0(arg),O_N,"newbytearray");          ac = argc(arg);
         len = QTOS((Q)ARG0(arg));          if ( ac == 1 ) {
         if ( len < 0 )                  if ( !OID((Obj)ARG0(arg)) ) error("newbytearray : invalid argument");
                 error("newbytearray : invalid size");                  switch ( OID((Obj)ARG0(arg)) ) {
         MKBYTEARRAY(array,len);                          case O_STR:
         if ( argc(arg) == 2 ) {                                  fname = BDY((STRING)ARG0(arg));
                                   fp = fopen(fname,"rb");
                                   if ( !fp ) error("newbytearray : fopen failed");
                                   if ( stat(fname,&sbuf) < 0 )
                                           error("newbytearray : stat failed");
                                   len = sbuf.st_size;
                                   MKBYTEARRAY(array,len);
                                   fread(BDY(array),len,sizeof(char),fp);
                                   break;
                           case O_N:
                                   if ( !RATN(ARG0(arg)) )
                                           error("newbytearray : invalid argument");
                                   len = QTOS((Q)ARG0(arg));
                                   if ( len < 0 )
                                           error("newbytearray : invalid size");
                                   MKBYTEARRAY(array,len);
                                   break;
                           default:
                                   error("newbytearray : invalid argument");
                   }
           } else if ( ac == 2 ) {
                   asir_assert(ARG0(arg),O_N,"newbytearray");
                   len = QTOS((Q)ARG0(arg));
                   if ( len < 0 )
                           error("newbytearray : invalid size");
                   MKBYTEARRAY(array,len);
                 if ( !ARG1(arg) )                  if ( !ARG1(arg) )
                         error("newbytearray : invalid initialization");                          error("newbytearray : invalid initialization");
                 switch ( OID((Obj)ARG1(arg)) ) {                  switch ( OID((Obj)ARG1(arg)) ) {
Line 474  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
Line 517  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
                                 if ( !ARG1(arg) )                                  if ( !ARG1(arg) )
                                         error("newbytearray : invalid initialization");                                          error("newbytearray : invalid initialization");
                 }                  }
         }          } else
                   error("newbytearray : invalid argument");
         *rp = array;          *rp = array;
 }  }
   
Line 483  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
Line 527  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
 void Pmemoryplot_to_coord(NODE arg,LIST *rp)  void Pmemoryplot_to_coord(NODE arg,LIST *rp)
 {  {
         int len,blen,y,i,j;          int len,blen,y,i,j;
         char *a;          unsigned char *a;
         NODE r0,r,n;          NODE r0,r,n;
         LIST l;          LIST l;
         BYTEARRAY ba;          BYTEARRAY ba;
Line 879  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
Line 923  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 1222  RESET:
Line 1266  RESET:
         }          }
 }  }
   
   /* XXX broken */
   int lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
   {
           Q **a0,**b;
           Q *aiq;
           N **a;
           N *ai;
           Q q,q1,dn2,a1,q0,bik;
           MAT m;
           unsigned int md;
           int n,ind,i,j,rank,t,inv,t1,ret,min,k;
           int **w;
           int *wi,*rinfo0,*rinfo;
           N m1,m2,m3,u,s;
   
           a0 = (Q **)mat->body;
           n = mat->row;
           if ( n != mat->col )
                   error("lu_dec_cr : non-square matrix");
           w = (int **)almat(n,n);
           MKMAT(m,n,n);
           a = (N **)m->body;
           UTON(1,m1);
           rinfo0 = 0;
           ind = 0;
           while ( 1 ) {
                   md = get_lprime(ind);
                   /* mat mod md */
                   for ( i = 0; i < n; i++ )
                           for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ )
                                   if ( q = aiq[j] ) {
                                           t = rem(NM(q),md);
                                           if ( t && SGN(q) < 0 )
                                                   t = (md - t) % md;
                                           wi[j] = t;
                                   } else
                                           wi[j] = 0;
   
                   if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue;
                   printf("."); fflush(stdout);
                   if ( !rinfo0 )
                           *perm = rinfo0 = rinfo;
                   else {
                           for ( i = 0; i < n; i++ )
                                   if ( rinfo[i] != rinfo0[i] ) break;
                           if ( i < n ) continue;
                   }
                   if ( UNIN(m1) ) {
                           for ( i = 0; i < n; i++ )
                                   for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) {
                                           UTON(wi[j],u); ai[j] = u;
                                   }
                           UTON(md,m1);
                   } else {
                           inv = invm(rem(m1,md),md);
                           UTON(md,m2); muln(m1,m2,&m3);
                           for ( i = 0; i < n; i++ )
                                   for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ )
                                           if ( ai[i] ) {
                                           /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                                                   t = rem(ai[j],md);
                                                   if ( wi[j] >= t )
                                                           t = wi[j]-t;
                                                   else
                                                           t = md-(t-wi[j]);
                                                   DMAR(t,inv,0,md,t1)
                                                   UTON(t1,u);
                                                   muln(m1,u,&s);
                                                   addn(ai[j],s,&u); ai[j] = u;
                                           } else if ( wi[j] ) {
                                                   /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                                                   DMAR(wi[j],inv,0,md,t)
                                                   UTON(t,u);
                                                   muln(m1,u,&s); ai[j] = s;
                                           }
                           m1 = m3;
                   }
                   if ( (++ind%8) == 0 ) {
                           ret = intmtoratm(m,m1,lu,dn);
                           if ( ret ) {
                                   b = (Q **)lu->body;
                                   mulq(*dn,*dn,&dn2);
                                   for ( i = 0; i < n; i++ ) {
                                           for ( j = 0; j < n; j++ ) {
                                                   q = 0;
                                                   min = MIN(i,j);
                                                   for ( k = 0; k <= min; k++ ) {
                                                           bik = k==i ? *dn : b[i][k];
                                                           mulq(bik,b[k][j],&q0);
                                                           addq(q,q0,&q1); q = q1;
                                                   }
                                                   mulq(a0[rinfo0[i]][j],dn2,&q1);
                                                   if ( cmpq(q,q1) ) break;
                                           }
                                           if ( j < n ) break;
                                   }
                                   if ( i == n )
                                           return;
                           }
                   }
           }
   }
   
   int nmat(N **m,int n)
   {
           int i,j;
   
           for ( i = 0; i < n; i++ ) {
                   for ( j = 0; j < n; j++ ) {
                           printn(m[i][j]); printf(" ");
                   }
                   printf("\n");
           }
   }
   
 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp)  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp)
 {  {
         MAT bmat,xmat;          MAT bmat,xmat;
Line 1262  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
Line 1421  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
                                 } else                                  } else
                                         wi[j] = 0;                                          wi[j] = 0;
   
                 if ( DP_Print ) {                  if ( DP_Print > 3 ) {
                         fprintf(asir_out,"LU decomposition.."); fflush(asir_out);                          fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
                 }                  }
                 rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);                  rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
                 if ( DP_Print ) {                  if ( DP_Print > 3 ) {
                         fprintf(asir_out,"done.\n"); fflush(asir_out);                          fprintf(asir_out,"done.\n"); fflush(asir_out);
                 }                  }
                 a = (Q **)almat_pointer(rank,rank); /* lhs mat */                  a = (Q **)almat_pointer(rank,rank); /* lhs mat */
Line 1305  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
Line 1464  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
                         wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));                          wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
                         for ( i = 0; i < wxsize; i++ ) wx[i] = 0;                          for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
                         for ( q = ONE, count = 0; ; ) {                          for ( q = ONE, count = 0; ; ) {
                                 if ( DP_Print )                                  if ( DP_Print > 3 )
                                         fprintf(stderr,"o");                                          fprintf(stderr,"o");
                                 /* wc = -b mod md */                                  /* wc = -b mod md */
                                 get_eg(&tmp0);                                  get_eg(&tmp0);
Line 1366  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
Line 1525  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
                                         ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);                                          ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);
                                         get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);                                          get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);
                                         if ( ret ) {                                          if ( ret ) {
                                                   rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                                                   cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                                                 for ( j = k = l = 0; j < col; j++ )                                                  for ( j = k = l = 0; j < col; j++ )
                                                         if ( cinfo[j] )                                                          if ( cinfo[j] )
                                                                 rind[k++] = j;                                                                  rind[k++] = j;
Line 1383  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
Line 1544  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
                                                                 print_eg("CHECK",&eg_check);                                                                  print_eg("CHECK",&eg_check);
                                                                 fflush(asir_out);                                                                  fflush(asir_out);
                                                         }                                                          }
                                                           *rindp = rind;
                                                           *cindp = cind;
                                                           for ( j = k = 0; j < col; j++ )
                                                                   if ( !cinfo[j] )
                                                                           cind[k++] = j;
                                                         return rank;                                                          return rank;
                                                 }                                                  }
                                         } else {                                          } else {
Line 1398  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
Line 1564  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
         }          }
 }  }
   
   int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT *nmmat,Q *dn,int **rindp,int **cindp)
   {
           MAT bmat,xmat;
           Q **a0,**a,**b,**x,**nm;
           Q *ai,*bi,*xi;
           int row,col;
           int **w;
           int *wi;
           int **wc;
           Q mdq,q,s,u;
           N tn;
           int ind,md,i,j,k,l,li,ri,rank;
           unsigned int t;
           int *cinfo,*rinfo;
           int *rind,*cind;
           int count;
           int ret;
           struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
           int period;
           int *wx,*ptr;
           int wxsize,nsize;
           N wn;
           Q wq;
           NumberField nf;
           DP m;
           int col1;
   
           a0 = (Q **)mat->body;
           row = mat->row; col = mat->col;
           w = (int **)almat(row,col);
           for ( ind = 0; ; ind++ ) {
                   md = get_lprime(ind);
                   STOQ(md,mdq);
                   for ( i = 0; i < row; i++ )
                           for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
                                   if ( q = (Q)ai[j] ) {
                                           t = rem(NM(q),md);
                                           if ( t && SGN(q) < 0 )
                                                   t = (md - t) % md;
                                           wi[j] = t;
                                   } else
                                           wi[j] = 0;
   
                   if ( DP_Print ) {
                           fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
                   }
                   rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
                   if ( DP_Print ) {
                           fprintf(asir_out,"done.\n"); fflush(asir_out);
                   }
                   for ( i = 0; i < col-1; i++ ) {
                           if ( !cinfo[i] ) {
                                   m = mb[i];
                                   for ( j = i+1; j < col-1; j++ )
                                           if ( dp_redble(mb[j],m) )
                                                   cinfo[j] = -1;
                           }
                   }
                   a = (Q **)almat_pointer(rank,rank); /* lhs mat */
                   MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
                   for ( j = li = ri = 0; j < col; j++ )
                           if ( cinfo[j] > 0 ) {
                                   /* the column is in lhs */
                                   for ( i = 0; i < rank; i++ ) {
                                           w[i][li] = w[i][j];
                                           a[i][li] = a0[rinfo[i]][j];
                                   }
                                   li++;
                           } else if ( !cinfo[j] ) {
                                   /* the column is in rhs */
                                   for ( i = 0; i < rank; i++ )
                                           b[i][ri] = a0[rinfo[i]][j];
                                   ri++;
                           }
   
                           /* solve Ax+B=0; A: rank x rank, B: rank x ri */
                           MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
                           MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
                           /* use the right part of w as work area */
                           wc = (int **)almat(rank,ri);
                           for ( i = 0; i < rank; i++ )
                                   wc[i] = w[i]+rank;
                           *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                           *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                           init_eg(&eg_mul); init_eg(&eg_inv);
                           init_eg(&eg_check); init_eg(&eg_intrat);
                           period = F4_INTRAT_PERIOD;
                           nsize = period;
                           wxsize = rank*ri*nsize;
                           wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
                           for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
                           for ( q = ONE, count = 0; ; ) {
                                   if ( DP_Print )
                                           fprintf(stderr,"o");
                                   /* wc = -b mod md */
                                   get_eg(&tmp0);
                                   for ( i = 0; i < rank; i++ )
                                           for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
                                                   if ( u = (Q)bi[j] ) {
                                                           t = rem(NM(u),md);
                                                           if ( t && SGN(u) > 0 )
                                                                   t = (md - t) % md;
                                                           wi[j] = t;
                                                   } else
                                                           wi[j] = 0;
                                   /* wc = A^(-1)wc; wc is not normalized */
                                   solve_by_lu_mod(w,rank,md,wc,ri,0);
                                   /* wx += q*wc */
                                   ptr = wx;
                                   for ( i = 0; i < rank; i++ )
                                           for ( j = 0, wi = wc[i]; j < ri; j++ ) {
                                                   if ( wi[j] )
                                                           muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);
                                                   ptr += nsize;
                                           }
                                   count++;
                                   get_eg(&tmp1);
                                   add_eg(&eg_inv,&tmp0,&tmp1);
                                   get_eg(&tmp0);
                                   for ( i = 0; i < rank; i++ )
                                           for ( j = 0; j < ri; j++ ) {
                                                   inner_product_mat_int_mod(a,wc,rank,i,j,&u);
                                                   addq(b[i][j],u,&s);
                                                   if ( s ) {
                                                           t = divin(NM(s),md,&tn);
                                                           if ( t )
                                                                   error("generic_gauss_elim_hensel:incosistent");
                                                           NTOQ(tn,SGN(s),b[i][j]);
                                                   } else
                                                           b[i][j] = 0;
                                           }
                                   get_eg(&tmp1);
                                   add_eg(&eg_mul,&tmp0,&tmp1);
                                   /* q = q*md */
                                   mulq(q,mdq,&u); q = u;
                                   if ( count == period ) {
                                           get_eg(&tmp0);
                                           ptr = wx;
                                           for ( i = 0; i < rank; i++ )
                                                   for ( j = 0, xi = x[i]; j < ri;
                                                           j++, ptr += nsize ) {
                                                           for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );
                                                           if ( k >= 0 ) {
                                                                   wn = NALLOC(k+1);
                                                                   PL(wn) = k+1;
                                                                   for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];
                                                                   NTOQ(wn,1,wq);
                                                                   subq(xi[j],wq,&u); xi[j] = u;
                                                           }
                                                   }
                                           ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);
                                           get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);
                                           if ( ret ) {
                                                   for ( j = k = l = 0; j < col; j++ )
                                                           if ( cinfo[j] > 0 )
                                                                   rind[k++] = j;
                                                           else if ( !cinfo[j] )
                                                                   cind[l++] = j;
                                                   get_eg(&tmp0);
                                                   ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
                                                   get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);
                                                   if ( ret ) {
                                                           if ( DP_Print > 3 ) {
                                                                   fprintf(stderr,"\n");
                                                                   print_eg("INV",&eg_inv);
                                                                   print_eg("MUL",&eg_mul);
                                                                   print_eg("INTRAT",&eg_intrat);
                                                                   print_eg("CHECK",&eg_check);
                                                                   fflush(asir_out);
                                                           }
                                                           return rank;
                                                   }
                                           } else {
                                                   period = period*3/2;
                                                   count = 0;
                                                   nsize += period;
                                                   wxsize += rank*ri*nsize;
                                                   wx = (int *)REALLOC(wx,wxsize*sizeof(int));
                                                   for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
                                           }
                                   }
                           }
           }
   }
   
 int f4_nocheck;  int f4_nocheck;
   
 int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind)  int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind)
Line 1731  void red_by_compress(int m,unsigned int *p,unsigned in
Line 2082  void red_by_compress(int m,unsigned int *p,unsigned in
   
 void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)  void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
 {  {
         register unsigned int up,lo;          unsigned int up,lo,dmy;
         unsigned int dmy;  
   
         *p++ = 0; r++; len--;          *p++ = 0; r++; len--;
         for ( ; len; len--, r++, p++ )          for ( ; len; len--, r++, p++ )
Line 1984  int find_lhs_and_lu_mod(unsigned int **a,int row,int c
Line 2334  int find_lhs_and_lu_mod(unsigned int **a,int row,int c
         return d;          return d;
 }  }
   
   int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo)
   {
           int i,j,k;
           int *rp;
           unsigned int *t,*pivot;
           unsigned int inv,m;
   
           *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int));
           for ( i = 0; i < n; i++ ) rp[i] = i;
           for ( k = 0; k < n; k++ ) {
                   for ( i = k; i < n && !a[i][k]; i++ );
                   if ( i == n ) return 0;
                   if ( i != k ) {
                           j = rp[i]; rp[i] = rp[k]; rp[k] = j;
                           t = a[i]; a[i] = a[k]; a[k] = t;
                   }
                   pivot = a[k];
                   inv = invm(pivot[k],md);
                   for ( i = k+1; i < n; i++ ) {
                           t = a[i];
                           if ( m = t[k] ) {
                                   DMAR(inv,m,0,md,t[k])
                                   for ( j = k+1, m = md - t[k]; j < n; j++ )
                                           if ( pivot[j] ) {
                                                   unsigned int tj;
                                                   DMAR(m,pivot[j],t[j],md,tj)
                                                   t[j] = tj;
                                           }
                           }
                   }
           }
           return 1;
   }
   
 /*  /*
   Input    Input
         a : n x n matrix; a result of LU-decomposition          a : n x n matrix; a result of LU-decomposition
Line 2247  void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
Line 2631  void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
         }          }
 }  }
   
   void Plu_mat(NODE arg,LIST *rp)
   {
           MAT m,lu;
           Q dn;
           Q *v;
           int n,i;
           int *iperm;
           VECT perm;
           NODE n0;
   
           asir_assert(ARG0(arg),O_MAT,"lu_mat");
           m = (MAT)ARG0(arg);
           n = m->row;
           MKMAT(lu,n,n);
           lu_dec_cr(m,lu,&dn,&iperm);
           MKVECT(perm,n);
           for ( i = 0, v = (Q *)perm->body; i < n; i++ )
                   STOQ(iperm[i],v[i]);
           n0 = mknode(3,lu,dn,perm);
           MKLIST(*rp,n0);
   }
   
 void Plu_gfmmat(NODE arg,LIST *rp)  void Plu_gfmmat(NODE arg,LIST *rp)
 {  {
         MAT m;          MAT m;
Line 2258  void Plu_gfmmat(NODE arg,LIST *rp)
Line 2664  void Plu_gfmmat(NODE arg,LIST *rp)
         VECT perm;          VECT perm;
         NODE n0;          NODE n0;
   
         asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");          asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");
         asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");          asir_assert(ARG1(arg),O_N,"lu_gfmmat");
         m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));          m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
         mat_to_gfmmat(m,md,&mm);          mat_to_gfmmat(m,md,&mm);
         row = m->row;          row = m->row;
Line 2557  void inner_product_int(Q *a,Q *b,int n,Q *r)
Line 2963  void inner_product_int(Q *a,Q *b,int n,Q *r)
                         t = wma; wma = sum; sum = t;                          t = wma; wma = sum; sum = t;
                 }                  }
         }          }
         GC_free(wm);          GCFREE(wm);
         GC_free(wma);          GCFREE(wma);
         if ( !sgn ) {          if ( !sgn ) {
                 GC_free(sum);                  GCFREE(sum);
                 *r = 0;                  *r = 0;
         } else          } else
                 NTOQ(sum,sgn,*r);                  NTOQ(sum,sgn,*r);
Line 2615  void inner_product_mat_int_mod(Q **a,int **b,int n,int
Line 3021  void inner_product_mat_int_mod(Q **a,int **b,int n,int
                         t = wma; wma = sum; sum = t;                          t = wma; wma = sum; sum = t;
                 }                  }
         }          }
         GC_free(wm);          GCFREE(wm);
         GC_free(wma);          GCFREE(wma);
         if ( !sgn ) {          if ( !sgn ) {
                 GC_free(sum);                  GCFREE(sum);
                 *r = 0;                  *r = 0;
         } else          } else
                 NTOQ(sum,sgn,*r);                  NTOQ(sum,sgn,*r);
Line 3050  void Pnd_det(NODE arg,P *rp)
Line 3456  void Pnd_det(NODE arg,P *rp)
                 nd_det(0,ARG0(arg),rp);                  nd_det(0,ARG0(arg),rp);
         else          else
                 nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);                  nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);
   }
   
   void Pmat_col(NODE arg,VECT *rp)
   {
           int i,j,n;
           MAT mat;
           VECT vect;
   
           asir_assert(ARG0(arg),O_MAT,"mat_col");
           asir_assert(ARG1(arg),O_N,"mat_col");
           mat = (MAT)ARG0(arg);
           j = QTOS((Q)ARG1(arg));
           if ( j < 0 || j >= mat->col) {
                   error("mat_col : Out of range");
           }
           n = mat->row;
           MKVECT(vect,n);
           for(i=0; i<n; i++) {
                   BDY(vect)[i] = BDY(mat)[i][j];
           }
           *rp = vect;
 }  }
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.62
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>