[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.10 and 1.11

version 1.10, 2018/12/10 22:24:42 version 1.11, 2019/03/03 05:21:17
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.9 2018/10/19 23:27:38 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.10 2018/12/10 22:24:42 noro Exp $ */
 #include "ca.h"  #include "ca.h"
 #include "gmp.h"  #include "gmp.h"
 #include "base.h"  #include "base.h"
Line 1229  int intmtoratm(MAT mat,Z md,MAT nm,Z *dn)
Line 1229  int intmtoratm(MAT mat,Z md,MAT nm,Z *dn)
     return 0;      return 0;
   row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
   bshiftz(md,1,&t);    bshiftz(md,1,&t);
   isqrtz(t,&s);    isqrt(t,&s);
   bshiftz(s,64,&b);    bshiftz(s,64,&b);
   if ( !b ) b = ONE;    if ( !b ) b = ONE;
   dn0 = ONE;    dn0 = ONE;
Line 1469  void isqrtz(Z a,Z *r)
Line 1469  void isqrtz(Z a,Z *r)
   Z two;    Z two;
   
   if ( !a ) *r = 0;    if ( !a ) *r = 0;
   else if ( UNIQ(a) ) *r = ONE;    else if ( UNIZ(a) ) *r = ONE;
   else {    else {
     k = z_bits((Q)a); /* a <= 2^k-1 */      k = z_bits((Q)a); /* a <= 2^k-1 */
     bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */      bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */
Line 1592  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
Line 1592  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn
   Z wq;    Z wq;
   
 #if SIZEOF_LONG == 8  #if SIZEOF_LONG == 8
   return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp);    return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,0);
 #endif  #endif
 init_eg(&eg_mul1); init_eg(&eg_mul2);  init_eg(&eg_mul1); init_eg(&eg_mul2);
   a0 = (Z **)mat->body;    a0 = (Z **)mat->body;
Line 1706  get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);
Line 1706  get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);
                 if ( !cinfo[j] )                  if ( !cinfo[j] )
                   cind[k++] = j;                    cind[k++] = j;
               return rank;                return rank;
             }              } else
                 goto reset;
           } else {            } else {
   reset:
             period = period*3/2;              period = period*3/2;
             count = 0;              count = 0;
           }            }
Line 1743  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1745  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
   Z wq;    Z wq;
   DP m;    DP m;
   
   #if SIZEOF_LONG == 8
     return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,mb);
   #endif
   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 1815  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1820  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
         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);
         /* wc = A^(-1)wc; wc is normalized */          /* wc = A^(-1)wc; wc is not normalized */
         solve_by_lu_mod(w,rank,md,wc,ri,1);          solve_by_lu_mod(w,rank,md,wc,ri,0);
         /* 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]);
Line 1829  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1834  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
               mpz_init_set(uz,BDY(b[i][j]));                mpz_init_set(uz,BDY(b[i][j]));
             else              else
               mpz_init_set_ui(uz,0);                mpz_init_set_ui(uz,0);
             for ( k = 0; k < rank; k++ ) {              for ( k = 0; k < rank; k++ )
               if ( a[i][k] && wc[k][j] ) {                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]);                    mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
               }  
             }  
             MPZTOZ(uz,u);              MPZTOZ(uz,u);
             divsz(u,mdq,&b[i][j]);              divsz(u,mdq,&b[i][j]);
           }            }
   
         count++;          count++;
         /* q = q*md */          /* q = q*md */
         mulz(q,mdq,&u); q = u;          mulz(q,mdq,&u); q = u;
Line 1859  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
Line 1860  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT 
                 if ( !cinfo[j] )                  if ( !cinfo[j] )
                   cind[k++] = j;                    cind[k++] = j;
               return rank;                return rank;
             }              } else
                 goto reset;
           } else {            } else {
   reset:
             period = period*3/2;              period = period*3/2;
             count = 0;              count = 0;
           }            }
Line 2131  RESET:
Line 2134  RESET:
 }  }
 #endif  #endif
   
 int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp)  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp,DP *mb)
 {  {
   MAT r;    MAT r;
   Z z;    Z z;
Line 2151  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
Line 2154  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
   int count;    int count;
   int ret;    int ret;
   int period;    int period;
     DP m;
   
   a0 = (Z **)mat->body;    a0 = (Z **)mat->body;
   row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
Line 2168  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
Line 2172  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
     if ( DP_Print > 3 ) {      if ( DP_Print > 3 ) {
       fprintf(asir_out,"done.\n"); fflush(asir_out);        fprintf(asir_out,"done.\n"); fflush(asir_out);
     }      }
   
       if ( mb ) {
         /* this part is added for inv_or_split_dalg */
         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 = (mpz_t **)mpz_allocmat(rank,rank); /* lhs mat */      a = (mpz_t **)mpz_allocmat(rank,rank); /* lhs mat */
     b = (mpz_t **)mpz_allocmat(rank,col-rank);      b = (mpz_t **)mpz_allocmat(rank,col-rank);
     for ( j = li = ri = 0; j < col; j++ )      for ( j = li = ri = 0; j < col; j++ )
       if ( cinfo[j] ) {        if ( cinfo[j] > 0 ) {
         /* the column is in lhs */          /* the column is in lhs */
         for ( i = 0; i < rank; i++ ) {          for ( i = 0; i < rank; i++ ) {
           w[i][li] = w[i][j];            w[i][li] = w[i][j];
Line 2181  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
Line 2198  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
             mpz_set_ui(a[i][li],0);              mpz_set_ui(a[i][li],0);
         }          }
         li++;          li++;
       } else {        } else if ( !cinfo[j] ) {
         /* the column is in rhs */          /* the column is in rhs */
         for ( i = 0; i < rank; i++ ) {          for ( i = 0; i < rank; i++ ) {
           if ( a0[rinfo[i]][j] )            if ( a0[rinfo[i]][j] )
Line 2249  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
Line 2266  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
         count++;          count++;
         /* q = q*md */          /* q = q*md */
         mpz_mul_ui(q,q,md);          mpz_mul_ui(q,q,md);
         fprintf(stderr,".");  
         if ( count == period ) {          if ( count == period ) {
           ret = mpz_intmtoratm(x,rank,ri,q,nm,den);            ret = mpz_intmtoratm(x,rank,ri,q,nm,den);
           if ( ret ) {            if ( ret ) {
             for ( j = k = l = 0; j < col; j++ )              for ( j = k = l = 0; j < col; j++ )
               if ( cinfo[j] )                if ( cinfo[j] > 0 )
                 rind[k++] = j;                  rind[k++] = j;
               else                else if ( !cinfo[j] )
                 cind[l++] = j;                  cind[l++] = j;
             ret = mpz_gensolve_check(mat,nm,den,rank,rind,cind);              ret = mpz_gensolve_check(mat,nm,den,rank,rind,cind);
             if ( ret ) {              if ( ret ) {
Line 2272  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
Line 2288  int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *
                 }                  }
               MPZTOZ(den,*dn);                MPZTOZ(den,*dn);
               return rank;                return rank;
             }              } else
                 goto reset;
           } else {            } else {
   reset:
             fprintf(stderr,"F");              fprintf(stderr,"F");
             period = period*3/2;              period = period*3/2;
             count = 0;              count = 0;

Legend:
Removed from v.1.10  
changed lines
  Added in v.1.11

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>