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

Diff for /OpenXM_contrib2/asir2018/engine/Q.c between version 1.2 and 1.3

version 1.2, 2018/09/21 07:06:51 version 1.3, 2018/09/24 22:26:43
Line 25  void gc_free(void *p,size_t size)
Line 25  void gc_free(void *p,size_t size)
   
 void init_gmpq()  void init_gmpq()
 {  {
   mp_set_memory_functions(Risa_GC_malloc,gc_realloc,gc_free);    mp_set_memory_functions(Risa_GC_malloc_atomic,gc_realloc,gc_free);
   
   mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOZ(ONEMPZ,ONE);    mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOZ(ONEMPZ,ONE);
   gmp_randinit_default(GMP_RAND);    gmp_randinit_default(GMP_RAND);
 }  }
   
   void pmat(Z **a,int row,int col)
   {
     int i,j;
   
     for ( i = 0; i < row; i++, printf("\n") )
       for ( j = 0; j < col; j++, printf(" ") )
         printexpr(CO,a[i][j]);
     printf("\n");
   }
   
 Z utoz(unsigned int u)  Z utoz(unsigned int u)
 {  {
   mpz_t z;    mpz_t z;
Line 161  void mulz(Z n1,Z n2,Z *nr)
Line 171  void mulz(Z n1,Z n2,Z *nr)
   
 void muladdtoz(Z n1,Z n2,Z *nr)  void muladdtoz(Z n1,Z n2,Z *nr)
 {  {
 #if 1  #if 0
   Z t;    Z t;
   
   if ( n1 && n2 ) {    if ( n1 && n2 ) {
Line 171  void muladdtoz(Z n1,Z n2,Z *nr)
Line 181  void muladdtoz(Z n1,Z n2,Z *nr)
         mpz_addmul(BDY(*nr),BDY(n1),BDY(n2));          mpz_addmul(BDY(*nr),BDY(n1),BDY(n2));
         if ( !mpz_sgn(BDY(*nr)) )          if ( !mpz_sgn(BDY(*nr)) )
           *nr = 0;            *nr = 0;
     }    }
 #else  #else
   Z t,s;    Z t,s;
   
Line 183  void muladdtoz(Z n1,Z n2,Z *nr)
Line 193  void muladdtoz(Z n1,Z n2,Z *nr)
   
 void mul1addtoz(Z n1,long u,Z *nr)  void mul1addtoz(Z n1,long u,Z *nr)
 {  {
   #if 0
   Z t;    Z t;
   
   if ( n1 && u ) {    if ( n1 && u ) {
Line 196  void mul1addtoz(Z n1,long u,Z *nr)
Line 207  void mul1addtoz(Z n1,long u,Z *nr)
         if ( !mpz_sgn(BDY(*nr)) )          if ( !mpz_sgn(BDY(*nr)) )
           *nr = 0;            *nr = 0;
     }      }
   #else
     Z t,s;
   
     mul1z(n1,u,&t); addz(*nr,t,&s); *nr = s;
   #endif
 }  }
   
 void mul1z(Z n1,long n2,Z *nr)  void mul1z(Z n1,long n2,Z *nr)
Line 1421  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
Line 1437  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
   int *rind,*cind;    int *rind,*cind;
   int count;    int count;
   int ret;    int ret;
   struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;    struct oEGT eg_mul1,eg_mul2,tmp0,tmp1,tmp2;
   int period;    int period;
   int *wx,*ptr;    int *wx,*ptr;
   int wxsize,nsize;    int wxsize,nsize;
   Z wn;    Z wn;
   Z wq;    Z wq;
   
   init_eg(&eg_mul1); init_eg(&eg_mul2);
   a0 = (Z **)mat->body;    a0 = (Z **)mat->body;
   row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
   w = (int **)almat(row,col);    w = (int **)almat(row,col);
Line 1478  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
Line 1495  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
       MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;        MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;
       MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;        MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;
       wc = (int **)almat(rank,ri);        wc = (int **)almat(rank,ri);
       for ( i = 0; i < rank; i++ )  
         wc[i] = w[i]+rank;  
       *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));        *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
       *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));        *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
   
       period = F4_INTRAT_PERIOD;        period = F4_INTRAT_PERIOD;
       for ( q = ONE, count = 0; ; ) {        for ( q = ONE, count = 0; ; ) {
           /* check Ax=B mod q */
         if ( DP_Print > 3 )          if ( DP_Print > 3 )
           fprintf(stderr,"o");            fprintf(stderr,"o");
         /* wc = b mod md */          /* wc = b mod md */
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) {            for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
             wi[j] = remqi((Q)bi[j],md);              wi[j] = remqi((Q)bi[j],md);
             if ( wi[j] && sgnz(bi[j]) < 0 )          /* wc = A^(-1)wc; wc is not normalized */
               wi[j] = md-wi[j];          solve_by_lu_mod(w,rank,md,wc,ri,0);
           }  
         /* wc = A^(-1)wc; wc is normalized */  
         solve_by_lu_mod(w,rank,md,wc,ri,1);  
         /* x += q*wc */          /* x += q*wc */
   get_eg(&tmp0);
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);            for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
         /* b =(A*wc+b)/md */          /* b =(b-A*wc)/md */
   get_eg(&tmp1); add_eg(&eg_mul1,&tmp0,&tmp1);
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0; j < ri; j++ ) {            for ( j = 0; j < ri; j++ ) {
             u = b[i][j];              mpz_t uz;
             for ( k = 0; k < rank; k++ ) mul1addtoz(a[i][k],wc[k][j],&u);  
               if ( b[i][j] )
                 mpz_init_set(uz,BDY(b[i][j]));
               else
                 mpz_init_set_ui(uz,0);
               for ( k = 0; k < rank; k++ ) {
                 if ( a[i][k] && wc[k][j] ) {
                   if ( wc[k][j] < 0 )
                     mpz_addmul_ui(uz,BDY(a[i][k]),-wc[k][j]);
                   else
                     mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
                 }
               }
               MPZTOZ(uz,u);
             divsz(u,mdq,&b[i][j]);              divsz(u,mdq,&b[i][j]);
           }            }
   get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);
         count++;          count++;
         /* q = q*md */          /* q = q*md */
         mulz(q,mdq,&u); q = u;          mulz(q,mdq,&u); q = u;
         if ( count == period ) {          if ( count == period ) {
           ret = intmtoratm(xmat,q,*nmmat,dn);            ret = intmtoratm(xmat,q,*nmmat,dn);
           if ( ret ) {            if ( ret ) {
               print_eg("MUL1",&eg_mul1);
               print_eg("MUL2",&eg_mul2);
             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 1623  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1654  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
       MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;        MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;
       MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;        MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;
       wc = (int **)almat(rank,ri);        wc = (int **)almat(rank,ri);
       for ( i = 0; i < rank; i++ )  
         wc[i] = w[i]+rank;  
       *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));        *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
       *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));        *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
   
Line 1634  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1663  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
           fprintf(stderr,"o");            fprintf(stderr,"o");
         /* wc = b mod md */          /* wc = b mod md */
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) {            for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
             wi[j] = remqi((Q)bi[j],md);              wi[j] = remqi((Q)bi[j],md);
             if ( wi[j] && sgnz(bi[j]) < 0 )  
               wi[j] = md-wi[j];  
           }  
         /* wc = A^(-1)wc; wc is normalized */          /* wc = A^(-1)wc; wc is normalized */
         solve_by_lu_mod(w,rank,md,wc,ri,1);          solve_by_lu_mod(w,rank,md,wc,ri,1);
         /* x += q*wc */          /* x += q*wc */
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);            for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
         /* b =(A*wc+b)/md */          /* b =(b-A*wc)/md */
         for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
           for ( j = 0; j < ri; j++ ) {            for ( j = 0; j < ri; j++ ) {
             u = b[i][j];              mpz_t uz;
             for ( k = 0; k < rank; k++ ) mul1addtoz(a[i][k],wc[k][j],&u);  
               if ( b[i][j] )
                 mpz_init_set(uz,BDY(b[i][j]));
               else
                 mpz_init_set_ui(uz,0);
               for ( k = 0; k < rank; k++ ) {
                 if ( a[i][k] && wc[k][j] ) {
                   if ( wc[k][j] < 0 )
                     mpz_addmul_ui(uz,BDY(a[i][k]),-wc[k][j]);
                   else
                     mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
                 }
               }
               MPZTOZ(uz,u);
             divsz(u,mdq,&b[i][j]);              divsz(u,mdq,&b[i][j]);
           }            }
         count++;          count++;
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>